the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Analysis of model error in forecast errors of extended atmospheric Lorenz 05 systems and the ECMWF system

### Holger Kantz

Forecast error growth as a function of lead time of atmospheric phenomena is caused by initial and model errors. When studying the initial error growth, it may turn out that small-scale phenomena, which contribute little to the forecast product, significantly affect the ability to predict this product. The question under investigation is whether omitting these atmospheric phenomena will improve the predictability of the resulting value. The topic is studied in the extended Lorenz (2005) system. This system shows that omitting small spatiotemporal scales that significantly affect prediction ability will reduce predictability more than modeling it. In other words, a system with model error (omitting phenomena) will not improve predictability. A hypothesis explaining and describing this behavior is developed, with the difference between systems (model error) produced at each time step seen as the error of the initial conditions. The resulting model error is then defined as the sum of the increments of the time evolution of the initial conditions so defined. The hypothesis is compared to the fit parameters that define the model error in certain approximations of the average forecast error growth. Parameters are interpreted in this context, and the approximations are used to estimate the errors described in the hypothesis. A method is proposed to distinguish increments of prediction error growth from small-spatiotemporal-scale phenomena and model errors. Results are presented for the error growth of the ECMWF system, where a 40 % reduction in model error between 1987 and 2011 is calculated based on the developed hypothesis, while over the same period the instability (error growth rate) of the system with respect to initial condition errors has grown.

- Article
(2423 KB) - Full-text XML
- BibTeX
- EndNote

Forecast errors in numerical weather prediction systems grow over time due to the inaccuracy of the initial state (initial error), which is amplified by the chaotic nature of the system itself and model imperfections (model error). In the setting of classical low-dimensional chaos, one would observe an exponential error growth of any tiny initial error whose exponent is given by the largest Lyapunov exponent of the system, with some saturation when the error reaches the magnitude of the standard deviation of the quantity to be predicted. In contrast to this, several authors have observed in the past (Toth and Kalnay, 1993; Lorenz, 1969; Aurell et al., 1996, 1997; Boffetta et al., 1998) that the proper Lyapunov exponent of a dynamical system might not be a relevant description of the initial error growth. Brisch and Kantz (2019) and Zhang et al. (2019) associated initial error growth with scale-dependent error growth, where tiny errors grow much faster than larger ones. Lorenz (1996) gave a sketch of such error growth: a typical quantity to be predicted is a superposition of the dynamics on different scales. After fast growth of small-scale errors with saturation at these very same small scales, the large-scale errors continue to grow at a slower rate until even these saturate. Therefore, Lyapunov exponents of structures of various spatiotemporal scales are taken as the previously mentioned scale-dependent quantity, and they determine the error growth on their respective scales. It is also evident that in practice initial errors are not infinitesimal in the mathematical sense, and therefore the exponential growth of infinitesimal errors might be irrelevant for the growth of forecast errors in operational weather forecasts.

In numerical weather predictions, the average forecast error as function of lead time is influenced by many deviations from a simple exponential growth: there is the saturation effect of larger errors, the potential scale dependence of the instability, and also the fact that real initial errors are not infinitesimal and might not point into the locally most unstable direction. Moreover, if the model error is due to neglecting small-scale and fast phenomena, it is highly fluctuating along the model trajectory. In order to describe such effects, the literature contains different phenomenological approximations for the time derivative of the average error magnitude (error growth rate or tendency) as a function of the error magnitude in numerical weather predictions. We will briefly recall these here since we will use suitable fits to the observed error growth as function of error magnitude for our atmospheric model simulations later. We will show that initial errors of magnitudes that are comparable to real weather forecasts do not play a dominant role in our studied model systems, while model errors do play such a role. Moreover, we will present an explanation of the observed error growth in terms of an averaged model error called the “drift”.

In low-dimensional-bounded chaotic systems with at least one positive Lyapunov exponent, the growth of infinitesimal errors is exponential for a finite time interval, given by a linear time derivative:

where *E*(*t*) is the error magnitude, *t* is time, and *λ* is the largest Lyapunov exponent of the system. Since *E*_{exp}(*t*) in Eq. (1) grows unboundedly, this can be true only as long as *E*(*t*) is small since every error has to saturate at the latest when it has grown to the order of magnitude of the diameter of the attractor (the invariant set). This saturation effect was considered by Lorenz (1982), who introduced the “quadratic hypothesis” *E*_{qv}(*t*):

where *E*_{lim} is the limit (saturation) value of the error magnitude. As a function of time, the error *E*_{qv}(*t*) shows a sigmoidal shape; see Fig. 5b.

For a scale-dependent error growth in the spirit of Lorenz (1969), Brisch and Kantz (2019) proposed using a power law divergence of the effective, scale-dependent Lyapunov exponent $\mathit{\lambda}\left(E\right)\propto {E}^{-b}$, which gives the following time evolution:

where the exponent *b* connects Lyapunov exponents and limit errors of the different scales (Brisch and Kantz, 2019), while the coefficient *a* determines the degree of the scales' coupling (Bednář and Kantz, 2022). The forecast error then grows as a power law in time, *E*_{p}(*t*), with a very fast growth rate when it is still small and a slow growth rate when it is large. Since *E*_{p}(*t*) (see Eq. 3) again grows unboundedly, Bednář and Kantz (2022) introduced the extended power law *E*_{ep}(*t*) that allows saturation using the same trick as Lorenz (1982):

Zhang et al. (2019) described scale-dependent error growth differently. They took a two-parametric hypothesis:

where *λ*_{r} is a synoptic-scale error growth rate and *β*_{r} is an upscale error growth rate from small-scale processes. If ${\mathit{\beta}}_{\mathrm{r}}/{\mathit{\lambda}}_{\mathrm{r}}$ is large, this leads to a super-exponential growth of small errors and to the classical exponential error growth when *E*_{r}(*t*) is large. We can and should again include saturation of the error by the factor $(\mathrm{1}-E/{E}_{\text{lim}})$:

The two-parametric model d*E*_{r} Eq. (5) was originally designed to describe initial and model error growth (Leith, 1978). In this interpretation, *λ*_{r} is the largest Lyapunov exponent of the system (similar to *λ*_{exp} in Eq. 1), while *β*_{r} is the model error source term due to the imperfect representation of the atmosphere. In addition, d*E*_{q}, Eq. (6) is called the quadratic hypothesis with model error for the same reason as for d*E*_{r} (Savijarvi, 1995; Dalcher and Kalnay, 1987), although in contrast to Eq. (2) it includes a constant term and therefore allows for some skewness in $\mathrm{d}E/\mathrm{d}t$ as a function of *E*.

When we present the results of our error growth numerical analysis, we particularly present the error growth rate as a function of error magnitude. This allows us to better distinguish between these different error growth models than studying the error magnitude as a function of time, even though error magnitude as function of time is relevant in predictions.

While the above-listed error growth approximations are supposed to approximate the effectively observed average error growth in operational forecasts, let us now focus on the model error. The model is given as a set of the first order in time differential equations of the form $\frac{\mathrm{d}\mathit{X}}{\mathrm{d}t}=\mathit{G}\left(\mathit{X}\right(t\left)\right)$, where ** X** is a (high-dimensional) phase space vector that describes the current state of the atmosphere and

**(**

*G***) is a vector-valued function that defines the rate of change of this vector at every possible state. In operational weather forecasts, the core of such a system in given by wind speed, pressure, density, and temperature, and the minimal setting for**

*X***is then called the “primitive equations” (Phillips, 1973). Following (part of) the meteorological literature, we will call the right-hand side**

*G***(**

*G***) the “model tendency”, while in the context of dynamical systems it is called the “vector field”.**

*X*Following, e.g., Orrell et al. (2001), the model error *G*_{e} at a model space point $\stackrel{\mathrm{\u0303}}{\mathit{X}}\left(t\right)$ is described as the difference between the model vector field (tendency) $\frac{\mathrm{d}\mathit{X}}{\mathrm{d}t}=\mathit{G}\left(\mathit{X}\right(t\left)\right)$ and the observed time derivative (tendency) $\frac{\mathrm{d}\stackrel{\mathrm{\u0303}}{\mathit{X}}}{\mathrm{d}t}$ of the projection $\stackrel{\mathrm{\u0303}}{\mathit{X}}\left(t\right):=\mathit{P}\left(\stackrel{\mathrm{\u0311}}{\mathit{X}}\right(t\left)\right)$ of the “reality” $\stackrel{\mathrm{\u0311}}{\mathit{X}}\left(t\right)$ into the model space:

Let us stress that in operational forecasts, since we do not know the perfect model, the true time derivative $\frac{\mathrm{d}\stackrel{\mathrm{\u0303}}{\mathit{X}}}{\mathrm{d}t}$ is only known by observation, while in our later model studies, we have a mathematical expression for the vector field of reality as well. In a very strong simplification, one could assume that the absolute value of *G*_{e} is, on average, the constant *β* in Eq. (5), which irrespective of initial condition errors will lead to a deviation of the model solution from reality. While it is evident how to define the model error in a single time step, we will later discuss how model errors propagate in time and how model errors at different positions along a trajectory accumulate, and we will also introduce the notion of drift for that purpose.

The main issue about forecasts is how far into the future they might be useful. The “prediction horizon” quantifies this as the time when the forecast error has grown to a certain percentage of the climatological uncertainty of the forecast target, where the latter is approximated by *E*_{lim} in the above error growth assumptions.

For exponential growth ${E}_{\text{exp}}\left(t\right)={E}_{\mathrm{0}}{e}^{{\mathit{\lambda}}_{\text{exp}}t}$ and for an initial error *E*_{0} going to zero, the time *t*_{lim} at which the error reaches a limiting value *E*_{lim} goes to infinity:

However, a strict predictability limit *t*_{lim} exists for scale-dependent error growth even when the initial error *E*_{0} vanishes (Palmer et al., 2014; Brisch and Kantz, 2019). For a description by a power law d*E*_{p} (see Eq. 3) the predictability limit *t*_{lim} is as follows:

For an exponential growth with a non-zero *β*_{r} parameter d*E*_{r} (see Eq. 5), the prediction horizon *t*_{lim} is as follows:

Scale-dependent error growth implies that both model assumptions *E*_{p} and *E*_{r} grow faster than exponentially when errors are small, thereby limiting the prediction horizon because further and further improvements of the precision of the initial condition are compensated by a faster initial error growth. In the context of weather prediction, this means that the influence of small-scale atmospheric phenomena, which contribute little to the final value, significantly affect the ability to predict this value. Figures 1–3 show such behavior simulated in the extended Lorenz (2005) system with one (L05-1), two (L05-2), and three (L05-3) spatiotemporal scales (see Appendix A for more information on these systems). Figure 1a shows the values of the L05-1 system variables at a given time. Because this is a single spatiotemporal scaled system, the average growth of the initially small error is exponential. The two initially nearby trajectories begin to diverge significantly in this setting after about 30 d (Fig. 1b). Adding a considerably smaller scale (L05-2 system) that does not significantly affect the overall value in sum (Fig. 2a) reduces the closeness of the two initially nearby trajectories of an overall variable by 10 d (Fig. 2b). By adding a third medium scale (Fig. 3a, L05-3 system), the two initially nearby trajectories of an overall variable start to diverge significantly after about 10 d (Fig. 3b), which is about 3 times earlier than for the L05-1 system. This is a consequence of the much faster growth of the small-scale errors.

