# Machine Learning for Stochastic Parameterization:
Generative Adversarial Networks in the Lorenz ’96 Model

###### Abstract

Stochastic parameterizations account for uncertainty in the representation of
unresolved sub-grid processes by sampling from the distribution of possible
sub-grid forcings. Some existing stochastic parameterizations utilize data-
driven approaches to characterize uncertainty, but these approaches require
significant structural assumptions that can limit their scalability. Machine
learning models, including neural networks, are able to represent a wide range
of distributions and build optimized mappings between a large number of inputs
and sub-grid forcings. Recent research on machine learning parameterizations
has focused only on deterministic parameterizations. In this study, we develop
a stochastic parameterization using the generative adversarial network (GAN)
machine learning framework. The GAN stochastic parameterization is trained and
evaluated on output from the Lorenz ’96 model, which is a common baseline
model for evaluating both parameterization and data assimilation techniques.
We evaluate different ways of characterizing the input noise for the model and
perform model runs with the GAN parameterization at weather and climate
timescales. Some of the GAN configurations perform better than a baseline
bespoke parameterization at both timescales, and the networks closely
reproduce the spatio-temporal correlations and regimes of the Lorenz ’96
system. We also find that in general those models which produce skillful
forecasts are also associated with the best climate simulations.

Journal of Advances in Modeling Earth Systems (JAMES)

National Center for Atmospheric Research, Boulder, Colorado, USA Atmospheric,
Oceanic and Planetary Physics, University of Oxford, UK University of Colorado
Boulder, Boulder, Colorado, USA School of Earth and Ocean Sciences, University
of Victoria, Victoria, British Columbia, Canada

David John Gagne<EMAIL_ADDRESS>
We have developed a generative adversarial network machine learning stochastic
parameterization of sub-grid forcing for the Lorenz ’96 dynamical model.

We test different configurations of the generative adversarial network
parameterization at weather and climate timescales and find that some have
representation errors better than a bespoke polynomial parameterization.

The generative adversarial networks closely reproduce the spatio-temporal
correlations and regimes of the Lorenz ’96 system.

## Plain Language Summary

Simulations of the atmosphere must approximate the effects of small-scale
processes with simplified functions called parameterizations. Standard
parameterizations only predict one output for a given input, but stochastic
parameterizations can sample from all the possible outcomes that can occur
under certain conditions. We have developed and evaluated a machine learning
stochastic parameterization, which builds a mapping between large-scale
current conditions and the range of small-scale outcomes from data about both.
We test the machine learning stochastic parameterization in a simplified
mathematical simulation that produces multi-scale chaotic waves like the
atmosphere. We find that some configurations of the machine learning
stochastic parameterization perform slightly better than a simpler baseline
stochastic parameterization over both weather- and climate-like time spans.

##  1 Introduction

A large source of weather and climate model uncertainty is the approximate
representation of unresolved sub-grid processes through parameterization
schemes. Traditional, deterministic parameterization schemes represent the
mean or most likely sub-grid scale forcing for a given resolved-scale state.
While model errors can be reduced to a certain degree through improvements to
such parameterizations, they cannot be eliminated. Irreducible uncertainties
result from a lack of scale separation between resolved and unresolved
processes. Uncertainty in weather forecasts also arises because the chaotic
nature of the atmosphere gives rise to sensitivity to uncertain initial
conditions. Practically, uncertainty is represented in forecasts using
ensembles of integrations of comprehensive weather and climate prediction
models, first suggested by Leith:1975jc. To produce reliable probabilistic
forecasts, the generation of the ensemble must include a representation of
both model and initial condition uncertainty.

Initial condition uncertainty is addressed by perturbing the initial
conditions of ensemble members, for example by selecting directions of optimal
perturbation growth using singular vectors [Buizza Palmer (1995)], or by
characterizing initial condition uncertainty during the data assimilation
cycle [Isaksen . (2010)]. One approach for representing irreducible model
uncertainty is stochastic parameterization of unresolved physical processes. A
stochastic parameterization represents the probability distribution of
possible sub-grid scale tendencies conditioned on the large scale. Each
ensemble member experiences one possible, equally likely realization of the
sub-grid-scale tendencies. A more detailed motivation for including stochastic
parameterizations in weather and climate models is presented in Palmer:2012fe.

Stochastic approaches for numerical weather prediction (NWP) were originally
proposed for use in the European Center for Medium-Range Weather Forecasts
(ECMWF) ensemble prediction system [Palmer . (1997), Buizza . (1999)]. They
were shown to substantially improve the quality of initialized ensemble
forecasts, and so became widely adopted by meteorological services around the
world, where they are used to produce operational ensemble weather, sub-
seasonal and seasonal forecasts [Teixeira Reynolds (2010), Reyes . (2009),
Berner . (2010), Stockdale . (2011), Palmer, Doblas-Reyes Weisheimer (2009),
Palmer (2012), Suselj . (2013), Sušelj . (2014), Weisheimer . (2014), Sanchez
. (2016), Leutbecher . (2017)].

Recent work has assessed the impact of stochastic parameterization schemes in
both idealized and state-of-the-art climate models for long term integration
[PD. Williams (2012), Ajayamohan . (2013), Juricke Jung (2014), Dawson Palmer
(2015), Wang . (2016), HM. Christensen . (2017), Davini . (2017), Strømmen .
(2018)]. These studies demonstrate that including a stochastic representation
of model uncertainty can go beyond improving initialized forecast reliability,
and can also lead to improvements in the model mean state [Palmer (2001),
Berner . (2012)], climate variability [Ajayamohan . (2013), Dawson Palmer
(2015), HM. Christensen . (2017)], and change a model’s climate sensitivity
[Seiffert von Storch (2010)]. These impacts occur through non-linear
rectification, noise enhanced variability, and noise-induced regime
transitions [Berner . (2017)]. In this way, small-scale variability
represented by stochastic physics can impact large spatio-temporal scales of
climate variability.

Despite the historical disconnect between the weather and climate prediction
communities, the boundaries between weather and climate prediction are
somewhat artificial [Palmer . (2008), Hurrell . (2009), Shapiro . (2010)].
This disconnect is highlighted by recent advances in prediction on timescales
from weather to sub-seasonal-to-seasonal and decadal by operational weather
forecasting centers around the world [Moncrieff . (2007), Vitart Robertson
(2012)]. Nonlinearities in the climate system lead to an upscale transfer of
energy (and therefore error) from smaller to larger scales [Lorenz (1969),
Palmer (2001), Tribbia Baumhefner (2004)]. At the same time, slowly evolving
modes of variability can produce predictable signals on shorter timescales
[Hoskins (2013), Vannitsem Lucarini (2016)]. Under the ‘seamless prediction’
paradigm, the weather and climate communities should work together to develop
Earth-system models [Brunet . (2010), HM. Christensen Berner (2019)], as
developments made in one community are expected to benefit the other. The
development and use of stochastic parameterizations is a good example of this
paradigm at work.

Recent years have seen substantial interest in the development of stochastic
parameterization schemes. Pragmatic approaches, such as the Stochastically
Perturbed Parameterization Tendencies (SPPT) scheme [Buizza . (1999), Palmer,
Buizza . (2009)] are widely used due to their ease of implementation and
beneficial impacts on the model [Sanchez . (2016), Leutbecher . (2017), HM.
Christensen . (2017)]. Other schemes predict the statistics of model
uncertainty using a theoretical understanding of the atmospheric processes
involved, such as the statistics of convection [Craig Cohen (2006), Khouider .
(2010), Sakradzija Klocke (2018), Bengtsson . (2019)]. A third approach is to
make use of observations or high-resolution simulations to characterize
variability that is unresolved in a low-resolution forecast model [Shutts
Palmer (2007)]. This last approach can be used to develop data-driven
stochastic schemes [Dorrestijn . (2015), Bessac . (2019)] or to constrain
tunable parameters in stochastic parameterizations [Shutts Pallares (2014),
HM. Christensen . (2015), HM. Christensen (2019)]. A drawback of these data-
driven approaches is that assumptions are made about the structure of the
stochastic parameterization (e.g. the physical process to focus on, or the
distribution of the stochastic term conditioned on the resolved state) in
order to make the analysis tractable using conventional methods.

Machine learning models offer an approach to parameterize complex nonlinear
sub-grid processes in a potentially computationally efficient manner from data
describing those processes. The family of machine learning models consist of
mathematical models whose structure and parameters (often denoted weights)
optimize the predictive performance of _a priori_ unknown relationships
between input (“predictor”) and output (“predictand”) variables. Commonly used
machine learning model frameworks range in complexity from simple linear
regression to decision trees and neural networks. More complex methods allow
modelling of broader classes of predictor-predictand relationships. Machine
learning models minimize overfitting through the use of regularization
techniques that impose constraints on the model structure and weights.
Regularization is critical for more complex machine learning models, so that
they can converge to optimal and robust configurations in large parameter
spaces. For the parameterization problem, regularization and physical
constraints are critical for ensuring that the machine learning model
predictions match the distribution of observed values. Machine learning for
parameterizations has been considered since [Krasnopolsky . (2005)] and
recently multiple groups have begun developing new parameterizations for a
variety of processes [Schneider . (2017), Gentine . (2018), Rasp . (2018),
Bolton Zanna (2019)]. However, these schemes have focused exclusively on
deterministic parameterization approaches, but the need for stochastic
perturbations is being recognized [Brenowitz Bretherton (2019)].

One active area in current machine learning research is generative modeling,
which focuses on models that create synthetic representative samples from a
distribution of arbitrary complexity without the need for a parametric
representation of the distribution. Generative adversarial networks, or GANs
[Goodfellow . (2014)], are a class of generative models that consist of two
neural networks in mutual competition. The generator network produces
synthetic samples similar to the original training data from a latent vector,
and the critic, or discriminator, network examines samples from the generator
and the training data in order to determine if a sample is real or synthetic.
The critic acts as an adaptable loss function for the generator by learning
features of the training data and teaching those features to the generator
through back-propagation. The original GAN formulation used a latent vector of
random numbers as the only input to the generator, but subsequent work on
conditional GANs [Mirza Osindero (2014)] introduced the ability to incorporate
a combination of fixed and random inputs into the generator, enabling sampling
from conditional distributions. Because the stochastic parameterization
problem can be framed as sampling from the distribution of sub-grid tendencies
conditioned on the resolved state, conditional GANs have the potential to
perform well on this task.

