Additive Model Perturbations Scaled by Physical Tendencies for Use in Ensemble Prediction

1.1 Model uncertainty

A numerical weather prediction model computes the forecast by time stepping. At each time step, the model has on input, typically, the current model state (defined on a spatial grid), and computes the model state at the next time step or, equivalently, computes the forecast tendency, the difference between the next-time-step and the present model states. The model uncertainty is, by definition, the uncertainty in the forecast tendency, e.g. Orrell et al. (2001). Being accumulated and transformed during the time stepping, the uncertainty in the tendency contributes to the uncertainty in the forecast fields.

The model uncertainty is caused by the following imperfections in the forecast model (listed in order of increasing importance).

1. Atmospheric model’s partial differential equations normally involve a number of simplifications like the neglect of variations (horizontal and vertical) in the gravitational force or the ideal gas law assumption. The model may also omit some processes like chemical reactions or development and impact of electric charges of hydrometeors.

2. A classical (i.e., based on laws of physics, not a neural network) meteorological forecast model solves a set of time and space-discretized differential equations. The discretization leads to truncation error.

3. Subgrid-scale processes are accounted for in atmospheric models in an approximate manner. On the one hand, these processes cannot be reproduced by the discretized model equations because the model grid is too coarse. On the other hand, subgrid-scale processes do impact grid-scale fields due to the nonlinearity of physical laws that govern the evolution of the atmosphere. In forecast models, this impact is assessed using simplified sub-models known as physical parametrization schemes. Simplifications (made in these schemes for computational reasons) along with the inherent uncertainty in the unresolved scales lead to uncertainties/errors in physical parametrizations, e.g., Palmer et al. (2005). Note that some processes in the atmosphere (such as turbulence, convection, and gravity waves) can become increasingly resolved by the model equations just by refining the computational grid, that is, by increasing the model’s spatial and temporal resolution. Therefore, these processes can, actually, be regarded as part of truncation error. The respective physical parametrization schemes, thus, attempt to reduce truncation errors while introducing their own (presumably, smaller) errors/uncertainties.

Our focus in this research is on uncertainties/errors in physical parametrizations of subgrid-scale processes.

1.2 Error and uncertainty

At each model time step, nonlinear interactions of the subgrid-scale field x_SGS with itself and with the grid-scale field x_GS produce a spectrum of combination wavenumbers some of which fall within the resolvable (grid-scale) range. This yields a contribution of the subgrid-scales to the (grid-scale) forecast tendency, see, e.g., Touil et al., 2007, for a formulation of the problem in turbulence theory. Having a perfect model, we could compute this contribution, let it be denoted π(x_SGS, x_GS). Deterministic physical parametrization schemes, having, by definition, access only to x_GS, attempt to assess this impact, producing a physical tendency, P(x_GS). The uncertainty in all physical parametrization schemes combined is the difference between P(x_GS) and the true physical tendency π(x_SGS, x_GS).

The subgrid-scale field x_SGS is not explicitly defined in the model, therefore it is unknown and even unknowable to the model (to the extent that x_SGS is not determined by x_GS). Therefore we assume that x_SGS is a random field with some conditional probability density p(x_SGS | x_GS) (we condition on x_GS because it is available to the model). This density induces the probability density p(π | x_GS). Then, the ideal (best in the mean square sense) deterministic physical tendency is the conditional expectation of the true tendency given the grid-scale field, e.g., Kwasniok (2012); Shutts and Pallarès (2014), $P_{\text{ideal}}^{\text{determ}}(x_{\text{GS}})=\mathbb{E}(\pi |x_{\text{GS}})$, where 𝔼 stands for expectation. It is $P_{\text{ideal}}^{\text{determ}}(x_{\text{GS}})$ that a deterministic physical parametrization seeks to approximate by P(x_GS). In these terms, it is meaningful to call

the error in the deterministic physical tendency P. However, the deterministic physical tendency is incapable of simulating the variability of π around its conditional mean value $\mathbb{E}(\pi |x_{\text{GS}})$,

(Palmer, 2012). This is the inherent (irreducible) uncertainty in π(x_SGS, x_GS) given x_GS. Then, the full uncertainty is the difference between the error ɛ_P and the irreducible uncertainty δ_P:

The irreducible uncertainty δ_P is especially important in so-called ‘gray zones’, where a process is partly resolved by the model, meaning that a length scale of x_SGS is comparable to the grid spacing. In this case, in each grid cell, the grid-scale impact of the subgrid-scale process is, effectively, a sum of a small (and random) number, ν, of random contributions: $\pi ={\Sigma }_{\ell =1}^{\nu }{\pi }_{\ell }$. With convection, the contributions ${\pi }_{\mathcal{l}}$ to the convective tendency are due to individual convective plumes (Plant and Craig, 2008). With boundary-layer turbulence, ${\pi }_{\mathcal{l}}$ are due to individual turbulent eddies (Kober and Craig, 2016). Assuming, for presentation purposes, that all ${\pi }_{\mathcal{l}}$ are the same in a grid cell and their number ν follows the Poisson distribution, we readily obtain that $\left|\mathbb{E}\pi \right|/\text{sd}\pi =\sqrt{\mathbb{E}\nu }$, where sd denotes standard deviation. This expression shows that the standard deviation (which measures randomness) and the mean value of the grid-scale impact of subgrid-scale processes are comparable to each other in magnitude if $\mathbb{E}\nu ~1$, which is the case in a gray zone as discussed above. This makes the irreducible uncertainty in π indeed substantial and implies that it needs to be properly accounted for in ensemble prediction.