Including small spatiotemporal scales, i.e., improving the model's spatial and temporal resolution, therefore enhances the instability (error growth rate) with respect to initial condition errors. The question under investigation in this paper is whether omitting small-scale atmospheric phenomena, which contribute little to the final value, will improve the predictability of the resulting value. In other words, how does the average forecast error growth change in a model where small-scale phenomena are omitted but where model errors are therefore introduced, compared to a model where all phenomena are present but the average forecast error growth is scale-dependent.

Buizza (2010), Magnusson and Kallen (2013), or Jacobson (2001) show that improving the model's spatial and temporal resolution will improve the prediction ability, especially for short forecast ranges (Buizza, 2010). However, the cited studies work with models that do not model small spatiotemporal phenomena (they are parameterized) and whose initial condition error magnitude is larger than the magnitude of these phenomena. We have verified the fact that the high-resolution model (that models small scales) is less stable than the low-resolution model (that does not model small scales) against initial condition errors (Bednář and Kantz, 2022; Budanur and Kantz, 2022) and that therefore the issue of omitting small scales has another facet. Our new approach models and omits small spatiotemporal scales using the one- and two-scale Lorenz (2005) system (L05-1 and L05-2) and its three-scale extension (L05-3) introduced before in Bednář and Kantz (2022). The omitted scale is the small scale for the L05-2 system, and the small and medium scale are for the L05-3 system. L05 system definition and further details can be found in Appendix A.

The protocol for how to measure the initial error growth for the L05 systems is defined in Sect. 2.1, and the results are presented and compared in Sect. 3.1. The model error scenario, where the L05-1 system is the model and the L05-2 and L05-3 systems are reality, is defined in Sect. 2.2, and the results are presented and compared in Sect. 3.2. The variant with initial and model error is defined in Sect. 2.3, and the results are presented and compared in Sect. 3.3. The results of different error growth scenarios are compared and discussed in Sect. 3.4. Section 4.1 includes the calculation of the model error (drift) defined in Sect. 2.4 and a hypothesis linking the error so defined with the growth of the average model error determined by the difference between model and reality. The meaning of the model error source *β* in d*E*_{r} of Eq. (5) and d*E*_{q} Eq. (6) and how to link the value of *β* with the value of the model error (drift) is discussed and explained in Sect. 4.2. Section 5 presents a similar analysis for the ECMWF forecast system data. Conclusions and discussions are then presented in the final section.

The average error magnitude for L05 systems is calculated numerically using the method introduced by Lorenz (1996, 2005). Generally, we define an “error” as the distance between two trajectories where one, the reference trajectory, is supposed to be the “truth”, while a second trajectory is generated under perturbation of the initial condition, under perturbation of the dynamical equations, or both. We measure the error magnitude *e*(*t*) after fixed time intervals. We then calculate the mean error magnitude *E*(*t*) after fixed times, calculate the average growth tendency $\frac{\mathrm{d}E}{\mathrm{d}t}$ during the last time interval, and report the mean error magnitudes vs. time and the mean growth tendencies (rates) vs. mean error magnitudes.

An alternative method for calculating scale-dependent error growth is called the “finite-size Lyapunov exponent” (Aurell et al., 1996, 1997; Boffetta et al., 1998; Cencini and Vulpiani, 2013). In brief, a finite-size Lyapunov exponent $\mathit{\lambda}=(\mathrm{1}/E)(\mathrm{d}E/\mathrm{d}t)$ or finite-size error growth tendency $\mathrm{d}E/\mathrm{d}t$ can be defined as the ergodic average over phase space of the growth rate of perturbations of a given magnitude *E*, where the growth rate is defined as the inverse of the time *t*_{f} needed for the error magnitude to increase by a pre-defined factor *f*, hence $\mathit{\lambda}\left(E\right)=(\mathrm{1}/{t}_{f})\mathrm{ln}f$. We choose the former method because it is closer to the process of calculating the average forecast error magnitude of numerical weather prediction systems (Lorenz, 1982; Savijarvi, 1995; Froude et al., 2013; Zhang et al., 2019) and because it is more consistent with the performance of forecasts.

For numerical weather prediction systems, errors in initial conditions and model errors (inaccurate representation of atmospheric processes by the model) are sources of prediction inaccuracy. For the L05 systems, we simulate the initial error growth (perfect model assumption), the model error growth (perfect initial conditions assumption), a combination of both (initial and model error assumption), and the model error growth as defined by Orrell et al. (2001) (drift assumption). To calculate the average error magnitude, a reference trajectory (considered the truth or verification) and a trajectory that is the numerical solution of the systems with a given error are repeatedly generated. For this scheme to be meaningful, we have to ensure that the reference trajectory is on the system's attractor and that the repetition of this scheme samples the whole attractor with correct weights (the invariant measure). We solve this issue in the following way: we first integrate the system over 10 years (175 200 steps), starting from arbitrary initial conditions, and assume that after discarding this transient that the trajectory is on the attractor. We continue to integrate this single trajectory and consider segments of it as reference trajectories for error growth, i.e., the many reference trajectories are simply segments of one very long trajectory, which ensures not only that all these segments are located on the attractor but also that they sample the attractor according to the invariant measure.

## 2.1 Initial error growth

When we use the term “initial error growth”, we denote the growth of errors in the initial conditions, which limit predictability if a system is chaotic. In order to numerically determine the largest Lyapunov exponent, we have to ensure that initial perturbations already point into the locally most unstable direction since otherwise errors might even shrink over short time periods (this is also a relevant issue in ensemble forecasts, and there its solution is found using bred vectors; see Toth and Kalnay, 1997). We solve these issues in the following way: we start with a random perturbation of the reference trajectory of very small amplitude and let this trajectory evolve over time before determining its distance toward the reference trajectory. In other words, we discard some initial time interval of error growth since this is affected by some transient behavior before it starts to grow with the maximum Lyapunov exponent.

We calculate the initial error growth in systems with one (L05-1), two (L05-2), and three (L05-3) scales to illustrate the behavior in systems with a different number of spatiotemporal scales. The three spatial scales *X*_{1}, *X*_{2}, and *X*_{3} for the L05-3 system and two spatial scales *X*_{1} and *X*_{2} for the L05-2 system cannot be separated in terms of a coordinate transform but are intrinsically coupled and superimposed in the variables *X*_{tot} of the system. The initial conditions of the reality for L05-3 and L05-2 systems are called ${X}_{\text{tot},\mathrm{0},n}$, from which one finds ${X}_{\mathrm{1},\mathrm{0},n}$, ${X}_{\mathrm{2},\mathrm{0},n}$, and ${X}_{\mathrm{3},\mathrm{0},n}$ through Eqs. (A10), (A11), and (A12), respectively, for the L05-3 system and ${X}_{\mathrm{1},\mathrm{0},n}$ and ${X}_{\mathrm{2},\mathrm{0},n}$ through Eqs. (A4) and (A5), respectively, for the L05-2 system. The initial conditions of the reality for the L05-1 system are called *X*_{0,n}. The initial values of the “prediction” are then for the L05-1 system ${X}_{\mathrm{0},n}^{\prime}={X}_{\mathrm{0},n}+{e}_{n}\left(\mathrm{0}\right)$, where *e*_{n}(0) are the initial errors randomly selected from the normal distribution $\text{ND}\phantom{\rule{0.25em}{0ex}}(\mathit{\mu}=\mathrm{0};\mathit{\sigma}=\mathrm{0.01})$. Since the system's state *X*_{tot} is the sum over all spatiotemporal components, for L05-3 and L05-2 systems any arbitrary but small error with spatially uncorrelated components affects only the smallest-scale component. Only a spatially correlated initial error would appear in another component. However, since this error would immediately propagate into the small-scale variables and then grow fastest in these, a perturbation with initial errors in the smallest-scale component is the only practical choice. The initial values of the prediction for the L05-3 system are then ${X}_{\text{tot},\mathrm{0},n}^{\prime}={X}_{\mathrm{1},\mathrm{0},n}+{X}_{\mathrm{2},\mathrm{0},n}+{X}_{\mathrm{3},\mathrm{0},n}+{e}_{\mathrm{3},n}={X}_{\text{tot},\mathrm{0},n}+{e}_{\mathrm{3},n}$, where *e*_{3,n}(0) is the initial error randomly selected from the normal distribution $\text{ND}\phantom{\rule{0.25em}{0ex}}(\mathit{\mu}=\mathrm{0};\mathit{\sigma}=\mathrm{0.001})$. The initial error *e*_{2,n}(0) of the L05-2 system is also randomly selected from the same normal distribution, where the initial value of the prediction is ${X}_{\text{tot},\mathrm{0},n}^{\prime}={X}_{\mathrm{1},\mathrm{0},n}+{X}_{\mathrm{2},\mathrm{0},n}+{e}_{\mathrm{2},n}={X}_{\text{tot},\mathrm{0},n}+{e}_{\mathrm{2},n}$.

From the initial values of reality and prediction, we integrate the L05 system equations (Eqs. A1, A8, and A9) for 41.7 d (*K*=2000 steps). In each time step, *k* of the numerical integration ${X}_{\mathit{\tau},k,n}$ and ${X}_{\mathit{\tau},k,n}^{\prime}$ are obtained. The size of the error at a given time *k*Δ*t* is ${e}_{\mathit{\tau},n}(k\cdot \mathrm{\Delta}t)={X}_{\mathit{\tau},k,n}^{\prime}-{X}_{\mathit{\tau},k,n}$, where $k=\mathrm{1},\mathrm{\dots},K$, $n=\mathrm{1},\mathrm{\dots},N$ (*N*=360 variables for all used systems). *τ* defines a scale or sum of scales ($\mathit{\tau}=\text{tot},\mathrm{1},\mathrm{2},\mathrm{3}$ (L05-3 system), $\mathit{\tau}=\text{tot},\mathrm{1},\mathrm{2}$ (L05-2 system)) and is therefore omitted for the L05-1 system. We perform *M*=400 runs to calculate the average error growth. In each new run, the initial values ${X}_{\mathit{\tau},\mathrm{0},n}$ are the last values ${X}_{\mathit{\tau},K,n}$ of the previous run. The average initial error growth *E*(*t*) is calculated as the geometric mean of the runs of the Euclidean distances between reality and prediction:

The geometric mean is chosen because of its suitability for comparison with growth governed by the largest Lyapunov exponent. For further information, see Bednář et al. (2014) or Ding and Li (2011). As a result, we have numerical averages for the error growth as a function of time steps after perturbing the reference trajectories in the full phase space and for each scale. We can convert these results into the error growth tendency (rate) as a function of the error magnitude.

## 2.2 Model error growth

By “model error growth”, we denote the growth of errors caused by the inaccurate description of reality by the “model”. This inaccuracy involves small-scale atmospheric processes unresolved by the model, which for numerical weather prediction systems are approximated to the resolved scale by a procedure called parameterization. It also denotes model biases that are either unknown or have not yet been addressed (Allen et al., 2006). It is a common expectation that model errors in numerical weather forecasts can be reduced by improving the spatial and temporal resolution of the forecast system.

