In his pioneering work (Lorenz 1969, hereafter L69), Lorenz first showed that a flow with many length scales, like the atmosphere, might have an intrinsic finite range of predictability.
Although Lorenz studied the simple 2D vorticity turbulence model in his paper, the conclusion of his study is profound and intriguing. Follow-up studies using more sophisticated models (e.g., Leith and Kraichnan 1972;Daley 1981;Foude et al. 2013;Sun and Zhang 2016;Judt 2018;Zhang et al. 2019) further supported Lorenz' results and the concept of "butterfly effect" has been widely accepted since then. Butterfly effect depicts that even the smallest unresolved errors by numerical models will propagate upscale and ruin our practical weather prediction at the synoptic-scale after a finite length of time (Palmer et al. 2014). Inspired by L69, estimations of this finite range of predictability has since been done extensively (e.g., Smagorinsky 1969;Lorenz 1982;Foude et al. 2013). For the synoptic weather system in mid-latitudes, more recent studies agree with Lorenz that this finite number should be around two weeks (L69; Reeves 2014; Zhang et al. 2019;Judt 2020). With this intrinsic predictability limit, current operational forecasts still have quite some room for improvement. In general, our operational weather forecast is skillful for less than 10 days in the mid-latitudes despite decades of "quiet revolution" (Bauer et al. 2015;Alley et al. 2019). To push our numerical weather prediction (NWP) skill closer to its intrinsic limit, we must understand further the error growth dynamics that limit NWP.
Over the years, conceptually and numerically simple turbulence frameworks, as used in L69, have contributed a lot to our understanding. In a turbulent fluid, the inverse cascade rate of the errors from small to large scales, which is the essence of the "butterfly effect," is noted to be intimately connected with the eddy turnover timescales that are determined by the slope of the Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
background energy spectrum of the fluid. For a flow with energy spectra of power-law behavior (k -p ), previous studies find that if the slope p < 3, the eddy doubling time decreases with scale and the upscale spreading of initially small-scale error provides an intrinsical limit to the predictability of such flows; if p ≥ 3, it is concluded that there is no such a limit (L69; Rotunno andSnyder 2008, hereafter RS2008).
Most of these studies mentioned above generally assume one single slope for the atmosphere. However, our real world is more complicated. Instead of one constant 𝑝 , observational studies (e.g., Nastrom and Gage 1985) indicate that the energy spectra in the atmosphere show a distinct transition from a slope of around -3 at synoptic scales (~1000s km) to a shallower -5/3 slope at mesoscales (~100s km) in the mid-latitudes. Numerous realistic simulations, using both regional (Skamarock 2004;Waite and Snyder 2013;Sun and Zhang 2016) and global high-resolution model (Skamarock et al. 2014), also successfully reproduce the transition of the slope, consistent with the observational estimates. The mechanism(s) that determine the slopes of the kinetic energy spectra are still under debate (Charney 1971;Tulloch and Smith 2006;Callies et al. 2014). Nevertheless, according to L69 and RS2008, we would expect an intrinsic predictability limit for our atmosphere due to this shallower slope at the small-scale end of the kinetic energy spectra.
Based on the observed kinetic energy spectra, we here propose a novel and simple theoretical framework for understanding error growth from minute perturbations in the real atmosphere. This framework features a "two-stage" error growth process, which connects to the two different slopes of the observed kinetic energy spectra. Figure 1 shows a conceptual schematic for the canonical atmospheric kinetic energy spectrum and the proposed error growth behavior linked to this spectrum. An initially minute error will, in the first stage, grow much faster at small Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
scales due to decreasing eddy turnover time within the -5/3 slope range. Within an inherently finite time, these small-scale errors within the -5/3 slope wavelength range will start to saturate while projecting to larger and eventually synoptic scales. In the second stage, the errors at synoptic scales (corresponding to wavelength range within the -3 slope) will grow quasi-exponentially until saturation due to near-constant eddy turnover time in this wavelength range.
With this conceptual picture, the next step is to quantify the growth of the errors under a simple analytical framework. As a tool to help our understanding of complex and chaotic nonlinear interaction, simple analytic equations have been used along with the earlier numeric studies on error growth dynamics. Lorenz (1982) showed that the growth of error variance E could be reasonably well parameterized by a simple exponential growth equation. Dalcher and Kalnay (1987) proposed a modified version based on Lorenz (1982) to describe the evolution of the error variance
by introducing an external error source S. This equation is adopted and widely used in studies of forecast uncertainty of operational weather prediction (e.g., Magnusson and Kallen 2013;Herrera et al. 2016;Žagar et al. 2017). However, very limited analytical work focused on the intrinsic predictability limit of weather systems where the external error source is eliminated. Selz and Craig (2015) fitted the errors in their "identical twin experiments" to an analytical equation they constructed. The reasonable agreement in their study between the full-physics model and simple analytical equations implies that we may also use analytical equations to investigate the intrinsic predictability limit. More recently, Zhang et al. (2019) found that Eq. ( 1) well captured the evolution of the intrinsic error dynamics in the full-physics model. However, both studies mentioned here did not provide detailed explanations behind this consistency between the results of complex full-physics atmospheric models and simple analytical equations considered. This paper serves as an extension of RS2008 and Zhang et al. (2019) and aims to provide a framework that helps us further understand the connection of error growth behavior, the background kinetic energy spectrum of the real atmosphere, and the detailed analytical equation proposed. In section 2, we first revisit L69's earlier model on error growth for different kinetic energy spectrum slopes.
Based on the results of the L69 model, we then propose our hybrid framework for the real atmosphere with hybrid kinetic energy spectra in section 3. A simple analytical equation is also derived in section 3 to further our understanding of the atmospheric predictability limit in the real atmosphere. A brief discussion is given in section 4.
The original model of L69 was devised to study the error growth and predictability of an atmospheric-like fluid system with homogeneous isotropic turbulence using a two-dimensional vorticity (2DV) equation. In this model, Lorenz assumed power-law behavior (k -p ) for the basicstate kinetic energy with specific considerations dedicated to the scenarios with p = 5/3, 7/3, and 3, respectively. While these calculations are robust, it is found that the downscale energy spectral slope of a large-scale forcing for the 2DV equation is -3 (Kraichnan 1967). Given that synopticscale forcing is the main driver for weather systems in the mid-latitudes, the physically consistent choice for the L69 model, therefore, is p=3, which raises concerns about his results for other scenarios. The model in L69 is elegantly generalized in RS2008 to include a surface quasigeostrophic (SQG) equation, which is known to have a -5/3 energy spectrum analogous to 3dimensional turbulence. Our study will adopt this generalized model in RS2008 and further illustrate different error growth scenarios in 2DV (-3 slope) and SQG (-5/3 slope).
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
The evolutions of the errors for the 2DV and SQG systems are detailed in L69 and RS2008.
We here briefly summarize their equations as a set of second-order initial-value problems:
where 𝐶 𝐾,𝐿 is a constant coefficient matrix derived to reflect the interactions between different length scales (K and L represent different spectral bands in the wavenumber space), N is the total number of spectral bands considered in the model 1 Despite similar forms in Eq. ( 2), very different error evolutions are found between the 2DV (-3 slope) case versus the SQG (-5/3 slope) case, which are rooted in their striking differences in 𝐶 𝐾,𝐿 (see Table 1 andTable 3 in RS2008) and the corresponding basic-state spectra. Through directly comparing their results with L69, RS2008 concluded that the basic-state energy spectrum was the determining factor in the error-energy evolution. They showed that a -5/3 spectrum would lead to limited predictability under varying dynamical models, while a -3 spectrum may have unlimited predictability when the initial perturbation becomes infinitesimally small.
Figure 2 depicts the error growth in different experiments using the 2DV equation under different initial condition errors. For each experiment, the initial error distribution is set so that the error field is limited to the small scales only. No initial error is added to the spectral bands that have larger length scales than the cutoff spectral band K (cutoff K in Fig. 2a)foot_1 . For length scales equal to or smaller than spectral band K, their initial error amplitudes are set to their saturation values. Increased 𝐾 means that the initial error is pushed to smaller scales, and thus its amplitude is exponentially reduced. We can find that, as the cutoff 𝐾 increases (initial error reduces exponentially), the time needed for the error to saturate at large scales increases linearly (Fig. 2b).
Therefore, if we could keep reducing the initial error to smaller and smaller scales, we could keep increasing the error saturation time at large scales without any limitation.
This linearity in Fig. 2b also implies that a similar amount of additional predictable time can be gained each time we increase K and therefore limit the initial condition errors to a smaller scale. In other words, error growth at different length scales can be characterized by a single growth rate in the 2DV case. Indeed, this uniform error growth rate agrees well with the turbulence assumption for a flow with a -3 spectrum. More specifically, if A is a measure of the amplitude of the total error energy, then the evolution of 𝑨( 𝑡) could be written as
assuming 𝛼 is the error growth rate. The error doubling time 𝜏 𝐷 can be then calculated to be 𝜏 𝐷 =
where 𝐷 is a constant on the order of unity. Therefore, we have,
which means that the error growth rate is constant for the 2DV case since 𝐸(𝑘) ∝ 𝑘 -3 . Given this constant 𝛼, the evolution of the total error energy in the 2DV case 𝑍 𝑡𝑜𝑡𝑎𝑙 2𝐷𝑉 can be simplified as
To include the error saturation effect at later times, we could also add an additional term as in Durran and Gingrich (2014) to force the time tendency of 𝑍 𝑡𝑜𝑡𝑎𝑙 2𝐷𝑉 to decrease smoothly to 0 as 𝑍 𝑡𝑜𝑡𝑎𝑙 2𝐷𝑉 approaches its saturation threshold 𝑍 𝑠𝑎𝑡 2𝐷𝑉 . With this adjustment, Eq. ( 5) becomes
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
While this additional term is ad hoc, Eq. 6) captures the error growth behavior reasonably well (Fig. 2c).
For the SQG scenario, with a -5/3 slope, the eddy turnover time in Eq. ( 3) will decrease exponentially with decreasing length scales. Hence, the gain in extra forecast lead time through further limiting the initial error to smaller scales will also become exponentially smaller. More specifically, by increasing cutoff K in the experiments (e.g., for experiments of K=12 and K=13
in Fig. 3a), the additional time we gain is simply the time it takes for the initial errors to propagate back and saturate larger scales (e.g., upscale growth from K=13 to K=12), which is on the order of the eddy turnover time at that scale (𝜏 𝐾=12 ). Given exponentially decreasing eddy turnover time under a -5/3 slope, Fig. 3b shows that the error saturation time at large scales can be extended at most by a few turnover cycles of the current smallest resolved scale, and it will eventually approach a near-constant value when the initial condition error approaches zero.
This limited predictability for the SQG (-5/3 slope) case could also be explained according to the turbulence energy cascade theory. After a finite time (on the order of eddy turnover time of the large-scale end if estimated using the turbulence assumption, more on this in Appendix B), errors will saturate no matter how small the initial error is. What we care about the most here is the characteristic finite timescale needed for the errors to saturate. For simplicity, we could write where 𝑍 𝑡𝑜𝑡𝑎𝑙 𝑆𝑄𝐺 is the total error for the SQG scenario. Assuming the saturation value of the total error for the SQG case is 𝑍 𝑠𝑎𝑡 𝑆𝑄𝐺 , then the time needed for the error to saturate, according to Eq. ( 7), is simply 𝑍 𝑠𝑎𝑡 𝑆𝑄𝐺 /𝛾. 𝛾 is the linear error growth rate that may vary with different base-state kinetic energy spectra and different initial condition errors. We acknowledge this linear error growth is not very realistic or physical. Yet, it is very simple and provides an estimation for the error saturation time if we know the value of 𝛾. Similar to Eq. ( 6) , we need to add a saturation term
to represent the saturation effect when 𝑍 𝑡𝑜𝑡𝑎𝑙 𝑆𝑄𝐺 approaches its saturation value 𝑍 𝑠𝑎𝑡 𝑆𝑄𝐺 .
Eq. ( 7) then becomes
Figure 3c further verifies that Eq. ( 8), which simply provides an estimation for the error saturation time, might not be a bad approximation for the original numerical solution of SQG-like error dynamics in Eq. ( 2).
Compared to Eqs. (2), Eq. ( 6) and Eq. ( 8) are more simplified with known analytical solutions that are much easier to understand. Moreover, we can estimate the parameters in both analytical error growth models from their respective basic-state spectrum. For example, 𝛼 can be estimated from Eq. ( 4), whereas 𝛾 is related to the eddy turnover time at the large-scale end of the -5/3 spectrum (more details will be discussed later). Next, we will combine and extend these simple analytical formulas to further explain the complex multiscale predictability of the real atmosphere.
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
In the real atmosphere, different from either 2DV or SQG, the observed kinetic energy spectra in the upper troposphere in the mid-latitudes usually consist of a -3 spectrum at the synoptic scales and a -5/3 spectrum at meso-and smaller scales (Nastrom and Gage 1985). Therefore, the error growth representative of the observed atmospheric energy spectra would have simultaneous contributions from both the 2DV-like spectrum at synoptic scales and the SQG-like spectrum at smaller scales. The ensemble means of the total error for any spectral bank 𝐾 can then be written as
which is a combination of the two ODEs in Eq. ( 2) and could be solved numerically as before. We should note here that the nonlinear saturation adjustment, as in Durran and Gingrich (2014), is also added to Eq. 9) when solving this equation. More details on this can be found in Appendix A. Due to this additional nonlinear saturation effect, the hybrid model of Eq. ( 9) cannot be linearly decoupled as the summation of a solution to the SQG-like system and a solution to the 2DV-like system.
Figure 4 shows an example of the error evolution solved from Eq. ( 9), with the saturation terms included. To solve this hybrid model, we first construct a hybrid basic-state energy spectrum similar to the observed spectrumfoot_2 and the 𝐶 𝐾,𝐿 2𝐷𝑉 and 𝐶 𝐾,𝐿 𝑆𝑄𝐺 are then computed based on the respective -3 and -5/3 parts of the kinetic energy spectrum (see Appendix A for more details).
Consistent with our schematic shown in Fig. 1, we can find that the errors first grow at small scales that are dominated by the -5/3 slope. These errors at the small scales then start to saturate at increasingly larger scales, and the total error growth will come predominantly from the -3 slope part of the kinetic energy spectrum after the smaller-scale errors saturate. Moreover, given that the evolution of small-scale errors is dominated by the SQG-like spectra, further reducing initial errors to infinitesimal scales does not help extend the predictability limit.
To derive an analytical expression for error growth under the observed hybrid spectra, we first separate the total errors ℰ 𝑡 into two parts according to different length scales. The total errors
represents the meso-small scales errors in the -5/3 slope regime, ℰ -3
represents the synoptic-scale errors in the -3 slope regime, Given the decreasing eddy turnover time within the -5/3 slope regime, the meso-small scales errors ℰ -5 3 feature SQG-like upscale growth. At the same time, this upscale growth process would transfer a small portion of these smaller-scale errors into the synoptic scales due to crossscale nonlinear interaction. While the physical mechanisms of the upscale error propagation in the real atmosphere are still under investigation (Zhang et al. 2007;Bierdel et al. 2018), this effect is included in 𝐶 𝐾,𝐿 terms in the numerical solution. In light of Eq. ( 7), the evolution of small-scale errors ℰ -5 3 could then be simplified as
where 𝛾 ′ represents the SQG-like upscale error growth as in Eq. ( 7) and Γ (ℰ -5 3 , ℰ -3 ) here represents the energy that is transited to the synoptic scales through interactions between the -3 slope and the -5/3 slope. An additional nonlinear saturation treatment as in Eq. ( 8) will be introduced later. Given that the -3 slope regime has much weaker cross-scale interaction compared to the -5/3 slope regime, it is reasonable that the small-scale errors are dominated by the SQG-like
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
error growth associated with the -5/3 slope only, which implies that the 𝛾 ′ term shall dominate the Γ (ℰ -5 3 , ℰ -3 ) term in Eq. ( 10). The numerical solution in Fig. 4 also suggests that the initial growth of the errors under a hybrid spectrum is mostly on the small-scale end. Therefore, assuming ℰ -5 3 𝑠𝑎𝑡 is the saturation error for the -5/3 slope regime, we can neglect the Γ (ℰ -5 3 , ℰ -3 ) term here and use
⁄ as a simple estimate of the characteristic time needed for the small scale errors ℰ -5 3 to saturate.
On the other hand, the errors at the synoptic scales will have both the 2DV-like exponential growth and the contributions from smaller scales. In light of Eq. ( 5) and the subtraction of the Γ (ℰ -5 3 , ℰ -3 ) term in Eq. ( 10), the evolution of ℰ -3 could also be approximately written as:
where 𝛼 is the corresponding error growth rate for the synoptic scales.
Combining Eq. ( 10) and Eq. ( 11), we can write the evolution for the total errors ℰ 𝑡 ,
Again, we can add (1 -ℰ 𝑡 ℰ 𝑡 𝑠𝑎𝑡 ) term to describe the saturation of ℰ 𝑡 , similar to Eq. ( 6) and Eq. ( 8).
The equation then becomes
Note that this equation is very similar to Eq. ( 1) that is used in Zhang et al. (2019) between Eq. ( 13) and Eq. ( 1) is that ℰ -3 in Eq. ( 13) is replaced with total error variance ℰ 𝑡 ( E in Eq. ( 1) ), which allows us to provide an analytical solution to the total error ℰ 𝑡 . Moreover, this change is a valid approximation of Eq. ( 13). When ℰ 𝑡 is small, the growth of the errors is dominated by SQG-like upscale process (the 𝛾 ′ term in Eq. ( 13) is much larger than 𝛼ℰ 𝑡 or 𝛼ℰ -3 ).
Changing ℰ -3 to ℰ 𝑡 only has minor impacts on the results. When ℰ 𝑡 becomes larger, then the 2DV-like growth dominates, ℰ -3 approximates to the value of ℰ 𝑡 due to the relatively small saturation value of ℰ -5 3 𝑠𝑎𝑡 . Therefore, we could approximately replace ℰ -3 with ℰ 𝑡 and define 𝜀 = ℰ 𝑡 ℰ 𝑡 𝑠𝑎𝑡 ⁄ , then Eq. ( 13) becomes
where ε(t) is the normalized error. ε=1 means error reaches a maximum or becomes saturated. 𝛼 is the error growth rate and 𝛽 = 𝛾 ′ ℰ 𝑡 𝑠𝑎𝑡 ⁄ . Figure 4b shows the evolution of normalized total errors derived by numerically solving saturation adjusted Eq. ( 9) versus the fitted curve using the analytical solution derived from Eq. ( 14). Both solutions agree with each other well. As mentioned earlier, Eq. ( 14) was also proposed in earlier studies and shown to be useful. Through simple derivation and approximation, our contribution here focuses on directly linking 𝛽 with the intrinsic upscale error growth (associated with the shallower -5/3 spectrum) under a nearly perfect model and nearly perfect initial condition scenario.
Indeed, similar to 𝛼, the parameter 𝛽 could also be estimated directly from the kinetic energy spectrum 𝐸(𝑘) under our framework. We use 𝛽 to represent upscale error growth processes from small convective scales to mesoscales within the -5/3 slope range (stage 1 of Fig. 1). After some time 𝑡 𝑘𝑡 , the mesoscale error will start to saturate, and the large-scale quasi-exponential error growth starts to dominate. In Eq. ( 14), the transition happens when 𝛼𝜀(𝑡 𝑘𝑡 ) = 𝛽, which implies Combining Eqs. ( 15) and ( 17), we have,
Recall Eqs. ( 7) and ( 10), the characteristic time needed for the smaller-scale errors ℰ -5 3 to saturate could also be estimated to be
Combining Eq. ( 18) and Eq. ( 19), we get wavenumber at the smallest scale resolved), 𝐸(𝑘) lies in the -5/3 regime. The kinetic energy in the -5/3 slope regime can then be written as ∫ 𝐸(𝑘) 𝑘 𝑠 𝑘 𝑡 𝑑𝑘. Similarly, we can also write the total kinetic energy as the sum of the kinetic energy in the -3 slope regime and the kinetic energy in the -5/3 slope regime, hence
where 𝑘 𝑙 is the wavenumber at the largest scale that a -3 slope might hold. Therefore, we have
Another easy way to estimate the value of 𝛽 is by utilizing the schematic shown in Fig. 1. At transition time 𝑡 𝑘𝑡 , the small-scale errors in the -5/3 regime start to saturate, and the large-scale errors in the -3 regime are still negligible. Therefore the normalized error could be estimated to be
Combining this with Eq. ( 15), once again, we have Eq. ( 22).
In light of Eq. ( 4) and Eq. ( 22), if the canonical atmospheric kinetic energy spectrum 𝐸(𝑘)
is known to us, then we can directly estimate the error growth behavior of the system using the analytical Eq. ( 14) proposed above, the parameter of this analytical equation can be calculated as follows:
To sum up, this simple analytical framework that we show is consistent with the error growth scenario described in Fig. 1. This framework is also well connected to the background Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
kinetic energy spectrum. All parameters in the analytical error growth model can be directly estimated from the energy spectrum of the background flow (Eq. 23).
It is natural to ask how well our proposed hybrid framework applies to the real atmosphere.
For a more direct comparison, simple dimensional results are used here. The largest length scale 𝐿 0 (corresponding to wavenumber 1) is chosen to be the circumference of a latitudinal cycle at the mid-latitude (~30000 km). The total kinetic energy of the background flow 𝐸 is estimated to be 150 m -2 s -2 , as in L69foot_3 . The units of distance and time are then 𝐿 * = 𝐿 0 2𝜋 ⁄ and 𝑇 * = 𝐿 * √𝐸 ⁄ ~ 4.25 𝑑𝑎𝑦, respectively. Therefore 𝑡 = 1 in the equation represents 4.25 days in the real atmosphere.
Table 1 shows the predictability limit derived from Eq. ( 9) of our hybrid framework and the results from L69. For L69, the predictability limit is simply the time when 𝑍 𝐾 (𝑡) = 𝑋 𝐾 . We use a 99% threshold for the calculation of the saturation time in our hybrid model under the Durran and Gingrich (2014) adjustment. Clearly, the predictability limits in our hybrid framework are much longer than L69. The reasons for this are twofold. On the one hand, L69 used a -5/3 slope across all the scales. By switching to the -3 slope at the synoptic scales as in the real atmosphere, our hybrid framework has less energy at smaller scales (consistent with observations), which leads to longer eddy turnover time, lower error growth rate, and hence longer predictability limit. On the other hand, the saturation approach we adopted from Durran and Gingrich (2014) will slow down the error growth rate as the errors approach their saturation threshold, which also will extend the predictability limit.
What is more intriguing is that the proposed hybrid framework shows an approximate twoweek limit for the synoptic scales at ~ 5000km. This limit agrees with our current understanding.
The same two-week limit for day-to-day weather predictability was first proposed by Lorenz in his early studies by analyzing the operational model products (Lorenz 1973(Lorenz , 1984;;Reeves 2014).
This limit is also found in today's most sophisticated numerical models (Foude et al. 2013;Judt 2018;Zhang et al. 2019). Moreover, the fitted 𝛼 and 𝛽 in Fig. 4b also agree well with the number estimated using full-physics convection-permitting global simulation (Fig. 3 in Zhang et al. 2019).
While all the predicted limits will vary proportionally with slightly different dimensional analysis (due to uncertainty in total kinetic energy E, for example), the ratio between the saturation times for different length scales will hold under dimensional process. Assuming this 2-week limit for the synoptic weather, then we learn from table 1 that the predictability limit for motions at ~1000 km is ~7 days, the predictability limit for ~500 km is ~ 5 days, and ~2 days for 100 km. All these numbers are generally consistent with the findings derived from complex, state-of-the-science modeling experiments in Zhang et al. (2019).
As we have mentioned before, we could also estimate the value of 𝛼 and 𝛽 from the kinetic energy spectrum directly. Utilizing the airplane data (Marenco et al. 2018), the observed atmospheric spectrum was fitted to a functional form in Lindborg (1999),
where 𝑑 1 = 9.1 × 10 -4 , 𝑑 2 = 3 × 10 -10 . Note 𝐸(𝑘) has units of m 3 ⋅ 𝑠 -2 , k has units of m -1 .
Therefore 𝑑1 has units of m 4/3 ⋅ 𝑠 -2 , 𝑑2 has units of 𝑠 -2 . The first term of Eq. ( 24) describes the shallower -5/3 wavelength range of the observed kinetic energy spectrum, while the second term Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
of Eq. ( 24) fits the synoptic -3 slope wavelength range. Figure 7 of Lindborg (1999) shows that this functional form of Eq. ( 24) works well for horizontal scales smaller than 1000 km when compared with the airplane observation data. Substituting Eq. ( 24) into Eq. ( 3) and assuming D equals 1.0, we find that the eddy turnover time for the synoptic-scale regime is around 16 hours that is independent of the wavenumber 𝑘 . The value of 𝛼 is then estimated to be around according to Eq. ( 23). Moreover, from Eq. ( 24) we can also tell the transition scale
can also be estimated to be on the order of 1/100 using Eq. ( 23). Hence, for the midlatitudes, an estimate of 𝛼 and 𝛽 from the energy spectrum can be given here,
These numbers, again, approximately match what we found in Fig. 4b and in earlier works done using full-physics models (Zhang et al. 2007(Zhang et al. , 2019)).
This proposed hybrid framework extends and complements earlier studies done by L69 and RS2008. By considering both the synoptic-scale 2DV-like dynamics (-3 slope) and the SQGlike motions (-5/3 slope) at smaller scales, our framework provides an improved understanding of the real atmosphere. The dimensional results also confirm that this hybrid framework gives more realistic estimations of the predictability limit compared to L69. To better understand the error growth process, we further derive a simple analytical equation for the evolution of the total error fields, which seems to work well with the idealized and full-physics simulation (Zhang et al. 2019).
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
We note here the canonical atmospheric kinetic energy spectra shown in Fig. 1 is really an idealization of the aircraft observations, which mainly reflect the upper levels in the mid-latitudes.
This canonical structure, however, is not always observed in the real atmosphere, which varies with different seasons, latitudes, and height levels. While the -3 slope and -5/3 slope are strictly assumed when building our new hybrid framework, it is easy to find that both numbers are likely not strictly required. On the smaller-scale end, the predictability limit, according to L69, will be limited as long as the slope is shallower than -3. L69 also showed that a -7/3 slope produced very close estimates of the predictability limit to the -5/3 scenario. Switching the slope of the smallscale spectrum from -5/3 to -7/3 or -2 in our hybrid model also has minor effects on the results and the predictability limits at large synoptic scales (not shown). The -3 slope for the synoptic scales plays a more significant role in determining the predictability limits in the sense that a steeper slope (e.g., -4 slope) will lead to longer predictability limits at the largest synoptic scales. As the 2DV and SQG dynamics do not support slopes other than -3 and -5/3, the experimental results of changing the slopes to different numbers based on the 2DV/SQG system are therefore less
convincing. Yet, based on the eddy turnover time argument, a slightly different slope shall not change the general picture ("two-stage" error growth process) shown in Fig. 1.
We shall acknowledge the fact that this newly proposed framework is based on L69.
Therefore all the assumptions made in L69 are still used in our current framework, which may pose limitations to the application of the framework. First, the statistical assumptions made in L69, such as homogeneity and isotropy, are not strictly valid for the real atmosphere. These assumptions do not allow any climatological mean motions and properties. It is also well known that systems like mountains and clouds are not randomly distributed. This heterogeneity could also be found in shows that the kinetic energy spectrum of the tropical region is different from the canonical spectrum shown in Fig. 1 (Judt 2020), indicating that the tropics may have very different error growth behavior than the mid-latitudes. Second, the specific dynamical equations adopted are also not strictly accurate for the real atmosphere. The two-dimensional vorticity equation is, at best, a very crude approximation for the large-scale dynamics. There is also no evidence showing that the small-scale dynamics could be described using surface quasi-geostrophic equations. Nonetheless, RS2008 has shown that, from the perspective of error growth, the kinetic energy spectrum slope might be more important compared to the dynamical equations used. The reason we choose these equations is also that their spectra are defensible on physical grounds so that we can combine them to construct a hybrid spectrum that is consistent with observational and full-physics modeling studies. With a realistic transition of the kinetic energy spectrum slope, we expect that our hybrid framework would capture the key components of the error growth behavior.
We also note here that we were not able to prove the causality between the kinetic energy spectrum and the error growth behavior in our framework. There are distinct differences between correlation and causation. It is also possible that the same physical processes/mechanisms lead to both the transition of the slopes and the error growth behavior simultaneously. For example, it has been hypothesized that moist convection and gravity waves generated by that might be responsible for the -5/3 slope at the small-scale end (Sun et al. 2017;Durran and Weyn 2016). Moist physics have also been shown to be the key for the upscale error propagation (Zhang et al. 2003(Zhang et al. , 2007;;Selz and Craig 2015). Hence, moist convection might be the actual source for both the shallower kinetic energy spectrum and the intrinsic predictability limit. If that is the case, slightly perturbing the moist physics scheme or the location of the convective grid will also lead to similar intrinsic predictability limit and error growth processes even in a coarse resolution model that is unable to Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
resolve the transition of the slope (e.g., Tribbia and Baumhefner 2004). Ongoing research is being done to study the underlying physical processes for error growth behavior and will be reported separately in the future.
Another limitation of this hybrid framework is the lack of vertical structure in the model.
Studies have shown that there is some degree of height dependence in the observed and simulated atmospheric spectra (Judt 2018). More critical differences may be found between the troposphere and the stratosphere (Skamarock et al. 2014). Given the connection between the background kinetic energy spectrum and the error growth behavior, it is very likely that differences in the spectrum could be associated with different predictability limits in the stratosphere and the troposphere. Their coupling will add another layer of complexity to the study (Butler et al. 2019).
With all these inadequacies aside, the most promising and encouraging finding of this study is that this simple new theoretical framework, which is built and based on the hybrid kinetic energy spectrum, could capture the error growth behavior found in the complex full-physics simulation.
This strong connection between the kinetic energy spectrum and the error growth process might also apply to other turbulent fluids, like the ocean, which is less understood now. Our simple model may provide a new perspective for the predictability of these turbulent fluids and beyond, as will be further examined in future studies.
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
for the convenience of dimensional analysis, the total energy ∑ 𝑋 𝑚 𝑛 𝑚=1
is set to be 1. As is shown in RS2008, computation of 𝐵 𝐾,𝐿 get increasingly difficult when 𝐾 = 𝐿 and they are both large, due to rapidly decreasing integration area. In this study, a total of 24 spectral bands are used whereas resolution factor 𝜌 is set to √2 , which means the smallest scale is ~10km after dimensional results.
Note we could also compute a single matrix 𝐶 and hence form an L69-type system of ODEs based on the hybrid spectrum in (A2) under 2DV dynamics (as done in Durran and Gingrich 2014), utilizing the fact that -5/3 is also an admissible spectral slope for 2DV if we add small-scale forcing. However, we believe this single matrix assumption is as unphysical as our current approach, if not more. If we examine the derivation process for matrix C in L69, the "inertial range" idea is implicitly adopted, where no energy source/sink is considered during the derivation.
Under this "inertial range" idea, it is unlikely for the 2DV system itself to present a hybrid spectrum automatically. It is also unphysical to assume the small-scale motions still obey the 2DV dynamics.
Therefore, we choose our current approach in the manuscript, which keeps the consistency between the "inertial range" assumption and the derivation of 𝐶 𝐾,𝐿 terms for different systems.
Nonetheless, this should have a minor effect on our results given the results shown in earlier studies (RS2008; Durran and Gingrich 2014)
The nonlinearity saturation effect introduced by Durran and Gingrich ( 2014) is also included in our study. The original set of n second-order differential equations in Eq. ( 2) can be rewritten to a set of 2n first-order differential equations. [𝐴𝑘 𝑙 (-𝑝+3) ] -1
2 ~ 𝜏(𝑘 𝑙 ) , 𝑝 < 3 (𝐵4)
Thus 𝑇 (the predictability limit), the time needed for errors at the smallest scales propagate to the largest scale, will grow larger and larger for a turbulent system with a steep slope 𝑝 ≥ 3. However, for 𝑝 < 3, predictability time remains finite no matter how we confine the initial error. And this finite predictability time has the same order of magnitude as the eddy turnover time at the largest scale 𝑘 𝑙 .
Accepted for publication in Journal of the Atmospheric Sciences. DOI10.1175/JAS-D-19-0271.1.
Adjacent spectral bands are differed by a constant resolution factor 𝜌 (𝜌 = √2 in this study). Assuming the length scales for all N spectral bands are 𝐷 0 , 𝐷 1 , … 𝐷 𝑁-1 , then we have 𝐷 0 = 𝜌 𝐾 𝐷 𝐾 for each spectral band K.
Given the power-law distribution of the base spectrum, the total initial error will decrease exponentially when we linearly increase K.
The observed spectrum transition happens at ~400km in the mid-latitudes, corresponding to zonal wavenumber ~70.
L69 use 148 m -2 s -2 for total energy. Density weighted total energy from reanalysis data give a strong seasonal variation, ranging from less than 100 m -2 s -2 in the summer and more than 200 m -2 s -2 in the winter.