Recalculation of error growth models’ parameters for the ECMWF forecast system

This article provides a new estimate of error growth models’ parameters approximating predictability curves and their differentials, calculated from data of the ECMWF forecast system over the 1986 to 2011 period. Estimates of the largest Lyapunov exponent are also provided, along with model error and the limit value of the predictability curve. The proposed correction is based on the ability of the Lorenz’s (2005) system to simulate predictability curve of the ECMWF forecasting 10 system and on comparing the parameters estimated for both these systems, as well as on comparison with the largest Lyapunov exponent (λ = 0.35 day) and limit value of the predictability curve (E∞ = 8.2) of the Lorenz’s system. Parameters are calculated from the Quadratic model with and without model error, as well as by the Logarithmic and General models and by the hyperbolic tangent model. The average value of the largest Lyapunov exponent is estimated to be in the <0.32; 0.41> day range for the ECMWF forecasting system, limit values of the predictability curves are estimated with lower theoretically 15 derived values and new approach of calculation of model error based on comparison of models is presented.


Introduction
Forecast errors in numerical weather prediction systems grow in time because of the inaccuracy of the initial state (initial error), chaotic nature of the weather system itself, and the model imperfections (model error). The growth of forecast error in weather prediction is exponential on average. As an error becomes larger, its growth slows down and then stops with the 20 magnitude saturating at about the average distance between two states chosen randomly from dynamically and statistically possible states (limit (saturated) error). This average growth of forecast error with increasing lead times is called the predictability curve.
Predictability curves (Froude et al., 2013) of the European Centre for Medium-Range Weather Forecasts (ECMWF) numerical weather prediction system are calculated by the approach developed by Lorenz (1982), where two types of error growth can 25 be obtained (Lorenz, 1982). The first type is calculated as the root mean square difference between forecast data of increasing lead times and analysis data valid for the same time. This error growth estimate consists of initial and model error and following Lorenz (1982) we will call it the lower bound predictability curve (L). The second type is calculated as the root-mean-square difference between pairs of forecasts, valid for the same time but with times differing by some fixed time interval (the 2 difference between two forecasts issued with 24-h lag but valid at the same time is used in this article). This type consists of 30 initial error and we will call it the upper bound predictability curve (U). Predictability curves of Lorenz's 05 system (L05; Lorenz, 2005) can be controlled by model parameters and by the size of the initial error and they are set to be as close to predictability curves of ECMWF forecasting system as possible.
Over the years several error growth models approximating predictability curves have been developed, aiming to quantify Lyapunov exponents, model errors (for the imperfect model case where the atmosphere is not perfectly modeled), and limit 35 (saturated) errors. The first, called Quadratic ( Km ), was designed by Lorenz (1969). Dalcher and Kalney (1987) added model error to the Quadratic model and Savijarvi (1995) changed it to the form ( Km  ), that is used today. An alternative, called Logarithmic model ( Lm ) was introduced by Trevisan et al. (1992;1993). General model ( Gm ) was introduced by Stroe and Royer (1993;1994). All these models approximate differences of predictability curves (error growth rate). Newer models approximate the predictability curve directly by the hyperbolic tangent ( Tm and Tm  ) (Žagar et al., 2017). 40 Values of parameters calculated from error growth models are used to evaluate the improvement of the ECMWF forecasting system (Magnusson and Kallen, 2013), to estimate the predictability or the limit error (Bengtsson et al., 2008), to quantify impacts of different model's resolutions (Buizza, 2010), to study chaos and model error in different spatial-temporal scales (Žagar et al., 2015; Žagar et al., 2017 ). They are also used by researchers when the need arises to estimate chaoticity, model error or predictability, but their validity can't be proved, because standard methods (Sprott, 2006) to calculate the largest 45 Lyapunov exponents for the ECMWF forecasting system can't be used due to a large number of variables. An independent value estimated from forecast and analysis anomalies can be calculated for the limit error (Simmons et al., 1995) and its validity will be discussed. The need for correct values of error growth models´ parameters increased these days because the Quadratic model with model error is used to describe multiscale weather (Zhang et al., 2019).
This article intends to provide a new estimate of parameters of error growth models in the ECMWF forecasting system 50 calculated from data over the 1986 to 2011 period. The correction is based on comparing the parameters calculated from the error growth models for the L05 system and the ECMWF forecasting system and on comparison with the largest Lyapunov exponent and the limit value of the predictability curve of the L05 system that can be calculated independently and with sufficient accuracy. To make the correction valid, predictability curves of the ECMWF forecasting system and the L05 systems are compared for two different methods (arithmetic and geometric averages) and the number of variables of the L05 system, 55 pertaining to the best match of the predictability curves is identified. As a result, a new approach to the calculation of model error based on a comparison of models is presented. This article is divided into six seven sections. The second describes the experimental setting. The third describes calculation of the predictability curves. The third fourth provides a comparison of predictability curves of the ECMWF forecasting system and the L05 system and the fourth fifth deals with the estimation of Lyapunov exponents, model, and limit errors of the 60 ECMWF forecasting system based on the correction. Discussion and conclusions are then presented in the final two sections.