To simulate this in the L05 systems environment (Appendix A), we use the L05-2 and L05-3 systems as the reality and the L05-1 system as the model. Thus, the unresolved or unknown scale is the small scale for the L05-2 system and the small and medium scale for the L05-3 system. This approach is justified by the fact that the L05-2 and L05-3 systems can be viewed as a variant of the L05-1 system:

where ${\stackrel{\mathrm{\u0303}}{F}}_{n}\left(t\right)={b}^{\mathrm{2}}[{X}_{\mathrm{2}},{X}_{\mathrm{2}}{]}_{\mathrm{1},n}+c[{X}_{\mathrm{2}},{X}_{\mathrm{1}}{]}_{\mathrm{1},n}-b{X}_{\mathrm{2},n}+F$ for the L05-2 system and ${\stackrel{\mathrm{\u0303}}{F}}_{n}\left(t\right)={b}_{\mathrm{1}}^{\mathrm{2}}[{X}_{\mathrm{2}},{X}_{\mathrm{2}}{]}_{\mathrm{1},n}+{b}_{\mathrm{2}}^{\mathrm{2}}[{X}_{\mathrm{3}},{X}_{\mathrm{3}}{]}_{\mathrm{1},n}+{c}_{\mathrm{1}}[{X}_{\mathrm{2}},{X}_{\mathrm{1}}{]}_{\mathrm{1},n}+{c}_{\mathrm{2}}[{X}_{\mathrm{3}},{X}_{\mathrm{2}}{]}_{\mathrm{1},n}-{b}_{\mathrm{1}}{X}_{\mathrm{2},n}-{b}_{\mathrm{2}}{X}_{\mathrm{3},n}+F$ for the L05-3 system are treated as a forcing, which varies in a complicated manner with time. We parameterize these small-scale phenomena contained in ${\stackrel{\mathrm{\u0303}}{F}}_{n}\left(t\right)$ by the average value of these phenomena, which is close to zero, and therefore we can write:

where 〈…〉 represents the mean calculated over a long orbit on the L05-2 and L05-3 system attractors.

To calculate the average model error growth, we first define initial conditions that are the same for model and reality (perfect initial conditions assumption) and are determined from the values ${X}_{\text{tot},\mathrm{0},n}$ of reality (L05-2 or L05-3 systems) at the end of the initial transient. Let us stress that we can use ${X}_{\text{tot},\mathrm{0},n}$ of our high-resolution L05-3 or L05-2 system without any projection as the initial state of the L05-1 system and that the lack of smaller scales is only expressed by the lack of feedback from the smaller scales in the equation of motion.

From these initial values, we integrate the L05-2 or L05-3 system equations (reality) and the L05-1 system equations (model) forward for 41.7 d (*K*=2000 steps). In each time step *k* of the numerical integration, ${X}_{\text{tot},k,n}$ (reality) and *X*_{k,n} (model) are obtained. The size of the error at a given time *k*Δ*t* is ${e}_{M,n}(k\cdot \mathrm{\Delta}t)={X}_{\text{tot},k,n}-{X}_{k,n}$, where $k=\mathrm{1},\mathrm{\dots},K$, $n=\mathrm{1},\mathrm{\dots},N$ (*N*=360 variables for all used systems). We perform *L*=400 runs to calculate the average error growth. In each new run, the initial values ${X}_{\text{tot},\mathrm{0},n}$ are the last values ${X}_{\text{tot},K,n}$ of the previous run. The average model error growth *E*_{M}(*t*) is calculated as the geometric mean of the runs of the Euclidean distances between reality and model:

As a result, we have numerical averages for the model error growth as a function of time steps. Note that in this framework, only ${X}_{\text{tot},k,n}$ (L05-2 and L05-3 systems) values compared to *X*_{k,n} (the L05-1 system) and not the individual scales. We can convert these results into the error growth tendency (rate) as a function of the error magnitude.

## 2.3 Initial and model error growth

By “initial and model error growth”, we denote the combination of the initial error growth defined in Sect. 2.1 and the model error growth defined in Sect. 2.2. We describe the L05-2 and L05-3 systems as reality and the L05-1 system with perturbations in the initial conditions of reality as “model prediction”.

In this setting, we do not discard the initial time interval of initial error growth because this transition period is negligible compared to the model error growth. The initial conditions of the reality for L05-3 and L05-2 systems are called ${X}_{\text{tot},\mathrm{0},n}$ and determined in the same way described above. The initial values of the model prediction for the L05-1 system are then ${X}_{\mathrm{0},n}^{\prime}={X}_{\text{tot},\mathrm{0},n}+{e}_{n}\left(\mathrm{0}\right)$, where *e*_{n}(0) values are the initial errors randomly selected from the normal distributions $\text{ND}(\mathit{\mu}=\mathrm{0};\mathit{\sigma}=\mathrm{0.01})$ and $\text{ND}(\mathit{\mu}=\mathrm{0};\mathit{\sigma}=\mathrm{0.2})$. From initial values, we integrate the L05-2 or L05-3 system equations (reality) and the L05-1 system equations (model prediction) forward for 41.7 d (*K*=2000 steps). In each time step *k* of the numerical integration, ${X}_{\text{tot},k,n}$ (reality) and ${X}_{k,n}^{\prime}$ (model prediction) are obtained. The size of the error at a given time *k*Δ*t* is ${e}_{M+\text{ie},n}(k\cdot \mathrm{\Delta}t)={X}_{\text{tot},k,n}-{X}_{k,n}^{\prime}$, where $k=\mathrm{1},\mathrm{\dots},K$, $n=\mathrm{1},\mathrm{\dots},N$ (*N*=360 variables for all used systems). We perform *L*=400 runs to calculate the average error growth. In each new run, the initial values ${X}_{\text{tot},\mathrm{0},n}$ are the last values ${X}_{\text{tot},K,n}$ of the previous run. The average initial and model error growth *E*_{M+ie}(*t*) is calculated as the geometric mean of the runs of the Euclidean distances between reality and model prediction:

As a result, we have numerical averages for the initial and model error growth as a function of time steps. Note that in this framework only ${X}_{\text{tot},k,n}$ (L05-2 and L05-3 systems) values are compared to *X*_{k,n} (the L05-1 system) and not the individual scales. We can convert these results into the error growth tendency (rate) as a function of the error magnitude.

## 2.4 Drift

Section 2.2 describes how we can numerically measure the effects of the model error on forecast accuracy. However, if we want to understand how the model error drives the model trajectory away from reality, we need an additional concept. The reason is that model errors at different positions along the trajectory are only weakly correlated. This is a consequence of the fact that the lack of small scales and fast degrees of freedom in the model equations dominates model errors. However, if model errors at different positions along a trajectory are uncorrelated, then they can partially compensate each other, and their effect is not the same as if we assume that model errors along a trajectory are about the same everywhere. Therefore, We will recall the concept of drift as discussed by Orrell et al. (2001). For these purposes, let us first generally define the model (L05-1 system in our case) as $\mathrm{d}\mathit{X}\left(t\right)/\mathrm{d}t=\mathit{G}\left(\mathit{X}\right(t\left)\right)$, where ** X**∈ℝ

^{n}is the model state space vector (

*n*=360 in our case) and the reality state space vector is $\stackrel{\mathrm{\u0303}}{\mathit{X}}\left(t\right)\in {\mathbb{R}}^{\stackrel{\mathrm{\u0303}}{n}}$. In general, $\stackrel{\mathrm{\u0303}}{n}\ne n$, and it is necessary to project $\stackrel{\mathrm{\u0303}}{\mathit{X}}$ from the state space of reality to the state space of model (Data Assimilation for Numerical Prediction Models). In our case, $\stackrel{\mathrm{\u0303}}{n}=n=\mathrm{360}$, $\stackrel{\mathrm{\u0303}}{\mathit{X}}={\mathit{X}}_{\text{tot}}$, and we use either the L05-2 system $\mathrm{d}{\mathit{X}}_{\text{tot}}\left(t\right)/\mathrm{d}t=\stackrel{\mathrm{\u0303}}{\mathit{G}}\left({\mathit{X}}_{\mathrm{1}}\right(t),{\mathit{X}}_{\mathrm{2}}(t\left)\right)$ or the L05-3 system $\mathrm{d}{\mathit{X}}_{\text{tot}}\left(t\right)/\mathrm{d}t=\stackrel{\mathrm{\u0303}}{\mathit{G}}\left({\mathit{X}}_{\mathrm{1}}\right(t),{\mathit{X}}_{\mathrm{2}}(t),{\mathit{X}}_{\mathrm{3}}(t\left)\right)$ as reality. The model error

*G*_{e}at the point

*X*_{tot}(

*t*) is then the difference between the model vector field (tendency) and the tendency of the projection of reality into the model space. In our case, we can write the following equation:

The drift vector ** d**(

*τ*) was introduced by Orrell et al. (2001) as

This is an accumulation of model errors along a piece of the model trajectory. As we will see in numerical simulations (Fig. 11), the absolute value of drift |** d**(

*τ*)| will not grow approximately linearly in time, i.e., it is not the same as accumulating the absolute value of the model error |

*G*_{e}| along the same piece of the trajectory.

This is a consequence of the here-considered case of neglected small-scale motion: since the ignored scales fluctuate quickly, the model errors at successive positions on the trajectory lose their correlations. We checked this for our L05 models explicitly by calculating the auto-correlation function of the drift vectors as a function of their time lag and found a very fast decay within a few time steps. Therefore, different from model errors in low-dimensional systems, which can be assumed to be spatially highly correlated, one here accumulates random vectors, and the drift therefore follows a path that resembles a Brownian path, as already suggested in Orrel et al. (2001). There and in Orrell (2002) it is also shown how to approximate the integral by summing a series of short-time model errors over finite time steps Δ*t*. The absolute value of drift |** d**(

*τ*)| as a function of

*τ*grows sublinearly, as will be demonstrated later, which gives a more realistic estimate of the role of model errors. What, however, is ignored here is that a model error in the first time step creates a kind of initial condition error for the second time step, which would then grow as an initial condition error. We will discuss this later.

To calculate the average drift *D* comparable to previous cases, we first calculate the time evolution of reality ${X}_{\text{tot},k,n}$ (L05-2 or L05-3 systems), calculated from the initial conditions after the transient period. From each time step *k* of the time evolution of ${X}_{\text{tot},k,n}$ reality (up to *K*=2000 steps), we calculate the one-step Δ*t* time evolution of the model. ${X}_{\text{tot},k,n}={X}_{k,n}$ values are therefore viewed as initial conditions for the one-step Δ*t* time evolution of the model. The size of the drift at a given time *k*Δ*t* is the sum of all previous and current error vectors: ${\mathit{d}}_{n}(k\cdot \mathrm{\Delta}t)={\sum}_{j=\mathrm{0}}^{k-\mathrm{1}}\left({\mathit{X}}_{j,n}\right((j+\mathrm{1})\cdot \mathrm{\Delta}t)-{\mathit{X}}_{\text{tot},j+\mathrm{1},n})$, where $k=\mathrm{1},\mathrm{\dots},K$, $n=\mathrm{1},\mathrm{\dots},N$ (*N*=360 variables for all used systems). Notice that it is not the absolute values of the Δ*t* errors that are accumulated but the vectors (see Fig. 4), meaning that in the summation there can be cancellation effects and hence a slower-than-linear growth of the drift with time.

