Nonlinear Processes in Geophysics Spatio-temporal error growth in the multi-scale Lorenz ’ 96 model

The influence of multiple spatio-temporal scales on the error growth and predictability of atmospheric flows is analyzed throughout the paper. To this aim, we consider the two-scale Lorenz’96 model and study the interplay of the slow and fast variables on the error growth dynamics. It is shown that when the coupling between slow and fast variables is weak the slow variables dominate the evolution of fluctuations whereas in the case of strong coupling the fast variables impose a non-trivial complex error growth pattern on the slow variables with two different regimes, before and after saturation of fast variables. This complex behavior is analyzed using the recently introduced Mean-Variance Logarithmic (MVL) diagram.


Introduction
The growth of perturbations, or errors, in weather or climate models has been discussed in the literature from different points of view using both simplified "toy" models (Gutiérrez et al., 2008) and real operational systems (Fernández et al., 2008).A common behavior has been identified in these cases.When plotting in a diagram the variance against the mean of the logarithms of perturbations (the so called Mean-Variance Logarithmic (MVL) diagram) a typical trajectory appears as the fingerprint of the spatio-temporal dynamics of the perturbations.These MVL trajectories usually show two distinct regimes.The first, with increasing variance, has been identified as a genuine infinitesimal fluctuation where the increasing variance appears as a result of the collapse to the main Liapunov vector.The second, showing a decrease of variance, appears by nonlinear effects (hence it is a genuine finite fluctuation) that destroy the spatial correlation and localization of the perturbations (Primo et al., 2007;Gutiérrez Correspondence to: S. Herrera (herreras@unican.es)et al., 2008).Hence, this methodology allows the characterization of both the temporal and spatial growth components, and their interplay, explaining the characteristic localized structure of error growth.
In all cases analyzed until now with this methodology, the system had a single temporal scale, and so these studies did not clearly establish the role of the multiple scales existing in the atmosphere.Moreover, little is known about the behavior with more temporal scales.This study is interesting from both fundamental and applied points of view.Many basic questions about the meaning and generation of infinitesimal and finite perturbations in models with multiple scales still remain open.Also, questions of practical interest such as the way of initializing ensembles with coupled ocean atmosphere models are unsolved.
This paper addresses the role of the multiple scales in the dynamics of the perturbation growth.A new complex growth pattern with several regimes of evolution is shown and analyzed.In Sect. 2 we describe the one-and two-scale Lorenz'96 models comparing their behavior.In Sect. 3 a quantitative description of the growth patterns obtained in one-scale L96 model is presented.Section 4 analyzes the growth of perturbations in the two-scale L96 model.Finally, some conclusions are given in Sect. 5.