Experimental setting
L05 model is based on the low-dimensional atmospheric system presented by Lorenz (1996) (Lorenz, 1963), it cannot be derived from any atmospheric dynamic 70 equations. The motivation was to formulate the simplest possible set of dissipative chaotically behaving differential equations that share some properties with the "real" atmosphere. One of the model´s properties is to have 5 to 7 main highs and lows that correspond to planetary waves (Rossby waves) and a number of smaller waves that correspond to synoptic-scale waves. For Eq. (1) this is only valid for 30 N = and that is, as it will be seen, not sufficient for the experimental setting. Therefore, spatial continuity modification of L05 system is used, where the Eq. (1) is rewritten to the form: 75 If L is even, ∑' denotes a modified summation, in which the first and last terms are to be divided by 2. If L is odd, ∑' denotes an ordinary summation. Generally, L is much smaller than N and J = L/2 if K is even and J = (L-1)/2 if L is odd. For comparison 80 with predictability curves of the ECMWF forecasting system, we choose N =30; 60; 90; 120; 150; 360. To keep a desirable number of main pressure highs and lows, Lorenz (2005) suggested to keep ratio N/L = 30 and therefore L = 1; 2; 3; 4; 5; 12.
For even values of L we have J = 1; 2; 6 and for odd values of L we have J = 0; 1; 2. Parameter F = 15 is selected as a compromise between too high Lyapunov exponent (smaller F) and undesirable shorter waves (larger F). For this setting and by the method of numerical calculation (Sprott, 2006), the global largest Lyapunov exponents are calculated (Table 2). By the 85 definition of Lorenz (1969): ""A bounded dynamical system with a positive Lyapunov exponent is chaotic ". Because the 4 value of the largest Lyapunov exponent is positive and the system under study is bounded, it is chaotic. Strictly speaking (Aligood et al., 1996), we also need to exclude the asymptotically periodic behavior, but such a task is impossible to fulfill for the numerical simulation. The choice of parameters F and time unit = 5 days is made to obtain a similar value of the largest 90 Lyapunov exponent as the ECMWF forecasting system.

Calculation of predictability curves
To calculate predictability curves (Lorenz, 1996), arbitrary values of the variables n X are chosen, and, using a fourth-order Runge-Kutta method with a time step ∆t = 0.05 or 6 hours, it is integrated forward for 14400 steps or 10 years. , Formulas to calculate prediction errors by geometric means (G) are: For an overview of the symbols see Table 1.
To calculate predictability curves for the ECMWF forecasting system (EFS) values of 500 hPa geopotential height are used.
Comparisons of model predictability curves are done through values normalized by the limit (saturated) errors ( . Because maximum forecast time for the ECMWF forecasting system is 10 days, presented predictability curves don't reach their limit value. Independent measure of limit error can be calculated as: is the time-averaged anomaly with respect to climate and ( ) ac − is the time-averaged analysis anomaly with 135 respect to climate. The climate is defined from ERA-Interim daily climatology.
,U E  and ,L E  differ if the ECMWF forecasting system does not sufficiently describe the variability of the atmosphere (model error). More information can be found in (Simmons et al., 1995). Because it will be shown that values of limit error calculated by this method aren´t correct, predictability curves of the ECMWF forecasting system are normalized by values calculated by Eq. (15).

Comparison of predictability curves 140
Predictability curves of the ECMWF and L05 systems are compared to find a setting of the L05 system (number of variables (N), the size of the initial errors, preference of arithmetic or geometric mean) that gives the most similar progress of systems' predictability curves.
Predictability curves of the L05 system show negative growth for the first time step (6 hours) but turn into an increase thereafter. At the second time step (12 hours) values of predictability curves reach approximately the same values as it had 145 initially. A possible explanation could be that initial errors set the initial state off the attractor and decrease occurs because the first tendency is to get on the attractor (Brisch and Kantz, 2019). With an increase of average errors, chaotic behavior becomes dominant. Predictability curves of the ECMWF forecasting system do not exhibit this type of behavior. This may be because of larger time steps or methods of objective analysis. We aim to get the most similar predictability curves of both models and therefore the first two time steps (up to 12 hours) of L05 model's predictability curves are filtered out. 150 Description of symbols that indicate the type of prediction error E in the text is provided in Table 1  Predictability curves calculated by arithmetic and geometric mean show a significant difference for the L05 system (all N), which agrees with Ruiqiang and Jianping (2011), and minor difference for the ECMWF forecasting system. For the L05 system and upper and lower bound predictability curves, the maximal difference is between 6.5 %  Table 2 (for a description of the symbols see Table 1). For the L05 system predictability curves are calculated with N = 60; 175 90; 120; 150 variables and by arithmetic and geometric mean. For the ECMWF forecasting system only arithmetic mean is used.
A comparison of lower bound predictability curves (Fig. 2) shows the most similar predictability curves of the ECMWF forecasting system and the L05 system for the L05 system calculated by arithmetic mean with N = 90. For upper bound predictability curves (Fig. 1), predictability curves for the L05 system with N = 90 are the most similar but to the year 1999 180 for predictability curves of the L05 system calculated by geometric mean and after 1999 by the arithmetic mean.

Estimation of parameters
Parameters of error growth models are the Lyapunov exponent, model error, and limit error. They are estimated from approximations of predictability curves or differences of predictability curves where t is time and 0.25 t = day (Figs. 3 and 4). Error growth models considered here are: 185  Figs. 3 and 4) and p, A, B, a, b are parameters.