It is worth remarking that the irreducible uncertainty δ_P caused by the randomness of π given x_GS is the aleatory (truly random) uncertainty, therefore it cannot be called error. Whereas ɛ_P is a kind of systematic error (reducible conditional bias), which is caused by imperfections in the physical parametrization schemes. This latter kind of uncertainty due to lack of knowledge is generically called epistemic. Both aleatory and epistemic uncertainties need to be taken into account in building an ensemble prediction scheme. The ideal random physical tendency $P_{\text{ideal}}^{\ast }(x_{\text{GS}})$ is a ‘possible true π consistent with the grid-scale field x_GS’, that is, a random draw from the probability distribution p(π | x_GS). It is reasonable to anticipate that the aleatory part of the uncertainty, being caused by random subgrid noise, should be associated with small spatial and time scales. On the contrary, the epistemic part of the uncertainty caused by systematic and, likely, flow-dependent deficiencies of the physical schemes, may be characterized by larger spatio-temporal scales.

In an ensemble, representing both kinds of uncertainty in physical tendency using pseudo-random spatial fields results in model perturbations introduced at each model time step during the forecast. Techniques that add stochasticity on top of existing deterministic physical parametrization schemes or replace the deterministic schemes with stochastic ones are known as ‘stochastic physics’.

1.3 Existing stochastic physics schemes

The first question is how to represent the aleatory uncertainty in physical tendency in ensemble prediction schemes. We believe that the most sensible way to address this question is to build intrinsically stochastic (rather than traditional deterministic) parametrization schemes. Research in this direction is underway, see Plant and Craig (2008); Dorrestijn et al. (2013); Sakradzija et al. (2015); Hirt et al. (2019); Machulskaya and Seifert (2019); Clark et al. (2021) and others. The second question is how to represent the epistemic uncertainty in physical tendency produced either by a deterministic or a stochastic physical parametrization scheme.

We state that, on the one hand, stochastic parametrizations have not yet replaced deterministic ones. So we still need techniques to represent uncertainties in deterministic physical parametrizations. Ad-hoc approaches are in wide use here. On the other hand, stochastic parametrizations are not going to be devoid of epistemic uncertainties either. To represent those uncertainties, we will, most likely, resort to ad-hoc schemes, too. Our focus in this research is on ad-hoc model perturbation schemes.

Currently, the two most popular ad-hoc techniques to represent uncertainties in physical parametrizations are Stochastically Perturbed Parametrization Tendencies (SPPT, Buizza et al., 1999; Leutbecher et al., 2017) and Stochastically Perturbed Parametrizations (SPP, Ollinaho et al., 2017). SPPT generates model perturbations relying on the assumption that the magnitude of the error in the physical tendency is proportional to the magnitude of the physical tendency itself. There is also a flavor of SPPT called iSPPT in which tendencies from different physical parametrizations are perturbed independently (Christensen et al., 2017). SPPT has proven to be very useful in practical ensemble prediction schemes despite its lack of physical consistency: it does not respect conservation laws because it perturbs tendencies without modifying fluxes, see Leutbecher et al. (2017), Lang et al. (2021), and references therein.

In SPP, selected parametric and structural elements of the model’s parametrization schemes are made spatio-temporal random fields rather than fixed numbers and fixed choices. Advantages of SPP include (i) reliance on expert knowledge to select perturbed elements and design probability distributions, (ii) internal consistency of the resulting numerical scheme and conservation properties, and (iii) capability of generating significant spread in the ensemble (Lang et al., 2021; McTaggart-Cowan et al., 2022). The disadvantages of SPP are more conceptual. First, it accounts only for uncertainties that can be captured by perturbing the specific parametrization schemes used in the forecast model in question. Given the inevitably simplified nature of many parametrization schemes, the resulting tendency perturbations may not explore some relevant directions in phase space. Second, it is hard to even suggest how, say, parameter-perturbation probability distributions can be objectively justified. The reason is that parameters (and structural elements) of parametrization schemes may have no counterparts in nature (there is no ‘diffusion coefficient’ in nature) and even in a high-resolution model used as a proxy to the truth. Besides, both SPPT and SPP can lead to biases, see, e.g., Leutbecher et al. (2017) and Bouttier et al. (2022), respectively.

We selected SPPT (described in more detail in section 2.1) as a starting point for our development in this study because it attempts to do exactly what is needed to represent uncertainty in physical parametrizations: it perturbs the physical tendency (see above section 1.2).

In this study, we analyze limitations of SPPT and build a technique that attempts to address those limitations. The new scheme termed Additive Model-uncertainty perturbations scaled by Physical Tendencies (AMPT) is tested in numerical experiments with a convective-scale ensemble prediction system.

2.1 Background on SPPT and notation

To facilitate the presentation of AMPT, we first outline SPPT (Leutbecher et al., 2017). By $\text{P}(x,y,\hspace{0.17em}\hspace{0.17em}\zeta ,\hspace{0.17em}\hspace{0.17em}t)=(P_1(x,y,\hspace{0.17em}\hspace{0.17em}\zeta ,\hspace{0.17em}\hspace{0.17em}t),\dots ,P_{n_{\text{fields}}}(x,y,\hspace{0.17em}\hspace{0.17em}\zeta ,\hspace{0.17em}\hspace{0.17em}t))$ we denote the vector-valued physical tendency (the net physical tendency, that is, generated by all physical parametrizations combined) at the spatial grid point with the Cartesian horizontal coordinates (x, y), the vertical coordinate ζ, and the forecast time t. Here P_i(x, y, ζ, t) is the component of P(x, y, ζ, t) in the i-th model field (variable) X_i, and n_fields is the number of model fields selected to be perturbed (most often, temperature, winds, and humidity, e.g., Christensen et al., 2017).

In SPPT, the perturbed physical tendency P* is postulated to be

where ξ(x, y, t) is the zero-mean and unit-variance spatio-temporal random field and κ the scalar (i.e. the same for all model variables) parameter that controls the magnitude of the perturbation. For stability reasons (to avoid sign reversal of the physical tendency, see, e.g., Leutbecher et al., 2017), the support of the probability distribution of κ ξ(x, y, t) is limited to the segment [–1,1] so that

This constraint limits the magnitude of perturbations in SPPT.