The purpose of this study is to evaluate how well GANs can parameterize the
sub-grid tendency component of an atmospheric model at weather and climate
timescales. A key question is whether a GAN can learn uncertainty
quantification within the parameterization framework, removing the need to
retrospectively develop separate stochastic representations of model
uncertainty. While the ultimate goal is to test these ideas in a full general
circulation model (GCM: left for a future work), as a proof of concept we will
use the two-timescale model proposed in lorenz1996, henceforth the L96 system,
as a testbed for assessing the use of GAN in atmospheric models. Simple
chaotic dynamical systems such as L96 are useful for testing methods in
atmospheric modeling due to their transparency and computational cheapness.
The L96 system has been widely used as a testbed in studies including
development of stochastic parameterization schemes [Wilks (2005), Crommelin
Vanden-Eijnden (2008), Kwasniok (2012), Arnold . (2013)], data assimilation
methodology [Fertig . (2007), Law . (2016), Hatfield . (2018)], as well as
using ML approaches to learn improved deterministic parameterization schemes
[Schneider . (2017), Dueben Bauer (2018), Watson (2019)].

The evaluation consists of four primary questions. First, given inputs from
the “true” L96 model run, how closely does the GAN approximate the true
distribution of sub-grid tendency? Second, when an ensemble of L96 models with
stochastic GAN parameterizations are integrated forward to medium range
weather prediction time scales, how quickly does the prediction error increase
and how well does the ensemble spread capture the error growth? Third, when
the L96 model with a stochastic GAN parameterization is integrated out to
climate time scales, how well does the GAN simulation approximate the
properties of the true climate? Fourth, how closely does the GAN represent
both different regimes within the system and the probability of switching
between them?

Details of the Lorenz ’96 model and the GAN are presented in Section 2, and
the results of the weather and climate analyses described above are presented
in Section 3. Section 4 provides a discussion of GAN parameterization ‘health
risks’, while an overall discussion of results is presented in Section 5.
Conclusions follow in Section 6.

##  2 Methods

###  2.1 Lorenz ’96 Model Configuration

The L96 system was designed as a ‘toy model’ of the extratropical atmosphere,
with simplified representations of advective nonlinearities and multi-scale
interactions [Lorenz (1996)]. It consists of two scales of variables arranged
around a latitude circle. The large scale, low-frequency $X$ variables are
coupled to a larger number of small-scale high-frequency $Y$ variables, with a
two-way interaction between the $X$s and $Y$s. It is the interaction between
variables of different scales that makes the L96 system ideal for evaluating
new ideas in parameterization development. The L96 system has proven useful in
assessing new techniques that were later developed for use in GCMs [Crommelin
Vanden-Eijnden (2008), Dorrestijn . (2013)].

The $X$ and $Y$ variables evolve following:

$\displaystyle\frac{\mathrm{d}X_{k}}{\mathrm{d}t}$
$\displaystyle=-X_{k-1}(X_{k-2}-X_{k+1})-X_{k}+F-\frac{hc}{b}\sum_{j=J(k-1)+1}^{kJ}Y_{j};\;$
$\displaystyle k=1,...,K$ (1a)
$\displaystyle\frac{\mathrm{d}Y_{j}}{\mathrm{d}t}$
$\displaystyle=-cbY_{j+1}(Y_{j+2}-Y_{j-1})-cY_{j}+\frac{hc}{b}X_{(\mathrm{int}[(j-1)/J]+1)};\;$
$\displaystyle j=1,...,JK,$ (1b)

where in the present study the number of $X$ variables, $K=8$ and the number
of $Y$ variables per $X$ variable, $J=32$. Further, we set the coupling
constant, $h=1$, the spatial-scale ratio, $b=10$ and the temporal-scale ratio
$c=10$. The forcing term $F=20$ is set large enough to ensure chaotic
behavior. The chosen parameter settings, which were used in [Arnold . (2013)],
are such that one model time unit (MTU) is approximately equivalent to five
atmospheric days, deduced by comparing error-doubling times in L96 and the
atmosphere [Lorenz (1996)].

In this study the full Lorenz ‘96 equations are treated as the ‘truth’ which
must be forecast or simulated. In the case of the atmosphere, the physical
equations of motion of the system are known. However, due to limited
computational resources, it is not possible to explicitly simulate the
smallest scales, which are instead parameterized. Motivated by this
requirement for weather and climate prediction, a forecast model for the L96
system is constructed by truncating the model equations, and parameterizing
the impact of the small $Y$ scales on the resolved $X$ scales:

$\frac{\mathrm{d}X^{*}_{k}}{\mathrm{d}t}=-X^{*}_{k-1}(X^{*}_{k-2}-X^{*}_{k+1})-X^{*}_{k}+F-\hat{U}(X^{*}_{k},t);\;k=1,...,K,$
(2)

where $X^{*}_{k}(t)$ is the forecast estimate of $X_{k}(t)$ and
$\hat{U}(X^{*}_{k},t)$ is the parameterized sub-grid tendency. The
parameterization $\hat{U}$ approximates the true sub-grid tendencies:

$U(X,Y)=\frac{hc}{b}\sum_{j=J(k-1)+1}^{kJ}Y_{j},$ (3)

which can be estimated from realizations of the “truth” time series as

$U_{k}(t)=[-X_{k-1}(t)\left(X_{k-2}(t)-X_{k+1}(t)\right)-X_{k}(t)+F]-\left(\frac{X_{k}(t+dt_{f})-X_{k}(t)}{dt_{f}}\right).$
(4)

following Arnold2013. The time step $dt_{f}$ equals the time step used in the
forecast model for consistency.

A long “truth” run of the L96 model is performed to generate both training
data for the machine learning models and a test period for both weather and
climate evaluations. The “truth” run is integrated for 20000 MTU using a
fourth-order Runge-Kutta timestepping scheme and a time step $dt=0.001$ MTU.
Output from the first 2000 MTU are used for training, and the remaining 18000
MTU are used for testing. A burn-in period of 2 MTU is discarded. All
parameterized forecast models of the L96 use a forecast timestep of
$dt_{f}=0.005$ MTU and a second order Runge-Kutta timestepping scheme. This
scheme is to represent the time discretization of the equations representing
the resolved dynamics in an atmospheric forecasting model.

###  2.2 GAN Parameterizations

The GAN parameterization developed for the Lorenz ’96 model in this study
utilizes a conditional dense GAN to predict the sub-grid tendency at the
current time step given information about the state at the previous time step.
We will investigate two classes of predictors of $U_{t}$: both $X$ and $U$ at
the previous forecast timestep, and $X$ alone. In the following discussion we
focus on GANs based on the first of these predictor sets; the construction of
those associated with the second set is analogous to that of the first. Note
that in this section we move to a discrete-time notation, with the time index
indicated by the subscript, $t$, where adjacent indices are separated by the
forecast time step, $dt_{f}$.

The GAN generator accepts $X_{t-1,k}$, $U_{t-1,k}$, and a latent random vector
$Z_{t-1,k}$ as input to estimate $\hat{U}_{t,k}$, or the predicted $U$ at time
$t$. The discriminator accepts $X_{t-1,k}$, $U_{t-1,k}$, and $V_{t,k}$ as
inputs (where $V_{t,k}$ may be either $U_{t,k}$ if from the training data or
$\hat{U}_{t,k}$ if from the generator) and outputs the probability that
$V_{t,k}$ comes from the training data. All inputs and outputs are re-scaled
to have a mean of 0 and standard deviation of 1 based on the training data
distributions. Note that we choose to develop local GANs, i.e. GANs which
accept $X$ and $U$ for a given spatial index, $k$, and that predict $\hat{U}$
for that index, $k$, as opposed to GANs that accept vector $X$ and $U$ and
thus include spatial information. This is to match the local nature of
parameterization schemes in weather and climate models.

Figure 1: (Top) A diagram of how the GAN networks are connected for training.
(Bottom) A diagram of the GAN network architectures used for the stochastic
parameterization.

Each GAN we consider consists of the same neural network architecture with
variations in the inputs and how noise is scaled and inserted into the
network. A diagram of the architecture of the GAN networks is shown in Fig. 1.
Both the generator and discriminator networks contain two hidden layers with
16 neurons in each layer. The weights of the hidden layers are regularized
with a L2, or Ridge, penalty [Hoerl Kennard (1970)] with scale factor
$\lambda$ of 0.001. Scaled exponential linear unit (SELU) activation functions
[Klambauer . (2017)] follow each hidden layer. SELU is a variation of the
common Rectified Linear Unit (ReLU) activation function with a scaled
exponential transform for the negative values that helps ensure the output
distribution retains a mean of 0 and standard deviation of 1. Larger numbers
of neurons per hidden layer were evaluated and did not result in improved
performance. Gaussian additive noise layers before each hidden layer and
optionally the output layer inject noise into the latent representations of
the input data. A batch normalization [Ioffe Szegedy (2015)] output layer
ensures that the output values retain a mean of 0 and and standard deviation
of 1, which helps the generator converge to the true distribution faster.