The two-scale Lorenz'96 model
The Lorenz'96 system (L96) is a simple conceptual model of atmosphere-like multi-scale dynamics (Lorenz, 1996).It consists of m slow variables x i coupled to m×n fast variables y j,i whose evolution is governed by the following nonlinear equations modeling advection, constant forcing, and linear damping: Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
x 1 x 2 y y 1,1 4,1 Fig. 1.Schematic illustrations of the L96 system with m=8 and n=4, with a total of m×n grid points in the inner circle.
where i = 1,...,m, j = 1,...,n; h is the coupling constant, b and c are the scaling constants and F is a constant forcing.
Both the x i 's and the y j,i 's are assumed to be periodic, i.e., x m+1 = x 1 ; x 0 = x m and y n+1,i = y 1,i+1 ; y 0,i = y n,i−1 .L96 was originally introduced to mimic multi-scale mid-latitude weather, considering an unspecified scalar meteorological quantity at m equidistant grid points along a latitude circle (see Fig. 1); moreover, each of these points was associated with n fast variables providing a higher resolution description of the system (m×n grid points).
In this paper we considered the conventional parameter values F =8 and m=40.This leads to a chaotic behavior with a time unit corresponding approximately to five days in an equivalent atmospheric model (Lorenz, 1996).This is the time unit used throughout the paper.We also selected n=4, b=10 and c=10, leading to a two-scale model where the fast variables fluctuate 10 times faster than the slow ones, but with a smaller amplitude.The coupling parameter h can be varied in the range (0, 1) to obtain different degrees of coupling between both scales.This model has been extensively used in the literature to study the influence of multiple spatio-temporal scales on the predictability of atmospheric flows.For instance, Wilks (2005) studied stochastic parametrization, Lorenz and Emanuel (1998) analyzed targeted observations and Orrell (2003) predictability over different timescales.
A simplified version of L96 has also been used to investigate the dynamics of the slow variables alone.In this case, the system is given by: Both systems were integrated using a fourth-order Runge-Kutta method with dt=0.005.
Figure 2a shows trajectories obtained by integrating (2) from a random initial condition z i (0).After a short transient these figures show the typical pattern of L96 model emulating atmospheric flows: waves with mean length of 5 units travel westward in a noisy background induced by chaos spending 7-9 time units to complete a latitude circle.Figure 2b and c shows the trajectory for the two-scale model (1) integrated from a random initial condition (x i (0),y j,i (0)) = (z i (0),y j,i (0)) and weak coupling (h=0.3)The z i and x i trajectories are initially very close, but they diverge very soon due to non-linear behavior and the effect of the fast variable dynamics.Contrarily, the spatial patterns of the X and Y trajectories are, after rescaling, very similar at all times.Figure 2f and g shows the trajectories of the two-scale model with the same initial conditions but with a stronger coupling (h=1).The main difference with the previous case can be seen in the fast variable, that now exhibits a high variability at small scales.In order to appreciate this fact we plot in Fig. 2d and h the trajectories of a reference point x 1 , y 1,1 and the coupling term S 1 = c n j =1 y j,1 .Note that this latter term is the true contribution, as an additive source, of the fast variables to the slow ones at a given point i = 1. Figure 2d shows that with weak coupling the dynamics of the fast variable is conditioned by the slow variable.It is a coupling dominated by the slow dynamics.Figure 2h shows the opposite effect, clearly the dynamics is dominated by fast variables.This effect is seen as producing high frequency spatial and temporal structures.
In all cases the last time step of this trajectory was used as the initial condition in the following section, in order to have a flow-assimilated initial condition.

Spatio-temporal growth of perturbations in one-scale systems
In this section we describe the growth of perturbations (e.g., errors in the initial conditions) in a one-scale system such as the one generated by Eq. ( 2).First, we show a qualitative description showing the phenomena of spatial localization and saturation.Then, we present a more quantitative description taking logarithms of the absolute perturbations and using the MVL diagram to obtain the characteristic curve for the spatio-temporal error growth.

