# Additive Model Perturbations Scaled by Physical Tendencies for Use in
Ensemble Prediction111 This is the authors’ final version of the paper
published in Tellus A: Dynamic Meteorology and Oceanography. DOI:
https://doi.org/10.16993/tellusa.3224

Michael Tsyrulnikov(1)222 Corresponding author<EMAIL_ADDRESS>ORCID 0000-0002-7357-334X  Elena Astakhova(1) and Dmitry Gayfulin(1)
(1)HydroMetCenter of Russia, Moscow, Russia

###### Abstract

Imperfections and uncertainties in forecast models are often represented in
ensemble prediction systems by stochastic perturbations of model equations. In
this article, we present a new technique to generate model perturbations. The
technique is termed Additive Model-uncertainty perturbations scaled by
Physical Tendencies (AMPT). The generated perturbations are independent
between different model variables and scaled by the local-area-averaged
modulus of physical tendency. The previously developed Stochastic Pattern
Generator is used to generate space and time-correlated pseudo-random fields.
AMPT attempts to address some weak points of the popular model perturbation
scheme known as Stochastically Perturbed Parametrization Tendencies (SPPT).
Specifically, AMPT can produce non-zero perturbations even at grid points
where the physical tendency is zero and avoids perfect correlations in the
perturbation fields in the vertical and between different variables. Due to a
non-local link from physical tendency to the local perturbation magnitude,
AMPT can generate significantly greater perturbations than SPPT without
causing instabilities. Relationships between the bias and the spread caused by
AMPT and SPPT were studied in an ensemble of forecasts. The non-hydrostatic,
convection-permitting forecast model COSMO was used. In ensemble prediction
experiments, AMPT perturbations led to statistically significant improvements
(compared to SPPT) in probabilistic performance scores such as spread-skill
relationship, CRPS, Brier Score, and ROC area for near-surface temperature.
AMPT had similar but weaker effects on near-surface wind speed and mixed
effects on precipitation.

[figure]style=plain,subcapbesideposition=top

Keywords:  Model uncertainty, SPPT, ensemble prediction, stochastic
perturbations.

##  1 Introduction

Forecasting natural phenomena such as weather cannot be perfect. Knowing the
degree of the imperfection (the forecast uncertainty) is always desirable and
sometimes vital, for example, in decision-making. A rough estimate of forecast
uncertainty can be obtained by comparing forecasts with verifying observations
and averaging the forecast-minus-observation statistics over time and space.
This simple approach often results in useful estimates. However, the
‘climatological’ estimates can be insufficient if we predict the state of a
nonlinear chaotic system like the Earth’s atmosphere, where the forecast
uncertainty can vary depending on the weather situation (on the local
structure of the atmospheric flow) and the observational coverage.

To allow for a situation-dependent assessment of the forecast uncertainty, a
forecast of forecast uncertainty is needed. In probabilistic terms, we have to
predict not just the state of the system (in our case, the atmosphere) but
also the probability distribution of the unknown true atmospheric state around
the forecast. A Monte-Carlo-based approach known as ensemble prediction, in
which the probability distribution of the truth is represented by a small
number (tens to hundreds) of points in state space (ensemble members), has
proven to be both feasible and useful in predicting major features of forecast
uncertainty, see, e.g., Leutbecher & Palmer (2008), Wilks (2011) and
references therein.

For ensemble prediction to be successful in predicting the forecast
uncertainty, ensemble members need to be pseudo-random draws from a
probability distribution that is reasonably close to the conditional
probability distribution of the truth given the data available prior to the
forecast. The most promising approach here is, arguably, to identify all
sources of uncertainty that affect the forecast and then stochastically model
each of those ‘input’ uncertainties individually. The forecast uncertainty is
caused by uncertainties in (i) meteorological observations, (ii) data
assimilation techniques, (iii) boundary conditions, and (iv) the forecast
model itself. In this study we are concerned with the latter source, the model
uncertainty.

###  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. 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. 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. 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
${\bf x}_{\rm SGS}$ with itself and with the grid-scale field ${\bf x}_{\rm
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
$\boldsymbol{\pi}({\bf x}_{\rm SGS},{\bf x}_{\rm GS})$. Deterministic physical
parametrization schemes, having, by definition, access only to ${\bf x}_{\rm
GS}$, attempt to assess this impact, producing a physical tendency, ${\bf
P}({\bf x}_{\rm GS})$. The uncertainty in all physical parametrization schemes
combined is the difference between ${\bf P}({\bf x}_{\rm GS})$ and the true
physical tendency $\boldsymbol{\pi}({\bf x}_{\rm SGS},{\bf x}_{\rm GS})$.

The subgrid-scale field ${\bf x}_{\rm SGS}$ is not explicitly defined in the
model, therefore it is unknown and even unknowable to the model (to the extent
that ${\bf x}_{\rm SGS}$ is not determined by ${\bf x}_{\rm GS}$). Therefore
we assume that ${\bf x}_{\rm SGS}$ is a random field with some conditional
probability density $p({\bf x}_{\rm SGS}\mid{\bf x}_{\rm GS})$ (we condition
on ${\bf x}_{\rm GS}$ because it is available to the model). This density
induces the probability density $p(\boldsymbol{\pi}\mid{\bf x}_{\rm 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 & Pallarès (2014), ${\bf P}_{\rm
ideal}^{\rm determ}({\bf x}_{\rm
GS})=\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf x}_{\rm GS})$, where
$\operatorname{\mathbb{E}}$ stands for expectation. It is ${\bf P}_{\rm
ideal}^{\rm determ}({\bf x}_{\rm GS})$ that a deterministic physical
parametrization seeks to approximate by ${\bf P}({\bf x}_{\rm GS})$. In these
terms, it is meaningful to call

$\varepsilon_{\bf P}={\bf P}-{\bf P}_{\rm ideal}^{\rm determ}={\bf
P}-\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf x}_{\rm GS})$ (1)

the error in the deterministic physical tendency ${\bf P}$. However, the
deterministic physical tendency is incapable of simulating the variability of
$\boldsymbol{\pi}$ around its conditional mean value
$\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf x}_{\rm GS})$,

$\delta_{\bf
P}=\boldsymbol{\pi}-\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf x}_{\rm
GS})$ (2)

(Palmer 2012). This is the inherent (irreducible) uncertainty in
$\boldsymbol{\pi}({\bf x}_{\rm SGS},{\bf x}_{\rm GS})$ given ${\bf x}_{\rm
GS}$. Then, the full uncertainty is the difference between the error
$\varepsilon_{\bf P}$ and the irreducible uncertainty $\delta_{\bf P}$:

$\boldsymbol{\varepsilon}={\bf P}-\boldsymbol{\pi}\equiv[{\bf
P}-\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf x}_{\rm
GS})]-[\boldsymbol{\pi}-\operatorname{\mathbb{E}}(\boldsymbol{\pi}\mid{\bf
x}_{\rm GS})]=\boldsymbol{\varepsilon}_{\bf P}-\boldsymbol{\delta}_{\bf P}.$
(3)