We perform *L*=400 runs in order to calculate the average error growth. In each new run, the initial values ${X}_{\text{tot},\mathrm{0},n}$ are the last values ${X}_{\text{tot},K,n}$ of the previous run. The average drift *D*(*t*) is defined as the geometric mean of the runs of the Euclidean distances between reality and model:

As a result, we have numerical averages for the drift as a function of time step. We can convert these results into the drift growth tendency as a function of the drift magnitude.

Based on the described methods, we calculate the average prediction error growth for L05 systems. We approximate the numerical error growth curves using the hypotheses Eqs. (1)–(6) and try to identify the most appropriate description. We use these results to determine how the average forecast error growth changes in a model where small-scale phenomena are omitted, but the model error is therefore created (perfect initial conditions assumption or initial and model error assumption) compared to a model where all phenomena are present, whereas the average forecast error growth is scale-dependent (perfect model assumption). The resulting behavior will be explained using the drift.

## 3.1 Initial error growth

Figure 5a shows the initial error growth rate (tendency) $\mathrm{d}E/\mathrm{d}t$ as a function of the error magnitude *E* for the L05-1 system, while Fig. 5b shows the error magnitude as a function of time. We also show the best-fit results of the error growth models represented by Eqs. (1)–(6). It turns out that the initial part of the growth rate is linear without any significant offset, i.e., we see a linear increase with a beginning at ($E=\mathrm{0},\mathrm{d}E/\mathrm{d}t=\mathrm{0}$). Therefore, constants in the error growth models that were included to represent the model error are consequently close to zero. In addition, the power law fit yields a power close to 1. Because of the saturation of the error after some time, the error growth rate decays to zero when the error is large, which can be represented well by the factor ($\mathrm{1}-E/{E}_{\text{lim}}$) in the error growth models. Hence, all models with this saturation term allow good fits to the error growth rate and the error magnitude as a function of time in the whole range and confirm that the L05-1 system is a classical chaotic system with the largest Lyapunov exponent of about *λ*≈0.33 d^{−1}.

The behavior is obviously different for the L05-2 system, which contains additional small-scale degrees of freedom, as shown in Fig. 6. The initial part of the error growth rate (for small *E*) is already curved, and hence the exponential growth model does not provide a good fit anymore. Introducing a non-vanishing error growth rate right from the beginning, i.e., starting from ($E=\mathrm{0},\mathrm{d}E/\mathrm{d}t={\mathit{\beta}}_{\mathrm{r}}$), which is the description by d*E*_{r}, the approximation moves closer to the data, but this is in clear contradiction to the initial error growth idea. Due to the lack of model errors, the growth rate starts from 0. In addition, the quadratic hypothesis is unable to reproduce this curvature well enough. Therefore, the data are best approximated by the power law in the initial part and by the extended power law with saturation on the whole range.

What we found for the initial error growth of the L05-2 system is even more pronounced in the L05-3 system with three spatiotemporal scales. The superiority of approximations d*E*_{p} and d*E*_{ep} over the other approximations is enhanced by the even faster growth of *E*_{tot}(*t*) compared to the exponential growth and *E*_{tot}(*t*) for the L05-2 system (Fig. 7). The reason for this behavior is described in Brisch and Kantz (2019) or in Bednář and Kantz (2022). Lorenz's (1969) statement can summarize it: a typical quantity to be predicted is a superposition of the dynamics on different scales. After a fast growth of the small-scale errors with saturation at these very same small scales, the large-scale errors continue to grow at a slower rate until even these saturate. This is the phenomenon of scale-dependent error growth. We also see that if we interpret the three systems L05-1 to L05-3 as low- and high-resolution models, the high-resolution model has larger instability and hence a shorter time until an ensemble of initial conditions has spread out on the attractor. If this were of relevance for the prediction horizon, then the high-resolution model would be less useful for forecasting than the low-resolution model.

## 3.2 Model error growth

We use the L05-2 system as reality and make forecasts using the L05-1 system. Their suitably averaged differences give rise to the model error as a function of lead time. Figure 8a (solid black curve) shows the model error growth rate $\mathrm{d}E/\mathrm{d}t$ as a function of the error magnitude *E*, while Fig. 8b shows the time evolution of this error. We see an initially very fast error growth caused by the differences in the equations of motion of reality and the model. After a short transient, we see in both panels of Fig. 8a behavior compatible with our error growth models. Those models with a constant term (i.e., the quadratic hypothesis with model error and the exponential growth with model error) provide the best fits, where to be good in the whole range the factor $(\mathrm{1}-E/{E}_{\text{lim}})$ of the quadratic hypothesis with model error is needed. In view of what will follow, we stress that based on the data both *E*_{r} and *E*_{q} provide good fits up to error magnitudes of about three units, with different values *λ*_{q}≈0.27 d^{−1} and *λ*_{r}≈0.17 d^{−1} of the largest Lyapunov exponent.

When we use L05-3 as reality and L05-1 as the model, the same conclusions are valid for the model error growth rate $\mathrm{d}E/\mathrm{d}t$ as a function of the error magnitude *E* (Fig. 9a) and *E*_{M}(*t*) (Fig. 9b). Note, however, that the rates $\mathrm{d}E/\mathrm{d}t\left(E\right)$ have much larger maximal values and that *E*_{M}(*t*) grows faster than when taking L05-2 as reality. However, if we ignore the very initial part of the error growth rate for small values *E*, which the error growth models cannot reproduce, we see that *E*_{r} and *E*_{q} provide the best fits.

## 3.3 Initial and model error growth

In both settings, we also show the results when we include a small initial condition error in addition to the model error. This initial condition error implies that the forecast error as a function of time starts with a non-zero value and correspondingly with a much lower growth rate than the model error alone, but apart from that there are no strong effects. Figures 8 and 9 show that it is not the net sum of the initial error growth *E*_{tot}(*t*) and the model error growth *E*_{M}(*t*). *E*_{M+ie}(*t*) goes from the initial value *E*_{M+ie}(0) through some transition period to the model error growth curve *E*_{M}(*t*). *E*_{M}(*t*) is therefore the limiting value to which *E*_{M+ie}(*t*) is attracted. Indeed, solid and dashed black curves in the insets of Figs. 8b and 9b show that already after time *t*=0.2 the model error alone has grown so much that there is no effect of even of the larger initial condition error of magnitude *E*(0)=0.2 anymore. The larger the initial error and the smaller the model error, the longer the transition period, but it is still short for realistic values of the initial condition error. For these reasons, it can be seen that the appropriate approximation for describing the variant with initial and model error remains the same as for describing the variant with model error only, which is the exponential growth with model error, d*E*_{r} Eq. (5), for the early growth phase and the quadratic hypothesis with model error d*E*_{q}, Eq. (6), for the entire length of the evolution (Figs. 8 and 9).

## 3.4 Comparison of initial and model error growth

We want to use the approximation formulae to construct the error curves for *E*_{tot}(0)=0, =0.1, and =0.2. The initial error magnitudes of 0.1 and 0.2 correspond to the relative values of the initial errors of current numerical weather prediction models for the L05 models. For the initial error growth *E*_{tot}(*t*) (for simplicity, let us redefine *E*_{tot}(*t*) to *E*_{ie}(*t*)), we use the extended power law solution and find the parameter values $\mathrm{d}{E}_{\text{ep}}=\mathrm{0.28}\cdot {E}^{\mathrm{0.66}}(\mathrm{1}-E/\mathrm{7})$ for the L05-2 system and $\mathrm{d}{E}_{\text{ep}}=\mathrm{0.38}\cdot {E}^{\mathrm{0.41}}(\mathrm{1}-E/\mathrm{7.1})$ for the L05-3 system, with initial values *E*_{ie(0)}(0)→0 (Fig. 10, solid red curve), *E*_{ie(0.1)}(0)=0.1 (Fig. 10a, dashed red curve), and *E*_{ie(0.2)}(0)=0.2 (Fig. 10b, dashed red curve).

For the model error growth *E*_{M}(*t*), we use the quadratic hypothesis with model error with the following best-fit parameters: $\mathrm{d}{E}_{\mathrm{q}}\left(t\right)/\mathrm{d}t=(\mathrm{0.27}\cdot E+\mathrm{0.34})(\mathrm{1}-E/\mathrm{7.6})$ and *E*_{M}(0)≈0.1 for the L05-2 system and $\mathrm{d}{E}_{\mathrm{q}}\left(t\right)/\mathrm{d}t=(\mathrm{0.38}\cdot E+\mathrm{1.47})(\mathrm{1}-E/\mathrm{7.8})$ with *E*_{M}(0)≈0.1 for the L05-3 system.

The reason why *E*_{0} is non-zero when using the approximation can be found in Sect. 4.2. We can use the same approximations for the initial and model error growth *E*_{M+ie}(*t*) as for the model error growth alone, *E*_{M}(*t*) with ${E}_{M+\text{ie}\left(\mathrm{0.1}\right)}\left(\mathrm{0}\right)=\mathrm{0.1}$ for the L05-2 system, and ${E}_{M+\text{ie}\left(\mathrm{0.2}\right)}\left(\mathrm{0}\right)=\mathrm{0.2}$ for the L05-3 system. The justification can be found in Sect. 4.3.

In addition to the graphical representation (Fig. 10), we compare the variants by expressing the times *t*_{95 %}, *t*_{71 %}, *t*_{50 %}, and *t*_{25 %} when the error magnitude *E* reaches 95 %, 71 %, 50 %, and 25 %, respectively, of its limiting (saturated) value *E*_{lim}. In the literature, *t*_{95 %} is understood as a practical predictability limit, while *t*_{71 %} corresponds to climatic variability, according to Savijarvi (1995). For Fig. 10, let us first comment on the difference between the limiting (saturation) values *E*_{lim} of the initial error curves *E*_{ie}(*t*) and the model error curves *E*_{M}(*t*) or *E*_{M+ie}(*t*). The pure initial error curves are produced in the perfect model scenario, meaning that forecast and reality are obtained with the same system. When model error is present, the variability of forecast and reality is different, and hence the limiting error is larger, in agreement with Simmons et al. (1995) or Li et al. (2018).

Figure 10 also provides insight into our question of whether high-resolution or low-resolution models will produce better forecasts. We recall that high-resolution models suffer from much faster initial condition error growth due to the stronger instability of small-scale motion. However, the small scales usually do not contribute to the forecast target. Figure 10 shows that predictability is significantly worsened by the model error rather than by the initial error growth. This means that the average forecast error in a model where small-scale phenomena are omitted but the model error is therefore present (black curves *E*_{M}(*t*) and *E*_{M+ie}(*t*)) grows significantly faster compared to a model where all scales are resolved (no model error) but the average forecast error growth is scale-dependent (red curves *E*_{ie}(*t*)). However, this phenomenon depends on the relative magnitudes of the model error term *β*, the (effective) Lyapunov exponent *λ*, and the initial condition error. If the growth rate *λ* were larger and the model error smaller, then the exponential error growth would overwhelm the contribution of model error in every time step. In our numerical experiments using the L05 models, the magnitude of the initial condition error is tuned to values corresponding to initial condition errors in real weather forecasts, so we are tempted to believe that in weather forecasts high-resolution models should also be superior.

