Buoyancy driven distributed chaos and ensemble weather forecasting

It is shown, using results of direct numerical simulations, that strong thermal convection in horizontal layer and on a hemisphere can be well described by the distributed chaos approach with the stretched exponential kinetic energy spectrum. Two relevant cases: the vorticity and helicity dominated distributed chaos, have been considered. The results obtained with the Weather Research and Forecast Model (moist convection with and without Coriolis effect), with NCEP/Stage IV precipitation products, with the Deterministic Weather Forecasting Systems based on ensemblevariational data assimilation, and with the Coupled Ocean-Atmosphere Mesoscale Prediction System (COAMPS) were also used in order to demonstrate applicability of the distributed chaos approach. In the later case the ensemble forecasting of the real snowstorms was considered for this purpose.

It is shown, using results of direct numerical simulations, that strong thermal convection in horizontal layer and on a hemisphere can be well described by the distributed chaos approach with the stretched exponential kinetic energy spectrum. Two relevant cases: the vorticity and helicity dominated distributed chaos, have been considered. The results obtained with the Weather Research and Forecast Model (moist convection with and without Coriolis effect), with the Deterministic Weather Forecasting Systems based on ensemble-variational data assimilation, and with the Coupled Ocean-Atmosphere Mesoscale Prediction System (COAMPS) were also used in order to demonstrate applicability of the distributed chaos approach. In the later case the ensemble forecasting of the real snowstorms was considered for this purpose.

DISTRIBUTED CHAOS
Power spectrum of most dynamical systems with the chaotic behaviour (including some atmospheric phenomena) decays exponentially with frequency [1]- [8]. For some so-called multiscale systems, described by equations with partial (time and spatial) derivatives the spectra decay exponentially with wavenumber as well. Figure  1, for instance, shows (in the semi-logarithmical scales) perturbation kinetic energy spectrum obtained in a direct numerical simulation (DNS) of so-called 'isotropic homogeneous turbulence' described by the Navier-Stokes equations at Reynolds number Re ≃ 2500 [9]. The dashed straight line indicates the exponential decay At this DNS a copy of the velocity field u 1 was perturbed using a slight perturbation of the forcing function f at one particular timestep. A new field u 2 was created by this perturbation. The forcing was applied at the low wavenumbers only. The perturbation field δu = u 1 − u 2 was then computed as well as its power spectrum, E d (k, t) after a steady state was reached This spectrum is shown in the Fig. 1.
As one can see from the insert the peak of the E d (k) spectrum is reached at the same value of k ≃ k 0 that appears in the exponential spectral decay Eq. with the characteristic scale k 0 and the chaotic dynamics at the higher wavenumbers (see also next Section). If we have an ensemble of the fields (obtained by varying the initial conditions, for instance) with the statistically varying parameters a and k 0 , then the ensemble averaged spectrum should be computed using equation where the P (a, k 0 ) is a joint probability distribution of the variables a and k 0 . For statistically independent parameters a and k 0 where the P (k 0 ) is a distribution of the parameter k 0 .
Let us consider an important particular case of the distributed chaos. If the characteristic velocity u 0 , corresponding to the scale k 0 , varies with k 0 in a scale in-non-peer reviewed EarthArxiv preprint (2019) DOI: 10.31223/osf.io/hej3w variant (scaling) form and the vorticity ω(x, t) correlation integral (where < ... > denotes an ensemble average) governs the scaling Eq. (7), then we obtain from the dimensional considerations Normal (Gaussian) distribution of the characteristic velocity (see Section 3 below) results in the chi-squared (χ 2 ) distribution of the variable k 0 where k β is a constant. Substitution of the Eq. (10) into Eq. (6) results in

THERMAL CONVECTION
At the Rayleigh-Bénard (thermal) convection a horizontal layer of the fluid (or gas) is cooled from the top and heated from below. Under Boussinesq approximations the nondimensional equations of the thermal convection can be written as 1 Pr where Ra is the Rayleigh number, P r is the Prandtl number, θ is deviation of temperature from the heat conduction state andẑ is the buoyancy direction [10]. There are also other ways of normalization of the equations. In particular the non-trivial case lim P r → 0 [11], [12] has been numerically studied in the Ref. [13] and exhibits a broadband kinetic energy spectrum already near onset of chaos (at Ra/Ra c = 1.1, where Ra c is a critical value of the Ra for the onset of the convection,). Figure 2 shows the kinetic energy spectrum, computed at Ra/Ra c = 1.1 and at the lim P r → 0, in the log-log scales (here and in all other figures log k ≡ log 10 k). The spectral data were taken from Fig. 1 of the Ref. [13]. The dashed curve indicates the exponential spectrum Eq. (3) with k 0 ≃ 0.96, that indicates a large-scale nature of the exponential decay.
In a recent paper Ref. [14] results of a direct numerical simulation (DNS) of the Rayleigh-Bénard (thermal) convection were reported and Figure 3 shows kinetic energy spectrum computed for this DNS at Ra = 10 7 and P r = 10 2 (the spectral data were taken from Fig.  10 of the Ref. [14]). The dashed curve indicates the stretched exponential spectrum Eq. (11) in the log-log scales ( log k ≡ log 10 k).
In another recent Ref. [15] the Weather Research and Forecast Model [16] was used for simulation of the atmospheric moist convection. Seven localized warm bubbles (with a positive temperature anomaly) were put into the initial condition in order to initiate convection. These bubbles, were interacting with each other under influence of a wind shear (see the Ref. [15] for more details). Figure 4 shows in the log-log scales kinetic energy spectrum averaged between 0 and 15 km of the vertical height and over 4-6 hours of the time evolution (the spectral data were taken from Fig. 10 of the Ref. [15]). The dashed curve indicates the stretched exponential spectrum Eq. (11) in the log-log scales.