From Eq. (4), the perturbation of the physical tendency in SPPT is seen to be multiplicative with respect to the model physical tendency:

Here and elsewhere Δ denotes a perturbation.

2.2 Motivation

The following deficiencies of SPPT led us to propose the new approach.

1. In SPPT, perturbations are large (small) when and where the physical parametrizations generate a large (small) physical tendency P. This formulation gives rise to a meaningful scaling of perturbations in situations when the model predicts a high or moderate intensity of subgrid-scale processes and produces a large or moderate physical tendency. However, it cannot cover situations in which the physical tendency appears to be small or even zero whereas the model uncertainty is, in fact, large. This may occur if, for example, in a grid cell, convection is initiated in nature whilst the convective parametrization fails to be activated (note that in this case switching to iSPPT would not help either).

2. Equation (6) implies that the multivariate perturbation vector ΔP(x, y, ζ, t) is strictly proportional to the physical tendency vector P(x, y, ζ, t). As noted by Leutbecher et al. (2017), this implies that SPPT tacitly assumes that only the magnitude of the vector P is in error and not its direction, which is highly unlikely. In other words, SPPT ‘assumes’ that the ratios of the physical tendencies in different variables i and j at the same point (x, y, ζ, t) are error-free. As a consequence, the SPPT perturbations (and thus the assumed model uncertainties) are perfectly (100%) correlated or perfectly (–100%) anticorrelated for any pair of model variables, which is not realistic.

   Indeed, from Eq. (6) written component-wise, we have ΔP_i = κ ξ P_i, where κ is a positive constant and 𝔼 ξ = 0. This implies that 𝔼 (ΔP_i | P_i) = 0, and Cov (ΔP_i, ΔP_j | P_i, P_j) = 𝔼 (ΔP_i ΔP_j | P_i, P_j) = κ²(sd ξ)² P_iP_j (where Cov stands for covariance). Taking into account that sd (ΔP_i | P_i)= κ |P_i| sd ξ and sd (ΔP_j | P_j) = κ |P_j| sd ξ, the correlation, that is, the covariance normalized by the product of the two standard deviations, becomes Corr (ΔP_i, ΔP_j | P_i, P_j) = ±1.

3. Similarly (and also noted by Leutbecher et al., 2017), since the SPPT random pattern ξ(x, y, t) does not depend on the vertical coordinate, the SPPT perturbations are perfectly coherent (correlated) for all variables at all levels in a vertical column, which again is unrealistic.

4. Moreover, it has appeared that for SPPT to give rise to a significant spread in the ensemble, the length and time scales need to be really large for convective-scale models. E.g., in Maurer et al. (2014), the tuned length and time scales in the 2.2-km resolution COSMO model were as large as 500 km and 6 h, respectively. This implies that the above unphysical ±100% correlation of SPPT perturbations approximately holds for all model variables in huge 4D volumes spanning the whole atmosphere in the vertical, hundreds of kilometers in the horizontal, and hours of forecast time.

5. In some models, with their specific configurations of model perturbation schemes, SPPT appeared to generate not enough spread in ensemble forecasts, preventing them from producing reliable probabilistic forecasts, see e.g. Christensen et al. (2017); Frogner et al. (2022). One possible reason for that is the limitation on the magnitude of perturbation, Eq. (5).

With iSPPT (Christensen et al., 2017), the above points 1–3 pertain to the tendency due to a single parametrization scheme. iSPPT alleviates weak points 4–5, allowing for smaller spatial scales of the random pattern and having the capability to generate a somewhat larger spread in the ensemble than SPPT (Wastl et al., 2019; Christensen et al., 2017).

For completeness, it is worth mentioning that in developing AMPT we did not aim at alleviating SPPT’s lack of physical consistency.

2.3 Approach

In AMPT, we propose to address the above deficiencies of SPPT as follows.

2.4 Treatment of specific model fields

In the atmosphere, we experimented with perturbations of the 3D fields of temperature T, pressure p, horizontal wind components u, v, and specific humidity q_v. We also tried perturbing cloud ice and cloud water but found that those perturbations had little overall impact, so we abandoned them. In the soil, we perturbed 3D fields of soil temperature and soil moisture.

2.5 Stochastic pattern generator (SPG)

In this study, we rely on the limited-area Stochastic Pattern Generator (SPG) developed by (Tsyrulnikov and Gayfulin, 2017) to generate independent four-dimensional spatio-temporal pseudo-random fields ξ_i needed by AMPT, see Eq. (10). Each ξ_i is computed by solving the stochastic pseudo-differential equation

Here x, y, z are the SPG-space spatial coordinates, ${\nabla }^2$ is the 3D Laplacian, α(t, x, y, z) is the Gaussian white noise, and λ, U, σ are the parameters. λ determines the spatial scale of ξ. Given λ, the characteristic velocity U determines the time scale of ξ. Given λ and U, the parameter σ determines the variance of ξ and is selected to ensure that sd ξ = 1. The computational domain is the cube of size 2π and periodic boundary conditions in all three dimensions.

The design of Eq. (20) (note the third order of the equation in time and the square root of the negated and shifted Laplacian) is dictated by two requirements:

1. The solution ξ satisfies the so-called proportionality of scales property (Tsyroulnikov, 2001). The meaning of this property is the following. For any t, let us expand ξ in Fourier series in space: $\xi (t,\hspace{0.17em}\hspace{0.17em}x,\hspace{0.17em}\hspace{0.17em}y,\hspace{0.17em}\hspace{0.17em}z)=\Sigma {\tilde{\xi }}_{mn\ell }(t){\text{e}}^{\text{i}(mx+ny+\ell z)}$, where $m,n,\mathcal{l}$ are the spatial wavenumbers and ${\tilde{\xi }}_{mn\ell }(t)$ are the Fourier coefficients. Then Eq. (20) decouples into a series of ordinary stochastic differential equations for the random processes ${\tilde{\xi }}_{mn\ell }(t)$. The proportionality of scales means that for a large total wavenumber $K=\sqrt{m^2+n^2+{\ell }^2}$, the time scale of the process ${\tilde{\xi }}_{mn\ell }(t)$ is proportional to its spatial scale $1/\mathcal{K}$, a property often possessed by natural spatio-temporal processes.