Specifically, for the L05-2 system (Fig. 10a), the forecast with model error (i.e., using L05-1 for the forecast), the time *t*_{25} is more than 3 times shorter than the limit without model and initial condition error (*E*_{ie(0)}(0)→0), being ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{4}$ d and ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{13}$ d, respectively. With increasing error magnitude, the ratio of forecast times decreases (${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{6}$ d vs. ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{19}$ d and ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{9}$ d vs. ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{24}$ d) until ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{16}$ d, which is approximately half as large as ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{37}$ d. Adding the error of the initial condition does not significantly change the error growth for the model error variant (${E}_{M}\left(t\right)\approx {E}_{M+\text{ie}\left(\mathrm{0.1}\right)}\left(t\right)$; see Sect. 4.3 for details). The error growth naturally increases for the variant with the initial error, and the ratio is reduced to twice the growth rate over the entire growth period (${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.1)}}=\mathrm{4}$ d vs. ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.1)}}=\mathrm{9}$ d, ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.1)}}=\mathrm{6}$ d vs. ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.1)}}=\mathrm{14}$ d, ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.1)}}=\mathrm{9}$ d vs. ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.1)}}=\mathrm{19}$ d, and ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.1)}}=\mathrm{16}$ d vs. ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.1)}}=\mathrm{32}$ d).

For the L05-3 system (Fig. 10b) taken as truth, the effect is even more dramatic. Without initial error (respectively for *E*_{ie(0)}(0)→0), ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{1}$ d is 7 times smaller than ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{7}$ d. Gradually, the ratio decreases (${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{2}$ d vs. ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{12}$ d and ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{4}$ d vs. ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{18}$ d) until ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M}=\mathrm{8}$ d, which is approximately 4 times smaller than ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0)}}=\mathrm{34}$ d. Adding the error of the initial condition does not significantly change the error growth for the model error variant (${E}_{M}\left(t\right)\approx {E}_{M+\text{ie(0.2)}}\left(t\right)$; see Sect. 4.3 for details). The growth changes for the variant with the initial error, and the ratio is reduced to 5 times the growth rate in the first half of growth (${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.2)}}=\mathrm{1}$ d vs. ${t}_{\mathrm{25}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.2)}}=\mathrm{5}$ d, ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.2)}}=\mathrm{2}$ d vs. ${t}_{\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.2)}}=\mathrm{10}$ d, ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.2)}}=\mathrm{4}$ d vs. ${t}_{\mathrm{71}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.2)}}=\mathrm{16}$ d, and ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},M+\text{ie(0.2)}}=\mathrm{8}$ d vs. ${t}_{\mathrm{95}\phantom{\rule{0.125em}{0ex}}\mathit{\%},\text{ie(0.2)}}=\mathrm{34}$ d).

In the numerical error growth study described so far, we used L05-1 as a model and compared its performance to two types of reality with different spatiotemporal resolutions. We complemented this study using as a single reality the L05-3 system and comparing the performance of the L05-1 and L05-2 models when using them for forecasts as two forecasts with different spatiotemporal resolutions. These findings (not shown) confirm the dominance of the model error over the initial condition error growth, i.e., the model with the higher resolution again provides better forecasts even for larger errors and therefore has a better prediction horizon. However, the effect in this setting is not as dramatic as in the results presented above .

In this section, we present considerations that explain the dominance of the model error in our forecast studies, where the notion of the drift (Sect. 2.4) plays a relevant role. We will consider the feedback of the drift as an initial error, which modifies the interpretation of the parameters in the exponential growth with model error d*E*_{r} (Eq. 5) and the quadratic hypothesis with model error d*E*_{q} (Eq. 6) and its relationship with the drift.

## 4.1 Drift and its role in explaining the model error growth

Figure 11 compares the model error growth *E*_{M}(*t*) (Eq. 14), the drift *D*(*t*) (Eq. 18), and the exponential growth approximation with model error *E*_{r}(*t*) (Eq. 5). The curves are determined from the difference between the L05-2 and the L05-1 systems (Fig. 11a) and between the L05-3 and the L05-1 systems (Fig. 11b) for error magnitudes where the saturation effect is not yet present (up to 2 and 3 d). We intend to understand the behavior of the model error since this is the issue in real forecasts. The drift *D*(*t*) (Fig. 11, blue curve) describes the early part of the model error evolution *E*_{M}(*t*) (Fig. 11, black curve) very well, while at longer lead times it is the exponential growth with model error contribution (Fig. 11, dashed black line) which can be fitted well to the model error curve. Notice, however, that the drift is calculated numerically from the simulation data, like the model error curve, and hence does not contain any free parameters, while the exponential growth with model error has two free parameters, which allow us to optimally achieve the agreement between the solid and dashed black curves.

In particular, in the L05-3 model, for the times where the model error growth *E*_{M}(*t*) can be fitted by *E*_{r}(*t*), the drift *D*(*t*) first overestimates the model error growth *E*_{M}. Following this, a significantly slower growth of the drift *D* relative to the model error growth *E*_{M} is observed, corresponding to a decrease in the growth rate (tendency) of the drift $\mathrm{d}D/\mathrm{d}t\left(D\right)$ compared to the increase in the model error growth rate (tendency) (insets in Fig. 11).

As already mentioned in Sect. 2.4, the drift can be viewed as the sum of the displacements from reality, created at each time step of the model, similar to how an error in the initial conditions will create an initial displacement at the initial time. If we interpret the drift increment $D\left({t}_{k}\right)-D\left({t}_{k-\mathrm{1}}\right)$ at each time step $\mathrm{\Delta}t={t}_{k}-{t}_{k-\mathrm{1}}$, $k=\mathrm{1},\mathrm{\dots},K$ as a new initial error at time *t*_{k}, then (similar to Orrell et al., 2001) we can model the model error growth *E*_{M} by applying the same time evolution assumption to $D\left({t}_{k}\right)-D\left({t}_{k-\mathrm{1}}\right)$ as the initial error growth, i.e., the exponential growth ${e}^{{\mathit{\lambda}}_{D}t}$ driven by the largest Lyapunov exponent *λ*_{D} of the model (L05-1 system). However, this growth should set in only at some time in the future since $\mathit{D}\left({t}_{k}\right)-\mathit{D}\left({t}_{k-\mathrm{1}}\right)$ does not point into the locally most unstable direction (see Sects. 2.1 and 3.1 for a description of the initial error growth for the L05-1 system). We approximate this behavior in two different ways and will explore which one gives better results. (Fig. 12). In the first, $D\left({t}_{k}\right)-D\left({t}_{k-\mathrm{1}}\right)$ evolves with time *t*_{i} in a constant approximation as follows:

while in a linearly decaying approximation it evolves as follows:

*M* and *σ* are parameters that we fix empirically. We propose the hypothesis *E*_{D}(*t*) as a description of the model error growth *E*_{M}(*t*) based on the sum of the individual increments of the drift:

where ap is the symbol for the constant (con) or linear (lin) approximation.

In comparison to the exponential error growth with model error (with or without saturation), this is a modification in two relevant details. First, the model error is not the same at all time steps into the future like the constant *β*_{r} or *β*_{q}, but it is the time-into-the-future-dependent increment of the drift *D*(*t*). Second, the new contribution in a given time step is not amplified exponentially in the next step, but because it does not point into the locally most unstable direction we let it either remain constant or even decay for a few time steps.

Figure 13 demonstrates how well the model error growth can be approximated by *E*_{D}(*t*), particularly by the approximation with an initial linear decay *E*_{D,lin}(*t*). As the numerical value of the Lyapunov exponent *λ*_{D} in the hypothesis *E*_{D}(*t*) we use the one determined as *λ*_{q} in the quadratic hypothesis with model error d*E*_{q} (Eq. 6), which is *λ*_{q,L05-2}=0.27 d^{−1} from the difference between the L05-1 and L05-2 systems (Fig. 8) and *λ*_{q,L05-3}=0.38 d^{−1} from the difference between the L05-1 and L05-3 systems (Fig. 9). The justification for using the parameter *λ*_{q} as *λ*_{D} can be found in Bednář et al. (2021). In short, we argue that *λ*_{q} is a better estimate of the Lyapunov exponent of the model system than, e.g., *λ*_{r}. The parameters *M* and *σ* of the hypotheses *E*_{D}(*t*) are determined empirically. It is, however, a bit puzzling that when using the value *λ*_{q} in *E*_{D,lin}(*t*), we can approximate that part of the model error growth curve *E*_{M}(*t*), which can also be well approximated by the exponential growth with model error d*E*_{r} (Fig. 13), leads to a value of the exponent that is different from *λ*_{q}. Therefore, in the next section, we discuss the relationship between the drift *D* and parameters *λ* and *β* in d*E*_{q} (Eq. 6) and d*E*_{r} (Eq. 5).

## 4.2 Understanding the drift through parameters of the quadratic hypothesis and exponential growth both with model error

We saw two meaningful approximations to the model error growth curve over lead time: the quadratic hypothesis with model error and the exponential growth with model error. The best-fit parameter values of these two approximations are listed in Table 1. This table shows that $\stackrel{\mathrm{\u203e}}{\mathrm{d}D/\mathrm{d}t}\approx {\mathit{\beta}}_{\mathrm{q}}$, but the most evident is the difference of the exponential growth exponent *λ*, where for both realities L05-2 and L05-3 the exponent *λ*_{q} of the quadratic hypothesis is larger than *λ*_{r} of the best-fit exponential error growth.

To understand the meaning of the parameters of the exponential growth with model error d*E*_{r} (Figs. 8 and 9), let us first define the model error growth *E*_{D,0}(*t*) based on the drift *D*, ignoring the initial decrease caused by ** D**(

*t*) not pointing into the locally most unstable direction:

where *λ*_{D} is the largest Lyapunov exponent of the model (L05-1 system). The time derivative (calculated from the difference at successive time steps) of *E*_{D,0}(*t*) is

where $\mathrm{d}D/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right({t}_{k}\left)\right)=\mathrm{d}D/\mathrm{d}t\left({t}_{k}\right)$, meaning that due to the monotonicity of *E*_{D,0}(*t*) in time, we can exchange the dependence on *t* by the dependence on *E*_{D,0}(*t*) (see Fig. 14). Equation (23) now claims that the dashed red curve, which is *λ*_{D}*E*, plus the values of the blue curve taken at corresponding times *t*_{k}, sum up to yield the red curve *E*_{D}(*E*(*t*)), where we approximate the slope of the linear increase of the red curve by the slope of the black curve, which describes the observed total error.