HELICITY DOMINATED DISTRIBUTED CHAOS
The above considered example of the vorticity dominated distributed chaos indicates the stretched exponential type of spectrum (see also Ref. [17]). Let us consider a generalization: If one denotes the distribution of the characteristic velocity u 0 as P(u 0 ), then P(u 0 )du 0 ∝ P (k 0 )dk 0 (16) or, taking into account the Eq. (7), On the other hand, for the stretched exponential spectrum Eq. (15) the P (k 0 ) asymptote at k 0 → ∞ is where b is a constant [18]. From the Eqs. (7), (17) and (18) one can conclude that the distribution P(u 0 ) is Gaussian (with zero mean), and the β is related to the α by the equation For the helicity dominated distributed chaos one should use the helicity correlation integral   instead of the above used vorticity correlation integral (let us recall that the helicity h = (ω · u) ). The I h is also known as the Levich-Tsinober invariant [19], and is usually related to the helical waves [20]. Substituting the helicity correlation integral into Eq. (7) and using the dimensional considerations one obtains: and then from the Eq. (19) β = 1/3, i.e. Figure 5 shows kinetic energy spectrum computed for the Weather Research and Forecast Model [16] numerical simulation of the atmospheric moist convection with the Coriolis effect (the spectral data were taken from Fig.  11a of the Ref. [15]). The dashed curve in the Fig. 5 indicates the stretched exponential spectrum Eq. (22) in the log-log scales (cf. previous Section, Fig. 4).
In a recent paper Ref. [21] results of a direct numerical simulation (DNS) of a Rayleigh-Bénard-like (thermal) convection on a hemisphere were reported. The system was heated at the equator and the gradient of temperature between the equator and the pole generates thermal plumes near the equator that move up from the equator zone toward the pole creating a thermal convection. Figure 6 shows kinetic energy spectrum computed for this DNS (the spectral data were taken from Fig.  18 of the Ref. [21] for the stationary state spectrum, Ra = 10 10 and P r = 7). The variable k is the spherical wavenumber. The dashed curve indicates the stretched exponential spectrum Eq. (22) in the log-log scales. Figure 7 shows the mean spectrum of kinetic energy in short-range weather forecasts (48-hours) experiment at 500 hPa. The spectral data were taken from the Fig. 7b of the Ref. [22] (obtained with the Deterministic Weather Forecasting Systems based on ensemble-variational data assimilation at Environment Canada). The dashed curve indicates the stretched exponential spectrum Eq. (22) in the log-log scales and covers Meso-and Synoptic scales up to the Planetary scales (indicated by the dotted vertical line).

ENSEMBLE WEATHER FORECASTING
The authors of a recent paper Ref. [23] have reported results of 100-member ensembles forecast for an East Coast winter snowstorm. An ensemble Kalman filter [24] was used in order to construct the 100-member ensembles which were integrated (with slightly altered initial conditions) for 36 hours ensemble forecast using COAMPS -Coupled Ocean-Atmosphere Mesoscale Prediction System [25].   Figure 8 shows in the log-log scales the mean background kinetic energy spectrum (ensemble-and meridional-averaged) at 500 hPa. The simulation was started at 1200 UTC 25 Dec 2010 (cyclogenesis for the storm) with real-atmospheric data. Figure 9 shows perturbation kinetic energy spectrum (ensemble-and meridional-averaged) at the 36 hours of the lead time (the spectral data were taken from the Fig. 6b of the Ref. [23]). In this case the perturbations were initially generated (initial error). The perturbation is the difference between the ensemble mean and one ensemble member. The dashed curves in the Figs. 8 and 9 indicate the stretched exponential decay Eq. (22). The Fig. 9 (see also Fig. 11) confirms the idea [23], [26]- [28] that the error growth can be mainly result of amplification of the errors at all wavenumbers.
Another winter storm in the Puget Sound region of the Pacific Northwest was studied by the same method in the Ref. [26]. Figure 10 shows average of the horizontal kinetic energy spectra at 700-hPa at 1200 UTC 17 Dec  2008 (the spectral data were taken from the Fig. 13 of the Ref. [26]). Figure 11 shows spectrum of the kinetic energy of the perturbations about the ensemble mean at 700 hPa at the 36 hours of the lead time (the spectral data were taken from the Fig. 14d of the Ref. [26]). The simulation was started at 0000 UTC 17 Dec 2008 with real-atmospheric data. The dashed curves in the Figs. 10 and 11 indicate the stretched exponential decay Eq. (11).

DISCUSSION
After the seminal paper Ref. [29] a vast amount of studies was devoted to the multiscale systems' predictability for the cases with power-law kinetic energy spectra like E(k) ∝ k −5/3 , E(k) ∝ k −3 etc. (see, for recent development Refs. [15], [30], [31] and references therein). One can see from the above considered examples that the stretched exponential spectra Eq. (15) (characteristic for the distributed chaos) should be taken into account as well, especially at the ensemble weather forecasting [32].
ACKNOWLEDGEMENT I thank A. Berera and R.D.J.G. Ho for sharing their data and discussions, and S. Vannitsem for comments.