The GAN training procedure iteratively updates the discriminator and generator
networks until the networks reach an adversarial equilibrium in which the
discriminator should not be able to distinguish “true” data from generator
samples. The inputs, outputs, and connections between networks are shown in
Fig. 1. A batch, or random subset of samples drawn without replacement from
the training data, of truth run output is split in half and one subset is fed
through the generator and then into the discriminator or is sent directly to
the discriminator. The discriminator is then updated based on errors in
predicting which samples have passed through the generator. Another batch of
samples are drawn and sent through the generator and then to a connected
discriminator with frozen weights. The gradient with respect to the generator
is calculated based on the reversed labeling that the generated samples
originate from the true training data. This reversal forces the discriminator
to teach the generator features that would worsen the discriminator’s own
performance. Gaussian noise is injected into the neural network at each
iteration through both the input and hidden layers. We consider hidden layer
Gaussian noise scaled to standard deviations at different orders of magnitude
in order to determine how the magnitude of the noise affects the forecast
spread and the representation of the model climate. In the training process,
neural network weights are updated based on subsamples of examples drawn
randomly from the training data without regard for temporal dependence. The
GAN is trained as a map between $X_{t-1,k}$ (and possibly $U_{t-1,k}$) and the
sub-grid-scale tendency $U_{t,k}$. In forecast mode, we test providing both
white, or uncorrelated noise, and red, or correlated noise to the GAN. The red
noise is generated using an AR(1) process with a correlation equal to the
lag-1 autocorrelation of the deterministic residuals of the GAN. The color of
the noise is not relevant during the training process: both white- and red-
noise representations are trained in the same way. The noise values are kept
constant through the integration of a single timestep. The difference between
the white- and red-noise representations only manifests when they are
incorporated as parameterizations in the full model (Eqn. 1).

The GANs are all trained with a consistent set of optimization parameters. The
GANs are updated through stochastic gradient descent with a batch size (number
of examples randomly drawn without replacement from the training data) of 1024
and a learning rate of 0.0001 with the Adam optimizer [Kingma Ba (2015)]. The
GANs are trained for 30 epochs, or passes through the training data. The model
weights are saved for analysis every epoch for the first 20 epochs and then
every 2 epochs between epochs 20 and 30. The GANs are developed with the Keras
v2.2 machine learning API coupled with Tensorflow v1.13.

The GAN configurations considered in this study are summarized in Table 1. A
short name of the format ‘predictors–noise magnitude–noise correlation’ is
introduced to simplify identification of different GANs. For example, ‘XU-
med-r’ refers to the GAN that takes X and U as predictors, and uses medium
(med) magnitude red (r) noise. While most GANs tested include the optional
additive noise layer before the output layer, the sensitivity to this choice
was also considered. GANs that do not include this final noise layer follow
the naming convention above, but are indicated by an asterisk.

Table 1: Summary of the GAN configurations tested. Short name | Input variables | Noise magnitude | Noise correlation | Output layer noise?
---|---|---|---|---
XU-lrg-w |  $X_{t-1,k}$, $U_{t-1,k}$ | 1 | white | yes
XU-med-w |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.1 | white | yes
XU-sml-w |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.01 | white | yes
XU-tny-w |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.001 | white | yes
X-med-w | $X_{t-1,k}$ | 0.1 | white | yes
X-sml-w | $X_{t-1,k}$ | 0.01 | white | yes
X-tny-w | $X_{t-1,k}$ | 0.001 | white | yes
XU-lrg-r |  $X_{t-1,k}$, $U_{t-1,k}$ | 1 | red | yes
XU-med-r |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.1 | red | yes
XU-sml-r |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.01 | red | yes
XU-tny-r |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.001 | red | yes
X-med-w | $X_{t-1,k}$ | 0.1 | red | yes
X-sml-r | $X_{t-1,k}$ | 0.01 | red | yes
X-tny-r | $X_{t-1,k}$ | 0.001 | red | yes
XU-lrg-w* |  $X_{t-1,k}$, $U_{t-1,k}$ | 1 | white | no
XU-med-w* |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.1 | white | no
XU-sml-w* |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.01 | white | no
XU-tny-w* |  $X_{t-1,k}$, $U_{t-1,k}$ | 0.001 | white | no
X-sml-w* | $X_{t-1,k}$ | 0.01 | white | no
X-tny-w* | $X_{t-1,k}$ | 0.001 | white | no

###  2.3 Polynomial Regression Parameterization

The GAN stochastic parameterization is evaluated against a cubic polynomial
regression parameterization, $\hat{U}_{t,k}$, similar to the model used in
Arnold2013.

$\displaystyle\hat{U}_{t,k}$ $\displaystyle=U^{d}_{t,k}+\epsilon_{t,k}$
$\displaystyle U^{d}_{t,k}$
$\displaystyle=aX_{t-1,k}^{3}+bX_{t-1,k}^{2}+cX_{t-1,k}+d$ (5)

The parameters $[a,b,c,d]$ are determined by a least squares fit to the
$(X,U)$ data from the L96 “truth” training run. It is known that the simple
cubic polynomial deterministic parameterization $U^{d}_{t,k}$ is a poor
forecast model for the L96 system [Wilks (2005), Arnold . (2013), H.
Christensen . (2015)], as $X$ does not uniquely determine $U$. The variability
in the $(X,U)$ relationship is accounted for using a temporally correlated
additive noise term:

$\displaystyle\epsilon_{t,k}$
$\displaystyle=\phi\epsilon_{t-1,k}+\sigma_{\epsilon}(1-\phi^{2})^{1/2}z_{t,k},$
(6)

where $z\sim\mathcal{N}(0,1)$, the first order autoregressive parameters
$(\phi,\sigma_{\epsilon})$ are estimated from the residual process
$r_{t}=U_{t}-U^{d}_{t}$, and the $\epsilon_{k}$ processes are independent for
different $X$ variables.

The polynomial parameterization has been specifically designed to represent
the impact of the $Y$ variables in this version of the L96 model, just as
traditional parameterization schemes are designed to represent a given process
in a GCM. Previous studies have demonstrated that the polynomial
parameterization with additive noise performs very well [Wilks (2005), Arnold
. (2013), H. Christensen . (2015)]. This ‘bespoke’ parameterization is
therefore a stringent benchmark against which to test GAN parameterizations.

##  3 Results

###  3.1 Metrics

The accuracy of ensemble weather forecasts can be summarized by evaluating the
root mean square error (RMSE) of the ensemble mean. The lower the RMSE, the
more accurate the forecast. The RMSE at lead time $\tau$ is defined as:

$RMSE(\tau)={\sqrt{\frac{1}{N}\sum\limits_{t=1}^{N}\left(X^{o}(t)-X^{m}(t_{\mathrm{init}},\tau)\right)^{2}}}\hskip
2.84526pt,$ (7)

where $N$ is the number of forecast-observation pairs, $X^{o}(t)$ is the
observed state at time $t$, and $X^{m}(t_{\mathrm{init}},\tau)$ is the
ensemble mean forecast at lead time $\tau$, initialized at
$t_{\mathrm{init}}$, such that $t=t_{\mathrm{init}}+\tau$.

If an ensemble forecast correctly accounts for all sources of uncertainty such
that the forecast of the spread of the ensemble and measured probabilities of
an event are statistically consistent, the forecast is said to be _reliable_
[Wilks (2011)]. In this study we assess the reliability of the ensemble using
the spread-error relationship [Leutbecher Palmer (2008), Leutbecher (2010)].
This states that, for an unbiased and reliable forecasting system, the root
mean square error in the ensemble mean is related to the average ensemble
variance:

$\frac{M}{M-1}\;\overline{\text{estimate ensemble
variance}}=\frac{M}{M+1}\;\overline{\text{mean square error}},$ (8)

where $M$ is the number of ensemble members, and the variance and mean error
have been estimated by averaging over many forecast-verification pairs. For
the large ensemble size used in this study, $M=40$, we can consider the
correction factor to be close to 1. A skilful probabilistic forecast will have
as small an RMSE as possible, while also demonstrating a statistical match
between RMSE and average ensemble spread.

The simplest definition of the ‘climate’ of the L96 system is the probability
density function (PDF) of the individual $X_{t,k}$ values. The climatological
skill can therefore be summarized by quantifying the difference between the
true and forecast PDF. The Hellinger distance $H$, is calculated for each
forecast model:

$H(p,q)=\frac{1}{2}\int\left(\sqrt{p(x)}-\sqrt{q(x)}\right)^{2}\mathrm{d}x,$
(9)

where $p(x)$ is the forecast PDF, and $q(x)$ is the verification PDF [Pollard
(2002)]. The smaller $H$, the closer the forecast pdf is to the true pdf. We
also considered the Kullback-Leibler ($KL$) divergence [Kullback Leibler
(1951)], motivated by information theory, but found it provided no additional
information over the Hellinger distance, so results for the $KL$ are not shown
for brevity.

Evidence that the L96 model displays distinct dynamical regimes of behavior
for the parameter set considered was presented in [H. Christensen . (2015)],
in which regime affiliations were determined using a metric based on the
temporally local covariance. [H. Christensen . (2015)] found that during the
more common regime (regime frequency $\sim$ 80%), the eight $X$ variables
exhibit wave-2 like behavior around the ring, while in the rarer regime, the
$X$ variables exhibit wave-1 type behavior. Another approach to characterizing
regime structure that makes use of both the instantaneous state and the recent
past of the system is Hidden Markov Model (HMM) analysis [Rabiner (1989),
Franzke . (2008), Monahan . (2015)]. In an HMM analysis, it is assumed that
underlying the observed state variables is an unobserved Markov chain taking
discrete values. The HMM algorithm provides maximum likelihood estimates of
the probability distributions of the state variables conditioned on the
instantaneous hidden state values, the stochastic matrix $Q$ of transition
probabilities for each time step, and an optimal estimate of the hidden state
sequence.

###  3.2 Offline assessment of GAN performance

Figure 2: Offline assessment of GAN performance. Hellinger distances between
the GAN sub-grid tendency distributions given input $X$ and $U$ values from
the truth run, and the truth run sub-grid forcing distribution as a function
of training epoch.