Qualitative description
Let us consider a perturbed trajectory around a control forecast z i (t) for the one-scale L96 model, which is computed integrating (2) from a flow-assimilated initial condition z i (0) (see previous section).The perturbed trajectory corresponds to the integration forward in time of the perturbed initial condition z i (0) = z i (0) + δz i (0), where δz i (0) are normallydistributed independent random perturbations (e.g.observational errors in the initial conditions) with zero mean and  The initial perturbations were randomly chosen according to a Gaussian distribution with standard deviation σ for δz and δx, and 0 for δy.Middle row: the same as in the upper row but in logarithmic scale.Lower row: the same as in the middle row but taking away the mean at each time.In order to ease their comparison, panels a, g, m are respectively identical to panels (d, j, p).Notice also the different scales in each case.
a given standard deviation σ , defined as a fraction of its stationary value, σ = 10 −4 σ st , with σ st ∼ 3.2.The spatiotemporal evolution of the ensemble is characterized by the non-infinitesimal fluctuations δz i (t) ≡ z i (t) − z i (t) between the perturbed member and the control trajectory.These fluctuations follow an equation that can be written separating linear and non-linear terms as: O(δz 2 ) representing terms of second, or larger, order.When fluctuations are small, nonlinear terms can be neglected and the evolution can be considered as infinitesimal.As the growth is exponential, due to the multiplicative character of the linear Eq. (3), in a short time the nonlinear terms become important and a new evolution emerges, which is not infinitesimal (finite fluctuations).Figure 3a shows an example of the growth of perturbations in one scale systems.This figure shows that, after some time (around t = 5) with this exponential growth, the perturbations reach locally the amplitude of the system (the scale of the plot is the same as the corresponding panel in Fig. 2) and saturate.Before that time, the spatio-temporal growth of errors follows a characteristic dynamics where temporal and spatial components interplay leading to characteristic localized growth patterns.In order to see these initial stages of the spatial evolution, Fig. 3g shows the same growth process, but in logarithmic scale, that is, we represent the log-perturbations In logarithmic scale the progressive growth in time is easily observed.Moreover, to see the localization phenomenon graphically Fig. 3m shows the deviation of the log-perturbation to its spatial mean, h i − 1 m h i .The initially Gaussian perturbations get quickly localized in a flowdependent structure which is kept until saturation due to nonlinear effects.
This behavior is represented in Fig. 4a-e, which shows the histograms of the perturbations δz i (t) at five different times: the initial (t = 0) and saturated (t = 7) times (with Gaussianlike distributions) and three intermediate stages indicated with dashed lines in Fig. 3m (with localized distributions).To get a higher statistical significance we used 20 random initial perturbations to build the histograms.The behavior shown is typical of the growth of a random initial perturbation in a spatio-temporal chaotic system with only one active scale.Initially the spatial distribution is Gaussian by construction (panel a).The further evolution shows a progressive localization (Fig. 4b, c, d) whose distributions exhibits a non-gaussian profile.This spatial localization is an effect of chaotic systems in their infinitesimal evolution (when perturbations are in the tangent manifold).It is well known that initial perturbations in chaotic systems tend to their main Liapunov vector (Baker and Gollub, 1996).Roughly speaking we can say that the initial random perturbation is going to the chaotic attractor adopting the form of the main Liapunov vector.The localization process ends when either the Liapunov vector is reached or the amplitude is large enough to induce the appearance of nonlinear effects.As said above, when nonlinear terms are important a second regime appears where the localization gained is destroyed (Fig. 4e).(see Gutiérrez et al., 2008, for more details).