195
The calculation is done for the ECMWF forecasting system and the L05 system ( 90 N = ), for arithmetic (A) and geometric (G) means, for upper bound predictability curves (U) and lower bound predictability curves (L). See Tables 3 and 4 for RMS values of parameters  , lim E ,  and p , that are calculated over all used initial errors for the L05 system and all calculated years for the ECMWF forecasting system.
The average values of parameters  , lim E are higher for the lower bound predictability curves than for the upper bound 200 predictability curves. Upper bound predictability curves should not include model error (theoretically 0  = ) but from Table   4 it can be seen that for the L05 system (arithmetic mean) the values are even higher than for the lower bound predictability curves. For the ECMWF forecasting system the values of  are higher for lower bound predictability curves which is theoretically more acceptable, but  is not zero for the upper bound predictability curves. A possible explanation can be the sensitivity to correct approximation (cases with higher  have lower  ), but this can not fully explain the discrepancy. For 205 p the values of upper and lower bound predictability curves are similar to each other (L05 system and ECMWF forecasting system).
There are significant differences of parameters  , lim E ,  and p between predictability curves calculated by arithmetic and geometric mean for the L05 system (for the ECMWF forecasting system only arithmetic mean is presented). The most 9 significant differences are detected for  and p , where for  values are closer to zero for geometric mean and values of 210 predictability curves calculated by arithmetic mean are two or three times higher. Values of parameter p are closer to 1 p = for geometric mean. This means that differences of predictability curves calculated by geometric mean have a shape that is close to a symmetric parabola (for example Fig. 3a) but for the arithmetic mean the parabolic shape is skewed to the left (for example Fig. 3c).
Note that the description of symbols that indicate the type of parameters of error growth models  ,  , p and lim E in the text 215 is provided in Table (Fig. 5a). For lower bound predictability curves (the L05 system with N = 90 calculated by arithmetic mean), the average value EFS U  over all error growth models is in the range 0.32; 0.41 day -1 (Fig. 5b). RMSEs of EFS L  are in the range 0.01; 0.02 day -1 .
For comparison, RMSEs of EFS L  are in the range 0.03; 0.07 day -1 (Fig. 5b). The average value EFS  over upper and lower 225 bound predictability curves is shown in Fig. 6 1987, 1988, 1995, 1997, 2003, and 2011 (15) are found for the L05 system with N = 90 (for lower bound predictability curves calculated by arithmetic mean and for upper bound predictability curves calculated by geometric 245 mean to 1999 and after by arithmetic mean). The most similar predictability curves of the L05 system and the ECMWF forecasting system with EFS E  calculated by Eq. (7) are found for the L05 system with N = 90 by the arithmetic mean for upper and lower bound predictability curves. It means that if the comparison is valid and model error is constant for the L05 system (same number of variables over years in the L05 system means constant model error over years), it must be constant also for the ECMWF forecasting system, but the calculation of parameters EFS L  shows a decreasing trend with increasing time (Fig.  250   8b). This can't help yet. But parameters EFS U  have non zero values (Fig. 8a) that are close to EFS L  for some years and that is inconsistent with the theoretical expectation that upper bound predictability curves should be without model error and therefore  should be 0 m/day. This inconsistency can be solved by the new definition of the model error. From Fig. 6  Gm has parameter p that defines skewness of the originally 255 parabolic shape of the difference of predictability curves. 1 p = pertains to symmetrical parabolic shape ( Gm becomes Km ) and 0 p = means the greatest skewness to the left ( Gm becomes Lm ). Parameters  also skew the originally parabolic shape ( Figs. 3 and 4). The model error can be seen as a difference between skewness of upper and lower bound predictability curves and the new definition of model error would be: Results (Fig. 9a)  There is also good agreement with trends of LU pp − (Fig. 9b). Because constant values of LU  − for predictability curves with EFS E  calculated by Eq. (7) are not theoretically possible, predictability curves with EFS E  calculated by Eq. (15) are 265 favored. The reason for the decreasing trend of 05 L LU  − , found for predictability curves of the L05 system with N = 90 that are the most similar with predictability curves of the ECMWF forecasting system normalized by EFS E  calculated by Eq. (15), is that they are partly calculated by geometric and partly by the arithmetic mean.
These arguments are taken as proof of the validity of EFS  , EFS E  calculated by Eq. (15). The reason for the overestimation of EFS E  calculated by Eq. (7) (Fig. 7) can be found in the multiscale behavior of weather. If some events are predictable on a 270 timescale longer than ten days (for example long-lived anomalies in sea surface temperature or soil moisture) than they wouldn´t be captured by medium-range weather forecast (Simmons et al., 1995;Brisch and Kantz, 2019). It is also possible that the overestimation is due to the different source of data used for calculation of EFS E  by Eqs. (7) and (15): For EFS E  calculated by Eq. (7) only data from ERA-Interim (Janoušek 2011) are used but for EFS E  calculated by Eq. (15) data from operational forecast are employed. 275 At the end of this section, it is important to remind the readers about the importance of the correct values of the parameters.