The GAN parameterizations are first evaluated on how closely their output sub-
grid forcing distributions match those of the truth run when the GANs are
supplied with input $X$ and $U$ values from the truth run. This is summarized
by the Hellinger distance in Figure 2. Most of the GANs show a trend of
decreasing Hellinger distance for the first few epochs followed by mostly
stable oscillations. GANs with both $X_{t-1,k}$ and $U_{t-1,k}$ as input tend
to perform better in the offline analysis than those with only $X_{t-1,k}$.
Larger input noise standard deviations seem to reduce the amount of
fluctuation in the Hellinger distance between epochs, but there does not
appear to be a consistent correlation with noise standard deviation and
Hellinger distance. Note that the weights as fitted at the end of epoch 30 are
used in the forecast networks, regardless of whether the GAN at this epoch
shows the minimum Hellinger distance.

###  3.3 GAN simulation of sub-grid-scale tendency distribution

The joint distributions of $X_{t-1}$ and $U_{t}$ from the different model runs
reveal how the noise standard deviation affects the model climate (Fig. 3).
Larger noise standard deviations increase the range of $X$ values appearing in
the run but do not appear to change the range of $U$ values output by the GAN.
The behavior of the XU-tny-r GAN devolved into oscillating between two
extremes. The X-only GANs did the best job in capturing the shape of the truth
distribution although they underestimated the variance at the extremes. While
the XU-w* and X-w* series captured the bulk of the distribution well, there
were spurious points outside the bounds of the truth distribution for all of
these models. The Polynomial model captured the conditional mean and variance
of the distribution well but smoothed over some of the subtle details of the
joint distribution that some of the GANs were able to capture.

Figure 3: Joint distributions (2D histograms) of $X_{t-1}$ and $U_{t}$ for
each GAN configuration. The truth joint distribution is overlaid in red
contours on each forecast model distribution.

###  3.4 Weather Evaluation

Figure 4: Summary of performance of different parameterized models (x-axis)
(a) Weather forecast performance. Ensemble spread (circles) and RMSE (crosses)
for experiments with white and red noise in GANs at timestep 201. The
horizontal dashed line indicates the RMSE for the polynomial forecast model.
Ideally, a forecast model will produce forecasts with small RMSE while
maintaining the match between spread and RMSE. (b) Climate performance. The
Hellinger distance between each forecast pdf and the ‘true’ pdf. The metric in
question is calculated for the best estimate of the climatological pdf,
averaging across all $X$ variables (circles), as well as for each X variable
in turn, e.g. comparing true and forecast $X_{1}$ pdfs, etc (crosses). The
latter gives an indication of the sampling uncertainty. The horizontal dashed
line indicates the mean value of $H$ for the polynomial model. Figure 5: RMSE
(lines with dots) and ensemble spread (dashed lines) for a subset of
experiments with white and red noise in GANs. Note that 400 forecast time
steps corresponds to 2 MTU, or 10 ‘atmospheric days’.

The parameterized models for the Lorenz ’96 system are evaluated in a weather
forecast framework. An extensive set of re-forecast experiments were produced
for 751 initial conditions selected from the attractor. An ensemble of 40
forecasts was produced from each initial condition (i.e. no initial condition
perturbations are used). Different random seeds are used for each ensemble
member to generate the stochastic perturbations used in the GAN or polynomial
parameterizations.

Figure 4 shows the RMSE and spread for all weather experiments at 1 MTU. A
reduction in RMSE indicates an ensemble forecast that more closely tracks the
observations. A good match between RMSE and ensemble spread indicates a
reliable forecast. The best performing GANs in terms of RMSE are X-tny-r,
X-tny-w, and X-tny-w*. All of these models performed slightly better than the
polynomial regression, which was competitive with most GANs in terms of both
RMSE and the ratio of RMSE to spread. The spread of the white noise GANs is
generally underdispersive. Most of the red noise GANs, on the other hand, are
somewhat overdispersive with X-med-r having the spread/error ratio closest to
1. Figure 5 shows the RMSE and ensemble spread for a subset of the ensemble
forecasts of the $X$ variables performed using the GAN parameterizations or
the bespoke polynomial parameterization. X-sml-w* demonstrates both low RMSE
and similar spread to RMSE throughout the forecast period. X-sml-r and X-sml-w
feature similar RMSE through 1 MTU, but X-sml-r has smaller RMSE after that
point and a better spread-to-error ratio throughout the period. The XU GANs
have higher RMSEs than their X counterparts because the XU models may have
overfit to the strong correlation between $U_{t-1,k}$ and $U_{t,k}$ in the
training data. Inspection of the input weights revealed that XU GANs generally
weigh $U_{t-1,k}$ more highly than $X_{t-1,k}$. At maxima and minima in the
waves, the XU models may be biased toward extending the current growth
forward, which can be a source of error in forecast runs.

###  3.5 Climate Evaluation

Figure 6: The skill of the forecast models at reproducing the climate of the
Lorenz ‘96 system defined as the pdf of the $X$ variables. (a) The pdf of the
Lorenz ‘96 system (black dashed) compared to a subset of the forecast models.
(b) The difference between the forecast and true pdfs shown in (a). Sampling
variability is indicated by shading the range in metrics for each of the 8 $X$
variables of the full Lorenz ‘96 system.

The GAN parameterizations are also tested on their ability to characterize the
climate of the Lorenz ‘96 system. First, the ability to reproduce the pdf of
the $X$ variables is evaluated. Each forecast model and the full L96 system
were used to produce a long simulation of length 10,000 MTU. Figure 6(a) shows
kernel density estimates of the marginal pdfs of $X_{t,k}$ from the full L96
system and from a sample of parameterized models. The pdf of the true L96
system is markedly non-Gaussian, with enhanced density forming a ‘hump’ at
around $X=8$. Compared to the true distribution, the XU-sml-w and XU-sml-r
models both perform poorly, producing very similar pdfs with too large a
standard deviation and that are too symmetric. However, the other models shown
skilfully reproduce the true pdf. Figure 6(b) shows the difference between
each forecast pdf and the true pdf. Several GANs perform as well if not better
than the benchmark bespoke polynomial parameterization.

Figure 4(b) shows the Hellinger distance $H$ evaluated for each parameterized
model. The filled circles indicate the value of the metric when the pdf is
evaluated across all $X$ variables in both parameterized and truth timeseries.
The crosses give an indication of sampling variability, and indicate the
metrics comparing pairs of $X$ variables, i.e. $X_{t,j}$ and $X_{t,k}$ for
$j\neq k$. Quantifying the parameterized model performance in this way allows
for easy ranking of the different models. While the AR(1) stochastic
polynomial parameterized forecast model is very skillful [Arnold . (2013)],
several GANs outperform the polynomial model.

Figure 7: The skill of the parameterized models (colors) at reproducing the
‘true’ (a) spatial correlation and (b) temporal correlation of the $X$
variables in the Lorenz ‘96 system (black), calculated from the climatological
simulation. The sampling variability in these metrics, as indicated by the
variability between the metric calculated for different $X$ variables, is
narrower than the plotted line width.

In addition to correctly capturing the distribution of the $X$ variables, it
is desirable that a parameterized model will capture the spatio-temporal
behavior of the system. This is assessed by considering the spatial and
temporal correlations in both the true system and parameterized models. The
diagnostic is shown for a subset of the tested paramererized models in Figure
7. It is evident that the parameterized models that skillfully capture the pdf
of $X$ also skillfully represent the spatio-temporal characteristics of the
system. The X-sml-w* scheme performs particularly well, improving over the
stochastic polynomial approach.

Figure 8: Joint distributions of the projections of $X$ on wavenumbers 1 and 2
($|Z_{1}|$ and $|Z_{2}|$ respectively) conditioned on HMM regime occupation
(red and blue contours) for the truth simulation subset of parameterized
simulations. In all but the upper left panel, thin black lines display the
unconditional joint distribution. In all panels the grey curves denote the
full joint distribution (without regime occupation conditioning) from the
truth simulation. The conditional distributions have been scaled by the
relative probabilities of each state. Inset: HMM stochastic matrix $Q$.

Following the regime results presented in [H. Christensen . (2015)], we use
HMM analysis to classify into two clusters the instantaneous states in the
4-dimensional space spanned by the norms of the projections of $X$ on
wavenumbers one through four (denoted $Z_{j},j=1,...,4$). Because of the
spatial homogeneity of the $X_{t,k}$ statistics, these wavenumber projections
correspond to the Empirical Orthogonal Functions.

Kernel density estimates of the joint pdfs of the projections of $X$ on
wavenumbers 1 and 2 are presented in Figure 8, along with estimates of the
joint distributions conditioned on the HMM regime sequence. The conditional
distributions have been scaled by the probability of each state so that the
full joint pdf is the sum of the conditional pdfs. For reference, the
unconditional joint pdf for truth is shown in gray contours in each panel. The
stochastic matrix $Q$ shows the probability of remaining in each regime
(diagonal values) or transitioning from one regime to the other (off-diagonal
values). As in Figures 6 and 7, only a subset of results are displayed.

The clear separation of the truth simulation into two distinct regimes is
modeled well by the polynomial parameterization and most of the GAN
parameterizations. With the exception of XU-sml-w and XU-sml-r, the regime
spatial structures and stochastic matrices are captured well. The GANs are
slightly more likely to transition between regimes than the truth run, while
the polynomial run is slightly more likely to stay in the same regime.
Consistent with the other climate performance results presented above, the
joint distribution of $|Z_{i}|$ produced by XU-sml-w and XU-sml-r are strongly
biased, with no evidence of a meaningful separation into two distinct regimes
of behavior.

To quantify the forecast model skill at reproducing the L96 behavior in
wavenumber space, Figure 9 shows $H$ calculated over the marginal pdfs of
$|Z_{1}|$ and $|Z_{2}|$ as upward and downward triangles respectively. Several
of the GAN evaluated are competitive with, or improve upon, the polynomial
parameterization scheme. The best performing GANs also performed the best in
terms of other climate metrics (e.g. see Figure 4).