We focus now on those parts of the three curves *E*>0.3 for L02-5 or for *E*>0.5 for L05-3. As has been said above, we observe that in this range of *E*, $\mathrm{d}{E}_{D,\mathrm{0}}/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right)$ (Fig. 14, red curve) has the same growth rate (tendency) as $\mathrm{d}{E}_{M}/\mathrm{d}t\left({E}_{M}\right)$ (Fig. 14, black curve), which is expressed by *λ*_{r}, fitting the *E*_{r}(*t*) behavior of Eq. (5) to the data (Fig. 14, dashed black curve, *λ*_{r,L05-2}=0.17 d^{−1}, *λ*_{r,L05-3}=0.25 d^{−1}) If we compare these parameter values to the fit using the quadratic hypothesis with model error, we see that *λ*_{r} is smaller than *λ*_{q} of d*E*_{q} (*λ*_{q,L05-2}=0.27 d^{−1}, *λ*_{q,L05-3}=0.38 d^{−1}). In our interpretation of the model errors involving the drift, this is due to the decrease in the drift growth rate $\mathrm{d}D/\mathrm{d}t\left(D\right({t}_{k}\left)\right)$ over time (Fig. 14, blue curve). Hence, *β*_{r} of the d*E*_{r} approximation of $\mathrm{d}{E}_{D,\mathrm{0}}/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right)$ is then an extrapolation of the linear decline of $\mathrm{d}D/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right({t}_{k}\left)\right)$ to ${E}_{D,\mathrm{0}}\left({t}_{\mathrm{0}}\right)=\mathrm{0}$. Therefore, if we solve Eq. (23) for $\mathrm{d}D/\mathrm{d}t$ and define $\mathrm{d}D/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right({t}_{k}\left)\right)\approx {\mathit{\beta}}_{D}-{\mathit{\alpha}}_{D}{E}_{D,\mathrm{0}}$, then *β*_{r}=*β*_{D}. Equation (23) can be used to determine the drift decrease rate *α*_{D}: ${\mathit{\lambda}}_{\mathrm{r}}{E}_{D,\mathrm{0}}+{\mathit{\beta}}_{\mathrm{r}}={\mathit{\lambda}}_{\mathrm{q}}{E}_{D,\mathrm{0}}-{\mathit{\alpha}}_{D}{E}_{D,\mathrm{0}}+{\mathit{\beta}}_{D}\to {\mathit{\alpha}}_{D}=({\mathit{\lambda}}_{\mathrm{q}}-{\mathit{\lambda}}_{\mathrm{r}})$. However, *β*_{r}=*β*_{D} is not the same as *β*_{r} in $\mathrm{d}{E}_{M}/\mathrm{d}t$ (Fig. 14, black curve) because it is reduced by the transition term expressed by Eqs. (20) and (21) and *α*_{D} is valid for $\mathrm{d}D/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right)$ and not for $\mathrm{d}D/\mathrm{d}t\left(D\right)$. In general, because $\mathrm{d}D/\mathrm{d}t\left({E}_{D,\mathrm{0}}\right({t}_{k}\left)\right)\approx {\mathit{\beta}}_{D}-{\mathit{\alpha}}_{D}{E}_{D,\mathrm{0}}$ tends to decrease, *λ*_{q}>*λ*_{r}. Since $\mathrm{d}{E}_{M}/\mathrm{d}t\left({E}_{M}\right)$ is almost identical to $\mathrm{d}{E}_{M+\text{ie}}/\mathrm{d}t\left({E}_{M}\right)$ and differs only in the early stage of development, the approximations of d*E*_{q} and d*E*_{r} are only marginally affected (Figs. 8a and 9a). Therefore, information about the drift *D* can be derived from these hypotheses also for the variant with initial and model error.

For the sole initial condition error, we found that *λ*_{q}≤*λ*_{r}, but this describes a setting very different from model error. In the initial condition error, we compare the forecast and reality of a given high-resolution model, which indeed has much larger error growth exponents for short times and small errors due to small-scale degrees of freedom. As soon as we talk about model errors (with or without initial error), we use the low-resolution L05-1 model for forecasts, and hence its parameters are relevant for the propagation of errors.

In summary, we propose a new interpretation of the growth of forecast errors due to model errors: model errors in successive time steps of the forecasts are only weakly correlated. Therefore, modeling them by a constant term in the error growth d*E* is inappropriate. The observed model forecast error growth can be modeled much more accurately if we use the accumulated model errors called drift, interpret the drift increments as additional initial condition errors, and propagate these forward in time. The decrease of the drift growth rate over forecast time can then explain the growth rate of the model errors. Depending on what data are observed, one can either use the drift to predict the forecast errors or use the forecast errors to infer the drift due to model error.

For the ECMWF forecasting system, we cannot perform error growth experiments, but we can check average forecast errors as a function of lead time. We therefore apply the new way of assessing the model error to the error growth *E*_{EFS}(*t*) of the 500 hPa geopotential height values (Bednář, 2023) calculated (Magnusson and Kallen, 2013) as 25 annual averages over the Northern Hemisphere (20–90°) obtained daily from 1 January 1987 to 31 December 2011. Over this period, we determine the decline in the average initial displacements of the model from reality per unit of time using the parameter *β*_{q} of the quadratic hypothesis with model error. Since the parameter *β* is also used to describe an upscale error growth rate from small-scale processes (Zhang et al., 2019), we check whether *λ*_{q}>*λ*_{r}, as defined in Sect. 5.2, for *β*_{q} determined by model error and *λ*_{q}≤*λ*_{r} for *β*_{q} determined from small-scale processes (see Sect. 4).

## 5.1 Methods of calculation

To eliminate the effects of model errors, the initial error growth curve *E*_{EFS,ie}(*t*) is calculated as the differences between two operational forecasts issued with 1 d lag for the same day. Specifically, we evaluate these for 27 different lead times and used the following pairs of lead times in hours: 0–24, 6–30, ..., 96–120, with 6 h shift, and from 108–132, 120–144, ..., 216–240 with 12 h shift. Detailed information about calculating the error growth of the ECMWF forecasting system can be found in Lorenz (1982). The error growth rate (tendency) is $\mathrm{d}{E}_{\text{EFS,ie}}/\mathrm{d}t\approx \left({E}_{\text{EFS,ie}}\right(t+\mathrm{\Delta}t)-{E}_{\text{EFS,ie}}(t\left)\right)/\mathrm{\Delta}t$ with Δ*t*=6 h for the first 17 time steps and Δ*t*=12 h for the rest. It is evident that this data analysis solely used data produced by the very same model. One can understand this 1 d offset between two forecasts in the following way. At day 0, we use some initial condition and propagate it forward in time. At day 1, when a new forecast starts with new initial conditions, these can be interpreted as perturbations to the day 1 forecast started at day 0. Therefore, comparing these two forecasts for the very same day as a function of lead time gives us the initial condition error growth. The only disadvantage of this procedure is that we cannot control the magnitude of the perturbation. What we interpret as perturbation is the deviation of the true forecast at day 1 from the new analysis, which is used to initialize the new forecast.

The initial and model error growth curve ${E}_{\text{EFS},M+\text{ie}}\left(t\right)$ is calculated as differences between operational forecasts and analyses from ERA-Interim for a given day. Forecasts range from 0.5 d ago relative to the given day to 10 d ago, with time step Δ*t*=12 h. The difference between operational analysis and analysis from ERA-Interim is taken as the initial error. The error growth rate (tendency) is $\mathrm{d}{E}_{\text{EFS},M+\text{ie}}/\mathrm{d}t\approx \left({E}_{\text{EFS},M+\text{ie}}\right(t+\mathrm{\Delta}t)-{E}_{\text{EFS},M+\text{ie}}(t\left)\right)/\mathrm{\Delta}t$ with Δ*t*=12 h.

## 5.2 Results and comparisons

From the data, we calculate 25 annual averages of the initial error growth curve *E*_{EFS,ie}(*t*) and 25 annual averages of the initial and model error growth curve ${E}_{\text{EFS},M+\text{ie}}\left(t\right)$ and their growth rates (tendencies) $\mathrm{d}{E}_{\text{EFS,ie}}/\mathrm{d}t$ and $\mathrm{d}{E}_{\text{EFS},M+\text{ie}}/\mathrm{d}t$. We approximate the growth rates by the exponential growth with model error d*E*_{r} and by the quadratic hypothesis with model error d*E*_{q}. Because the data are only up to 10 d and therefore do not cover the entire growth curve and because d*E*_{q} is a three-parameter approximation, we first discuss the error in parameter estimation. Magnusson and Kallen (2013) showed that the error saturation parameter *E*_{lim} estimated from the d*E*_{q} hypothesis underestimates the true limiting value. Bednář et al. (2021) showed that the deviations of the values of *λ*_{q} and *β*_{q} of d*E*_{q} from the true values are anti-correlated, meaning that when one is overestimated, the other is underestimated. The average value of *λ*_{q} over 25 annual averages of *E*_{EFS,ie}(*t*) and ${E}_{\text{EFS},M+\text{ie}}\left(t\right)$ has been determined by Bednář et al. (2021) to be ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{q}}=\mathrm{0.35}$ d^{−1}, and to approximate the data, we fix *λ*_{q} to this value. Therefore, we decrease the oscillation of *β*_{q} and bring *E*_{lim} of d*E*_{q} closer to the values determined by Magnusson and Kallen (2013). The fact that *E*_{lim} is closer to the theoretical limit values estimated by Magnusson and Kallen (2013) justifies this approach.

The resulting values are shown in Fig. 15. We find that *λ*_{r,ie}≥*λ*_{q,ie} and ${\mathit{\lambda}}_{\mathrm{r},M+\text{ie}}<{\mathit{\lambda}}_{\mathrm{q},M+\text{ie}}$, which satisfies the hypothesis presented in Sect. 4.2. Here, ${\mathit{\lambda}}_{\mathrm{r},M+\text{ie}}$ (Fig. 15a, solid red curve) has an approximately constant value, showing approximately the same decrease in *α*_{D} = the drift rate $\mathrm{d}D/\mathrm{d}t\left(D\right)$ over the years because it is shown in Sect. 4.2 that ${\mathit{\alpha}}_{D}=({\mathit{\lambda}}_{\mathrm{q},M+\text{ie}}-{\mathit{\lambda}}_{\mathrm{r},M+\text{ie}})$.

For *λ*_{r,ie} (Fig. 15a, solid blue curve), an increase over time is observed, indicating an increase in the resolution of the ECMWF system with smaller spatiotemporal-scale phenomena with a larger error growth rate. The parameter ${\mathit{\beta}}_{\mathrm{q},M+\text{ie}}$ (Fig. 15b, solid blue curve) of the quadratic hypothesis with model error (Eq. 6) shows an approximately linear decrease over the years. These values indicate that the average displacement value per unit of time $\mathit{\beta}\approx \stackrel{\mathrm{\u203e}}{\mathrm{d}D/\mathrm{d}t}$ decreased by 40 % from 1987 to 2011, from 5.6 to 3.4 m d^{−1}. The parameter *β*_{q,ie} (Fig. 15b, solid red curve) for the variant with an initial error shows an approximately constant value. Together with a constant value of *λ*_{q,ie} (Fig. 15a, black curve), this means that the shape of the error growth rate (tendency) $\mathrm{d}{E}_{\text{ie}}/\mathrm{d}t\left({E}_{\text{ie}}\right)$ changes over the years only by adding a part for smaller *E*_{ie} as the model better describes smaller spatiotemporal-scale phenomena, but this does not change the overall approximation of d*E*_{q}, as can be seen in Fig. 16. These figures also show the similarity to the error growth rates $\mathrm{d}E/\mathrm{d}t$ of the L05-2 and L05-3 systems (Figs. 6–9) and the relevance of fixing *λ*_{q} to approximate the data using d*E*_{q}.