2. The spatial spectra of ξ should be convergent, that is, the process variance should be finite in the space-continuous case and thus should be bounded above as the spatial resolution increases in the space-discrete case.

The solution to Eq. (20) is a zero-mean and unit-variance homogeneous (stationary in time and space) 4D random field, which has ‘nice’ spatio-temporal correlations: the shape of the correlation function is the same along any direction in the 4D space (thus, including spatial and temporal correlations). This universal correlation function belongs to the very popular in spatial statistics Matérn class, see (Tsyrulnikov and Gayfulin, 2017) for details.

Since the SPG computational domain is the cube whereas the limited-area-model domain can be approximated by a rectangular parallelepiped, we use the anisotropic Laplacian $\nabla {\hspace{0.17em}}^2$ instead of ${\nabla }^2$ in Eq. (20):

Here γ is the aspect ratio of the model domain in the horizontal and δ controls the length scale along the z axis.

The SPG internal parameters (for the model variable labeled by i) λ_i and δ_i are computed from the respective horizontal and vertical length scales, L_i and H_i, which are SPG external parameters (this is done by taking advantage of the known correlation functions of the 4D SPG field in space and time). The time scale is ${\mathcal{T}}_i=L_i/U$, where U is the above characteristic velocity (another external parameter).

With each realization of the driving white noise α, the solution to Eq. (20) is computed in 3D Fourier space using the backward Euler finite difference scheme in time. Once per ${\mathcal{T}}^{\text{FFT}}$ model time, the solution is converted to physical space (using the 3D inverse fast Fourier transform, FFT). The resulting physical-space random pattern ξ(t,x,y,z) is then mapped onto the model domain and used to compute model perturbations (according to Eq. (10)) for the next ${\mathcal{T}}^{\text{FFT}}$ model time interval. Avoiding the application of the inverse Fourier transform every time step significantly reduces the computational cost of the numerical scheme. ${\mathcal{T}}^{\text{FFT}}$ is selected to be a few times less than the time scale of the process, ${\mathcal{T}}_i$.

SPG can also generate 3D fields ξ(t, x, y) and ξ(x, y, z) as well as 2D fields ξ(x, y).

2.6 Mapping fields from SPG grid to model grid

In each of the three spatial dimensions, to avoid significant unphysical correlations in ξ_i between the opposite sides of the model domain due to the periodicity of the SPG domain, a segment of length 2L_i in the SPG domain is, first, discarded. Then, the rest of the SPG domain (the ‘working domain’) is mapped onto the model domain, where the field is interpolated to the model grid points. In the horizontal, the mapping is piecewise linear. In the vertical, two options were explored.

With Option 1, the mapping z ⟼ ζ (where ζ is the model’s vertical coordinate) is piecewise linear (just as in the horizontal) so that in the model space, the resulting field is stationary (homogeneous) in the vertical as a function of ζ.

With Option 2, a non-stationarity (inhomogeneity) in the vertical is introduced. Specifically, we assume that the model vertical levels are unevenly positioned by the model designer to account for the variable vertical length scale of the model fields. Say, in the planetary boundary layer, the model levels are dense, reflecting shallower meteorological structures (and thus shorter vertical length scales) than in the troposphere and stratosphere, where the vertical grid is much coarser because the vertical length scales are greater.

To obtain the same vertical inhomogeneity in AMPT perturbations interpolated to the model grid, we employ a mapping that is linear from the SPG coordinate z to the model vertical (continuous) ‘computational’ coordinate defined to be equal to $\mathcal{l}$ at the model level $\mathcal{l}$. With this mapping, the correlations between all adjacent model grid points in the vertical are the same (due to stationarity of the SPG field as a function of z). As a result, in model space, the vertical correlation falls to the same value at short vertical distances in the boundary layer (because the model grid is dense there) and at longer vertical distances above the boundary layer (where the model grid is coarser) thus implying the desired inhomogeneity in the vertical.

2.7 Properties of AMPT

4.1 Role of supersaturation limiter

The question was how the supersaturation limiter impacts the spread of the ensemble and the bias of the ensemble mean forecast.

First, we found that the impact on the near-surface temperature was rather small, so we focused on the impact on precipitation. Figure 2 shows the bias-spread scatterplots for two configurations of AMPT: without and with the supersaturation limiter. One can see that the bias-to-spread ratio was almost the same for the two configurations whilst the spread was larger without the supersaturation limiter (a qualitatively similar result was obtained in the setting without humidity AMPT perturbations, not shown). Therefore, we decided that the supersaturation limiter is not worth being activated in AMPT. So, for the experiments described in the rest of the article, the supersaturation limiter was off in both SPPT and AMPT.

Figure 2

4.2 Roles of humidity and temperature perturbations

Here we compare three configurations of AMPT and the basic configuration of SPPT in terms of their impacts on the bias and the spread of the ensemble. In the AMPT-TUVQv configuration, T, u, v, q_v fields were perturbed. In the AMPT-TUV configuration, only T, u, v fields were perturbed (that is, without humidity model perturbations). In the AMPT-UV configuration, only u,v fields were perturbed. With AMPT-UV, we found that the spread in the ensemble was somewhat too small. This led us to increase the amplitude multiplier κ in this configuration from its default value of 0.75 (see section 3) to 1.

Figure 3 shows bias-spread scatterplots for near-surface temperature. One can see that perturbing humidity on top of T, u, v did not change much neither the bias nor the spread (AMPT-TUV and AMPT-TUVQv performed similarly in these experiments). The other pair of schemes, SPPT and AMPT-UV, had significantly lower (in modulus) bias and spread. In dry conditions (the top left plot), SPPT had a smaller spread and bias-to-spread ratio than AMPT-UV. In wet conditions (the top right plot), SPPT and AMPT-UV performed similarly in terms of both spread and the bias-to-spread ratio, but the signs of the biases were different: negative for SPPT and positive for AMPT. In very wet conditions (the bottom plot), it was AMPT-UV that had the smallest bias-to-spread ratio while having nearly the same spread as SPPT.