We note that in particular, models which accurately capture the regime
behavior will also show good correlation statistics when averaged over a long
time series. The regime analysis can help diagnose why a model shows poor
correlation statistics. For example, X-sml-w* accurately captures the marginal
distributions of the two regimes, but the frequency of occurrence of the red
regime, dominated by wavenumber-1 behavior, is slightly too high (at 54%
compared to 51% in the ‘truth’ run). This reduces the magnitude of the
negative (positive) correlation at a lag of 0.75 (1.5) MTU observed in figure
11.

Figure 9: The Hellinger distance between each forecast pdf and the ‘true’ pdf,
considering the projection of the $X$ variables onto (upward triangles)
wavenumber 1, and (upward triangles) wavenumber 2.

###  3.6 Wavelet Analysis

Figure 10: (Top) Wavelet power spectra for the white and red noise GAN climate
runs as well as the polynomial and truth runs. (Bottom) Mean absolute percent
differences between the truth wavelet power spectrum and each model.

To further investigate why the white noise and red noise GANs differ in
operational performance, a wavelet analysis is performed on time series of
outputs from the climate runs. A discrete wavelet transform decomposes the
time series into contributions from different periods. The total energy $E$
for a given period is represented as

$E=\frac{1}{T}\sum_{t=1}^{T}w_{t}^{2}$ (10)

where $w$ is the wavelet magnitude at a given timestep $t$. The total power
for each period from each model is shown in Fig. 10. All sml GANs except the
XU-sml-r and XU-sml-w follow the truth power curve closely. follow The
polynomial regression follows the truth closely although it tends to
underestimate the power slightly for each period. The GANs peak in power at
the same period when the temporal correlation in Fig. 7 crosses 0. The GANs
with poor Hellinger distances also contained more energy for longer periods.

A clearer evaluation of the wavelet differences can be found by calculating
the mean absolute percentage difference from the truth run at different
wavelengths. The absolute difference between the truth and parameterized runs
increases with increasing wavelengths, so the percentage difference is
employed to control for this trend:

$MAPD=\frac{1}{T}\sum_{t=1}^{T}\frac{|E_{g,t}-E_{u,t}|}{E_{u,t}}$ (11)

The MAPD scores with wavelength in Fig. 10 shows that none of the GANs
consistently perform better at all periods, but some do provide closer matches
to the truth spectrum for the longer periods. In the peak energy period, the
different GANs have minimum error for slightly different periods before
increasing in error again. The X-sml-r GAN uniquely shows decreasing MAPD with
increasing period, while the white noise GANs generally have similar
differences across the range of evaluated periods.

##  4 GAN Health Risks

The disparity between the offline verification statistics and those from the
climate and weather runs highlights the challenges in training GANs for
parameterization tasks. Neither the values of the generator loss or the
offline evaluation of the GAN samples correlated with their performance in the
forecast model integrations. The generator and discriminator networks optimize
until they reach an equilibrium, but there is no guarantee that the
equilibrium is stable or optimal. Some of the differences in the results may
be due to particular GANs converging on a poor equilibrium state as opposed to
other factors being tested. GANs and other machine learning parameterization
models are trained under the assumption that the data are independent and
identically distributed, but in practice are applied to spatially- and
temporally-correlated fields sequentially, potentially introducing nonlinear
feedback effects. GANs are more complex to train than other machine learning
methods because they require two neural networks and do not output an
informative loss function. Larger magnitude additive noise appears to help
prevent runaway feedbacks from model biases at the expense of increasing
weather prediction errors. The inclusion of the batch normalization output
layer appeared to assist both training and prediction by limiting the possible
extremes reached during integration.

##  5 Discussion

In this study, we chose to focus on GANs for stochastic parameterization
primarily because the framework offers a way to directly embed stochasticity
into the model instead of adding it a posteriori to a deterministic
parameterization. The use of the discriminator network as an adaptive loss
function is also attractive because it reduces the need for developing a hand-
crafted loss function and can be scaled to higher dimensional and more complex
outputs.

Several of the GANs tested show a weather and climate skill that is
competitive with a bespoke polynomial parameterization scheme. For
climatological skill, several different metrics were considered including the
distribution of the $X$ variable, spatio-temporal correlation statistics, and
regime behavior. We found that forecast models that performed well according
to one metric also performed consistently well for all metrics. The good
performance of the GAN is encouraging, demonstrating that GANs can indeed be
used as explicit stochastic parameterizations of uncertain sub-grid processes
directly from data. Furthermore, a small number of GANs improve upon the
bespoke polynomial approach, indicating the potential for such machine-learned
approaches to improve on our current hand-designed parameterization schemes,
given suitable training data. While this has been proposed and demonstrated
for deterministic parameterizations [Schneider . (2017), Rasp . (2018)], this
is a first demonstration for the case of an explicitly stochastic
parameterization.

Figure 11: Correlation between weather forecast skill and climate performance.
(a) Weather forecast RMSE versus climate Hellinger distance. (b) Weather
forecast spread-error ratio versus climate Hellinger distance. Colors indicate
forecast model, as in Figure 4

Comparison of Figure 4(a) and (b) indicates a relationship between forecast
models that perform well on weather and climate timescales. To quantify this
further, Figure 11 shows the correlation between weather forecast skill and
climate performance for each model considered. Models which produce weather
forecasts with a lower RMSE also show good statistics on climate timescales.
In contrast, the _reliability_ of weather forecasts, i.e., the statistical
match between spread and error, is a poor predictor of climate performance.
This is reflected by the competitive performance of the white noise GAN for
producing a realistic climate, whereas on weather timescales, red noise
increases the spread and can thereby substantially improve the forecast
reliability (e.g. consider the X-med, X-sml and X-tny GAN for white and red
noise respectively: Figure 4). This relationship between initialized and
climatological performance has been discussed in the context of global models
[KD. Williams . (2013)]. It provides further evidence that parameterizations
can first be tested in weather forecasts before being used in climate models,
as promoted by the ‘seamless prediction’ framework [Brunet . (2010)].

Wavelet analysis helped uncover differences in model performance across
different time scales. While the other evaluation metrics focused on
distributional or error metrics in the time domain, the wavelet power spectrum
separated the time series into different periods and enabled comparisons of
the energy embedded in different scales. In particular, the wavelet analysis
revealed that some of the GANs added energy to the system either at long
periods or across all periods in some cases.

The standard deviation of the noise does impact both the training of the GANs
and the resulting weather and climate model runs. Using too large a standard
deviation limits the ability of the GAN to discover the structure of the true
distribution of the data. Standard deviations that are too small may result in
either the generator or discriminator becoming overly good at their tasks
during training, which results in the GAN equilibrium being broken. During
simulations, noise standard deviations that are too small can result in the
system becoming trapped within one regime and never escaping.

In addition to the GAN configurations evaluated in the paper, other GAN
settings were tested and were found to have similar or worse performance.
Given the relative simplicity of the Lorenz ’96 system, adding neurons in each
hidden layer did not improve performance. Using the SELU activation function
generally resulted in faster equilibrium convergence than the ReLU. Varying
the scaling factor for the L2 regularization on each hidden layer did affect
model performance. Using a larger value greatly reduced the variance of the
predictions, but using a smaller value resulted in peaks of the final
distributions being too far apart, especially when both $X_{t-1,k}$ and
$U_{t-1,k}$ were used as inputs. We also tested deriving $U$ from a 1D
convolutional GAN that reproduced the set of $Y$ values associated with each
$X$. That approach did produce somewhat realistic $Y$ values but contained
”checkerboard” artifacts from the convolutions and upscaling, especially at
the edges. The sum of the $Y$ values was also not equal to $U$ derived from
Eq. 5.

The L96 system is commonly used as a testbed for new ideas in
parameterization, and ideas tested using the system can be readily developed
further for use in higher complexity Earth system models. However, the L96
system has many fewer dimensions than an Earth system model and a relatively
simple target distribution. The relative simplicity of the L96 system may have
also led to the more complex GAN overfitting to the data compared with the
simpler polynomial parameterization. For more complex, higher dimensional
systems, the extra representational capacity of the GAN may provide more
benefit than can be realized in L96. The computational simplicity of L96 also
allows for the production of extremely large training data sets with little
compute resources. Higher complexity Earth system model output can provide
training set coverage spatially but will be limited temporally by the amount
of computational resources available.

##  6 Conclusions

In this study, we have developed an explicitly stochastic GAN framework for
the parameterization of sub-grid processes in the Lorenz ’96 dynamical system.
After testing a wide range of GANs with different noise characteristics, we
identified a subset of models that outperform a benchmark bespoke cubic
polynomial parameterization. Returning to the questions posed in the
Introduction, we found that this subset of GANs approximates well the joint
distribution of the resolved state and the sub-grid-scale tendency. This model
subset also produces the most accurate weather forecasts (i.e. lowest RMSE in
the ensemble mean). Some GANs produce reliable forecasts, in which the
ensemble spread is a good indicator of the error in the ensemble mean.
However, these are not necessarily the GANs that produce the most accurate
forecasts. The subset of models with the most accurate weather forecasts
produce the most accurate climate simulations, as characterized by probability
distributions, space and time correlations, and regime dynamics. However, we
note that the GANs which produce skilful weather and climate forecasts were
different to those which performed well in “offline” mode.

Although the GANs required an iterative development process to maximize model
performance and were very sensitive to the noise magnitude and other
hyperparameter settings, they do show promise as an approach for stochastic
parameterization of physical processes in more complex weather and climate
models. Applying other recently developed GAN losses and regularizers [Kurach
. (2019)] could help reduce the chance for the GAN to experience a failure
mode.

The experiments presented here demonstrate that GANs are a promising machine
learning approach for developing stochastic parameterization in complex GCMs.
Key lessons learned and unanswered questions include:

  * •

While including the tendency from the previous timestep provides a natural
approach for building temporal dependence into the parameterization, it can
lead to accumulation of error in the forecast, such that local-in-time
parameterization should also be considered.

  * •

Autocorrelated noise is important for a skillful weather forecast, but appears
less important for capturing the climatological distribution.

  * •