Nowadays,
Km  is used in the ECMWF forecasting system to estimate the influence of different spatiotemporal scales where parameter  newly represents the intrinsic upscale error growth and propagation from small scales and  represents synopticscale error growth (Zhang et al., 2019). The results of our analysis well support this approach by the new definition of model error (Eq. (16)) and by showing the errors of approximations for individual error growth models. 280

Conclusion
The values of error growth models' (Eqs. (8) -(13)) parameters that approximate predictability curves and their differences (Figs. 3 and 4) in the ECMWF forecast system (Tables 3 and 4) were recalculated. It is based on similarities of normalized upper and lower bound predictability curves (Figs. 1 and 2) of the ECMWF forecasting system (annual arithmetic mean of geopotential heights of 500 hPa from years 1986 -2011) and the L05 system (N = 90, arithmetic mean for lower bound 285 predictability curves; geometric mean up to 1999 and arithmetic mean after 1999 for upper bound predictability curves). It is also based on knowledge of the largest Lyapunov exponent (λ = 0.35 day -1 ) and the limit value of the predictability curve (E∞ = 8.2) of the L05 system.
Lyapunov exponents of the ECMWF forecasting system were recalculated by Eq. (14). The average value over all error growth models for upper bound predictability is in the range 0.33; 0.41 day -1 (Fig. 5a) and RMSEs are mostly about 0.01 day -1 . For 290 lower bound predictability curves average value over all error growth models is in the range 0.32; 0.41 day -1 (Fig. 5b).
RMSEs are in the range 0.01; 0.02 day -1 . The average value over upper and lower bound predictability curves is shown in with increasing years for predictability curves with limit values calculated by Eq. (15), and almost constant trend with increasing time (slight decrease can be due to the error of approximations) for predictability curves with limit values calculated by Eq. (7), which is theoretically impossible (Fig. 9a). This new model error calculated as a difference of model error parameters between the upper (Fig. 8a) and lower (Fig. 8b) bound predictability curves well support model error parameters calculated for upper bound predictability curves that are used to represents the intrinsic upscale error growth and propagation 305 from small scales (Zhang et al., 2019).

Code and data availability
The ECMWF forecasting system dataset was obtained from the personal repository of Linus Magnusson (Magnusson, 2013).
L05 system dataset, products from the ECMWF forecasting system dataset, codes, and figures were conducted in Wolfram Mathematica and they are permanently stored at http://www.doi.org/10.17605/OSF.IO/CEK32. 310

Author contributions
HB proposed the idea, carried out the experiments, and wrote the paper. AR and JM supervised the study and co-authored the paper.     ) and over all years for the ECMWF

380
forecasting system of parameters  and p (for description see Table 1).