Figure 3

Figure 4 shows bias-spread scatterplots for accumulated precipitation. In both wet (the left plot) and very wet (the right plot) conditions, the differences between the four schemes were more systematic than in Figure 3. Like in Figure 3, the magnitudes of precipitation bias and spread of SPPT and AMPT-UV were quite similar. Again, like for T_2m, AMPT-TUV and AMPT-TUVQv had significantly greater (than SPPT and AMPT-UV) bias and spread, but unlike Figure 3, the bias-to-spread ratio for AMPT-TUV was significantly lower than for AMPT-TUVQv.

Figure 4

It is also worth noting that the precipitation biases in Figure 4 were negative for SPPT and positive for AMPT. Investigating mechanisms that lead a complex nonlinear model to generate biases of different signs under different zero-mean stochastic model perturbations would be interesting but is beyond the scope of this research.

We summarize findings from Figures 3 and 4 as follows.

1. The magnitudes of forecast bias and spread due to SPPT were similar to those due to AMPT-UV (i.e., AMPT without temperature and humidity perturbations).

2. The addition of temperature perturbations in AMPT (in AMPT-TUV) led to significantly greater spread but somewhat higher (i.e. worse) bias-to-spread ratio.

3. Perturbing humidity on top of T, u, v (in AMPT-TUVQv) led to a significant growth in spread for precipitation but at the expense of a significant degradation in the bias-to-spread ratio.

Thus, without temperature perturbations in AMPT (i.e. perturbing only u, v), we could not outperform SPPT in terms of spread. However, with humidity perturbations in AMPT (i.e., perturbing T, u, v, q_v), forecast biases for precipitation were too large. Therefore, in the ensemble prediction experiments described in the next subsection, AMPT temperature (and winds) perturbations were switched on whereas humidity perturbations were switched off.

5.1 Setup

The general experimental setup described in section 3 applies to the experiments presented in this section with the caveat that there was no unperturbed member in the ensemble here. The results were verified against near-surface observations (about 40 stations) using the VERification System Unified Survey (VERSUS) developed within the COSMO consortium (Gofa et al., 2010).

The list of experiments and their basic features are presented in Table 1. In AMPT-SOIL, initial soil perturbations generated following section 2.4.4 were added to the respective fields of the initial-conditions ensemble members.

Besides the four model perturbation schemes listed in Table 1, we also tested a hybrid of AMPT and SPPT. We expected that such a hybrid might lead to an improvement because SPPT and AMPT perturbations explore different volumes in phase space. Indeed, the direction of the SPPT perturbation vector (comprised of all model variables at a grid point or a grid column) exactly coincides with the direction of the physical tendency vector, Eq. (6). In contrast, the AMPT perturbation vector can have any direction because ξ_i are uncorrelated, see Eq. (10) and section 2.3.3, so it has nothing to do with the direction of the model’s physical tendency vector at all. Besides, the magnitude of the SPPT perturbation is proportional to the magnitude of the local physical tendency, which is never the case with AMPT. Somewhat surprisingly, we could not see any benefit from the hybridization of the two schemes (not shown), so we abandoned the hybrid. As suggested by an anonymous reviewer, this outcome might be caused by the smaller magnitude of SPPT perturbations, which therefore could not make a difference.

During the two-month trial, there were two cases (out of 59) in which one of the forecasts exhibited catastrophic instability. On March 15, three out of the four tested configurations ‘exploded’: NOPERT (in which, we recall, there were no model perturbations, only initial and lateral-boundary conditions were perturbed), SPPT, and AMPT-NOSOIL, whereas the most successful AMPT configuration, AMPT-SOIL, performed normally. On March 20, both AMPT configurations failed whereas NOPERT and SPPT were OK. These two cases were excluded from the forecast performance statistics. So, our experiments showed that numerical instability can happen with AMPT, but its frequency is nearly as low as with SPPT. A longer test period is needed to make a more precise statement.

In this section, we show verification scores for T_2m forecasts averaged over the whole two-month period. Similar results were obtained for each of the two months separately (not shown). For near-surface wind speed V_10m, verification results were, mostly, similar to those for T_2m, except that soil perturbations had a significantly smaller impact on the wind speed than on temperature. This could be expected because perturbations of soil temperature and soil moisture directly impact surface fluxes of sensible and latent heat but not the flux of momentum. Due to the similarity of verification scores for V_10m compared to T_2m, we only display one typical plot for V_10m below. We also present a comparison of AMPT and SPPT for upper-level fields.

With precipitation, the only statistically significant effect of AMPT was an increase in ensemble spread. This can be seen in the model-perturbations-only experiments reported above in section 4 (see Figure 4) and this was observed in the experiments described in this section (not shown). Other verification scores for precipitation forecasts were ambiguous because, on the one hand, the time period was rather dry and on the other hand, the observation network was too scarce to detect relatively rare and localized precipitation events (not shown).

In contrast to section 4, where the results were obtained by averaging over large samples (the whole spatial grid and two months of data), the ensemble forecasts examined in this section were verified against a relatively scarce observation network. This caused a larger statistical uncertainty and therefore required statistical significance testing, see Appendix B for methodological details. The default settings of the statistical significance calculations were as follows. We examined statistical significance of improvements due to AMPT-NOSOIL relative to SPPT because both schemes involved no soil perturbations and thus were comparable. We averaged the scores over all lead times up to 48 h to reduce statistical uncertainty. The significance level was set at α = 0.05. The number of bootstrap samples was 10⁵. Along with the p-value of the above bootstrap test, p_bootstrap, we also computed the p-value of the classical one-sided Student’s t-test, p_Student.