It is possible that spatial correlations are also important in a higher
complexity Earth system model, which could not be assessed here due to the
simplicity of the Lorenz ‘96 system.

  * •

It is possible that the noise characteristics could also be learned by the GAN
framework to automate the tuning of the stochasticity.

Future work will use these lessons to develop machine-learned stochastic
parameterization schemes for use in higher complexity Earth system models.
GANs of a similar level of complexity to those used for L96 could emulate
local effects, such as some warm rain formation processes. Other generative
neural network frameworks, such as variational autoencoders, should also be
investigated to determine if they can provide similar or better performance
with a less sensitive training process.

###### Acknowledgements.

This research started in a working group supported by the Statistical and
Applied Mathematical Sciences Institute (SAMSI). DJG and HMC were funded by
National Center for Atmospheric Research Advanced Study Program Postdoctoral
Fellowships and by the National Science Foundation through NCAR’s Cooperative
Agreement AGS-1852977. HMC was funded by Natural Environment Research Council
grant number NE/P018238/1. AHM acknowledges support from the Natural Sciences
and Engineering Research Council of Canada (NSERC) and thanks SAMSI for
hosting him in the autumn of 2017. The NCAR Cheyenne supercomputer was used to
generate the model runs analyzed in this paper. Sue Ellen Haupt and Joseph
Tribbia provided helpful reviews of the paper prior to submission. The source
code for this paper can be viewed and downloaded at
https://github.com/djgagne/lorenz_gan/.