Based on the fact that scale-dependent error growth implies an intrinsic predictability limit, we examined whether omitting atmospheric phenomena, which contribute little to the final value, will improve the predictability of the resulting value. In other words, how does the average forecast error growth change in a model where small-scale phenomena are omitted, but the model error is therefore larger compared to a model where all phenomena are present and the average forecast error growth is scale-dependent. For this, we used the L05 systems defined by Lorenz (2005) and Bednář and Kantz (2022) and the ECMWF systems with data from Magnusson and Kallen (2013) stored in Bednář (2023).

We confirmed that for the multi-scale systems L05-2 and L05-3, the initial error growth *E*_{ie}(*t*) can be described well by the power law d*E*_{p} Eq. (3) or the extended power law d*E*_{w} Eq. (4), respectively, while a simple exponential growth with model error d*E*_{r} (Eq. 5) or the quadratic hypothesis with model error d*E*_{q} (Eq. 6) are less appropriate. However, the non-zero parameter *β* in d*E*_{r} and d*E*_{q} describing the model error also generally relates the multi-scale nature of the system. We showed that in the L05 and ECMWF data (in contrast to the model error scenario) *λ*_{q}≤*λ*_{r}, meaning that the approximation of $\mathrm{d}{E}_{\text{ie}}/\mathrm{d}t\left({E}_{\text{ie}}\right)$ in the early stage grows faster than the approximation of the whole curve due to the presence of only smaller spatiotemporal scales in this part.

For the scenarios of model error growth *E*_{M}(*t*) and both initial and model error growth *E*_{M+ie}(*t*), we showed the appropriateness of the description using exponential growth with model error d*E*_{r} and a quadratic hypothesis with model error d*E*_{q} (Figs. 8 and 5). For $\mathrm{d}{E}_{M}/\mathrm{d}t\left({E}_{M}\right)$, we explained the initial decline and subsequent growth described by d*E*_{r} using the drift *D* (Fig. 11) defined by Orrell et al. (2001), which we extended using a hypothesis that views the drift *D* as a succession of initial errors followed by an exponential time evolution driven by the largest Lyapunov exponent *λ* of the model after a transition period (Fig. 12). We identified that *λ*≈*λ*_{q}, and the validity of *λ*_{q}>*λ*_{r} based on the drift evolution was verified in the L05 and ECMWF systems (Fig. 14 and Eq. 21). For the L05 systems, we have demonstrated that ${\mathit{\beta}}_{\mathrm{q}}\approx \stackrel{\mathrm{\u203e}}{\mathrm{d}D/\mathrm{d}t}$, i.e., that from the *β*_{q} value of the quadratic hypothesis with model error d*E*_{q} the average displacement per unit time (average drift rate) can be determined.

For ECMWF systems (forecast of 500 hPa geopotential height), this means that from 1987 to 2011 we observe a decrease in the model error by approximately 40 % from an average displacement value of 5.6 to 3.4 m d^{−1}. Note that while in 1987 the error in initial displacement (initial error) was approximately 16 m (Magnusson and Kallen, 2013; Bednář et al., 2021), for the variant with initial and model error, a displacement of 5.6 m is produced every day in addition to this value. In 2011, the initial displacement was 6 m; for the variant with initial and model error, an average displacement of 3.4 m is produced daily. Thus, we observe a significant contribution of the model error to the error growth. This is also why the findings from the model error growth scenario can be applied to the initial and model error growth scenario, where the initial and model error growth goes asymptotically to the model error growth.

It is also why omitting atmospheric phenomena, which contribute little to the final value, will not improve the predictability of the resulting value. The average prediction error grows faster in a model where small-scale phenomena are omitted, but the model error is therefore created, compared to a model where all phenomena are present, but the average forecast error growth is scale-dependent (Fig. 10).

We now discuss the possibility that the growth of the displacement produced by the model error may also be scale-dependent. In our case, the model was an L05-1 system, i.e., a system with one scale and exponential error growth (*E*_{D}(*t*) hypothesis). A variant where the L05-2 system was used as the model and the L05-3 system as the reality was also tested. The resulting model error growth is approximately identical to the previous variant (L05-1 system as the model and L05-3 system as the reality), i.e., adding a small scale did not affect the exponential growth of the drift *D* increment. However, it should be noted that the magnitude of the average drift per unit of time ${\mathit{\beta}}_{\mathrm{q}}\approx \stackrel{\mathrm{\u203e}}{\mathrm{d}D/\mathrm{d}t}$ is much greater than the limit (saturation) value of small-scale error magnitude *E*_{2,lim}, and thus we are already in the region of exponential error growth of large-scale variables. For the ECMWF system, it can be seen (Fig. 15) that over the years that the values of a parameter *β*_{q} of the d*E*_{q} approximation of average initial error growth *E*_{ie}(*t*) and initial and model error growth *E*_{ie+M}(*t*) converge, and the growth curves of the two variants are similar for the later analyzed years, as also confirmed by Froude et al. (2013). This means that, in contrast to the presented results of the L05 system, the drift rate of these years is low, and the issue of scale-dependent growth of the drift *D* increment is relevant and should be further investigated.

Another topic for further research is extending the *E*_{D}(*t*) hypothesis (Eq. 21) to describe the model error growth *E*_{M}(*t*) over the entire range up to saturation rather than just the early part where exponential growth is valid. For the part of the time evolution of the model error where the growth slows down and reaches saturation, it can be seen that the drift must reach its limiting value and then no longer contributes to the model error growth and that the exponential growth of the drift increment must slow down and reach its limiting value. However, the specific form needs to be investigated.

The L05-1 system was introduced by Lorenz (2005) as a spatial continuity modification of the Lorenz (1996) system with *N* variables connected by governing equations:

where

$n=\mathrm{1},\mathrm{\dots},N$. *X*_{n} are unspecified (i.e., unrelated to actual physical variables) scalar meteorological quantities (units), *F* is a constant representing external forcing, and *t* is time. The index is cyclic so that ${X}_{n-N}={X}_{n+N}={X}_{n}$ and variables can be viewed as existing around a latitude circle. If *L* is even, $\sum {}^{\prime}$ denotes a modified summation in which the first and last terms are to be divided by 2. If *L* is odd, $\sum {}^{\prime}$ denotes an ordinary summation. Generally, *L* is much smaller than *N* and $J=L/\mathrm{2}$ if *L* is even and $J=(L-\mathrm{1})/\mathrm{2}$ if *L* is odd. To a certain extent, the model quantitatively describes weather systems, but unlike the well-known Lorenz model of atmospheric convection (Lorenz, 1963), it cannot be derived from any atmospheric dynamic equations. The motivation was to formulate the simplest possible set of dissipative chaotically behaving differential equations that share some properties with the “real” atmosphere. Although mechanisms such as potential vorticity generation are lacking in the equations, the model generates five to seven main highs and lows corresponding to planetary waves (Rossby waves). To keep five to seven main highs and lows, Lorenz (2005) suggested a ratio $N/L=\mathrm{30}$ and *F*=15. The choice of parameters *F* and the setting of time unit=5 d is also made to obtain a similar value of the largest Lyapunov exponent to that of the ECMWF forecasting system (Lorenz, 2005).

The L05-2 system (extension to two spatiotemporal scales) was also introduced by Lorenz (2005) as a modification of the two-scale Lorenz (1996) system, where scales were coupled by linear terms that together do not alter the large-scale plus small-scale energy and where small-scale variables were driven entirely by the coupling. Rewriting the equations of the L05-1 system, we would get

where *c* sets the rapidness of small scale compared to large scale and *b* sets the small-scale amplitude size compared to large scale. Equations (A2) and (A3) have an unrealistic property compared to the numerical weather prediction systems. The large-scale and small-scale features are represented by separate sets of variables *X*_{1} and *X*_{2} instead of appearing as superimposed features of a single set *X*_{tot}. Lorenz (2005) wanted to keep the system as simple as possible, so instead of, for example, Fourier analysis, a procedure for expressing variables *X*_{tot,n} as sums of *X*_{1,n} and *X*_{2,n} was introduced:

Parameters *α*, *ω*, and *I* are chosen so that *X*_{1} is a low-pass-filtered version of *X*_{tot}, and *X*_{2} represents the difference between the full signal *X*_{tot} and the filtered signal. By this procedure, *X*_{2} has a much smaller amplitude than *X*_{1}, and also its time evolution should be faster since the temporal derivative is related to the spatial derivative via the difference $({X}_{\mathrm{1},n+\mathrm{1}}-{X}_{\mathrm{1},n-\mathrm{2}})$, which for the low-pass-filtered signal *X*_{1} is typically smaller than for the signal *X*_{2}.

More precisely, Lorenz's (2005) idea is that the parameters *α* and *ω* are chosen so that *X*_{1} equals *X*_{tot} whenever *X*_{tot} changes quadratically over the longitudes (variables) *n*−*I* through *n*+*I*. It is when ${\sum}_{i=-I}^{I}{}^{\prime}(\mathit{\alpha}-\mathit{\omega}|i\left|\right)=\mathrm{1}$ and ${\sum}_{i=-I}^{I}{}^{\prime}{i}^{\mathrm{2}}(\mathit{\alpha}-\mathit{\beta}|i\left|\right)=\mathrm{0}$. By solving these equations, we get

The procedures (Eqs. A4 and A5) are functions of the interval length $[-I,I]$.

When creating a system $\mathrm{d}{X}_{\text{tot}}/\mathrm{d}t$ as the sum of $\mathrm{d}{X}_{\mathrm{1}}/\mathrm{d}t$ and $\mathrm{d}{X}_{\mathrm{2}}/\mathrm{d}t$ (sum of Eqs. A2 and A3), the coupling term *c**X*_{1,n} in Eq. (A3), which enables short waves to develop, is combined with the dissipation term $-{X}_{\mathrm{1},n}$ in Eq. (A2). Therefore, the coupling term can be canceled entirely, or it can appear in *X*_{1} rather than *X*_{2} when *X*_{tot} is analyzed, and there might be nothing to enable the short waves in *X*_{2} to grow. Lorenz (2005) reformulated the coupling process by adding a small fraction of *X*_{1} to *X*_{2} so small waves in *X*_{2} can amplify. It is done by replacing ${b}^{\mathrm{2}}[{X}_{\mathrm{2}},{X}_{\mathrm{2}}{]}_{\mathrm{1},n}+c{X}_{\mathrm{1},n}$ by $[{X}_{\mathrm{2}},{X}_{\mathrm{2}}+{c}^{\prime}{X}_{\mathrm{1}}{]}_{\mathrm{1},n}$ in Eq. (A3), and the L05-2 system would be

where $c={c}^{\prime}\cdot {b}^{\mathrm{2}}$.

Based on the L05-2 system (Eqs. A4–A8), Bednář and Kantz (2022) designed a three-level (scales) system (L05-3):