5.2 Accuracy of ensemble mean forecasts and reliability of probabilistic forecasts

Figure 5(left) shows root-mean-square errors (RMSE) of the T_2m ensemble mean forecast (the upper bunch of curves). Each curve corresponds to a model perturbation scheme from the list in Table 1. One can see in Figure 5(left) that, excluding the initial transient (spinup) period of some 3 h (likely, due to imbalances in initial conditions), the RMSE had a prominent diurnal cycle with a broad minimum at night and a narrower maximum shortly after midday.

Figure 5

To highlight the barely seen in Figure 5(left) differences between the curves in the upper bunch, Figure 5(right) displays the normalized reduction in the T_2m ensemble-mean RMSE with respect to NOPERT, that is, (RMSE_NOPERT – RMSE)/RMSE_NOPERT. SPPT perturbations led to a slight overall reduction in RMSE. AMPT-NOSOIL perturbations gave rise to a more significant decrease in the T_2m RMSE most of the time except for the rather short time period in the afternoon local time when the ensemble-mean forecast deteriorated. The effect of soil perturbations on the ensemble mean forecast looks positive: the RMSE reduction is seen to be nearly uniformly higher in the experiment AMPT-SOIL than in AMPT-NOSOIL. However, all these differences were not statistically significant at the 0.05 level. For near-surface wind speed, the results were similar, not shown.

The lower bunch of curves in Figure 5(left) shows the ensemble spread in T_2m for the four schemes examined. One can see that the spread was somewhat too low with both SPPT and AMPT, though AMPT perturbations induced substantially greater spread than SPPT perturbations (this conclusion was valid for the spread in V_10m as well, not shown). The advantage of AMPT-NOSOIL over SPPT in terms of spread was highly statistically significant: the p-values of both the bootstrap test and the Student’s t-test were less than 0.001. For all 8 weeks and all lead times, the forecast spread with AMPT-NOSOIL was substantially greater than the spread with SPPT. The advantage of AMPT-SOIL over SPPT is seen to be even bigger.

Note that in an ideal ensemble, its members are indistinguishable from the truth. Both members of such an ensemble and the truth can be viewed as independent draws from the same probability distribution. In particular, the expectation of any member equals the expectation of the truth. The standard deviation of any member at some point in space and time (let it be denoted by σ) is the same as the respective standard deviation of the truth. Then, the expected square of the ensemble spread (defined as the square root of the unbiased sample variance) equals σ² by construction. At the same time, the expected square of the error in the ensemble mean (the mean-square error, MSE) equals $\frac{N+1}{N}{\sigma }^2$, where N is the ensemble size (Fortin et al., 2014). Being averaged in space and time, the expected square of the ensemble spread becomes equal to $\overline{{\sigma }^2}$. Correspondingly, the averaged in space and time MSE becomes $\frac{N+1}{N}\overline{{\sigma }^2}$. Thus, it is clear that for large ensembles, the root-mean-square spread should be very close to RMSE. For an ensemble of size N = 10 used in this study, RMSE should, on average, exceed the spread by $\sqrt{\frac{N+1}{N}}-1\approx 0.05$, which is of course much smaller than the differences between the RMSE and the spread seen in Figure 5(left).

Another reason why forecast RMSE can be greater than ensemble spread is observation error. In complex terrain, the dominant source of error can be representativeness uncertainty, which accounts for the fact that observations can poorly represent grid-cell averaged fields provided by the model. We did not account for the contribution of observation error to RMSE.

It is worth mentioning that the variations in the RMSE as functions of lead time seen in Figure 5(left) are not accompanied by corresponding variations in the spread. The reason is that the diurnal-cycle variations in RMSE were largely caused by systematic forecast errors (biases). Indeed, verification results of COSMO model forecasts on different domains presented in (Rieger et al., 2021) show that the magnitude of diurnal variations in the forecast bias is between 1K and 2K, with a broad maximum at night and a narrow minimum in the afternoon. Stochastic zero-mean model perturbations, cannot, in general, represent model biases, implying that the forecast spread cannot reflect the contribution of the bias to the RMSE. Either the forecasts are to be debiased or a multi-physics/multi-model ensemble (in which different members have different biases) is to be used, see, e.g., Berner et al. (2015).

Higher in the atmosphere, the effects of AMPT were similar to what we found for the near-surface fields. Figure 6 shows the RMSE and spread for upper-level temperature and wind speed (averaged over four lead times: 12, 24, 36, and 48 h). RMSE was computed here with respect to the ECMWF analyses interpolated to the model grid. It is seen that AMPT produced a higher spread and nearly the same RMSE as SPPT. The superiority of AMPT over SPPT in terms of spread was statistically significant for both fields at all levels shown in Figure 6, with p < 0.001 for both the bootstrap test and Student’s t-test. It is interesting to note that for temperature, the RMSE and spread were maximal near the ground (where turbulence and small-scale disturbances due to orography are the strongest). For wind speed, on the contrary, both the RMSE and spread were the smallest at the lowest model levels. Their downward decrease was most pronounced in the surface layer (perhaps, because the wind speed itself is relatively low near the ground due to turbulent friction).

Figure 6

If the spread is systematically and significantly different (say, lower, as in our experiments) than RMSE, then the conditional cumulative distribution function of the verification (truth, observation) x_v given the ensemble cumulative distribution function F_e(x) of the model variable x, i.e., P(x_v ⩽x | F_e) (where P stands for probability), is systematically different from F_e(x). This kind of inconsistency is known as lack of reliability of the probabilistic forecast (Toth et al., 2003). In these terms, Figures 5(left) and 6 demonstrate that AMPT perturbations are capable of significantly improving reliability of ensemble-based probabilistic forecasts compared to SPPT perturbations.

