Higher-order statistics based multifractal predictability measures for anisotropic turbulence and the theoretical limits of aviation weather forecasting

Theoretical predictability measures of turbulent atmospheric flows are essential in estimating how realistic the current storm-scale strategic forecast skill expectations are. Atmospheric predictability studies in the past have usually neglected intermittency and anisotropy, which are typical features of atmospheric flows, rendering their application to the storm-scale weather regime ineffective. Furthermore, these studies are frequently limited to second-order statistical measures, which do not contain information about the rarer, more severe, and, therefore, more important (from a forecasting and mitigation perspective) weather events. Here we overcome these rather severe limitations by proposing an analytical expression for the theoretical predictability limits of anisotropic multifractal fields based on higher-order autocorrelation functions. The predictability limits are dependent on the order of statistical moment (q) and are smaller for larger q. Since higher-order statistical measures take into account rarer events, such more extreme phenomena are less predictable. While spatial anisotropy of the fields seems to increase their predictability limits (making them larger than the commonly expected eddy turnover times), the ratio of anisotropic to isotropic predictability limits is independent of q. Our results indicate that reliable storm-scale weather forecasting with around 3 to 5 hours lead time is theoretically possible.


Results
Since the concepts of both prediction and predictability of complex systems are now widely accepted to be probabilistic [35][36][37][38][39][40] , we utilize a stochastic Generalized Scale-invariant (GSI) multifractal 27,30,[41][42][43][44] based methodology here. Atmospheric space-time scaling laws are of the general form 30 where |Δr| is the isotropic spatial scale function, L is the integral length scale (usually taken as the size of the largest eddy), and T is the eddy turnover time corresponding to L. The main difference between real space and Fourier space scale functions is that they are symmetric with respect to different generators: G s and G s T (this superscript 'T' indicates the transpose of a matrix, not to be confused with the integral time scale T or the scale transformation operator T λ ). For spatially isotropic, self-affine cases or GSI cases with no differential rotation of structures (e = 0), G s is symmetric so that = G G s s T . The Fourier space scale function (indicated by the subscript FS) corresponding to the real space scale function can, therefore, be taken as where |k| is the isotropic spatial scale function in Fourier space, k = (k x , k y , k z ) is the wavevector, K i and Ω i are the (angular) wavenumber and (angular) frequency corresponding to the integral length and time scales (L and T), respectively. Even though this study prefers the theoretical physics convention of using angular frequency ω and angular wavenumber k (although the word angular is often dropped in the manuscript for convenience) in the Fourier space instead of the alternate spectroscopy convention of using frequency and spectroscopic wavenumber, both conventions result in the same final outcome as long as they are used consistently. Although real and Fourier space scale functions are equivalent in a scaling sense, they are generally not identical. The canonical scale functions in Eqs. (3 and 4) are simply convenient approximations.
Semi-fourier space scale functions. Taking into consideration, the scaling anisotropy between space and time (based on the Kolmogorov-Obukhov law 45,46 ) suggested by earlier studies 19,47 , ω | | ∝ | | k H t can be squared on both sides and non-dimensionalized using Ω i and K i to get ∝ . ω Since a spatial-scale dependent but position independent expression for predictability limit is what is needed, it is advantageous to work in semi-Fourier space, fully exploiting this spatial position independence property. The purpose of squaring, is to obtain an expression for ω H 1/ t in www.nature.com/scientificreports www.nature.com/scientificreports/ terms of |Δt|, that can be raised to power 2 and used in Eq. (4) to replace the ω , which when using |ω|= 2π/|Δt| only on the right hand side becomes which when squared and substituted in Eq. (4) as discussed above gives as usual, and the subscript 'SFS' denotes the semi-Fourier space. Using theoretical predictability limits of q-th order statistical moments. The polyspectrum is defined as a generalized spectrum of order q (for q = 2, 4, etc. the polyspectrum is the spectrum, tri-spectrum, etc.), and due to statistical stationarity, the total polyspectrum are the decorrelated and correlated polyspectra, respectively. At Following the pioneering work of Lorenz 20 , that defines the predictability limit as the time until which errors in prediction have not exceeded a prechosen magnitude which for practical purposes should be considerably greater than typical observational errors but less than the magnitude of difference between randomly chosen states of the system, here the theoretical predictability limit Δt p is taken as the time when the correlated polyspectrum equals decorrelated polyspectrum (i.e. . Applying the outcome of the Methods section along with Eq. (6) in this definition results in . Equation (7) when solved for Δt p (q) results in the analytical expression This predictability limit obtained in a spatially isotropic framework can be directly translated to that of a spatially anisotropic framework by simply replacing the spatially isotropic function |k| by a GSI scale function ||k|| FS and the integral length, time scales (L,T) by the sphero-scale l s (the spatial scale where structures are roundish), sphero-time l st (the turnover time corresponding to l s ) respectively resulting in the analytical expression for the predictability limit Δt p (q) of the q-th order statistical moment of the atmospheric field considered: having linear GSI parameters c = e = f = 0, d = 1, horizontal trivial anisotropy parameter a and vertical scaling anisotropy parameter H z . The sphero-scale l s is the scale at which the structures of the field are approximately roundish, the sphero-time l st is the turnover time of eddies of size l s , whereas k s and k st are the sphero-wavenumber and sphero-frequency, respectively. Since (ε is the spatially averaged energy flux) and |k| is independent of L, Eq. (8) is independent of L. Equation (9) on the other hand depends on l s since ||k|| FS depends on l s (although ). The critical ratio μ of the correlated polyspectra to that of the total polyspectra (the polyspectrum is a generalized spectrum of order q as described earlier) specifies how much error in prediction is acceptable, a and H z are the trivial (non-scaling) horizontal and scaling vertical anisotropy parameters respectively, whereas H t is the space-time anisotropy parameter. The scaling moment function K(q) is given as where C 1 is the codimension of the mean, and α is the index of multifractality. The basic cascade equation for a scale by scale conserved multifractal field is given as: , where the angular brackets denote ensemble averaging, q the order of the statistical moment, and the scale ratio λ= Largest scale/intermediate scale. The scaling moment function K(q) that describes how the statistical properties of each moment behave under scale transformations 48 is also the (base λ, Laplace) second characteristic function (SCF). The smoothness parameter H is used to obtain non-conservative observed fields from conservative multifractal cascade processes, whereas η is the exponent of the conservative turbulent flux (not to be confused with the Kolmogorov scale usually denoted by η in turbulence literature). The significance of η is that for observed non-conservative fields such as the velocity shear across a scale l, based on the physical notion of eddy turnover time or purely dimensional considerations are directly proportional to the η-th power of the conservative fields such as the energy flux density (for the case of horizontal wind shear η = 1/3). Spatial GSI scale functions are denoted by || ||, and in Fourier space (k x , k y , k z ) as ||k|| FS , where the wave vector k = (k x , k y , k z ) has the Euclidean

empirical-parameter based estimate of predictability limits. Empirical estimates of multifractal
parameters α = 1.5,C 1 = 0.15 used by earlier works 26 and 45,46 are used along with horizontal, vertical anisotropy parameters a = 1.6 (an ECMWF interim flux based estimate 49 ), H z = 5/9 (following Schertzer and Lovejoy 50 ), a typical storm-scale sphero-scale of 100 m (determined from CloudSat data by Lovejoy et al. 51 ), μ = 0.5 (following the critical ratio used by Schertzer and Lovejoy 19,22 ) in Eq. (9) for assessing the predictability limits of horizontal wind fields. In this study, the isotropic wind fields have a = 1, H z = 1; horizontally anisotropic fields have a = 1.6, H z = 1; whereas the horizontally and vertically anisotropic wind fields have a = 1.6, H z = 5/9. For such horizontally and vertically anisotropic wind fields in the convective regime which is typically 100 km in the horizontal 30 and 10 km in the vertical, the maximum predictability limits occur at the largest scales of the regime (100 km,100 km,10 km) as can be inferred from Eq. (9) and are about 5 hrs and 4 hrs for q = 2 and q = 4 respectively. The limits derived using larger q values are smaller as expected (as illustrated by Fig. 1), since rarer events are less predictable 19,22 . Fields with a larger sphero-scale of 1000 m have predictability limits that are smaller (the maximum values are about 4hrs and 3hrs for q = 2 and q = 4 respectively) than those of fields with 100 m spheroscale, as anticipated (as illustrated by Fig. 2). This is a consequence of the sub-and super-spheroscales being dominated by the buoyancy variance flux and energy flux respectively 50 , and systems with stratiform dynamics having lesser buoyancy variance flux than those with convective dynamics 30 . Finally, these figures also show that anisotropy improves predictability in accordance with earlier spectra based assessments. probability of occurrence. The scaling moment function, K(q) is related to the codimension of the order is denoted by γ q and is the solution of the equation c′(γ q ) = q, where c′(γ) is the first derivative of the codimension of order of singularities in γ. For q = 2, 4 the corresponding γ 2 ,γ 4 are computed using this equation with empirical estimates of α,C 1 (as discussed earlier). By substituting these order of singularity values, the corresponding codimension of order of singularities are then computed. The probability of occurrence of events above a scaling threshold 30 is then obtained for different scales but with the largest scale being 10000 km (since in the atmosphere, the scales typically vary from 1mm to 10000 km). A straightforward mathematical simplification of the above discussion shows that γ = α c C q ( ) q 1 (where both α and C 1 are positive), implying that larger the statistical order q larger the codimension of singularities γ c( ) q corresponding to that order and therefore smaller the probability of occurrence of flux events with order of singularities exceeding γ q . This means that higher-order statistical moments are more representative of less probable or extreme events. Figures 1 and 2, although informative do not directly show how the predictability limits are affected due to anisotropy at super-spheroscales (scales larger than the sphero-scale) and sub-spheroscales (scales smaller than the spheroscale). Since this sphero-scale l s is the same in all three directional planes, it is independent of the direction. Therefore, it is necessary to get the predictability limit into a similar direction independent format. To do this Δt p has to become independent of both the azimuthal and polar angles, and angular averaging seems to be the simplest way of doing this in a more generalized manner (this loss of directional information along with the scale being limited by the smallest of the three scales are the drawbacks of doing this). Angular averaging Eq. (9) for further investigation, results in the angular averaged predictability limit [Δt p (q)] AA given by = .

Discussion
= . C H 1 5, 015, 033, 1 and anisotropy parameters as discussed in the text. The wavenumbers k x , k y , k z have units of km −1 and are in the x, y, z directions respectively, whereas the predictability limits are in hours. The horizontal wavenumbers represent a scale range from 1 to 100 km, whereas the vertical wavenumber represents a scale range from 1 to 10 km (scales smaller than 1 km are not shown here as their predictability limits are less than 1 hour). The two rows differ only by the order of autocorrelation q used in deriving the predictability limits. (a) spectra based predictability limits (q = 2) of an isotropic horizontal wind field as a function of logarithmic wavenumber. (b) same as a, but for vertically anisotropic cases. (c) same as a, but for horizontally and vertically anisotropic cases. (d-f) are same as (a-c) but for q = 4. Comparing figures (a-c) with (d-f) shows that higher-order predictability limits are smaller than lower-order ones. Since higher-order statistical moments represent more extreme events, these figures indicate that such events are less predictable. Comparing figures (a-c) with each other and figures (d-f) with each other illustrates that fields that are more anisotropic are more predictable. = .
= . C H 1 5, 015, 033, 1 and anisotropy parameters as discussed in the text. The wavenumbers k x , k y , k z have units of km −1 and are in the x, y, z directions respectively, whereas the predictability limits are in hours. The horizontal wavenumbers represent a scale range from 1 to 100 km, whereas the vertical wavenumber represents a scale range from 1 to 10 km (scales smaller than 1 km are not shown here as their predictability limits are less than 1 hour). The two rows differ only by the order of autocorrelation q used in deriving the predictability limits. (a) spectra based predictability limits (q = 2) of an isotropic horizontal wind field as a function of logarithmic wavenumber. (b) same as a, but for vertically anisotropic cases. (c) same as a, but for horizontally and vertically anisotropic cases. (d-f) are same as a, b, c but for q = 4. Comparing figures (a-c) with (d-f) shows that higher-order predictability limits are smaller than lower-order ones. Since higher-order statistical moments represent more extreme events, these figures indicate that such events are less predictable. Comparing figures (a-c) with each other and figures (d-f) with each other illustrates that fields that are more anisotropic are more predictable. Comparing Fig. 1   where θ and φ are the polar and azimuthal angles, respectively. It is sufficient to consider only non-negative wavenumbers k x , k y and k z since the Fourier space scale function ||k|| FS we use is an even function, as can be seen from Eq. (9), due to which the results are symmetric in the negative part. Therefore, we consider only the first quadrant π ( ) 0, 2 where (sinφ cosθ), (sinφ sinθ) and (cosφ) are all non-negative. Figure 3a) shows the spectra based angular averaged predictability limits (q = 2) of a horizontal wind field with l s = 100 m, as a function of logarithmic wavenumber. The probability of occurrence shows the probability distribution of energy fluxes (ε) above the scaling threshold 30 Intermediate scale is the scale ratio and γ is the order of singularity) given by is the codimension of the order of singularities and is related to K(q)), corresponding to wavenumbers 10 0 ,10 0.5 ,10 1.0 ,10 1.5 ,10 2.0 ,10 2.5 . The angular averaged predictability limits of the isotropic, horizontally isotropic, vertically anisotropic, horizontally and vertically anisotropic wind fields are denoted by , , , respectively. At scales around 10 km, about 5% of the energy fluxes contribute to the second-order statistical moment, and the horizontally and vertically anisotropic horizontal wind field corresponding to these fluxes has an angular averaged predictability limit of 3.1hrs. Figure 3b (although they are not exactly equal in these figures, they are both related through l). At scales around 10 km, about 0.025% of the energy fluxes contribute to the fourth-order statistical moment, and the horizontally and vertically anisotropic horizontal wind field corresponding to these fluxes has an angular averaged predictability limit of 2.7 hrs. Vertical stratification seems to improve predictability only in the subsphero-wavenumbers. Figure 4 is the same as Fig. 3 but is for a 1000 m sphero-scale. At scales around 10 km, about 5% of the energy fluxes contribute to the second-order statistical moment, and the horizontally and vertically anisotropic horizontal wind field corresponding to these fluxes has an angular averaged predictability limit of 2.5 hrs. Figure 4b) same as 4a), but for q = 4. At scales around 10 km, about 0.025% of the energy fluxes contribute to the fourth-order statistical moment, and the horizontally and vertically anisotropic horizontal wind field corresponding to these fluxes has an angular averaged predictability limit of 2.2 hrs.
Finally, Fig. 5 shows the importance of sphero-scales. The two panels differ only by the sphero-scale l s , and the ratios of the angular averaged predictability limits are independent of the order of statistical moment q. In Fig. 5a) l s = 100 m, at scales around 10 km, about 5% and 0.025% of the energy fluxes contribute to the second and fourth-order statistical moments and the horizontal wind field corresponding to these fluxes when subject to both horizontally and vertically anisotropic has an angular averaged predictability limit that is 2.15 times the isotropic limit. Figure 5b) is the same as a), but has l s = 1000 m. At scales around 10 km, about 5% and 0.025% of the energy fluxes contribute to the second and fourth-order statistical moments and the horizontal wind field corresponding to these fluxes when subject to both horizontally and vertically anisotropic has an angular averaged predictability limit that is 1.75 times the isotropic limit. Vertical stratification enhances and diminishes predictability in the subsphero and supersphero-wavenumbers. Horizontal stratification improves predictability over all scales, although the improvement is not very significant. Figure 5a,b) show that wind fields with smaller spheroscales are better predictable as they are less dominated by convective dynamics.
In conclusion, (i) the super and subsphero-scale predictability is enhanced and diminished correspondingly for scaling anisotropy, (ii) for trivial anisotropy, predictability over the entire scale range is improved in accordance with spectra based estimates, (iii) reliably forecasting convectively less dominant systems that are more probable to occur with around 5 hours lead time seems to be theoretically possible (in case of less probable events this lead time is reduced to 4 hours), iv) whereas reliably forecasting convectively more active systems that are more probable to occur with around 4 hours lead time seems to be theoretically possible (in case of less probable events this lead time is reduced to 3 hours). Although convective scale numerical models are capable of reliably forecasting events with stronger large-scale forcing (organized convection) sometimes even out to 4 hours (around the theoretical limit proposed for convectively active and more probable events) in some cases (when initialized with high-resolution Doppler radar observations), their ability to predict air-mass type storms (unorganized convection) is still quite low [53][54][55][56] (not even close to the theoretical limit proposed here for convectively active and more extreme or less probable events). Even recent predictability studies 57,58 using storm-scale ensemble forecasting systems conclude that these convective-allowing sophisticated NWP models perform poorly (the predictability of scales smaller than 100 km is totally lost at around 1hr -which is quite low compared to the theoretical limits proposed in this present study) especially for quantitative precipitation forecasting, and that further effort is therefore needed in improving the basic understanding storm-scale weather predictability. The results of this (2019) 9:19829 | https://doi.org/10.1038/s41598-019-56304-2 www.nature.com/scientificreports www.nature.com/scientificreports/ study, demonstrates that the current expectations of reliable aviation weather forecasts with 2 to 6 hrs lead times (i.e., strategic aviation weather forecasting) are not totally unrealistic subject to the incorporation of multifractal cascade dynamics based modeling strategies (especially for unorganized convective weather phenomena which are difficult to forecast using conventional convective scale NWP models).