The MVL approach
Spatially localized structures follow logarithmic statistical distributions (log-normal, log-Poissonian,...) and thus they are better analyzed in their logarithmic representation.Moreover, it is known that the logarithmic transformation of infinitesimal perturbations behave as growing rough surfaces, which allows a precise analysis by using techniques borrowed from statistical physics (López et al., 2004).The Mean-Variance Logarithmic (MVL) diagram (Gutiérrez et al., 2008) was introduced with these premises in mind.The MVL diagram is based on the log-perturbations, defined as the logarithm of the absolute value of the perturbations defined above, i.e., h i (t) = log|δz i (t)|.It is a two-dimensional representation displaying in the abscissa the spatial average of the log-perturbations (the M-value): and in the ordinate the spatial variance (the V-value): The evolution (growth) of the initial perturbation is given by the parametric curve (M(t), V (t)) in the diagram.It can be shown that the horizontal axis represents the spatiallyaveraged error growth in time (related to the leading Liapunov exponent λ as M(t) ∼ λt), whereas the vertical evolution represents the spatial dynamics given by a changing spatial correlation length (characteristic spatial scale) of the perturbation field.The resulting spatio-temporal growth is a combined effect of both interplaying components (Gutiérrez et al., 2008), and its dynamics can be characterized with the MVL diagram.
As an example, Fig. 5c represents the MVL diagram for a perturbation growth similar to Fig. 3a.The curves shown correspond to an average of 10 different initial conditions (flow-assimilated by integrating a random initial condition during 10 time units and 20 random initial perturbations for each initial condition.This figure shows two typical patterns of systems with a single scale and with random Gaussian initial perturbations.Note that Gaussianity imposes a fix value of V = 1.26 at the initial time-step (Gutiérrez et al., 2008).The temporal evolution of the M-and V-values are displayed in Fig. 5 panels a and b, respectively.It can be shown that whereas the first magnitude grows until saturation, the later grows and decays indicating a gain and decay www.nonlin-processes-geophys.net/17/329/2010/ Nonlin.Processes Geophys., 17, 329-337, 2010 of spatial structure, or localization, in the original representation.In order to see the effect of the amplitude of the initial perturbation, two different values of σ were considered: 10 −4 (as in Fig. 3a) and 10 −6 .When the initial perturbation is small enough (the case σ = 10 −6 ) V saturates due to the finite size of the system, before the action of non-linear terms.
In other words, the initial perturbation, that tends to the main Liapunov vector, reaches this limit in an infinitesimal evolution before being perturbed by the non linear terms.
As we can see in this example, the MVL Diagram is a powerful tool to analyze the dynamics of initial perturbations.Its implementation and physical interpretation are very simple which increases their practical value: An increase (decrease) of variance (V-value) means an approach to (separation from) the tangent manifold of the chaotic attractor, which can be used as a quantifier of dynamic assimilation.Finally note that its use is not limited to the analysis of initial perturbations but, in general, it is valid for the analysis of any dynamical state as those appearing in data assimilation and synchronization (Szendro et al., 2009).

Growth of perturbations with multiple scales
The fluctuation dynamics with multiple scales is very dependent of the intensity of the coupling between both scales.In the case of weak coupling the fast and slow variables follow a similar dynamics of fluctuations: The saturation time of both variables is close (Fig. 3h and i), and the patterns of lo-calization are also similar (Fig. 3n and o).This result is not surprising since, as seen in Sect.2, the fast scale is deeply conditioned by the slow scale.The spatial and temporal variations of their fluctuations are then similar.
The strong coupling case is completely different.The growth of perturbations follows a more complex pattern.Firstly, the saturation time of the fast variable is shorter than that of the corresponding slow variable (Fig. 3l and k).The fast variable saturates at t=1, while the slow variable keeps a pattern similar to the simplified L96 case (Fig. 3l, k, and j).Secondly, the fast variable in the case of strong coupling shows a finer spatial resolution (Fig. 3l and r) than the corresponding ones in the weak coupling case (Fig. 3i and o).In fact, the characteristic pattern of traveling waves described before is more clearly observed in the case of strong coupling.Finally, the localization patterns shown by panels in the bottom row of Fig. 3 show a different spatial behavior, since the saturation of the fast variables produces a cross-over in the slow ones with two separate regimes.Before the saturation of the fast variables, the growth of the slow ones follow the characteristic regime where perturbations get progressively localized.This is better appreciated in the histograms (see Fig. 4f and g).When nonlinear saturation effects appear in the fast variables (yielding the Gaussian-like spatial distribution of perturbations in Fig. 4m), the slow variables are also forced to lose the localized structure (Fig. 4h).Afterwards, when the fast variables are fully saturated (cross-over point), the slow ones recover the localized spatial structure (Fig. 4i) until they also saturate (Fig. 4j).These two characteristic regimes have not been reported elsewhere and are a fundamental knowledge to understand the different dynamical behavior of perturbations in multi-scale spatio-temporal models, including weather models; for instance, atmosphere-ocean coupled models show different time scales for the different components: fast atmosphere and slow ocean.

MVL diagrams
In order to analyse these regimes in detail, Fig. 6 shows the MVL diagram for the one-and two-scale L96 models in three cases of weak (Fig. 6c), intermediate (Fig. 6g) and strong (Fig. 6j) coupling parameters.In the weak coupling case the diagram of the slow variable is very similar to the one of the one-scale system.The dynamics of the fast variable is conditioned to the dynamics of the slow variable, which drives the whole system.The slope of the logarithmic amplitude in the growth process, which is closely related with the first Liapunov exponent, is similar in the cases of one or two-scale systems (Fig. 6a).In the strong coupling case, the MVL Diagram shows a complicated behavior of the slow variable (Fig. 6j).Since now the fast scale dominates the evolution of the slow one there is an initial evolution regime, until the saturation of the fast scale.In this initial regime, the slow variable follows the dynamics of the fast variable with a greater slope in the logarithmic amplitude (Fig. 6h) than in the weak coupling configuration.There is a second regime, after the saturation of the fast variable, where the dynamics of the slow variables becomes independent and the slope of the logarithmic amplitude is smaller.In terms of the variance (Fig. 6b, f, and i) there are two effects on the fast variable when increasing the coupling parameter.The first is the earlier decay of the variance with the increasing coupling.This is related to the fast arrival to the non-linear barrier caused by the faster growth of the amplitude.The second effect is that increased coupling leads to a larger maximum variance (i.e.localization of the perturbation) reached by the fast variable.
That is, the stronger the coupling, the more relevant the fast variable becomes and, thus, the finer the spatial detail it can develop.
It is worth remarking the complex patterns of the slow variable in this case.Four regimes are observed in the MVL diagram (Fig. 6j).The first two are induced by the behavior of the fast variables.After the saturation of this variable the slow one continues increasing the structure until its own saturation.In Fig. 6h and i, it is apparent that when the amplitude of the fast variable is saturated (around t = 2), the variance continues its evolution in a slowing down process that lasts until the saturation time of the slow variable (around t = 7−8).
For the sake of completeness, we represent in Fig. 7 the behavior of fluctuations in the whole interval of the coupling parameter.We plot the maximum value of V for the fast scale, which gives the strength of spatial localization, as a function of parameter h.Three regions for the coupling parameter (labeled as weak, medium or strong) exhibiting different behaviors can be observed.As expected, in the weak coupling zone, h ∈ (0,0.3),V has a value close to the one of the single scale and no interaction between slow and fast dynamics takes place; in the intermediated zone, h ∈ (0.3,0.7), V grows linearly with h, and interactions start taking place.Finally, the variance V remaining saturated in the strong coupling zone h ∈ (0.7,1) and two different regimes appear in the slow variables as a consequence of the interplay with fast ones.

Fast dynamics and effective noise
It is interesting to investigate the role played by the fast variables in the dynamics of a multi-scale system.The possibility of parametrization of fast degrees of freedom by an equivalent noise is important to reduce the complexity of a problem (Just et al., 2001).For the sake of illustration we explore this possibility in our case.
In a case with weak coupling, where the slow variable is dominant, the fast variables play a secondary role and they could be substituted by a weak additive noise without altering the main dynamics.However, when the coupling is strong, the existence of two distinct regimes impede the existence of an equivalent effective noise, even to reproduce the second regime, when fast variables are saturated and the slow ones become dominant.We show in Fig. 8 the effect of adding a random Gaussian noise in a one-scale model with an intensity equivalent to those of the fast source term (i.e. with the same mean and standard deviation).The noise is added at the saturation time but we checked that the effect is similar when it is placed in any other previous time.We can see that the equivalent noise does not reproduce the dynamics of the two-scale model.The added noise is not dynamically assimilated, V decreases until the noise is assimilated, and a very different pattern is observed in the MVL diagram (Fig. 8).Note that starting in a previous time is not a solution since at these times the fast variable is dominant and cannot be substituted by a passive component.How to find a dynamically compatible noise to reduce the two-scale model in this case is a difficult problem that exceeds the scope of this paper.Figure 8 is a good illustration of this problem.Note that this may have important implications in order to define appropriate initial perturbations for ensemble prediction in coupled atmosphere-ocean models.

Conclusions
The spatio-temporal dynamics of the error growth in the twoscale model proposed by Lorenz (1996) has been characterized by means of the MVL diagram.For comparison purposes, a simplified version of the model with a single scale was also considered.This latter model was already characterized by Gutiérrez et al. (2008) and evolves from a small random initial perturbation by increasing the mean error and the localization until the finite-size and non-linear effects become large and destroy the localization gained reaching a stationary state.
Depending on the strength of the coupling between the fast and slow variables the two-scale model presents a more or less complex picture.The weak coupling does not introduce appreciable changes with respect to the one-scale system.Contrarily, strong coupling introduces a complex behavior, where the errors of the fast variables rapidly gain and lose local structure, affecting the slow variables in a non-trivial manner.The loss of structure arising from the saturation of the fast variables induces a loss also in the slow ones.When the fast variables saturate, the slow ones gain structure again forming a secondary hump, which ends as expected when the non-linear barrier is reached.On the other hand, the slow variables also have an effect on the fast ones.Namely, the localization of the fast variables does not reach a constant value (full saturation) until the slow ones saturate.
Finally, we tried to reproduce the second hump experienced by the slow variables in the MVL diagram by substituting the fast variables (once saturated) by an equivalent noise term in the governing equation for the slow variables, but our attempt was unsuccessful.
This complex behavior may be relevant for other systems such as the atmosphere-ocean system where there are fast varying variables (in the atmosphere) gaining and losing the structure much faster than others (in the ocean, for instance).This issue could be made evident by analyzing atmosphereocean coupled model output from existing ensemble prediction systems.Of course, this is out of the scope of the present work and will be dealt with in a forthcoming paper.

Fig. 2 .
Fig. 2.Comparison of the trajectories of the simplified (Eq.2) and two-scale (Eq. 1) L96 models with weak and strong coupling parameters; all integrated forward in time from the same initial conditions.The panels show the trajectories of the one-scale model (a, e), the slow (b) and fast (c) variables of the two-scale model with weak coupling h=0.3 and the time evolution (d) of a single point in space (i = 1) in the slow variable (black line), in one of the fast variables associated to this point (gray line) and the contribution of the fast variables to this point at each time (dotted line).Panels (f-h) are equivalent to (b-d) but for the case of strong coupling h=1.Note that panels (a) and (e) are identical in order to facilitate a visual comparison between patterns in each row.Notice also the different scales used to depict slow and fast variables.

Fig. 3 .
Fig. 3. Upper row: dynamics of a single realization of perturbations of the initial condition for (a, d) the simplified L96 model,and (b, c) the slow-fast variables, δx i and δy j,i of the L96 model with weak coupling h=0.3 and (e, f) with strong coupling h=1.The initial perturbations were randomly chosen according to a Gaussian distribution with standard deviation σ for δz and δx, and 0 for δy.Middle row: the same as in the upper row but in logarithmic scale.Lower row: the same as in the middle row but taking away the mean at each time.In order to ease their comparison, panels a, g, m are respectively identical to panels (d, j, p).Notice also the different scales in each case.

Fig. 4 .
Fig. 4. Histograms of 20 perturbations of an initial condition for (a)-(e) the simplified L96 model, (f)-(j) the slow variables x i of the L96 model with strong coupling (h=1), and (k)-(o) the fast variable y j,i of the L96 model.For each column, the different rows from top to bottom represent the initial, three intermediate and the final times, corresponding to the horizontal lines depicted in Fig. 3.
Fig. 5. (a) Temporal evolution of the M-component of the growth of perturbations of the example shown in Fig. 3a, (b) temporal evolution of the V-component, (c) MVL diagram.Two different values of σ for the standard deviation of the initial perturbations have been considered.

Fig. 6 .
Fig. 6.Time series of the mean (a, e, h) and variance (b, f, i) of the log-perturbations of a simplified L96 model (black) and the slow (x i , dark gray) and fast (y j,i , light gray) components of a two-scale L96 model.(c, g, j) MVL diagrams corresponding to the previous components.The upper row(a, b, c) is for weak coupling h=0.26, the medium row (e, f, g) is for medium coupling h=0.6 and the lower row is for strong coupling h=1.All curves are averages of 200 realizations (see text).

Fig. 7 .
Fig. 7. Maximum value of the variance of the log-perturbations of the fast components y j,i of the two-scale L96 model (V component of the MVL-diagram) as a function of the coupling parameter h.

Fig. 8 .
Fig. 8. Growth of perturbations in the two-scale model (x i ;y j,i ) with strong coupling (h=1) versus the single-scale model (z i ) and the single-scale model with an equivalent noise (x * i ): (a) temporal evolution of the M-component of the growth of perturbations, (b) temporal evolution of the V-component, (c) MVL diagram.