Figure 7 presents a more direct illustration of improvements in reliability brought about by AMPT-NOSOIL and AMPT-SOIL with respect to NOPERT and SPPT. It shows the reliability component of the Brier score² for the (most populated) event T_2m > 0°C. The improvement of AMPT-NOSOIL with respect to SPPT was statistically significant, p_bootstrap = 0.05 (though p_Student = 0.07 did not reach the significance level of 0.05).

Figure 7

5.3 Resolution and discrimination of probabilistic forecasts

A reliable uncalibrated ensemble-based probabilistic forecast tells the user that if the ensemble variance (i.e., the spread squared) at some point in space and time equals some number d, then the expected MSE (RMSE squared) also approximately equals d. This is an important property expected from an ensemble. However, reliability alone does not fully characterize the quality of a probabilistic forecast. Indeed, the constant spread equal to the ‘climatological’ RMSE would imply a perfectly reliable forecast, which is, however, not much useful. The ensemble needs to generate sufficiently variable spread. If it does and if variations in the spread correspond to variations in the RMSE (the property known as resolution), then the ensemble can provide relevant information about the spatially and temporally variable RMSE.

One useful measure of resolution combined with reliability is the continuous ranked probability score (CRPS). Figure 8 displays CRPS for T_2m. One can see that both SPPT and AMPT perturbations did improve CRPS (compared with NOPERT) while AMPT led to substantially greater improvements, especially with additional soil perturbations. The improvements of AMPT-NOSOIL compared with SPPT were statistically significant. The p-value of the bootstrap test computed following Appendix B was p_bootstrap = 0.03. The p-value of the Student’s t-test was p_Student = 0.02. AMPT-SOIL performed a fortiori statistically significantly better than AMPT-NOSOIL.

Figure 8

Figure 9 displays CRPS for near-surface wind speed V_10m. The impact of AMPT on the forecast performance in terms of V_10m was weaker than on T_2m but more stable, so the statistical significance of the advantage of AMPT-NOSOIL vs SPPT was stronger: p_bootstrap < 0.001 and p_Student = 0.004. Soil perturbations contributed little to the forecast CRPS. We mentioned a reason for that in section 5.1. Why the impact of AMPT on wind was weaker (but more stable) than the impact on temperature in the settings without soil perturbations remains unclear.

Figure 9

Another popular measure of reliability and resolution is the Brier score, which, in contrast to CRPS, is computed for specific events. If the event is defined as x < θ (or x > θ), where x is the continuous model variable in question (say, temperature) and θ is a threshold, then CRPS is the integral of the Brier score over all possible θ (Hersbach, 2000). Therefore, Figure 8 implies that the Brier score integrated over all temperature thresholds is better with AMPT than with SPPT. To find out what happens for specific distinct thresholds, we selected several meteorologically relevant events ($T_{\text{2m}}\text{ <}-5^{°}C,T_{\text{2m}}\text{ >}{0}^{°}C,T_{\text{2m}}\text{ >}{5}^{°}C$) and verified the ensemble forecasts using the Brier score. AMPT did outperform SPPT for all these events, we show the results for the most populated event T_2m > 0°C in Figure 10 (for near-surface wind speed, the results were similar, not shown). The improvements of AMPT-NOSOIL (and AMPT-SOIL as well) compared with SPPT were again statistically significant: p_bootstrap = 0.03 and p_Student = 0.01. Similar results were obtained for near-surface wind speed, not shown.

Figure 10

We also examined the resolution component of the Brier Score. With the reliability component substantially improved, see Figure 7, the resolution component was barely changed or slightly reduced (i.e., degraded, not shown). As noted by Candille and Talagrand (2005), an increase in spread may lead to a degradation in resolution because a larger spread is more akin to the climatological spread, which yields no resolution. So, we state that AMPT substantially improved reliability without significantly degrading resolution.

As discussed above, the ability of a probabilistic forecast to predict the uncertainty in the forecast can be measured by reliability and resolution through the conditional distribution of the verification x_v given the ensemble cumulative distribution function F_e(x), i.e., p(x_v | F_e) (where p stands for probability density). The reciprocal and complementary view on the accuracy of a probabilistic forecast is through p(F_e | x_v), i.e., conditional on verification. The capability of a probabilistic forecast to concentrate the probability mass close to the observed x_v (thus, discriminating between different outcomes x_v) is called discrimination (e.g., Wilks, 2011).

Discrimination is commonly assessed using the Relative Operating Characteristic (ROC) defined for a specific meteorological event and, more succinctly, using the area under the ROC curve (ROC area). Figure 11 shows the ROC area for the event T_2m > 0°C (similar results were obtained for the two other events, T_2m < –5°C and T_2m > 5°C, not shown). From this figure, we see that AMPT did outperform SPPT. Our statistical significance tests showed that the advantage of AMPT-NOSOIL over SPPT in terms of ROC area was statistically significant: p_bootstrap = p_Student = 0.01. As AMPT-SOIL is seen in Figure 11 to be almost uniformly better than AMPT-NOSOIL, we conclude that AMPT-SOIL was statistically significantly more skilful for near-surface temperature than SPPT, too. For near-surface wind speed, the results were mixed (not shown).

Figure 11

A few remarks regarding Figure 11 are in order. First, the overall level of the ROC area score was quite high (its perfect value is 1 whereas a value less than 0.5 indicates no skill) so the probabilistic forecasts in question were quite skillful. Second, as with many other plots above, the ensemble forecasts were the worst at zero lead time. This was, most likely, caused by the poor quality of the initial ensemble. Finally, we note that the ROC area is known to be insensitive to the degree of reliability and cannot be improved by calibration of a probabilistic forecast (e.g., Wilks, 2011). Therefore, the superiority of AMPT-SOIL over SPPT and NOPERT in terms of ROC area implies that AMPT did not just inflate spread (which could be done by calibration), it can yield ‘good spread’ in the sense of Eckel and Mass (2005), meaning greater spread when and where it is relevant.

5.4 Summary of ensemble prediction results