The irreducible uncertainty $\delta_{\bf P}$ is especially important in so-
called ‘gray zones’, where a process is partly resolved by the model, meaning
that a length scale of ${\bf x}_{\rm 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, $\nu$, of
random contributions:
$\boldsymbol{\pi}=\sum_{\ell=1}^{\nu}\boldsymbol{\pi}_{\ell}$. With
convection, the contributions $\boldsymbol{\pi}_{\ell}$ to the convective
tendency are due to individual convective plumes (Plant & Craig 2008). With
boundary-layer turbulence, $\boldsymbol{\pi}_{\ell}$ are due to individual
turbulent eddies (Kober & Craig 2016). Assuming, for presentation purposes,
that all $\boldsymbol{\pi}_{\ell}$ are the same in a grid cell and their
number $\nu$ follows the Poisson distribution, we readily obtain that
$|\operatorname{\mathbb{E}}\boldsymbol{\pi}|/{\mathsf{sd}}\,\boldsymbol{\pi}=\sqrt{\operatorname{\mathbb{E}}\nu}$,
where ${\mathsf{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 $\operatorname{\mathbb{E}}\nu\sim 1$, which is the case in a gray
zone as discussed above. This makes the irreducible uncertainty in
$\boldsymbol{\pi}$ indeed substantial and implies that it needs to be properly
accounted for in ensemble prediction.

It is worth remarking that the irreducible uncertainty $\delta_{\bf P}$ caused
by the randomness of $\boldsymbol{\pi}$ given ${\bf x}_{\rm GS}$ is the
aleatory (truly random) uncertainty, therefore it cannot be called error.
Whereas $\varepsilon_{\bf 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 ${\bf P}_{\rm ideal}^{*}({\bf
x}_{\rm GS})$ is a ‘possible true $\boldsymbol{\pi}$ consistent with the grid-
scale field ${\bf x}_{\rm GS}$’, that is, a random draw from the probability
distribution $p(\boldsymbol{\pi}\mid{\bf x}_{\rm 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 & Craig (2008), Dorrestijn et al. (2013), Sakradzija et
al. (2015), Hirt et al. (2019), Machulskaya & 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 Methodology

In this section, we review SPPT and introduce a new approach to generation of
model-uncertainty perturbations. The new AMPT scheme builds on SPPT and
attempts to avoid/relax some of the deficiencies of SPPT. AMPT is applied to
both atmosphere and soil.

###  2.1 Background on SPPT and notation

To facilitate the presentation of AMPT, we first outline SPPT (Leutbecher et
al. 2017). By ${\bf P}(x,y,\zeta,t)=(P_{1}(x,y,\zeta,t),\dots,P_{n_{\rm
fields}}(x,y,\zeta,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 $\zeta$, and the forecast time $t$. Here
$P_{i}(x,y,\zeta,t)$ is the component of ${\bf P}(x,y,\zeta,t)$ in the $i$-th
model field (variable) $X_{i}$, and $n_{\rm 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 ${\bf P}^{*}$ is postulated to be

${\bf P}^{*}(x,y,\zeta,t)=(1+\kappa\,\xi(x,y,t))\,{\bf P}(x,y,\zeta,t),$ (4)

where $\xi(x,y,t)$ is the zero-mean and unit-variance spatio-temporal random
field and $\kappa$ 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
${\kappa}\,\xi(x,y,t)$ is limited to the segment $[-1,1]$ so that

$|\xi(x,y,t)|<{\kappa}^{-1}.$ (5)

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:

$\Delta{\bf P}={\bf P}^{*}(x,y,\zeta,t)-{\bf
P}(x,y,\zeta,t)=\kappa\,\xi(x,y,t)\,{\bf P}(x,y,\zeta,t).$ (6)

Here and elsewhere $\Delta$ denotes a perturbation.

###  2.2 Motivation

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

  1. 1.

In SPPT, perturbations are large (small) when and where the physical
parametrizations generate a large (small) physical tendency ${\bf 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. 2.

Equation (6) implies that the multivariate perturbation vector $\Delta{\bf
P}(x,y,\zeta,t)$ is strictly proportional to the physical tendency vector
${\bf P}(x,y,\zeta,t)$. As noted by Leutbecher et al. (2017), this implies
that SPPT tacitly assumes that only the magnitude of the vector ${\bf 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,\zeta,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 $\Delta P_{i}=\kappa\,\xi
P_{i}$, where $\kappa$ is a positive constant and
$\operatorname{\mathbb{E}}\xi=0$. This implies that
$\operatorname{\mathbb{E}}(\Delta P_{i}\mid P_{i})=0$, and
${\mathsf{Cov}}\,(\Delta P_{i},\Delta P_{j}\mid
P_{i},P_{j})=\operatorname{\mathbb{E}}(\Delta P_{i}\Delta P_{j}\mid
P_{i},P_{j})=\kappa^{2}({\mathsf{sd}}\,\xi)^{2}P_{i}P_{j}$ (where
${\mathsf{Cov}}\,$ stands for covariance). Taking into account that
${\mathsf{sd}}\,(\Delta P_{i}\mid P_{i})=\kappa\,|P_{i}|\,{\mathsf{sd}}\,\xi$
and ${\mathsf{sd}}\,(\Delta P_{j}\mid
P_{j})=\kappa\,|P_{j}|\,{\mathsf{sd}}\,\xi$, the correlation, that is, the
covariance normalized by the product of the two standard deviations, becomes
${\mathsf{Corr}}\,(\Delta P_{i},\Delta P_{j}\mid P_{i},P_{j})=\pm 1$.

  3. 3.

Similarly (and also noted by Leutbecher et al. 2017), since the SPPT random
pattern $\xi(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. 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 $\pm 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. 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.3.1 Univariate AMPT design

We rely on the SPPT’s assumption that the standard deviation of the model-
uncertainty perturbation $\Delta X_{i}=\Delta P_{i}$ is proportional to the
modulus of its respective physical tendency, ${\mathsf{sd}}\,(\Delta X_{i}\mid
P_{i})=\kappa|P_{i}|$. In AMPT, we define this dependency to be more general
than just point-wise, postulating ${\mathsf{sd}}\,(\Delta X_{i}\mid P_{i})$ to
be proportional to an area-averaged $|P_{i}|$. This allows AMPT to generate
non-zero perturbations even at grid points with zero physical tendency — if
there are nearby points with non-zero physical tendency. Theoretically, this
approach can be justified as follows.

Consider the unknown true model-uncertainty field $\varepsilon_{i}({\bf r})$
(where ${\bf r}$ is the spatio-temporal coordinate vector $(x,y,\zeta,t)$).
Assume that $\varepsilon_{i}({\bf r})$ can be modeled as a random field with
zero mean and the unknown spatially variable standard deviation
$\sigma_{i}({\bf r})$. In these settings, SPPT, effectively, estimates
$\sigma_{i}({\bf r})$ as $\kappa|P_{i}({\bf r})|$, see Eq.(6), that is, the
only predictor to estimate $\sigma_{i}$ at the point ${\bf r}$ in space and
time is the physical tendency at the same point. We propose to acknowledge
that $|P_{i}({\bf r})|$ is a noisy ‘observation’ of the true model-uncertainty
standard deviation $\sigma_{i}({\bf r})$ and therefore it is worth looking for
other predictors, i.e., other data that can contain information on
$\sigma_{i}({\bf r})$. In AMPT, we hypothesize that these additional relevant
data are values of the absolute physical tendency in the vicinity of the
spatio-temporal point in question ${\bf r}$. Linearly combining these noisy
data, we obtain the AMPT’s estimate of the unknown $\sigma_{i}({\bf r})$ given
the known field $|P_{i}(.)|$:

$\widehat{\sigma}_{i}({\bf r})=\int W_{i}({\bf r,r}^{\prime})\,|P_{i}({\bf
r}^{\prime})|\,{\rm d}{\bf r}^{\prime},$ (7)

where $W_{i}({\bf r,r}^{\prime})$ is a weighting function that determines the
contribution of the absolute physical tendency, $|P_{i}|$, evaluated at the
grid point ${\bf r}^{\prime}$ to the estimate of the model-error standard
deviation, $\widehat{\sigma}_{i}$, evaluated at the grid point ${\bf r}$ and
the integral is over the model domain (and, possibly, over model time as
well). With the simplifying assumption that the weighting function is
homogeneous, $W_{i}({\bf r,r}^{\prime})=w_{i}({\bf r-r}^{\prime})$, we rewrite
Eq.(7) as

$\widehat{\sigma}_{i}({\bf r})=\kappa_{i}{\cal P}_{i}({\bf r}),$ (8)

where $\kappa_{i}=\int w_{i}({\bf r}){\rm d}{\bf r}$,

${\cal P}_{i}({\bf r})=\int K_{i}({\bf r-r}^{\prime})\;|P_{i}({\bf
r}^{\prime})|\,{\rm d}{\bf r}^{\prime}$ (9)

is the local-area averaged absolute physical tendency (we will call it the
scaling physical tendency), and $K_{i}({\bf r})=\frac{1}{\kappa_{i}}w_{i}({\bf
r})$ is the averaging kernel such that $\int K_{i}({\bf r})\,{\rm d}{\bf
r}=1$.

Having estimated $\sigma_{i}({\bf r})$, we can simulate the model uncertainty
as

$\Delta X_{i}({\bf r})=\widehat{\sigma}_{i}({\bf r})\,\xi_{i}({\bf
r})=\kappa_{i}\,{\cal P}_{i}({\bf r})\,\xi_{i}({\bf r}),$ (10)

where $\xi_{i}({\bf r})$ is a zero-mean and unit-variance Gaussian random
field postulated to be stationary in space and time. Its spatio-temporal
correlations are discussed below in section 2.3.4. Spatial (and temporal) non-
stationarity of the model-uncertainty field $\varepsilon_{i}$ defined by
Eq.(10) comes from variability in the scaling physical tendency ${\cal
P}_{i}({\bf r})$. Non-Gaussianity of the model-uncertainty field comes from
the randomness of ${\cal P}_{i}({\bf r})$.

####  2.3.2 Scaling physical tendency

Technically, with the gridded fields, the integral in Eq.(9) is replaced with
a sum. Given the layered structure of the atmosphere, we perform the averaging
in the horizontal only:

${\cal
P}_{i}(x,y,\zeta,t)=\sum_{q}\sum_{r}w_{qr}\left|P_{i}(x_{q},y_{r},\zeta,t)\right|,$
(11)

where $(x_{q},y_{r})$ is the horizontal grid point, $i$ (we recall) labels the
model variable, and $w_{qr}$ are the averaging weights. The latter are
specified to be non-zero only inside the averaging area
$|x_{q}-x|<A_{i},|y_{q}-y|<A_{i}$, where $A_{i}$ is the half-size of the
averaging area (the averaging length scale) in both $x$ and $y$ directions.
For simplicity and due to a lack of knowledge on the spatial structure of
uncertainties associated with $\left|P_{i}(x_{q},y_{r},\zeta,t)\right|$ as
‘observations’ on $\sigma_{i}(x,y,\zeta,t)$ we adopt the simplest design: the
weights $w_{qr}$ are equal to each other within the averaging area and
normalized so that for any $x,y$, we have $\sum_{qr}w_{qr}=1$. In the context
of limited area modeling, if the point $(x,y)$, where ${\cal P}_{i}$ is
evaluated, is near the model’s boundary, $A_{i}$ is reduced so that the
averaging area is still a square and is within the model domain.

If $A_{i}$ is greater than the size of the domain, the averaging is performed
over the whole domain so that ${\cal P}_{i}(x,y,\zeta,t)={\cal
P}_{i}(\zeta,t)$ is the same for all grid points in the horizontal. This is
how we computed the scaling physical tendency for atmospheric temperature,
winds, and soil temperature. For less Gaussian fields such as humidity and
soil moisture, ${\cal P}_{i}$ is computed by averaging over a significantly
smaller moving window centered at the grid point in question (see section 3
for details). It is worth noting that dependence of ${\cal P}_{i}(\zeta,t)$ on
model time $t$ is essential. The reason is that during the forecast, the mean
magnitude of physical tendency may undergo significant variations due to a
passage of a front or a convective system, changes in convection, etc.

####  2.3.3 Multivariate and 3D aspects

First, we allow in AMPT not only for errors in the modulus of the vector ${\bf
P}$ but also for errors in the direction of ${\bf P}$. We do so by introducing
independent driving random fields $\xi_{i}$ for different model variables
$X_{i}$ (see Eq.(10)). Though purely uncorrelated perturbations are
unphysical, we rely on the forecast model to introduce physically meaningful
relationships between the variables during its adaptation to the model
perturbations. Since the magnitude of the perturbation at each model time step
is small compared with the natural variability333In the experiments described
below, the mean absolute magnitudes of AMPT perturbations per time step were
as follows (at 8 a.m. local time). Temperature: somewhat less than 0.001 K in
the lower 3-km tropospheric layer and much less above it. A horizontal wind
component: from 0.015 m/s near the ground to 0.0001 m/s at 3 km height and
even less above 3 km., the adaptation is expected to go smoothly.

Second, to get rid of the perfect coherency of the perturbations in the
vertical, we switch from the 3D random pattern $\xi(x,y,t)$ in SPPT to the 4D
random fields $\xi_{i}(x,y,\zeta,t)$ in AMPT.

Third, the perturbation-magnitude multiplier $\kappa_{i}$ is variable-specific
in AMPT.

####  2.3.4 Length scales and perturbation magnitudes

Let us assess the horizontal length scale of the SPPT perturbation field, see
Eq.(6). To this end, with fixed $\zeta$ and $t$, suppose that the physical
tendency $P_{i}(x,y,\zeta,t)$ is a stationary (homogeneous) and isotropic
random field as a function of $x,y$. Let its correlation function be denoted
as $C_{P}(r)$, where $r$ is the horizontal distance. Since the random pattern
$\xi(x,y,t)$ is stationary, isotropic, and independent of $P_{i}$ by
construction, the horizontal correlation function of the SPPT perturbation
field $\Delta^{\rm SPPT}$ is $C^{\rm SPPT}_{\Delta}(r)=C_{\xi}(r)\,C_{P}(r)$,
where $C_{\xi}(r)$ is the correlation function of $\xi$. Assuming, further,
that $C_{\xi}(r)$ and $C_{P}(r)$ are twice differentiable, we obtain the
differential length scale (e.g., Monin & Yaglom 2013) $L^{\rm SPPT}_{\Delta}$
of the perturbation field as $(L^{\rm SPPT}_{\Delta})^{-2}=(C^{\rm
SPPT}_{\Delta})^{\prime\prime}(0)$. Since $C_{\xi}(0)=C_{P}(0)=1$ and
$C_{\xi}^{\prime}(0)=C_{P}^{\prime}(0)=0$, we obtain

$(L^{\rm SPPT}_{\Delta})^{-2}=(L^{\rm SPPT}_{\xi})^{-2}+L_{P}^{-2},$ (12)

where $L^{\rm SPPT}_{\xi}$ and $L_{P}$ are the differential length scales of
the random pattern $\xi$ and the physical tendency $P$, respectively.

As we noted in item 4 of section 2.2, the optimally tuned length scale of the
SPPT random pattern in convective-scale applications appears to be much
greater than the respective scales of natural variability. In terms of the
differential length scales, we may therefore conclude that $L^{\rm
SPPT}_{\xi}\gg L_{P}$. Then, Eq.(12) implies that $L^{\rm
SPPT}_{\Delta}\approx L_{P}$. Since it is plausible that the optimally tuned
$L^{\rm SPPT}_{\Delta}$ approximates the unknown model-uncertainty length
scale $L_{\varepsilon}$, we obtain

$L_{\varepsilon}\approx L_{P}.$ (13)

In AMPT, the random pattern $\xi$ multiplies the scaling physical tendency
${\cal P}$ so that

$(L^{\rm AMPT}_{\Delta})^{-2}=(L^{\rm AMPT}_{\xi})^{-2}+L_{\cal P}^{-2},$ (14)

where $L_{\cal P}$ is the differential length scale of ${\cal P}$. Since
${\cal P}$ is a spatially smoothed version of $|P|$, $L_{\cal P}\gg L_{P}$.
Therefore, for the length scale of the AMPT perturbation to be close to
$L_{P}$, we have to make $L^{\rm AMPT}_{\xi}\approx L_{P}$, that is, much
smaller than in SPPT. The optimally tuned length scales of the AMPT random
patterns are given in section 3 (they, indeed, appeared to be an order of
magnitude smaller than their SPPT counterparts).

As for the AMPT magnitude multipliers $\kappa_{i}$, the non-local dependence
of the scaling physical tendency on the modulus of the unperturbed physical
tendency (Eqs. (9) and (11)) suggests that $\kappa_{i}$ do not need to obey
the SPPT’s strict upper limit on the magnitude multiplier, Eq.(5). This allows
AMPT to cause greater spread than SPPT in the ensemble if needed.

####  2.3.5 Further details

The rest of the Methodology section is organized as follows. Details on AMPT
perturbations for specific model fields are given below in section 2.4. The
Stochastic Pattern Generator (SPG, Tsyrulnikov & Gayfulin 2017) is outlined in
section 2.5. Section 2.6 explains how SPG fields are mapped from the SPG
domain onto the model domain. In section 2.7 we briefly discuss stability,
conservation properties of the new scheme, possible biases due to nonlinearity
of the forecast model, and explain the terminology according to which we call
SPPT perturbations multiplicative and AMPT perturbations additive.

###  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_{\rm 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.4.1 Atmospheric temperature, pressure, and winds

Independent perturbations of $T,u,v$ are computed following Eq.(10). The
pressure perturbation $\Delta p$ is computed from $\Delta T$ by integrating
the hydrostatic equation (in which $\Delta q_{\rm v}$ is neglected as a small
contribution to a small perturbation) assuming zero pressure perturbation at
the top of the model domain.

As for wind perturbations, we note that theoretically, it is ‘better’ to rely
on mutually uncorrelated random stream function and velocity potential
perturbation fields — rather than on mutually uncorrelated $u$ and $v$ (i.e.,
zonal and meridional wind) perturbation fields. The reason is that the former
approach allows for isotropic vector-wind perturbations (Monin & Yaglom 2013,
section 12.3), unlike the latter approach. However, in practical terms, we
were not able to identify any practically significant flaw in the vector field
composed of two independent $u$ and $v$ perturbation fields. For this reason
and due to the lack of evidence on the actual structure of model
uncertainties, we stick to the simpler formulation of AMPT with mutually
independent $\Delta u(x,y,\zeta,t)$ and $\Delta v(x,y,\zeta,t)$ in this study.

####  2.4.2 Humidity

The salient difference of humidity $q_{\rm v}$ from $T,u,v$ is that $q_{\rm
v}$ has a narrow range of values (from zero to saturation or somewhat higher
than saturation). The range of $q_{\rm v}$ is narrow in the sense that it is
comparable to the standard deviation of the natural variability in $q_{\rm
v}$. To make sure that the AMPT-perturbed $q_{\rm v}$ is within this range
(i.e., from zero to saturation) and does not directly introduce any bias into
the model, we modify the above formulation of AMPT. Specifically, we employ a
kind of ‘perturbation symmetrization’ as follows. Consider a grid point ${\bf
s}$ and model time $t$, where and when the unperturbed specific humidity is
$q_{\rm v}({\bf s},t)$. We compute the tentative perturbation $\Delta q_{\rm
v}({\bf s},t)$ following Eq.(10) and then symmetrically truncate it at $\pm
c$, where $c=\min(q_{\rm v}({\bf s},t),\,q_{\rm sat}({\bf s},t)-q_{\rm v}({\bf
s},t))$ and $q_{\rm sat}({\bf s},t)$ is the saturated specific humidity. Due
to the truncation, the perturbed $q_{\rm v}^{*}({\bf s},t)=q_{\rm v}({\bf
s},t)+\Delta q_{\rm v}({\bf s},t)$ is, first, within the admissible range from
0 to $q_{\rm sat}({\bf s},t)$ and second, as the truncation is symmetric, no
bias is directly introduced: $\operatorname{\mathbb{E}}(q_{\rm v}^{*}({\bf
s},t)\mid q_{\rm v}({\bf s},t))=q_{\rm v}({\bf s},t)$ (which would not be the
case if we just truncated $q_{\rm v}^{*}({\bf s},t)$ at 0 and $q_{\rm
sat}({\bf s},t)$).

####  2.4.3 Soil fields

In the land (soil) model, tendencies of two model fields are perturbed: soil
temperature $T_{\rm so}$ and soil moisture (more specifically, soil water
content $W_{\rm so}$ per unit area within the soil layer in question).
Compared with the atmospheric AMPT, the differences in the treatment of the
soil fields are the following.

  1. 1.

The scaling tendencies ${\cal P}_{T_{\rm so}}$ and ${\cal P}_{W_{\rm so}}$ are
computed by averaging the total tendency, assuming that all processes in the
soil are modeled with substantial uncertainty.

  2. 2.

The averaging in Eq.(11) is performed over land only.

  3. 3.

The perturbation patterns in the soil, $\xi_{T_{\rm so}}(x,y,t)$ and
$\xi_{W_{\rm so}}(x,y,t)$, used to generate perturbations following Eq.(10)
are three-dimensional (not four-dimensional as in the atmosphere).

Whole-domain averaging is used to compute the scaling physical tendency ${\cal
P}_{T_{\rm so}}$. Local-area averaging (i.e., with the averaging length scale
$A_{W_{\rm so}}$ much smaller than the domain size) is employed to compute the
scaling physical tendency ${\cal P}_{W_{\rm so}}$. This choice is motivated by
higher variability/non-Gaussianity of soil moisture compared to soil
temperature (not shown).

The perturbed $W_{\rm so}$ is truncated to ensure that the volumetric soil
water content $\eta_{\rm so}$ is between the wilting-point $\eta_{\rm wp}$ and
field capacity $\eta_{\rm fc}$:

$\eta_{\rm wp}\leq\eta_{\rm so}\equiv\frac{W_{\rm
so}(x,y,Z,t)}{\Delta_{Z}}\leq\eta_{\rm fc},$ (15)

where $Z=1,2,\dots,n_{Z}$ labels the soil layer, $\Delta_{Z}$ is the thickness
of the $Z$-th layer, and $n_{\rm Z}$ is the number of layers.

####  2.4.4 Initial soil perturbations

In the soil, processes have much longer time scales than in the atmosphere,
therefore the role of model perturbations can be revealed only in long-range
forecasts or in cycled systems. Dealing in this study with short-range
forecasts without cycling and having an under-dispersive ensemble of initial
conditions, we had to develop a generator of initial $T_{\rm so}$ and $W_{\rm
so}$ perturbations.

The initial soil temperature perturbation is specified as

$\Delta T_{\rm so}^{\rm ini}(x,y,Z)=\kappa^{\rm ini}_{T_{\rm so}}\,c_{T_{\rm
so}}^{-Z}\,\xi(x,y),$ (16)

where $\kappa^{\rm ini}_{T_{\rm so}}$ is the external magnitude parameter,
$c_{T_{\rm so}}>1$ is the vertical-decay external parameter,
$Z=0,1,\dots,n_{\rm Z}$ (with $Z=0$ denoting the so-called surface temperature
and $Z=1,\dots,n_{\rm Z}$ labeling the soil layers), and $\xi(x,y)$ is the 2D
pseudo-random field. Note that the random pattern $\xi$ is the same for all
soil layers, whilst the magnitude of the perturbation exponentially decreases
downwards.

With the soil moisture $W_{\rm so}$, the technique (inspired by Schraff et al.
2016) is to perturb the Soil Moisture Index

$S=\frac{\eta_{\rm so}-\eta_{\rm wp}}{\eta_{\rm fc}-\eta_{\rm wp}}$ (17)

as follows:

$\Delta S(x,y,Z)=\kappa^{\rm ini}_{S}\,c_{W_{\rm so}}^{1-Z}\,\xi(x,y),$ (18)

where $\kappa^{\rm ini}_{S}$ is the magnitude parameter, $c_{W_{\rm so}}$ is
the vertical-decay parameter, and $Z=1,2,\dots,n_{\rm Z}$. If at a grid point,
the perturbed $S$ appears to lie outside the meaningful range $[0,1]$, the
perturbation $\Delta S$ is truncated accordingly. Using Eqs. (17) and (15), we
finally convert the perturbation $\Delta S(x,y,Z)$ into the perturbation of
$W_{\rm so}$:

$\Delta W_{\rm so}^{\rm ini}(x,y,Z)=\Delta_{Z}\,[\eta_{\rm wp}+\Delta
S(x,y,Z)(\eta_{\rm fc}-\eta_{\rm wp})].$ (19)

###  2.5 Stochastic pattern generator (SPG)

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

$\left(\frac{\partial}{\partial
t}+\frac{U}{\lambda}\sqrt{1-\lambda^{2}\nabla^{2}}\right)^{3}\xi(t,x,y,z)=\sigma\alpha(t,x,y,z).$
(20)

Here $x,y,z$ are the SPG-space spatial coordinates, $\nabla^{2}$ is the 3D
Laplacian, $\alpha(t,x,y,z)$ is the Gaussian white noise, and
$\lambda,U,\sigma$ are the parameters. $\lambda$ determines the spatial scale
of $\xi$. Given $\lambda$, the characteristic velocity $U$ determines the time
scale of $\xi$. Given $\lambda$ and $U$, the parameter $\sigma$ determines the
variance of $\xi$ and is selected to ensure that ${\mathsf{sd}}\,\xi=1$. The
computational domain is the cube of size $2\pi$ 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. 1.

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

  2. 2.

The spatial spectra of $\xi$ 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 & 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^{2}_{*}$ instead of $\nabla^{2}$ in Eq.(20):

$\nabla^{2}_{*}=\frac{\partial^{2}}{\partial
x^{2}}+\gamma^{2}\frac{\partial^{2}}{\partial
y^{2}}+\delta^{2}\frac{\partial^{2}}{\partial z^{2}}.$ (21)

Here $\gamma$ is the aspect ratio of the model domain in the horizontal and
$\delta$ controls the length scale along the $z$ axis.

The SPG internal parameters (for the model variable labeled by $i$)
$\lambda_{i}$ and $\delta_{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 ${\cal T}_{i}=L_{i}/U$, where
$U$ is the above characteristic velocity (another external parameter).

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

SPG can also generate 3D fields $\xi(t,x,y)$ and $\xi(x,y,z)$ as well as 2D
fields $\xi(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 $\xi_{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\mapsto\zeta$ (where $\zeta$ 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 $\zeta$.

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 $\ell$ at the model level $\ell$. 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

####  2.7.1 Stability

Stochastic-dynamic systems in which the magnitude of random forcing depends on
the current state of the system and the forcing (perturbation) is time-
correlated, can be unstable due to a positive feedback loop. Indeed, a
deviation of the model state from its mean value may result in a greater
forcing, which may lead to an even greater deviation of the state from the
mean, and so on until ‘explosion’. To break this vicious circle, we considered
two strategies.

A technically simpler one, which we adopted in this study, is to update the
scaling tendency ${\cal P}_{i}$ not at every time step but less frequently,
once per ${\cal T}_{{\cal P}_{i}}^{\rm update}$ model time during the
forecast, where ${\cal T}_{{\cal P}_{i}}^{\rm update}$ is an external
parameter defined below in section 3.

A somewhat more involved but potentially more powerful approach (which we left
for future research) is to calculate the scaling physical tendency from the
unperturbed (control) model run, ${\bf x}^{\rm control}$, let us denote it
${\cal P}_{i}({\bf x}^{\rm control}({\bf s},t))$ (where $i$ labels the model
variable, ${\bf s}$ is the spatial grid point, and $t$ is time). The scaling
physical tendency, being a moving average of the modulus of the physical
tendency, is a smooth field, so ${\cal P}_{i}({\bf x}^{\rm control}({\bf
s},t))$ can be stored during the control run in a file on a coarse spatio-
temporal grid (meaning that this should be feasible in terms of required
computer resources). Then, ${\cal P}_{i}({\bf x}^{\rm control}({\bf s},t))$
can be used to compute AMPT perturbations for all ensemble members. By
construction, this device will destroy the harmful positive feedback loop. We
remark that taking fields from the control run to make model perturbations
state-dependent can be used to prevent instabilities not only with AMPT but
with other techniques as well.

However, as noted by an anonymous reviewer, this technique cannot be applied
to long-range predictions where ensemble members may diverge from the
unperturbed member to the extent that the physical tendency produced by the
latter becomes irrelevant for the former. Another caveat is that an
application of that technique to an operational ensemble prediction system
will require that the unperturbed member be computed somewhat earlier than the
rest of the ensemble.

Another promising approach is to perturb fluxes instead of prognostic fields,
as discussed in section 2.7.2.

####  2.7.2 Conservation properties

Neither SPPT nor AMPT respects local conservation laws. The reason is that in
both SPPT and AMPT, the model fields are perturbed whereas the fluxes are not.
Switching from perturbing state variables to perturbing fluxes (as done by Van
Ginderachter et al. 2020, for deep convection), including boundary fluxes,
would solve the problem. However, a purely flux-based model perturbation
scheme aimed to represent uncertainties in multiple physical parametrizations
remains to be built and tested, which is beyond the scope of this study.

####  2.7.3 Biases

It is well known and easy to understand that feeding a nonlinear system with
an unbiased signal can lead to a biased system output. Biases can appear even
in a linear system if it is perturbed in a multiplicative way, like in SPPT.
Bouttier et al. (2012, Figs. 3 and 5) and Leutbecher et al. (2017, Fig. A1(a))
found that in their systems, SPPT led to the drying of the atmosphere and
reduced precipitation. In experiments by Romine et al. (2014) the atmosphere
also became too dry due to SPPT (their Figs. 5 and 8) but that was associated
with an increase in precipitation (their Fig. 9). One possible reason that led
to different outcomes in those studies is the impact of the supersaturation
limiter (which nullifies temperature and humidity perturbations at the grid
points where the water vapor is saturated or super-saturated). We
experimentally study forecast biases induced by various configurations of AMPT
and compare them with forecast biases induced by SPPT in section 4.

####  2.7.4 Additive vs. multiplicative model perturbations

In Appendix A, we discuss why SPPT perturbations are multiplicative and why we
call AMPT perturbations additive. We show that AMPT perturbations are truly
additive in the setting when the scaling physical tendency is taken from the
control run (this option was discussed in section 2.7.1).

In numerical experiments presented below, the model perturbation scheme is an
approximation to the setting in which the scaling physical tendency is
completely decoupled from the current model state. The coupling is only
relaxed by updating the scaling physical tendency much less frequently than
every model time step (see sections 2.7.1 and 3). In this regime, AMPT
perturbations are not multiplicative but not yet fully additive. Nevertheless,
with some abuse of terminology, we call them additive bearing in mind the
above scheme with the scaling physical tendency taken from the unperturbed
model run.

##  3 Experimental Settings

The program code of AMPT was built into the limited-area non-hydrostatic COSMO
model (Baldauf et al. 2011), which has both atmospheric and soil prediction
modules. The model was, in turn, embedded in a limited-area ensemble
prediction system. The AMPT-generated model perturbations along with an
ensemble of initial and lateral-boundary conditions allowed us to compute (and
verify) ensemble forecasts. The goal was to assess the effect of AMPT on
deterministic and probabilistic forecasts and compare it with the effect of
SPPT perturbations.

The COSMO model (version 5.01) was used in the convection-permitting
configuration with a horizontal grid spacing of 2.2 km, 172*132 grid points in
the horizontal, and 50 levels in the vertical. The model integration time step
was 20 s. The model’s vertical coordinate was the height-based hybrid (Gal-
Chen) coordinate (Gal-Chen & Somerville 1975).

The ensemble prediction system used in this study was COSMO-Ru2-EPS (Montani
et al. 2014, Astakhova et al. 2015), which was developed within the FROST-2014
international project (Kiktev et al. 2017) and the CORSO priority project of
the COnsortium for Small-scale MOdeling (COSMO) (Rivin et al. 2018). The
ensemble size was 10. COSMO-Ru2-EPS performed a dynamical downscaling of the
forecasts of the driving COSMO-S14-EPS system developed by the Italian
meteorological service ARPAE-SIMC (Montani et al. 2014). Thus, both initial
and lateral boundary conditions for the control forecast and ensemble members
were provided by COSMO-S14-EPS, which had a horizontal grid spacing of 7 km
and 40 vertical levels. COSMO-S14-EPS was a clone of the consortium ensemble
prediction system COSMO-LEPS (Montani et al. 2011) with a smaller ensemble
size.

The model domain is shown in Fig. 1. Note the complexity of the area, which
contains high mountains along with the adjoining valleys and sea. The center
of the domain is located nearly at $44^{\circ}$N, $40^{\circ}$E. The climate
at the sea level is humid subtropical. Numerical experiments were carried out
in this study mostly for the winter-spring season: in February and March, for
which we had access to all data needed to run and verify ensemble forecasts.
Some sensitivity experiments were also conducted for May cases.

Figure 1: Model domain and orography.

The following AMPT parameters were selected for numerical experiments. The
horizontal length scale $L_{i}$ of the SPG driving random fields $\xi_{i}$
(where $i$ labels the perturbed fields) was tuned to be 50 km (i.e., 20–25
horizontal mesh sizes) for $T,u,v$ and 35 km for $q_{\rm v}$. The vertical
length scale $H_{i}$ was about 3 km for $T,u,v$ and about 2 km for humidity.

The characteristic velocity $U=15$ m/s was selected from the physically
meaningful interval of 10–20 m/s.

After some experimentation, the (dimensionless) model-uncertainty magnitude
multipliers $\kappa_{i}$ were specified as 0.75 for $T,u,v$, $T_{\rm so}$,
$W_{\rm so}$ and 0.5 for $q_{\rm v}$.

The length scale $A_{i}$ of averaging the absolute physical tendency (see
section 2.3.2) was specified as equal to the length scale $L_{i}$ of the
respective SPG random field $\xi_{i}$ for $q_{\rm v}$ and $W_{\rm so}$. For
$T,u,v$ and $T_{\rm so}$, the respective $A_{i}$ were selected large enough to
ensure the whole-domain averaging of $|P_{i}|$.

The time update interval ${\cal T}_{{\cal P}_{i}}^{\rm update}$ of the scaling
physical tendency (see section 2.7.1) for the model field $X_{i}$ was set
equal to the time scale ${\cal T}_{i}$ of the respective random pattern
$\xi_{i}$. Note that with ${\cal T}_{{\cal P}_{i}}^{\rm update}$ much less
than ${\cal T}_{i}$, the perturbed model may become unstable (see section
2.7.1), whereas with ${\cal T}_{{\cal P}_{i}}^{\rm update}$ much greater than
${\cal T}_{i}$, the resulting area-averaged ${\cal P}_{i}$ may become
irrelevant in a rapidly developing meteorological situation.

It is worth reiterating at this point that (i) we updated the scaling physical
tendency once per ${\cal T}_{{\cal P}_{i}}^{\rm update}\approx 1$ h model
time, (ii) we computed the physical-space random patterns by performing the
inverse Fourier transform of the spectral-space SPG fields once per ${\cal
T}^{\rm FFT}=20$ min model time, and (iii) we added model perturbations every
model time step.

For the soil fields $T_{\rm so}$ and $W_{\rm so}$, the common time scale
${\cal T}_{\rm so}$ was specified 12 times as large as the atmospheric-
temperature time scale ${\cal T}_{\rm T}$ (so the time scales of the
perturbation fields were, roughly, 1 h in the atmosphere and 12 h in the
soil).

In initial soil perturbations (see section 2.4.4), the magnitude parameter
$\kappa^{\rm ini}_{T_{\rm so}}$, was 1 K. The magnitude parameter of the
initial soil moisture index perturbation, $\kappa^{\rm ini}_{S}$ was only 0.01
(larger values led to unrealistically large model tendencies in $T_{\rm so}$).
The vertical-decay parameters were $c_{T_{\rm so}}=1.75$ and $c_{W_{\rm
so}}=2$.

The mapping from the SPG space to the model space in the vertical was
performed using Option 1 described in section 2.6 (for technical reasons we
couldn’t perform enough numerical experiments with the more physical Option 2
to judge which option is better).

In SPPT, the spatial scale was about 500 km and the time scale was 6 h. The
random multiplier $\kappa\xi$ had the standard deviation 1 and then was
truncated at the absolute value 0.8. This SPPT setup implied greater
perturbations than those explored with the same-resolution COSMO model by
Maurer et al. (2014). We opted for stronger SPPT perturbations because
otherwise, they generated too little spread in the ensemble forecasts. The
supersaturation limiter was off in SPPT. In AMPT, the impact of the
supersaturation limiter is discussed below in section 4.1.

Tapering (i.e., gradual reduction) of perturbations (i) in the lower
troposphere towards the surface and (ii) in the stratosphere from the
tropopause upwards was handled in SPPT as follows. The stratospheric tapering
was always active because it is believed that the radiation tendency, which is
dominant in the stratosphere, is quite accurate in clear-sky conditions (e.g.,
Leutbecher et al. 2017). As for the lower-tropospheric tapering, which is
intended to prevent instabilities due to inconsistencies of perturbed physical
tendencies and unperturbed surface fluxes (Wastl et al. 2019), we found that
SPPT worked better without it. Specifically, our experiments showed that, on
the one hand, SPPT without tapering was stable in the boundary layer. On the
other hand, with tapering, SPPT led to an unacceptably small ensemble spread
in the near-surface fields, so we switched off the lower-tropospheric tapering
in SPPT. In AMPT, tapering was off everywhere.

Note that the above spatial and time scales employed in AMPT were an order of
magnitude less than those in SPPT (50 km vs. 500 km and 1 h vs. 6 h). As
discussed in section 2.3.4, this is reasonable. If in SPPT, $\xi(x,y,t)$ were
small-scale, the product $\xi\,{\bf P}$ would become too patchy, reducing the
effect of the perturbation on the forecast. To find out if this argument is
reasonable, we ran SPPT with a smaller time scale of 1 h and spatial scales of
50 km and 100 km. The resulting spread in the ensemble forecast was indeed
very small, confirming the conclusions of Maurer et al. (2014) (made for a
different domain, orography, physiography, etc.) and justifying the choice of
the SPPT parameters for our domain (the setup was also consistent with that
employed in AROME-EPS by Bouttier et al. (2012)).

Ensemble forecasts were initialized every day at 00 UTC during the two months
of February–March, 2014. The local time was UTC+4h.

Our focus in this study was on near-surface fields: temperature,
precipitation, and wind speed.

##  4 Bias and spread induced by AMPT and SPPT

In preliminary experiments, we found that humidity perturbations did
substantially increase ensemble spread but at the expense of introducing a
significant bias into the forecast. This led us to explore the impact of
humidity perturbations (and the supersaturation limiter) on the spread of
ensemble forecasts and the bias of the ensemble-mean forecast (in terms of
precipitation and near-surface temperature). The aim was to decide whether it
is worth perturbing humidity in the AMPT scheme at all. The (informal)
criterion was twofold: AMPT should generate significantly more spread in the
ensemble forecast than SPPT while having the bias-to-spread ratio as low as
possible. Note that by the bias-to-spread ratio, we mean the absolute value of
the bias divided by the spread.

In the experiments described in this section, to goal was to isolate the roles
of different model perturbations. To this end, we set the same initial and
lateral-boundary conditions for all ensemble members and switched off soil
perturbations. Additionally, hydrostatically balanced pressure perturbations
in AMPT were deactivated. The bias (caused solely by model perturbations) is
defined in this section as the domain-averaged difference between the ensemble
mean and the unperturbed deterministic forecast.

Results are presented in terms of bias-spread scatterplots for near-surface
temperature and accumulated total precipitation. Bias and spread are combined
for lead times from 3 h to 24 h (every three hours for precipitation and every
hour for near-surface temperature) and shown on a single scatterplot. Each
symbol on the plot represents the bias (the value on the $x$-axis) and the
spread (the corresponding value on the $y$-axis) for a single lead time.

###  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: Impact of the supersaturation limiter on spread and bias of
accumulated precipitation forecasts for lead times from 3 to 24 h. Averaging
over cases with the mean (over the model domain) daily accumulated
precipitation $\geq 1$ mm (left) and $\geq 8$ mm (right).

###  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_{\rm 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 $\kappa$ 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: Spread vs. bias of $T_{\rm 2m}$ forecasts for lead times from 3 to
24 h. Averaging over cases with the mean (over the model domain) daily
accumulated precipitation $<1$ mm (top, left), $\geq 1$ mm (top, right), and
$\geq 8$ mm (bottom).

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 Fig. 3. Like in Fig. 3,
the magnitudes of precipitation bias and spread of `SPPT` and `AMPT-UV` were
quite similar. Again, like for $T_{\rm 2m}$, `AMPT-TUV` and `AMPT-TUVQv` had
significantly greater (than `SPPT` and `AMPT-UV`) bias and spread, but unlike
Fig. 3, the bias-to-spread ratio for `AMPT-TUV` was significantly lower than
for `AMPT-TUVQv`.

It is also worth noting that the precipitation biases in Fig. 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.

Figure 4: Spread vs. bias of accumulated precipitation forecasts for lead
times from 3 to 24 h. Averaging over cases with the mean (over the model
domain) daily accumulated precipitation $\geq 1$ mm (left) and $\geq 8$ mm
(right).

We summarize findings from Figs. 3 and 4 as follows.

  1. 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. 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. 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_{\rm 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 Testing AMPT in the ensemble prediction system

###  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.

Table 1: List of ensemble prediction experiments Experiment | Model perturbations
---|---
`NOPERT` | None
`SPPT` | Atmospheric SPPT perturbations
`AMPT-NOSOIL` | Atmospheric AMPT perturbations
`AMPT-SOIL` | Atmospheric and soil AMPT perturbations

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 $\xi_{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_{\rm 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_{\rm 10m}$, verification results were, mostly, similar to those for $T_{\rm
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_{\rm 10m}$ compared to $T_{\rm 2m}$, we only
display one typical plot for $V_{\rm 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 Fig. 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 $\alpha=0.05$. The number of bootstrap samples
was $10^{5}$. Along with the $p$-value of the above bootstrap test, $p_{\rm
bootstrap}$, we also computed the $p$-value of the classical one-sided
Student’s t-test, $p_{\rm 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_{\rm 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 Fig.
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.

To highlight the barely seen in Fig. 5(left) differences between the curves in
the upper bunch, Fig. 5(right) displays the normalized reduction in the
$T_{\rm 2m}$ ensemble-mean RMSE with respect to `NOPERT`, that is,
${\rm(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_{\rm 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.

Figure 5: $T_{\rm 2m}$. Left: RMSE of ensemble mean (the upper bunch of
curves) and ensemble spread (the lower bunch of curves). Right: The normalized
reduction in ensemble-mean RMSE, that is,
${\rm(RMSE_{NOPERT}-RMSE)/RMSE_{NOPERT}}$, the higher the better.

The lower bunch of curves in Fig. 5(left) shows the ensemble spread in $T_{\rm
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_{\rm 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
$\sigma$) 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 $\sigma^{2}$ 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 Fig.
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 Fig. 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
Fig.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: RMSE and spread for AMPT-NOSOIL vs. SPPT in free atmosphere.
Left: Temperature. Right: Wind speed.

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_{\rm v}$ given the
ensemble cumulative distribution function $F_{\rm e}(x)$ of the model variable
$x$, i.e., $P(x_{\rm v}{\leq}x\mid F_{\rm e})$ (where $P$ stands for
probability), is systematically different from $F_{\rm e}(x)$. This kind of
inconsistency is known as lack of reliability of the probabilistic forecast
(Toth et al. 2003). In these terms, Figs. 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 score444 The
reliability component of the Brier score is a weighted mean-square difference
between the predicted event probabilities $p^{\rm f}_{k}$ (where $k$ labels
the discrete values of the forecast probability obtainable from a finite-size
ensemble) and the respective observed frequencies of the event for all cases
when the probability $p^{\rm f}_{k}$ was forecast. for the (most populated)
event $T_{\rm 2m}>0^{\circ}C$. The improvement of `AMPT-NOSOIL` with respect
to `SPPT` was statistically significant, $p_{\rm bootstrap}=0.05$ (though
$p_{\rm Student}=0.07$ did not reach the significance level of 0.05).

Figure 7: The reliability component of the Brier score for the event $T_{\rm
2m}>0^{\circ}C$. The lower the better.

###  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_{\rm 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_{\rm bootstrap}=0.03$. The $p$-value
of the Student’s t-test was $p_{\rm Student}=0.02$. `AMPT-SOIL` performed a
fortiori statistically significantly better than `AMPT-NOSOIL`.

Figure 8: CRPS for $T_{\rm 2m}$. The lower the better. Note that the y-axis
does not start at 0.

Figure 9 displays CRPS for near-surface wind speed $V_{\rm 10m}$. The impact
of AMPT on the forecast performance in terms of $V_{\rm 10m}$ was weaker than
on $T_{\rm 2m}$ but more stable, so the statistical significance of the
advantage of `AMPT-NOSOIL` vs. `SPPT` was stronger: $p_{\rm bootstrap}<0.001$
and $p_{\rm 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: CRPS for $V_{\rm 10m}$. The lower the better. Note that the y-axis
does not start at 0.

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<\theta$ (or $x>\theta$), where $x$ is the continuous model
variable in question (say, temperature) and $\theta$ is a threshold, then CRPS
is the integral of the Brier score over all possible $\theta$ (Hersbach 2000).
Therefore, Fig. 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_{\rm 2m}<-5^{\circ}C$, $T_{\rm 2m}>0^{\circ}C$, $T_{\rm
2m}>5^{\circ}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_{\rm 2m}>0^{\circ}C$ in Fig. 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_{\rm bootstrap}=0.03$ and $p_{\rm Student}=0.01$. Similar
results were obtained for near-surface wind speed, not shown.

Figure 10: Brier score for the event $T_{\rm 2m}>0^{\circ}C$. The lower the
better.

We also examined the resolution component of the Brier Score. With the
reliability component substantially improved, see Fig. 7, the resolution
component was barely changed or slightly reduced (i.e., degraded, not shown).
As noted by Candille & 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_{\rm v}$ given the
ensemble cumulative distribution function $F_{\rm e}(x)$, i.e., $p(x_{\rm
v}\mid F_{\rm e})$ (where $p$ stands for probability density). The reciprocal
and complementary view on the accuracy of a probabilistic forecast is through
$p(F_{\rm e}\mid x_{\rm v})$, i.e., conditional on verification. The
capability of a probabilistic forecast to concentrate the probability mass
close to the observed $x_{\rm v}$ (thus, discriminating between different
outcomes $x_{\rm 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_{\rm 2m}>0^{\circ}C$ (similar results were obtained
for the two other events, $T_{\rm 2m}<-5^{\circ}C$ and $T_{\rm
2m}>5^{\circ}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_{\rm bootstrap}=p_{\rm Student}=0.01$. As `AMPT-SOIL` is seen in Fig. 11 to
be almost uniformly better than `AMPT-NOSOIL`, we conclude that `AMPT-SOIL`
was statistically significantly more skillful for near-surface temperature
than SPPT, too. For near-surface wind speed, the results were mixed (not
shown).

A few remarks regarding Fig. 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 & Mass (2005),
meaning greater spread when and where it is relevant.

Figure 11: ROC area for the event $T_{\rm 2m}>0^{\circ}C$. The higher the
better. Note that the y-axis does not start at 0.

###  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.

##  6 Conclusions

A new technique called Additive Model-uncertainty perturbations scaled by
Physical Tendency (AMPT) has been proposed. AMPT addresses some issues of the
wide-spread SPPT scheme. First, in contrast to SPPT, AMPT can generate
perturbations with significant magnitude even at grid points where the model
physical tendency happens to be spuriously small or is small due to
cancellation of contributions from different physical parametrizations.
Second, AMPT improves on the SPPT’s unphysical perfect correlations of
perturbations between different model variables and in the vertical. Third,
due to the non-local link from the model fields to the scaling physical
tendency, which determines the magnitude of AMPT perturbations, the SPPT’s ban
on sign reversal of the physical tendency, Eq.(5), is removed in AMPT. This
enables AMPT to generate significantly greater perturbations than SPPT without
causing instabilities.

Taking the unperturbed model run to introduce state-dependence in model
perturbations for ensemble members is discussed. This approach breaks the
positive feedback loop from the model state to perturbations and then back to
the model state at a later lead time, therefore, it can be used to maintain
stability of ensemble forecasts perturbed in a state-dependent manner. We
found in this research that just relaxing the pointwise dependence of the
model perturbation magnitude on the model state (by using a spatial moving
average of the magnitude of physical tendency and by updating it less
frequently than each model time step) helped much in preventing instabilities.

AMPT employs the Stochastic Pattern Generator (SPG, Tsyrulnikov & Gayfulin
2017) to generate four-dimensional random fields with tunable spatio-temporal
correlations (but can be used with any pattern generator). The random fields
generated by SPG to perturb different model fields are mutually independent.
They are scaled by the area-averaged modulus of physical tendency and added to
the model fields at every model time step. AMPT perturbations of three-
dimensional atmospheric model fields of temperature, pressure (computed from
temperature perturbations via hydrostatics), wind, humidity, and three-
dimensional soil fields (temperature and moisture) were imposed and their
effects on convection-permitting ensemble forecasts (based on the COSMO
forecast model) were examined.

Practically, AMPT performed better than SPPT in our ensemble prediction
experiments but at the expense of being computationally more expensive than
SPPT. The bigger cost of AMPT was mostly because the AMPT’s random patterns
are three-dimensional (in the atmosphere), as opposed to two-dimensional
random patterns in SPPT. In two dimensions (in the horizontal, soil
perturbations), the cost of running AMPT was negligibly small. The overhead of
running AMPT with the perturbed three-dimensional $T,u,v$ fields (as well as
$p$, $T_{\rm SO}$, and $W_{\rm SO}$) was about 10% of the time of running the
COSMO model itself. There were three main computer-intensive parts of AMPT.
First, the spectral SPG solver, which was optimized and parallelized, was
responsible for about 25% of the total burden of AMPT. Second, the three-
dimensional fast Fourier transform (FFT) took another 25% of computer time.
FFT was only partly optimized and done with an old legacy code. The third
computer intensive part was the computation of the scaling physical tendencies
${\cal P}_{i}$ and spatial interpolations. This part was not parallelized at
the time of writing this paper. As a result, AMPT was about 10 times more
expensive than SPPT, but the gap can be reduced by a factor of three after the
AMPT code is fully optimized in preparation for an operational application.

Another difference between AMPT and SPPT is that AMPT has many more
parameters. On the one hand, this can be viewed as a disadvantage because AMPT
requires more effort to tune them. On the other hand, all the AMPT parameters
have clear physical meaning and are set to meaningful default values. This
offers the AMPT user a choice: either to set the parameters to their default
values indicated in section 3 or to tune some of the parameters taking into
account specifics of a particular application.

It is also worth noting that AMPT suffers from a lack of physical consistency
even to a greater extent than SPPT since AMPT perturbations are not balanced.
We argue that this might not be a problem because of a small magnitude of
model perturbations introduced every model time step so that the model manages
to balance them itself. Experimentally, we found no indications of problems
related to the lack of balance in AMPT perturbations.

The main findings of the study are the following.

  * •

Humidity perturbations generated significant ensemble spread in precipitation
forecasts but led to a high bias-to-spread ratio. For this reason, AMPT was
systematically tested without humidity perturbations.

  * •

Withholding temperature perturbations led to a substantial reduction in the
bias and in the bias-to-spread ratio, but at the expense of a substantial
reduction in spread as well. For this reason, temperature (and pressure) were
included in the list of perturbed model fields in the final ensemble
prediction experiments.

  * •

In ensemble prediction experiments, it was found that AMPT was much more
effective in generating spread in the ensemble than SPPT, thus considerably
improving reliability of the ensemble. This was the case for the upper-air
fields as well as for the near-surface fields.

  * •

Most probabilistic ensemble verification scores for near-surface temperature
and wind speed were improved due to atmospheric AMPT perturbations as compared
with SPPT (at the statistical significance level of 0.05).

  * •

Impacts of AMPT on the root-mean-square errors of the ensemble-mean
temperature and wind speed forecasts were mixed (compared with SPPT). This
held for free-atmosphere as well as near-surface variables.

  * •

AMPT perturbations of soil moisture and soil temperature significantly
improved deterministic and probabilistic verification scores for near-surface
temperature and had little effect on near-surface wind speed.

  * •

Probabilistic verification of ensemble forecasts in terms of precipitation
gave mixed results, perhaps, due to an insufficient number of precipitating
events during the time period examined and scarcity of the observation
network.

The technique can be further developed in the following directions. First, it
looks reasonable to introduce state-dependence not only in the magnitude of
perturbation fields but also in their spatial and temporal scales. This will
require a non-stationary (in space and time) stochastic random field generator
(e.g., based on the technique proposed in Tsyrulnikov & Rakitko (2019)) and a
method to specify the non-stationarity patterns. Second, in its current
formulation, AMPT takes the net physical tendency as input but it can be used
at the process level as well, i.e., for the output of each parametrization
scheme separately. Third, switching from perturbations of tendencies to
perturbations of fluxes looks as a promising way to achieve physical
consistency. Finally, we note that AMPT can be used in data assimilation as
well as in ensemble prediction systems: in the atmosphere, in the soil, and
also in the ocean.

## Acknowledgments

We are grateful to T.Paccagnella and A.Montani, who provided us with initial
and boundary conditions for our forecast ensemble. We would like to thank
members of the COSMO Working Group on Predictability and Ensemble Methods for
many fruitful discussions. V.Volkova assisted with the ensemble forecast
verification in the free atmosphere. A.Kirsanov and A.Bundel kindly helped us
to use the VERSUS verification tool. We also thank D.Blinov for assistance in
working with the COSMO model. Valuable comments of two anonymous reviewers led
to significant improvements of the manuscript. The study has been conducted
under Tasks 1.1.1 and 1.1.4 of the Scientific Research and Technology Plan of
the Russian Federal Service for Hydrometeorology and Environmental Monitoring.

## Data availability statement

The dataset on which this paper is based is too large to be retained or
publicly archived with available resources. Documentation and methods used to
support this study are available from the first author. The Fortran source
code of SPG is available from https://github.com/gayfulin/SPG. Reproducible R
code that realizes the methodology of statistical significance testing
(presented in Appendix B), along with raw data, can be found at
https://github.com/cyrulnic/AMPT.

## Disclosure statement

No potential conflict of interest was reported by the authors.

## Authors’ contributions

All authors contributed to the study conception and design. Michael
Tsyrulnikov wrote the manuscript. Elena Astakhova performed ensemble
prediction experiments and verification. Dmitry Gayfulin wrote the program
code and conducted preliminary sensitivity tests. All authors commented on
previous versions of the manuscript. All authors read and approved the final
manuscript.

## Appendices

##  Appendix 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 & Gayfulin (2017, Eq.(1.3)): ${{\rm
d}{\bf x}}/{{\rm d}t}={\bf F}({\bf x})-\boldsymbol{\xi}$ (where ${\bf x}$ is
the truth, ${\bf F}$ is the forecast-model operator, and $\boldsymbol{\xi}$ is
the random model perturbation). So, the model perturbation is multiplicative
if it is a function of ${\bf x}$ (the state of the ensemble member on which
the model perturbation is imposed) and additive if it does not depend on ${\bf
x}$.

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

$\Delta P_{i}^{\rm SPPT}({\bf s},t)=\kappa\,P_{i}({\bf x}({\bf
s},t))\,\xi({\bf s},t).$ (22)

With AMPT, if the scaling tendency is taken from the control (unperturbed)
run, so that ${\cal P}_{i}({\bf s},t)={\cal P}_{i}({\bf x}^{\rm control}({\bf
s},t))$, the perturbation ceases to depend on ${\bf x}$:

$\Delta P_{i}^{\rm AMPT}({\bf s},t)=\kappa_{i}\,{\cal P}_{i}({\bf x}^{\rm
control}({\bf s},t))\,\xi_{i}({\bf s},t),$ (23)

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.

##  Appendix 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 $\psi$ is greater for
scheme 1 than for scheme 2. More precisely, we are going to compare the scores
$\psi_{1}$ and $\psi_{2}$ 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\operatorname{\mathbb{E}}\widehat{\psi}_{1}>\operatorname{\mathbb{E}}\widehat{\psi}_{2}\equiv\psi_{2}$.
The default ‘null’ hypothesis states that there is no difference between the
two schemes: $H_{0}:\psi_{1}=\psi_{2}$.

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 $\psi$ 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,\dots,n$) and their differences $\delta_{j}^{\rm
obs}=\widehat{\psi}_{1j}-\widehat{\psi}_{2j}$. We treat $\delta_{j}^{\rm obs}$
as (observed) realizations of the respective random variables $\delta_{j}$.
Assuming that temporal dependencies in $\delta_{j}$ decay on a time scale
shorter than one week, we regard $\delta_{j}$ as a sample of statistically
independent and identically distributed random numbers.

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

To compute $p$ we need to know the probability distribution of $\delta_{j}$
under $H_{0}$. 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 $\delta_{j}$) equals the
empirical distribution (which is concentrated at the observed data, in our
case, $\delta_{j}^{\rm obs}$). However, this setup cannot be directly applied
here because $H_{0}$ assumes that $\operatorname{\mathbb{E}}\delta_{j}=0$
whereas the empirical distribution is, likely, biased ($\bar{\delta}^{\rm
obs}\neq 0$). Following Efron & 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}^{\rm obs}=\delta_{j}^{\rm obs}-\bar{\delta}^{\rm obs}$), compute
the mean in each bootstrap sample, $\bar{\delta}^{*}$, and calculate the
fraction of bootstrap samples in which $\bar{\delta}^{*}>\bar{\delta}^{\rm
obs}$. This fraction is the estimate of the $p$-value.

## References

  * (1)
  * Astakhova et al. (2015) Astakhova, E., Montani, A. & Alferov, D. Y. (2015), ‘Ensemble forecasts for the Sochi-2014 Olympic Games’, Russian Meteorology and Hydrology 40(8), 531–539.
  * Baldauf et al. (2011) Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M. & Reinhardt, T. (2011), ‘Operational convective-scale numerical weather prediction with the COSMO model: description and sensitivities’, Monthly Weather Review 139(12), 3887–3905.
  * Berner et al. (2015) Berner, J., Fossell, K., Ha, S.-Y., Hacker, J. & Snyder, C. (2015), ‘Increasing the skill of probabilistic forecasts: Understanding performance improvements from model-error representations’, Monthly Weather Review 143(4), 1295–1320.
  * Bouttier et al. (2022) Bouttier, F., Fleury, A., Bergot, T. & Riette, S. (2022), ‘A single-column comparison of model-error representations for ensemble prediction’, Boundary-Layer Meteorology 183(2), 167–197.
  * Bouttier et al. (2012) Bouttier, F., Vié, B., Nuissier, O. & Raynaud, L. (2012), ‘Impact of stochastic physics in a convection-permitting ensemble’, Monthly Weather Review 140(11), 3706–3721.
  * Buizza et al. (1999) Buizza, R., Miller, M. & Palmer, T. (1999), ‘Stochastic representation of model uncertainties in the ECMWF ensemble prediction system’, Quarterly Journal of the Royal Meteorological Society 125(560), 2887–2908.
  * Candille & Talagrand (2005) Candille, G. & Talagrand, O. (2005), ‘Evaluation of probabilistic prediction systems for a scalar variable’, Quarterly Journal of the Royal Meteorological Society 131(609), 2131–2150.
  * Christensen et al. (2017) Christensen, H., Lock, S.-J., Moroz, I. & Palmer, T. (2017), ‘Introducing independent patterns into the stochastically perturbed parametrization tendencies (SPPT) scheme’, Quarterly Journal of the Royal Meteorological Society 143(706), 2168–2181.
  * Clark et al. (2021) Clark, P., Halliwell, C. & Flack, D. (2021), ‘A physically based stochastic boundary layer perturbation scheme. Part I: Formulation and evaluation in a convection-permitting model’, Journal of the Atmospheric Sciences 78(3), 727–746.
  * Dorrestijn et al. (2013) Dorrestijn, J., Crommelin, D. T., Siebesma, A. P. & Jonker, H. J. (2013), ‘Stochastic parameterization of shallow cumulus convection estimated from high-resolution model data’, Theoretical and Computational Fluid Dynamics 27(1), 133–148.
  * Eckel & Mass (2005) Eckel, F. A. & Mass, C. F. (2005), ‘Aspects of effective mesoscale, short-range ensemble forecasting’, Weather and Forecasting 20(3), 328–350.
  * Efron & Tibshirani (1994) Efron, B. & Tibshirani, R. J. (1994), An introduction to the bootstrap, CRC press.
  * Fortin et al. (2014) Fortin, V., Abaza, M., Anctil, F. & Turcotte, R. (2014), ‘Why should ensemble spread match the RMSE of the ensemble mean?’, Journal of Hydrometeorology 15(4), 1708–1713.
  * Frogner et al. (2022) Frogner, I.-L., Andrae, U., Ollinaho, P., Hally, A., Hämäläinen, K., Kauhanen, J., Ivarsson, K.-I. & Yazgi, D. (2022), ‘Model uncertainty representation in a convection-permitting ensemble — SPP and SPPT in HarmonEPS’, Monthly Weather Review 150(4), 775–795.
  * Fuchs (2014) Fuchs, A. (2014), Nonlinear dynamics in complex systems, Springer.
  * Gal-Chen & Somerville (1975) Gal-Chen, T. & Somerville, R. C. (1975), ‘On the use of a coordinate transformation for the solution of the Navier-Stokes equations’, Journal of Computational Physics 17(2), 209–228.
  * Gofa et al. (2010) Gofa, F., Raspanti, A. & Charantonis, T. (2010), Assessment of the performance of the operational forecast model COSMOGR using a conditional verification tool, in ‘Proc. 10th COMECAP Conference of Meteorology, Patras. Greece, Hellenic Meteorological Society’.
  * Hersbach (2000) Hersbach, H. (2000), ‘Decomposition of the continuous ranked probability score for ensemble prediction systems’, Weather and Forecasting 15(5), 559–570.
  * Hirt et al. (2019) Hirt, M., Rasp, S., Blahak, U. & Craig, G. C. (2019), ‘Stochastic parameterization of processes leading to convective initiation in kilometer-scale models’, Monthly Weather Review 147(11), 3917–3934.
  * Kiktev et al. (2017) Kiktev, D., Joe, P., Isaac, G. A., Montani, A., Frogner, I.-L., Nurmi, P., Bica, B., Milbrandt, J., Tsyrulnikov, M., Astakhova, E., Bundel, A., Belair, S., Pyle, M., Muravyev, A., Rivin, G., Rozinkina, I., Paccagnella, T., Wang, Y., Reid, J., Nipen, T. & Ahn, K.-D. (2017), ‘FROST-2014: The Sochi winter Olympics international project’, Bulletin of the American Meteorological Society 98(9), 1908–1929.
  * Kober & Craig (2016) Kober, K. & Craig, G. C. (2016), ‘Physically based stochastic perturbations (PSP) in the boundary layer to represent uncertainty in convective initiation’, Journal of Atmospheric Sciences 73(7), 2893–2911.
  * Kwasniok (2012) Kwasniok, F. (2012), ‘Data-based stochastic subgrid-scale parametrization: an approach using cluster-weighted modelling’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370(1962), 1061–1086.
  * Lang et al. (2021) Lang, S. T., Lock, S.-J., Leutbecher, M., Bechtold, P. & Forbes, R. M. (2021), ‘Revision of the stochastically perturbed parametrisations model uncertainty scheme in the integrated forecasting system’, Quarterly Journal of the Royal Meteorological Society 147(735), 1364–1381.
  * Leutbecher et al. (2017) Leutbecher, M., Lock, S.-J., Ollinaho, P., Lang, S. T., Balsamo, G., Bechtold, P., Bonavita, M., Christensen, H. M., Diamantakis, M., Dutra, E. et al. (2017), ‘Stochastic representations of model uncertainties at ECMWF: State of the art and future vision’, Quarterly Journal of the Royal Meteorological Society 143(707), 2315–2339.
  * Leutbecher & Palmer (2008) Leutbecher, M. & Palmer, T. N. (2008), ‘Ensemble forecasting’, Journal of Computational Physics 227(7), 3515–3539.
  * Machulskaya & Seifert (2019) Machulskaya, E. & Seifert, A. (2019), ‘Stochastic differential equations for the variability of atmospheric convection fluctuating around the equilibrium’, Journal of Advances in Modeling Earth Systems 11(8), 2708–2727.
  * Maurer et al. (2014) Maurer, D., Walser, A. & Arpagaus, M. (2014), ‘First COSMO-E experiments with the stochastically perturbed parametrization tendencies (SPPT) scheme’, COSMO Newsletter 14, 19–27.
  * McTaggart-Cowan et al. (2022) McTaggart-Cowan, R., Separovic, L., Charron, M., Deng, X., Gagnon, N., Houtekamer, P. L. & Patoine, A. (2022), ‘Using stochastically perturbed parameterizations to represent model uncertainty. Part II: Comparison with existing techniques in an operational ensemble’, Monthly Weather Review 150(11), 2859–2882.
  * Monin & Yaglom (2013) Monin, A. S. & Yaglom, A. M. (2013), Statistical fluid mechanics, volume II: Mechanics of turbulence, Vol. 2, Courier Corporation.
  * Montani et al. (2014) Montani, A., Alferov, D., Astakhova, E., Marsigli, C. & Paccagnella, T. (2014), ‘Ensemble forecasting for Sochi-2014 Olympics: The COSMO-based ensemble prediction systems’, COSMO Newsletter (14), 88–94.
  * Montani et al. (2011) Montani, A., Cesari, D., Marsigli, C. & Paccagnella, T. (2011), ‘Seven years of activity in the field of mesoscale ensemble forecasting by the COSMO-LEPS system: main achievements and open challenges’, Tellus A: Dynamic Meteorology and Oceanography 63(3), 605–624.
  * Ollinaho et al. (2017) Ollinaho, P., Lock, S.-J., Leutbecher, M., Bechtold, P., Beljaars, A., Bozzo, A., Forbes, R. M., Haiden, T., Hogan, R. J. & Sandu, I. (2017), ‘Towards process-level representation of model uncertainties: Stochastically perturbed parametrisations in the ECMWF ensemble’, Quarterly Journal of the Royal Meteorological Society 143(702), 408–422.
  * Orrell et al. (2001) Orrell, D., Smith, L., Barkmeijer, J. & Palmer, T. N. (2001), ‘Model error in weather forecasting’, Nonlinear processes in geophysics 8(6), 357–371.
  * Palmer (2012) Palmer, T. (2012), ‘Towards the probabilistic earth-system simulator: A vision for the future of climate and weather prediction’, Quarterly Journal of the Royal Meteorological Society 138(665), 841–861.
  * Palmer et al. (2005) Palmer, T., Shutts, G., Hagedorn, R., Doblas-Reyes, F., Jung, T. & Leutbecher, M. (2005), ‘Representing model uncertainty in weather and climate prediction’, Annu. Rev. Earth Planet. Sci. 33, 163–193.
  * Plant & Craig (2008) Plant, R. & Craig, G. C. (2008), ‘A stochastic parameterization for deep convection based on equilibrium statistics’, Journal of the Atmospheric Sciences 65(1), 87–105.
  * Rieger et al. (2021) Rieger, D., Milelli, M., Boucouvala, D., Gofa, F., Iriza-Burca, A., Khain, P., Kirsanov, A., Linkowska, J. & Marcucci, F. (2021), ‘Verification of ICON in limited area mode at COSMO national meteorological services’, Reports on ICON (6).
  * Rivin et al. (2018) Rivin, G., Rozinkina, I., Astakhova, E. & Coauthors (2018), ‘The COSMO Priority Project CORSO (Consolidation of Operational and Research Results for the Sochi Olympic Games). Final Report’, COSMO Technical report (35), 249–262.
  * Romine et al. (2014) Romine, G. S., Schwartz, C. S., Berner, J., Fossell, K. R., Snyder, C., Anderson, J. L. & Weisman, M. L. (2014), ‘Representing forecast error in a convection-permitting ensemble system’, Monthly Weather Review 142(12), 4519–4541.
  * Sakradzija et al. (2015) Sakradzija, M., Seifert, A. & Heus, T. (2015), ‘Fluctuations in a quasi-stationary shallow cumulus cloud ensemble’, Nonlinear Processes in Geophysics 22(1), 65–85.
  * Schraff et al. (2016) Schraff, C., Reich, H., Rhodin, A., Schomburg, A., Stephan, K., Perianez, A. & Potthast, R. (2016), ‘Kilometre-scale ensemble data assimilation for the COSMO model (KENDA)’, Quarterly Journal of the Royal Meteorological Society 142(696), 1453–1472.
  * Shutts & Pallarès (2014) Shutts, G. & Pallarès, A. C. (2014), ‘Assessing parametrization uncertainty associated with horizontal resolution in numerical weather prediction models’, Phil. Trans. R. Soc. A 372(2018), 20130284.
  * Toth et al. (2003) Toth, Z., Talagrand, O., Candille, G. & Zhu, Y. (2003), Probability and ensemble forecasts, in ‘Jolliffe, I.T. and Stephenson, D.B. Forecast verification: A practitioner’s guide in atmospheric science’, John Wiley and Sons, pp. 137–163.
  * Touil et al. (2007) Touil, H., Hussaini, M., Gotoh, T., Rubinstein, R. & Woodruff, S. (2007), ‘Development of stochastic models for turbulence’, New Journal of Physics 9(7), 215.
  * Tsyroulnikov (2001) Tsyroulnikov, M. D. (2001), ‘Proportionality of scales: An isotropy-like property of geophysical fields’, Quarterly Journal of the Royal Meteorological Society 127(578), 2741–2760.
  * Tsyrulnikov & Gayfulin (2017) Tsyrulnikov, M. & Gayfulin, D. (2017), ‘A limited-area spatio-temporal stochastic pattern generator for simulation of uncertainties in ensemble applications’, Meteorol. Zeitschrift 26(5), 549–566.
  * Tsyrulnikov & Rakitko (2019) Tsyrulnikov, M. & Rakitko, A. (2019), ‘Impact of non-stationarity on hybrid ensemble filters: A study with a doubly stochastic advection-diffusion-decay model’, Quarterly Journal of the Royal Meteorological Society 145(722), 2255–2271.
  * Van Ginderachter et al. (2020) Van Ginderachter, M., Degrauwe, D., Vannitsem, S. & Termonia, P. (2020), ‘Simulating model uncertainty of subgrid-scale processes by sampling model errors at convective scales’, Nonlinear Processes in Geophysics 27(2), 187–207.
  * Wastl et al. (2019) Wastl, C., Wang, Y., Atencia, A. & Wittmann, C. (2019), ‘A hybrid stochastically perturbed parametrization scheme in a convection-permitting ensemble’, Monthly Weather Review 147(6), 2217–2230.
  * Wilks (2011) Wilks, D. S. (2011), Statistical methods in the atmospheric sciences, Academic press.