## References

  * Ajayamohan . (2013) Ajayamohan2013Ajayamohan, RS., Khouider, B. Majda, AJ.  2013\.  Realistic initiation and dynamics of the Madden-Julian Oscillation in a coarse resolution aquaplanet GCM Realistic initiation and dynamics of the Madden-Julian Oscillation in a coarse resolution aquaplanet GCM.  Geophysical Research Letters406252–6257.  https://doi.org/10.1002/2013GL058187
  * Arnold . (2013) Arnold2013Arnold, HM., Moroz, IM. Palmer, TN.  2013\.  Stochastic parameterizations and model uncertainty in the Lorenz ’96 system Stochastic parameterizations and model uncertainty in the lorenz ’96 system.  Phil. Trans. R. Soc. A371.  http://dx.doi.org/10.1098/rsta.2011.0479
  * Bengtsson . (2019) bengtsson2019Bengtsson, L., Bao, J., Pegion, P., Penland, C., Michelson, S. Whitaker, J.  2019\.  A model framework for stochastic representation of uncertainties associated with physical processes in NOAA’s Next Generation Global Prediction System (NGGPS) A model framework for stochastic representation of uncertainties associated with physical processes in NOAA’s next generation global prediction system (NGGPS).  Mon. Wea. Rev..  in press  doi:10.1175/MWR-D-18-0238.1
  * Berner . (2017) berner2017Berner, J., Achatz, U., Batte, L., De La Cámara, A., Christensen, H., Colangeli, M.Yano, YI.  2017\.  Stochastic Parameterization: Towards a new view of Weather and Climate Models Stochastic parameterization: Towards a new view of weather and climate models.  Bull. Amer. Meteor. Soc.98565–588.
  * Berner . (2012) Berner:2012gxBerner, J., Jung, T. Palmer, TN.  201207\.  Systematic Model Error: The Impact of Increased Horizontal Resolution versus Improved Stochastic and Deterministic Parameterizations Systematic Model Error: The Impact of Increased Horizontal Resolution versus Improved Stochastic and Deterministic Parameterizations.  J. Clim.25144946–4962.
  * Berner . (2010) Berner:2010fnBerner, J., Shutts, GJ., Leutbecher, M. Palmer, TN.  201002\.  A Spectral Stochastic Kinetic Energy Backscatter Scheme and Its Impact on Flow-Dependent Predictability in the ECMWF Ensemble Prediction System A Spectral Stochastic Kinetic Energy Backscatter Scheme and Its Impact on Flow-Dependent Predictability in the ECMWF Ensemble Prediction System.  Journal of the Atmospheric Sciences663603–626.
  * Bessac . (2019) bessac2019Bessac, J., Monahan, AH., Christensen, HM. Weitzel, N.  2019\.  Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes Stochastic parameterization of subgrid-scale velocity enhancement of sea surface fluxes.  Mon. Wea. Rev.1471447–1469.
  * Bolton Zanna (2019) Bolton2019Bolton, T. Zanna, L.  2019\.  Applications of Deep Learning to Ocean Data Inference and Subgrid Parameterization Applications of deep learning to ocean data inference and subgrid parameterization.  Journal of Advances in Modeling Earth Systems110376-399.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018MS001472 10.1029/2018MS001472
  * Brenowitz Bretherton (2019) Brenowitz2019Brenowitz, N. Bretherton, C.  2019\.  Spatially Extended Tests of a Neural Network Parameterization Trained by Coarse-Graining Spatially extended tests of a neural network parameterization trained by coarse-graining.  arXiv preprint.  https://arxiv.org/abs/1904.03327
  * Brunet . (2010) brunet2010Brunet, G., Shapiro, M., Hoskins, B., Moncrieff, M., Dole, R., Kiladis, GN.Shukla, J.  2010\.  Collaboration of the Weather and Climate Communities to Advance Subseasonal-to-Seasonal Prediction Collaboration of the weather and climate communities to advance subseasonal-to-seasonal prediction.  Bull. Amer. Meteor. Soc.91101397-1406.  10.1175/2010BAMS3013.1
  * Buizza . (1999) Buizza1999Buizza, R., Milleer, M. Palmer, TN.  1999\.  Stochastic representation of model uncertainties in the ECMWF ensemble prediction system Stochastic representation of model uncertainties in the ECMWF ensemble prediction system.  Q. J. R. Meteor. Soc.1252887-2908.  https://doi.org/10.1002/qj.49712556006
  * Buizza Palmer (1995) buizza1995Buizza, R. Palmer, TN.  1995\.  The Singular-Vector Structure of the Atmospheric Global Circulation The singular-vector structure of the atmospheric global circulation.  J. Atmos. Sci.5291434–1456.
  * H. Christensen . (2015) cmp15Christensen, H., Moroz, I. Palmer, T.  2015\.  Simulating weather regimes: impact of stochastic and perturbed parameter schemes in a simple atmospheric model Simulating weather regimes: impact of stochastic and perturbed parameter schemes in a simple atmospheric model.  Climate Dynamics447-82195-2214.
  * HM. Christensen (2019) christensen2019bChristensen, HM.  2019\.  Constraining stochastic parametrisation schemes using high-resolution simulations Constraining stochastic parametrisation schemes using high-resolution simulations.  Q. J. R. Meteor. Soc..  submitted
  * HM. Christensen Berner (2019) christensen2019Christensen, HM. Berner, J.  2019\.  From reliable weather forecasts to skilful climate response: a dynamical systems approach From reliable weather forecasts to skilful climate response: a dynamical systems approach.  Quarterly Journal of the Royal Meteorological Society.  in press
  * HM. Christensen . (2017) christensen2017Christensen, HM., Berner, J., Coleman, D. Palmer, TN.  2017\.  Stochastic parametrisation and the El Niño-Southern Oscillation Stochastic parametrisation and the El Niño-Southern Oscillation.  J. Clim.30117–38.
  * HM. Christensen . (2015) christensen2015dChristensen, HM., Moroz, IM. Palmer, TN.  2015\.  Stochastic and Perturbed Parameter Representations of Model Uncertainty in Convection Parameterization Stochastic and perturbed parameter representations of model uncertainty in convection parameterization.  J. Atmos. Sci.7262525–2544.
  * Craig Cohen (2006) craig2006Craig, GC. Cohen, BG.  2006\.  Fluctuations in an Equilibrium Convective Ensemble. Part I: Theoretical Formulation Fluctuations in an equilibrium convective ensemble. part i: Theoretical formulation.  J. Atmos. Sci.6381996-2004.
  * Crommelin Vanden-Eijnden (2008) crommelin2008Crommelin, D. Vanden-Eijnden, E.  2008\.  Subgrid-Scale Parametrisation with Conditional Markov Chains Subgrid-scale parametrisation with conditional markov chains.  J. Atmos. Sci.6582661–2675.
  * Davini . (2017) davini2017Davini, P., von Hardenberg, J., Corti, S., Christensen, HM., Juricke, S., Subramanian, A.Palmer, TN.  2017\.  Climate SPHINX: evaluating the impact of resolution and stochastic physics parameterisations in the EC-Earth global climate model Climate sphinx: evaluating the impact of resolution and stochastic physics parameterisations in the ec-earth global climate model.  Geosci. Model Dev.1031383–1402.
  * Dawson Palmer (2015) Dawson2015Dawson, A. Palmer, TN.  2015\.  Simulating weather regimes: impact of model resolution and stochastic parameterization Simulating weather regimes: impact of model resolution and stochastic parameterization.  Climate Dynamics442177-2193.  https://doi.org/10.1007/s00382-014-2238-x
  * Dorrestijn . (2013) dorrestijn2013Dorrestijn, J., Crommelin, DT., Biello, JA. Böing, SJ.  2013\.  A data-driven multi-cloud model for stochastic parametrization of deep convection A data-driven multi-cloud model for stochastic parametrization of deep convection.  Phil. Trans. R. Soc. A3711991.
  * Dorrestijn . (2015) dorrestijn2015Dorrestijn, J., Crommelin, DT., Siebesma, AP., Jonker, HJJ. Jakob, C.  2015\.  Stochastic parameterization of convective area fractions with a multicloud model inferred from observational data Stochastic parameterization of convective area fractions with a multicloud model inferred from observational data.  J. Atmos. Sci.72854–869.
  * Dueben Bauer (2018) dueben2018Dueben, PD. Bauer, P.  2018\.  Challenges and design choices for global weather and climate models based on machine learning Challenges and design choices for global weather and climate models based on machine learning.  Geosci. Model Dev113999–4009.  doi: 10.5194/gmd-2018-148
  * Fertig . (2007) fertig2007Fertig, EJ., Harlim, J. Hunt, BR.  2007\.  A comparative study of 4D-VAR and a 4D Ensemble Kalman Filter: Perfect model simulations with Lorenz-96 A comparative study of 4D-VAR and a 4D Ensemble Kalman Filter: Perfect model simulations with Lorenz-96.  Tellus, Series A: Dynamic Meteorology and Oceanography59196–100.
  * Franzke . (2008) fcfm08Franzke, C., Crommelin, D., Fischer, A. Majda, AJ.  2008\.  A hidden Markov model perspective on regimes and metastability in atmospheric flows A hidden Markov model perspective on regimes and metastability in atmospheric flows.  J. Clim.211740-1757.
  * Gentine . (2018) Gentine2018Gentine, P., Pritchard, M., Rasp, S., Reinaudi, G. Yacalis, G.  2018\.  Could Machine Learning Break the Convection Parameterization Deadlock? Could machine learning break the convection parameterization deadlock?  Geophysical Research Letters45115742-5751.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018GL078202 10.1029/2018GL078202
  * Goodfellow . (2014) Goodfellow2014Goodfellow, IJ., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S.Bengio, Y.  20146\.  Generative Adversarial Networks Generative Adversarial Networks.  Advances in Neural Information Processing Systems 27 Advances in neural information processing systems 27 ( 2672–2680).  http://papers.nips.cc/paper/5423-generative-adversarial-nets.pdf
  * Hatfield . (2018) hatfield2018Hatfield, S., Subramanian, A., Palmer, T. Düben, P.  2018\.  Improving weather forecast skill through reduced precision data assimilation Improving weather forecast skill through reduced precision data assimilation.  Mon. Wea. Rev.14649–62.  https://doi.org/10.1175/MWR-D-17-0132.1
  * Hoerl Kennard (1970) Hoerl1970Hoerl, AE. Kennard, RW.  1970\.  Ridge regression: Biased estimation for nonorthogonal problems Ridge regression: Biased estimation for nonorthogonal problems.  Techometrics1255–67.  https://doi.org/10.1080/00401706.1970.10488634
  * Hoskins (2013) hoskins2013Hoskins, B.  2013\.  The potential for skill across the range of the seamless weather-climate prediction problem: a stimulus for our science The potential for skill across the range of the seamless weather-climate prediction problem: a stimulus for our science.  Quarterly Journal of the Royal Meteorological Society139672573–584.  http://dx.doi.org/10.1002/qj.1991 10.1002/qj.1991
  * Hurrell . (2009) Hurrell2009Hurrell, J., Meehl, GA., Bader, D., Delworth, TL., Kirtman, B. Wielicki, B.  2009\.  A Unified Modeling Approach to Climate System Prediction A unified modeling approach to climate system prediction.  Bull. Amer. Meteor. Soc.901819–1832.  https://doi.org/10.1175/2009BAMS2752.1
  * Ioffe Szegedy (2015) Ioffe2015Ioffe, S. Szegedy, C.  2015\.  Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift Batch normalization: Accelerating deep network training by reducing internal covariate shift.  arXiv Preprint.  https://arxiv.org/abs/1502.03167
  * Isaksen . (2010) isaksen2010Isaksen, L., Bonavita, M., Buizza, R., Fisher, M., Haseler, J., Leutbecher, M. Raynaud, L.  2010\.  Ensemble of data assimilations at ECMWF Ensemble of data assimilations at ECMWF 636.  Shinfield park, ReadingEuropean Centre for Medium-Range Weather Forecasts.
  * Juricke Jung (2014) Juricke2014Juricke, S. Jung, T.  2014\.  Influence of stochastic sea ice parametrization on climate and the role of atmosphere–sea ice–ocean interaction Influence of stochastic sea ice parametrization on climate and the role of atmosphere–sea ice–ocean interaction.  Phil. Trans. R. Soc. A3722018.  https://www.doi.org/10.1098/rsta.2013.0283
  * Khouider . (2010) khouider2010Khouider, B., Biello, J. Majda, AJ.  2010\.  A Stochastic Multicloud Model for Tropical Convection A stochastic multicloud model for tropical convection.  Commun. Math. Sci.81187–216.
  * Kingma Ba (2015) Kingma2014Kingma, DP. Ba, J.  2015\.  Adam: A Method for Stochastic Optimization Adam: A method for stochastic optimization.  Proceedings of the 3rd International Conference for Learning Representations. Proceedings of the 3rd international conference for learning representations.  San Diego, CA.  https://arxiv.org/abs/1412.6980
  * Klambauer . (2017) Klambauer2017Klambauer, G., Unterthiner, T., Mayr, A. Hochreiter, S.  2017\.  Self-Normalizing Neural Networks Self-normalizing neural networks.  Advances in Neural Information Processing Systems (NIPS). Advances in neural information processing systems (nips).  https://arxiv.org/abs/1706.02515
  * Krasnopolsky . (2005) Krasnopolsky2005Krasnopolsky, VM., Fox-Rabinovitz, MS. Chalikov, DV.  2005\.  New Approach to Calculation of Atmospheric Model Physics: Accurate and Fast Neural Network Emulation of Longwave Radiation in a Climate Model New approach to calculation of atmospheric model physics: Accurate and fast neural network emulation of longwave radiation in a climate model.  Mon. Wea. Rev.13351370-1383.  https://doi.org/10.1175/MWR2923.1 10.1175/MWR2923.1
  * Kullback Leibler (1951) kullback1951Kullback, S. Leibler, RA.  1951\.  On Information and Sufficiency On information and sufficiency.  The Annals of Mathematical Statistics2279–86.  http://www.jstor.org/stable/2236703
  * Kurach . (2019) Kurach2019Kurach, K., Lučić, M., Zhai, X., Michalski, M. Gelly, S.  2019\.  A Large-Scale Study on Regularization and Normalization in GANs A large-scale study on regularization and normalization in gans.  International Conference on Machine Learning. International conference on machine learning.  https://arxiv.org/abs/1807.04720
  * Kwasniok (2012) kwasniok2012Kwasniok, F.  2012\.  Data-based stochastic subgrid-scale parametrisation: an approach using cluster-weighted modelling Data-based stochastic subgrid-scale parametrisation: an approach using cluster-weighted modelling.  Phil. Trans. R. Soc. A37019621061–1086.
  * Law . (2016) law2016Law, KJ., Sanz-Alonso, D., Shukla, A. Stuart, AM.  2016\.  Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators.  Physica D: Nonlinear Phenomena3251–13.
  * Leith (1975) Leith:1975jcLeith, CE.  1975\.  Numerical weather prediction Numerical weather prediction.  Reviews of Geophysics133681.
  * Leutbecher (2010) leutbecher2010Leutbecher, M.  2010\.  Diagnosis of Ensemble Forecasting Systems Diagnosis of ensemble forecasting systems.  Seminar on Diagnosis of Forecasting and Data Assimilation Systems, 7 - 10 September 2009 Seminar on diagnosis of forecasting and data assimilation systems, 7 - 10 september 2009 ( 235–266).  Shinfield Park, ReadingECMWF.
  * Leutbecher . (2017) leutbecher2017Leutbecher, M., Lock, SJ., Ollinaho, P., Lang, STK., Balsamo, G., Bechtold, P.Weisheimer, A.  2017\.  Stochastic representations of model uncertainties at ECMWF: State of the art and future vision Stochastic representations of model uncertainties at ecmwf: State of the art and future vision.  Quarterly Journal of the Royal Meteorological Society1437072315–2339.
  * Leutbecher Palmer (2008) leutbecher2008Leutbecher, M. Palmer, TN.  2008\.  Ensemble Forecasting Ensemble forecasting.  Journal of Computational Physics2273515–3539.  https://doi.org/10.1016/j.jcp.2007.02.014
  * Lorenz (1969) lorenz1969Lorenz, EN.  1969\.  The Predictability of a Flow which Possesses Many Scales of Motion The predictability of a flow which possesses many scales of motion.  Tellus213289–307.
  * Lorenz (1996) lorenz1996Lorenz, EN.  1996\.  Predictability– A Problem Partly Solved Predictability– a problem partly solved.  Proceedings, Seminar on Predictability, 4-8 September 1995 Proceedings, seminar on predictability, 4-8 september 1995 ( 1, 1–18).  Shinfield Park, ReadingECMWF.
  * Mirza Osindero (2014) Mirza2014Mirza, M. Osindero, S.  2014\.  Conditional Generative Adversarial Nets Conditional generative adversarial nets.  arXiv preprint.  https://arxiv.org/abs/1411.1784
  * Monahan . (2015) mrhm15Monahan, AH., Rees, T., He, Y. McFarlane, N.  2015\.  Multiple regimes of wind, stratification, and turbulence in the stable boundary layer Multiple regimes of wind, stratification, and turbulence in the stable boundary layer.  J. Atmos. Sci.723178-3198.
  * Moncrieff . (2007) Moncrieff2007Moncrieff, MW., Shapiro, MA., Slingo, JM. Molteni, F.  2007\.  Collaborative research at the intersection of weather and climate Collaborative research at the intersection of weather and climate.  WMO Bulletin56204–211.
  * Palmer (2001) Palmer:2001gwPalmer, TN.  200101\.  A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parametrization in weather and climate prediction models A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parametrization in weather and climate prediction models.  Q. J. R. Meteor. Soc.127572279–304.
  * Palmer (2012) Palmer:2012fePalmer, TN.  201204\.  Towards the probabilistic Earth-system simulator: a vision for the future of climate and weather prediction Towards the probabilistic Earth-system simulator: a vision for the future of climate and weather prediction.  Quarterly Journal of the Royal Meteorological Society138665841–861.
  * Palmer . (1997) Palmer:1997vaPalmer, TN., Barkmeijer, J., Buizza, R. Petroliagis, T.  199712\.  The ECMWF Ensemble Prediction System The ECMWF Ensemble Prediction System.  Meteorological Applications404301–304.
  * Palmer, Buizza . (2009) Palmer:2009uuPalmer, TN., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G.Weisheimer, A.  2009\.  Stochastic parametrization and model uncertainty. Stochastic parametrization and model uncertainty.
  * Palmer, Doblas-Reyes Weisheimer (2009) Palmer:2009jiPalmer, TN., Doblas-Reyes, FJ. Weisheimer, A.  2009\.  Toward seamless prediction: calibration of climate change projections using seasonal forecasts reply. Toward seamless prediction: calibration of climate change projections using seasonal forecasts reply.
  * Palmer . (2008) Palmer:2008ckPalmer, TN., Doblas-Reyes, FJ., Weisheimer, A. Rodwell, MJ.  2008\.  Toward seamless prediction: Calibration of climate change projections using seasonal forecasts Toward seamless prediction: Calibration of climate change projections using seasonal forecasts.  Bull. Amer. Meteor. Soc.894459–470.
  * Pollard (2002) pollard2002Pollard, D.  2002\.  A User’s Guide to Measure Theoretic Probability A user’s guide to measure theoretic probability.  Cambridge University Press.
  * Rabiner (1989) rabiner89Rabiner, LR.  1989\.  A tutorial on hidden Markov Models and selected applications in speech recognition A tutorial on hidden Markov models and selected applications in speech recognition.  Proc. IEEE77257-286.
  * Rasp . (2018) Rasp2018Rasp, S., Pritchard, MS. Gentine, P.  2018\.  Deep learning to represent subgrid processes in climate models Deep learning to represent subgrid processes in climate models.  Proceedings of the National Academy of Sciences115399684–9689.  https://www.pnas.org/content/115/39/9684 10.1073/pnas.1810286115
  * Reyes . (2009) Reyes:2009iyReyes, FD., Weisheimer, A., Déqué, M., Keenlyside, N., McVean, M., Murphy, JM.Palmer, TN.  2009\.  Addressing model uncertainty in seasonal and annual dynamical ensemble forecasts Addressing model uncertainty in seasonal and annual dynamical ensemble forecasts.  Q. J. R. Meteor. Soc.1351538–1559.
  * Sakradzija Klocke (2018) sakradzija2018Sakradzija, M. Klocke, D.  2018\.  Physically Constrained Stochastic Shallow Convection in Realistic Kilometer‐Scale Simulations Physically constrained stochastic shallow convection in realistic kilometer‐scale simulations.  JAMES10112755–2776.
  * Sanchez . (2016) Sanchez2016Sanchez, C., Williams, KD. Collins, M.  2016\.  Improved stochastic physics schemes for global weather and climate models Improved stochastic physics schemes for global weather and climate models.  Quarterly Journal of the Royal Meteorological Society142147–159.  https://doi.org/10.1002/qj.2640
  * Schneider . (2017) Schneider2017Schneider, T., Lan, S., Stuart, A. Teixeira, J.  2017\.  Earth System Modeling 2.0: A Blueprint for Models That Learn From Observations and Targeted High-Resolution Simulations Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations.  Geophysical Research Letters442412,396-12,417.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL076101 10.1002/2017GL076101
  * Seiffert von Storch (2010) Seiffert2010Seiffert, R. von Storch, JS.  2010\.  A Stochastic Analysis of the Impact of Small-Scale Fluctuations on the Tropospheric Temperature Response to CO2 Doubling A stochastic analysis of the impact of small-scale fluctuations on the tropospheric temperature response to co2 doubling.  J. Clim.232307–2319.  https://doi.org/10.1175/2009JCLI3043.1
  * Shapiro . (2010) Shapiro2010Shapiro, M., Shukla, J., Brunet, G., Nobre, C., Béland, M., Dole, R.Wallace, JM.  2010\.  An Earth-System Prediction Initiative for the Twenty-First Century An earth-system prediction initiative for the twenty-first century.  Bull. Amer. Meteor. Soc.91101377-1388.  https://doi.org/10.1175/2010BAMS2944.1
  * Shutts Pallares (2014) shutts2014Shutts, GJ. Pallares, AC.  2014\.  Assessing parametrization uncertainty associated with horizontal resolution in numerical weather prediction models Assessing parametrization uncertainty associated with horizontal resolution in numerical weather prediction models.  Phil. Trans. R. Soc. A3722018.
  * Shutts Palmer (2007) shutts2007Shutts, GJ. Palmer, TN.  2007\.  Convective Forcing Fluctuations in a Cloud-Resolving Model: Relevance to the Stochastic Parameterization Problem Convective forcing fluctuations in a cloud-resolving model: Relevance to the stochastic parameterization problem.  J. Clim.202187–202.
  * Stockdale . (2011) Stockdale:2011kmStockdale, TN., Anderson, DLT., Balmaseda, MA., Doblas-Reyes, F., Ferranti, L., Mogensen, K.Vitart, F.  201108\.  ECMWF seasonal forecast system 3 and its prediction of sea surface temperature ECMWF seasonal forecast system 3 and its prediction of sea surface temperature.  Climate Dynamics373-4455–471.
  * Strømmen . (2018) strommen2018Strømmen, K., Christensen, HM., Berner, J. Palmer, TN.  2018\.  The impact of stochastic parametrisations on the representation of the Asian summer monsoon The impact of stochastic parametrisations on the representation of the asian summer monsoon.  J. Clim.505–62269–2282.
  * Sušelj . (2014) Suselj:2014thSušelj, K., Hogan, TF. Teixeira, J.  2014\.  Implementation of a stochastic eddy-diffusivity/mass-flux parameterization into the Navy Global Environmental Model Implementation of a stochastic eddy-diffusivity/mass-flux parameterization into the Navy Global Environmental Model.  Weather and forecasting2961374–1390.
  * Suselj . (2013) Suselj:2013jdSuselj, K., Teixeira, J. Chung, D.  201307\.  A Unified Model for Moist Convective Boundary Layers Based on a Stochastic Eddy-Diffusivity/Mass-Flux Parameterization A Unified Model for Moist Convective Boundary Layers Based on a Stochastic Eddy-Diffusivity/Mass-Flux Parameterization.  Journal of Atmospheric Sciences7071929–1953.
  * Teixeira Reynolds (2010) Teixeira:2010hdTeixeira, J. Reynolds, CA.  201002\.  Stochastic Nature of Physical Parameterizations in Ensemble Prediction: A Stochastic Convection Approach Stochastic Nature of Physical Parameterizations in Ensemble Prediction: A Stochastic Convection Approach.  dx.doi.org1362483–496.
  * Tribbia Baumhefner (2004) tribbia2004Tribbia, JJ. Baumhefner, DP.  2004\.  Scale Interactions and Atmospheric Predictability: An Updated Perspective Scale interactions and atmospheric predictability: An updated perspective.  Mon. Wea. Rev.1323703–713.
  * Vannitsem Lucarini (2016) vannitsem2016Vannitsem, SS. Lucarini, V.  2016\.  Statistical and Dynamical Properties of Covariant Lyapunov Vectors in a Coupled Atmosphere-Ocean Model - Multiscale Effects , Geometric Degeneracy , and Error Dynamics Statistical and Dynamical Properties of Covariant Lyapunov Vectors in a Coupled Atmosphere-Ocean Model \- Multiscale Effects , Geometric Degeneracy , and Error Dynamics.  Journal of Physics A: Mathematical and Theoretical4922.
  * Vitart Robertson (2012) Vitart:2012ipVitart, F. Robertson, AW.  2012\.  Subseasonal to Seasonal Prediction Project: bridging the gap between weather and climate Subseasonal to Seasonal Prediction Project: bridging the gap between weather and climate.  World Meteorological Organisation Bulletin612.
  * Wang . (2016) wang2016Wang, Y., Zhang, GJ. Craig, GC.  2016\.  Stochastic convective parameterization improving the simulation of tropical precipitation variability in the NCAR CAM5 Stochastic convective parameterization improving the simulation of tropical precipitation variability in the NCAR CAM5.  Geophys. Res. Let.43126612–6619.  10.1002/2016GL069818
  * Watson (2019) watson2019Watson, PAG.  2019\.  Applying machine learning to improve simulations of dynamical systems using empirical error correction Applying machine learning to improve simulations of dynamical systems using empirical error correction.  JAMES.  Submitted
  * Weisheimer . (2014) Weisheimer:2014dfWeisheimer, A., Corti, S. Palmer, T.  2014\.  Addressing model error through atmospheric stochastic physical parametrizations: impact on the coupled ECMWF seasonal forecasting system Addressing model error through atmospheric stochastic physical parametrizations: impact on the coupled ECMWF seasonal forecasting system.  Phil. Trans. R. Soc. A372.
  * Wilks (2005) wilks2005Wilks, DS.  2005\.  Effects of Stochastic Parametrizations in the Lorenz ’96 System Effects of stochastic parametrizations in the Lorenz ’96 system.  Quarterly Journal of the Royal Meteorological Society131606389–407.
  * Wilks (2011) Wilks2011Wilks, DS.  2011\.  Statistical Methods in the Atmospheric Sciences Statistical Methods in the Atmospheric Sciences (3rd ).  Academic Press.
  * KD. Williams . (2013) williams2013Williams, KD., Bodas-Salcedo, A., Déqué, M., Fermepin, S., Medeiros, B., Watanabe, M.Williamson, DL.  2013\.  The transpose-AMIP II experiment and its application to the understanding of southern ocean cloud biases in climate models The transpose-AMIP II experiment and its application to the understanding of southern ocean cloud biases in climate models.  Journal of Climate26103258–3274.  10.1175/JCLI-D-12-00429.1
  * PD. Williams (2012) Williams2012Williams, PD.  2012\.  Climatic impacts of stochastic fluctuations in air-sea fluxes Climatic impacts of stochastic fluctuations in air-sea fluxes.  Geophysical Research Letters39.  https://doi.org/10.1029/2012GL051813