For near-surface temperature forecasts, atmospheric AMPT perturbations gave rise to statistically significant improvements (compared to SPPT) in the performance of the ensemble prediction system — in terms of spread-skill relationship, CRPS, Brier score, and ROC area. In terms of RMSE of the ensemble-mean forecast of near-surface temperature, the effects of AMPT were mixed. For near-surface wind speed forecasts, the effects of the atmospheric AMPT were similar to but weaker than the effects on near-surface temperature forecasts. Soil AMPT perturbations imposed in addition to atmospheric AMPT perturbations led to nearly uniform further improvements in the ensemble performance for near-surface temperature and smaller improvements in near-surface wind speed. The impact of AMPT on precipitation forecast scores was mixed.

In the free atmosphere, AMPT consistently generated bigger spread than SPPT, without degrading RMSE. The positive effect of AMPT, though, decreased away from the surface. We recall that our focus in this study was on the near-surface fields, so we did nothing to increase the spread induced by AMPT in the free atmosphere.

A Additive and multiplicative perturbations

This Appendix contains remarks on the suitability of the term ‘additive’ in the acronym AMPT.

In stochastic systems theory, the term ‘additive noise’ denotes a perturbation that does not depend on the system state. If the perturbation is a function of the system state, it is called multiplicative, e.g., Fuchs (2014, sec. 6.5). We recall that in ensemble prediction, an ensemble member can be viewed as a numerical solution to the stochastic dynamic equation of the unknown, random truth (system state), e.g., Tsyrulnikov and Gayfulin (2017, Eq.(1.3)): dx/dt = F(x) – ξ (where x is the truth, F is the forecast-model operator, and ξ is the random model perturbation). So, the model perturbation is multiplicative if it is a function of x (the state of the ensemble member on which the model perturbation is imposed) and additive if it does not depend on x.

With SPPT, the model perturbation at some point in space and time, (s, t), is proportional to the physical tendency P(x(s, t)), where x is the state of the ensemble member, Eq. (6). Therefore, the SPPT perturbation in each model variable i is a function of x and thus is multiplicative:

With AMPT, if the scaling tendency is taken from the control (unperturbed) run, so that $P_i(\text{s},\hspace{0.17em}\hspace{0.17em}t)=P_i({\text{x}}^{\text{control}}(\text{s},\hspace{0.17em}\hspace{0.17em}t))$, the perturbation ceases to depend on x:

it only depends on the known initial conditions of the control run. Thus, the AMPT perturbation is indeed additive in the setting in which the scaling physical tendency is taken from the control model run.

B Methodology of testing statistical significance

Here we describe how we calculated statistical significance of improvements due to AMPT in ensemble prediction experiments presented in section 5.

Let us test the hypothesis that a verification score ψ is greater for scheme 1 than for scheme 2. More precisely, we are going to compare the scores ψ₁ and ψ₂ that could be obtained from a very large sample of forecasts, whilst having their unbiased estimates ${\widehat{\psi }}_1$ and ${\widehat{\psi }}_2$ obtained from smaller samples. That is, the hypothesis to be tested is $H:{\psi }_1\equiv �{\widehat{\psi }}_1>�{\widehat{\psi }}_2\equiv {\psi }_2$. The default ‘null’ hypothesis states that there is no difference between the two schemes: H₀: ψ₁ = ψ₂.

Since atmospheric states are correlated in time and forecast errors depend on the atmospheric flow, forecast errors are also temporally correlated (dependent). Accounting for these time correlations would complicate statistical hypothesis testing, so for simplicity, we employ the following approach. We divide the whole two-month period into n=8 weeks and compute the score ψ for each of the two schemes for each week separately, getting the respective sequences ${\widehat{\psi }}_{1j}$ and ${\widehat{\psi }}_{2j}$ (for j = 1, 2, …, n) and their differences ${\delta }_j^{\text{obs}}={\widehat{\psi }}_{1j}-{\widehat{\psi }}_{2j}$. We treat ${\delta }_j^{\text{obs}}$ as (observed) realizations of the respective random variables δ_j. Assuming that temporal dependencies in δ_j decay on a time scale shorter than one week, we regard δ_j as a sample of statistically independent and identically distributed random numbers.

Informally, if ${\delta }_j^{\text{obs}}$ are mostly positive, we tend to believe that $\mathbb{E}{\delta }_j>0$, meaning that the hypothesis H we are testing is true. Formally, we reject H₀ (and thus accept H) if $p=P(\overline{\delta }\ge {\overline{\delta }}^{\text{obs}}|H_0)$, where the overbar denotes sample averaging, is sufficiently small. This means that if, under H₀, the probability of the observed (and more extreme) deviation of $\overline{\delta }$ away from zero is small, then the data we have are, likely, not consistent with H₀. Specifically, we select a significance level, α, and reject H₀ if the above probability p (known as the p-value) is less than α.

To compute p we need to know the probability distribution of δ_j under H₀. This can be easily done using an approach to statistical inference known as bootstrap. In its simplest flavor, bootstrap postulates that the data distribution (i.e., the distribution of δ_j) equals the empirical distribution (which is concentrated at the observed data, in our case, ${\delta }_j^{\text{obs}}$). However, this setup cannot be directly applied here because H₀ assumes that $�{\delta }_j=0$ whereas the empirical distribution is, likely, biased (${\overline{\delta }}^{\text{obs}}\ne 0$). Following Efron and Tibshirani (1994, ch.16), to define the bootstrap distribution, we shift the empirical distribution so that it has zero mean. Then, we proceed as usual, sampling (with replacement) from this discrete distribution (concentrated with equal probabilities at the points ${\Delta }_j^{\text{obs}}={\delta }_j^{\text{obs}}-{\overline{\delta }}^{\text{obs}}$), compute the mean in each bootstrap sample, $\overline{\delta }$, and calculate the fraction of bootstrap samples in which $\overline{\delta }>{\overline{\delta }}^{\text{obs}}$. This fraction is the estimate of the p-value.