Methods estimation of the scaling exponent of correlation polyspectra.
Following earlier studies 27,29 and using Eq. (1), the q-th order structure function, for even integer q can be written using the binomial theorem as . Equation (11) implies that the total statistical moment of order q  = .
= . C H 1 5, 015, 033, 1 and anisotropy parameters, as discussed in the text. The subscripts ip,hap,vap,hvap denote isotropic, horizontally anisotropic, vertically anisotropic, horizontally and vertically anisotropic cases, whereas the subscript AA indicates angular averaging. While the angular averaged predictability limits is in hours, the wavenumber k is in km −1 and corresponds to scales ranging from 10 m to 10 km. The two panels differ only by the order of autocorrelation q used in deriving the predictability limits. (although λ and k are not exactly equal in these figures, they are both related through l). (b) same as a, but for polyspectra based angular averaged predictability limits (q = 4). The curves in both figures (a,b) show that fields which are more anisotropic are more predictable (over scales larger than the sphero-scale l s ) and that larger events are more probable and predictable, whereas comparison between these two figures indicates that more extreme events (less probable) are less predictable. www.nature.com/scientificreports www.nature.com/scientificreports/ responding generalized correlation polyspectral ((q − 1)th order) density 61 (the subscript c q of the polyspectral density means that it is the correlated polyspectral density and is dependent on the order of the statistical moment q, but it does not mean that it is of the q-th order) are related to each other via the (q − 1)d dimensional inverse Fourier transform where the whole domain of integration is over the region in D (q−1)d (since here we deal with three-dimensional space and one-dimensional time, d = 4) as there are (q − 1) of these d-dimensional integrals. Due to the assumption of statistical translational invariance the right-hand side is independent of R and is dependent only on the lags Δ Δ … − R R , , q 1 1 . When (q − n) of these (q − 1) lags are equal to ΔR, and the remaining (n − 1) lags are zero Eq. (12) becomes  (although λ and k are not exactly equal in these figures, they are both related through l). (b) same as a, but for polyspectra based angular averaged predictability limits (q = 4). The curves in both figures (a,b) show that fields which are more anisotropic are more predictable (over scales larger than the sphero-scale l s ) and that larger events are more probable and predictable, whereas comparison between these two figures indicates that more extreme events (less probable) are less predictable. Comparing Fig. 3 with this figure shows that fields with smaller sphero-scales are more predictable.
Following a procedure similar to that used in the derivation of spectra 37 based predictability limits, the scaling exponent of P c q is derived, as shown here. From the generalized scaling law Eq. (1), it follows that Using the anisotropic functional scale equation 30 (i.e., Eq. 2) in Eq. (15) gives G is the Fourier space scaling operator with the Fourier space generator matrix  G, K = (k, ω) is the Fourier space wave vector and D el is the elliptical space-time dimension. The scale invariance of the scalar product K.ΔR implies that =  G G T . From Eq. (17) Figure 5. Ratio of the angular averaged theoretical predictability limits of anisotropic horizontal wind fields to that of isotropic horizontal wind fields. The wind fields have multifractal parameters α = . = .
= . C H 1 5, 015, 033, 1 and anisotropy parameters, as discussed in the text. The two panels differ only by the sphero-scale l s . The subscripts ip,hap,vap,hvap denote isotropic, horizontally anisotropic, vertically anisotropic, horizontally and vertically anisotropic cases, whereas the subscript AA indicates angular averaging. The wavenumber k = |k| is in km −1 and corresponds to scales ranging from 10 m to 10 km. The ratios of the angular averaged predictability limits are independent of the order of statistical moment q and, therefore, of how probable the occurrence of the event is. (a) the spatially anisotropic fields have l s = 100 m. (b) same as a, but the spatially anisotropic fields have l s = 1000 m. The probability of occurrence values in Fig. (a,b) are the same as those in Figs. 3 and 4. Vertical stratification enhances and diminishes predictability in the subsphero and supersphero-wavenumbers. Horizontal stratification improves predictability over all scales. Wind fields with smaller spheroscales are better predictable. Larger events are more probable and have larger anisotropic to isotropic angular averaged predictability limit ratios. ( The general solution of this functional equation (Eq. (19)) [found by adopting a procedure following chapter 7 of Lovejoy and Schertzer 30 , similar to that subsequently used in Appendix A of Ramanathan et al. 43 , and using it along with the anisotropic functional scale equation] is ∝ The q-th order correlated polyspectrum is therefore is the polyspectral exponent (H and η are the conservation or fluctuation exponent and exponent of the conservative turbulent flux respectively, whereas is the moment scaling function [not to be confused with the Fourier space wave vector K = (k, ω)]). By repeating the above steps but for spatial scaling laws instead of space-time scaling laws and assuming statistical stationarity, it follows that the total polyspectrum ∝ β − k E T q q .