where *c*_{1}, *c*_{2}, *b*_{1}, and *b*_{2} are parameters, and the procedures for expressing the variables are

where *I*_{1} and *I*_{2} set the length of the intervals $[-I,I]$.

The parameters of L05 systems (L05-1, L05-2, L05-3) should be set so that all scales behave chaotically (the largest Lyapunov exponent of each scale is positive) and that all scales have a significant difference in amplitudes and fluctuation rates. For the L05-1 system (Eq. A1), the chaotic behavior is determined by the value of *F* and the number of variables *N*. For Eqs. (A2) and (A3), where the forcing *F* acts only on a large scale, the chaotic behavior of the small scale is created by coupling. The coupling size is cascaded from a large scale to a small one. Because the values of large-scale variables are determined by the forcing *F*, the *F* value indirectly affects the small-scale chaotic behavior and must be chosen large enough to ensure chaotic behavior through coupling for all scales (levels). This fact must also apply to L05-2 and L05-3 systems, but procedures (A4) and (A5) for the L05-2 system and (A10)–(A12) for the L05-3 system also affect the scales' chaotic behavior, amplitude, and fluctuation rate through the choice of *I* (Lorenz, 2005).

To maintain the required properties *F*=15, *N*=360, *L*=12, and *J*=6 is chosen for the L05-1 system (Fig. 1a). To have the small scale 100 times smaller than the large scale, *F*=15, *N*=360, *L*=12, *J*=6, *b*=10, *c*=1, and *I*=10 are selected for the L05-2 system (Fig. 2a). For the L05-3 system, with requirements for the medium-scale amplitude to be about 10 times smaller than the large-scale amplitude, the small-scale amplitude to be about 10 times smaller than the medium-scale amplitude, and for the scales to have different oscillation rates (Fig. 3a), *F*=15, *N*=360, *L*=12, *J*=6, *b*_{1}=1, *b*_{2}=10, *c*_{1}=1, *c*_{2}=1, *I*_{1}=20, and *I*_{2}=10. The calculation is done using a fourth-order Runge–Kutta method with a time step $\mathrm{\Delta}t=\mathrm{1}/\mathrm{240}$ or 0.5 h.

The ECMWF forecasting system dataset was obtained from the personal repository of Linus Magnusson (Bednář, 2023). The L05-3 system dataset and products from the ECMWF forecasting system dataset, codes, and figures were conducted in Wolfram Mathematica, and they are permanently stored at https://doi.org/10.17605/OSF.IO/2EWXB (Bednář, 2023).

HB proposed the idea, carried out the experiments, and wrote the paper. HK supervised the study and co-authored the paper.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

The authors are grateful to Linus Magnusson for offering the use of the ECMWF forecasting system dataset from his personal repository.

This research has been supported by the Grantová Agentura České Republiky (grant no. 19-16066S).

This paper was edited by Chanh Kieu and reviewed by Quan Wang and three anonymous referees.

Allen, M., Frame, D., Kettleborough, J., and Stainforth, D.: Model error in weather and climate forecasting, in: Predictability of Weather and Climate, edited by: Palmer, T. and Hagedorn, R., Cambridge University Press, Cambridge, UK, 391–427, https://doi.org/10.1017/CBO9780511617652.004, 2006.

Aurell, E., Boffetta, G., Crisanti, A., Paladin, G., and Vulpiani, A.: Growth of noninfinitesimal perturbations in turbulence, Phys. Rev. Lett., 77, 1262, https://doi.org/10.1103/PhysRevLett.77.1262, 1996.

Aurell, E., Boffetta, G., Crisanti, A., Paladin, G., and Vulpiani, A.: Predictability in the large: an extension of the concept of Lyapunov exponent, J. Phys. A-Math. Gen., 30, 1–26, https://doi.org/10.1088/0305-4470/30/1/003, 1997.

Bednář, H.: Analysis of model error in forecast errors of Extended Atmospheric Lorenz' 05 Systems and the ECMWF system, OSF [code and data set], https://doi.org/10.17605/OSF.IO/2EWXB, 2023.

Bednář, H. and Kantz, H.: Prediction error growth in a more realistic atmospheric toy model with three spatiotemporal scales, Geosci. Model Dev., 15, 4147–4161, https://doi.org/10.5194/gmd-15-4147-2022, 2022.

Bednář, H., Raidl, A., and Mikšovský, J.: Initial Error Growth and Predictability of Chaotic Low-dimensional Atmospheric Model, International Journal of Automation and Computing, 11, 256–264, https://doi.org/10.1007/s11633-014-0788-3, 2014.

Bednář, H., Raidl, A., and Mikšovský, J.: Recalculation of error growth models' parameters for the ECMWF forecast system, Geosci. Model Dev., 14, 7377–7389, https://doi.org/10.5194/gmd-14-7377-2021, 2021.

Boffetta, G., Giuliani, P., Paladin, G., and Vulpiani, A.: An Extension of the Lyapunov Analysis for the Predictability Problem, J. Atmos. Sci., 23, 3409–3416, https://doi.org/10.1175/1520-0469(1998)055<3409:AEOTLA>2.0.CO;2, 1998.

Brisch, J. and Kantz, H.: Power law error growth in multi-hierarchical chaotic system-a dynamical mechanism for finite prediction horizon, New J. Phys., 21, 1–7, https://doi.org/10.1088/1367-2630/ab3b4c, 2019.

Budanur, N. R. and Kantz, H.: Scale-dependent error growth in Navier-Stokes simulations, Phys. Rev. E, 106, 1–7, https://doi.org/10.1103/PhysRevE.106.045102, 2022.

Buizza, R.: Horizontal resolution impact on short- and long-range forecast error, Q. J. Roy. Meteor. Soc., 136, 1020–1035, https://doi.org/10.1002/qj.613, 2010.

Cencini, M. and Vulpiani, A.: Finite Size Lyapunov Exponent: Review on Applications, J. Phys. A-Math. Theor., 46, 254019, https://doi.org/10.1088/1751-8113/46/25/254019, 2013.

Dalcher, A. and Kalnay, E.: Error growth and predictability in operational ECMWF analyses, Tellus A, 39, 474–491, https://doi.org/10.1111/j.1600-0870.1987.tb00322.x, 1987.

Ding, R. and Li, J.: Comparisons of two ensemble mean methods in measuring the average error growth and the predictability, Acta Meteorol. Sin., 25, 395–404, https://doi.org/10.1007/s13351-011-0401-4, 2011.

Froude, L. S., Bengtsson, L., and Hodges, K. I.: Atmospheric Predictability Revised, Tellus A, 63, 1–13, https://doi.org/10.3402/tellusa.v65i0.19022, 2013.

Jacobson, M. Z.: GATOR-GCMM: 2. A study of day- and nighttime ozone layers aloft, ozone in national parks, and weather during the SARMAP field campaign, J. Geophys. Res., 106, 5403–5420, https://doi.org/10.1029/2000JD900559, 2001.

Leith, C. E.: Objective methods for weather prediction, Annu. Rev. Fluid Mech., 10, 107–128, https://doi.org/10.1146/annurev.fl.10.010178.000543, 1978.

Li, J., Feng, J., and Ding, R.: Attractor radius and global attractor radius and their application to the quantification of predictability limits, Clim. Dynam., 51, 2359–2374, https://doi.org/10.1007/s00382-017-4017-y, 2018.

Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20, 130–141, https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2, 1963.

Lorenz, E. N.: The predictability of a flow which possesses many scales of motion, Tellus, 21, 289–307, https://doi.org/10.1111/j.2153-3490.1969.tb00444.x, 1969.

Lorenz, E. N.: Atmospheric predictability experiments with a large numerical model, Tellus, 34, 505–513, https://doi.org/10.1111/j.2153-3490.1982.tb01839.x, 1982.

Lorenz, E. N.: Predictability: a problem partly solved, in: Predictability of Weather and Climate, edited by: Palmer, T. and Hagedorn, R., Cambridge University Press, Cambridge, UK, 1–18, https://doi.org/10.1017/CBO9780511617652.004, 1996.

Lorenz, E. N.: Designing chaotic models, J. Atmos. Sci., 62, 1574–1587, https://doi.org/10.1175/JAS3430.1, 2005.

Magnusson, L. and Kallen, E.: Factors Influencing Skill Improvements in the ECMWF Forecasting System, Mon. Weather Rev., 141, 3142–3153, https://doi.org/10.1175/MWR-D-12-00318.1, 2013.

Orrell, D.: Role of the metric in forecast error growth: how chaotic is the weather?, Tellus, 54, 350–362, https://doi.org/10.1034/j.1600-0870.2002.01389.x, 2002.

Orrell, D., Smith, L., Barkmeijer, J., and Palmer, T. N.: Model error in weather forecasting, Nonlin. Processes Geophys., 8, 357–371, https://doi.org/10.5194/npg-8-357-2001, 2001.

Palmer, T. N., Döring, A., and Seregin, G.: The real butterfly effect, Nonlinearity, 27, 9, https://doi.org/10.1088/0951-7715/27/9/R123, 2014.

Phillips, N. A.: Principles of Large Scale Numerical Weather Prediction, in: Dynamic Meteorology, edited by: Morel, P., Springer, Dordrecht, https://doi.org/10.1007/978-94-010-2599-7_1, 1973.

Savijarvi, H.: Error Growth in a Large Numerical Forecast System, Mon. Weather Rev., 123, 212–221, https://doi.org/10.1175/1520-0493(1995)123<0212:EGIALN>2.0.CO;2, 1995.

Simmons, A. J., Mureau, R., and Petroliagis, T.: Error growth and estimates of predictability from the ECMWF forecasting system, Q. J. Roy. Meteor. Soc., 121, 1739–1771, https://doi.org/10.1002/qj.49712152711, 1995.

Toth, Z. and Kalnay, E.: Ensemble forecasting at NMC: The generation of perturbations, B. Am. Meteorol. Soc., 74, 2317–2330, https://doi.org/10.1175/1520-0477(1993)074<2317:EFANTG>2.0.CO;2, 1993.

Toth, Z. and Kalnay, E.: Ensemble forecasting at NCEP and the breeding method, Mon. Weather Rev., 12, 3297–3319, https://doi.org/10.1175/1520-0493(1997)125<3297:EFANAT>2.0.CO;2, 1997.

Zhang, F., Sun, Q., Magnusson, L., Buizza, R., Lin, S. H., Chen J. H., and Emanuel, K.: What is the Predictability Limit of Multilatitude Weather, J. Atmos. Sci., 76, 1077–1091, https://doi.org/10.1175/JAS-D-18-0269.1, 2019.

- Abstract
- Introduction
- Error growth in the L05 systems – types and methods of calculation
- Error growth in the L05 systems – results and comparisons
- Error growth in the L05 systems – discussion and explanation of results
- Error growth in the ECMWF systems
- Conclusion and discussion
- Appendix A: Lorenz's L05 systems
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Error growth in the L05 systems – types and methods of calculation
- Error growth in the L05 systems – results and comparisons
- Error growth in the L05 systems – discussion and explanation of results
- Error growth in the ECMWF systems
- Conclusion and discussion
- Appendix A: Lorenz's L05 systems
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References