# Echo State Network for two-dimensional turbulent moist Rayleigh-Bénard
convection

Florian Heyder  Institut für Thermo- und Fluiddynamik, Technische Universität
Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany  Jörg Schumacher  Institut
für Thermo- und Fluiddynamik, Technische Universität Ilmenau, Postfach 100565,
D-98684 Ilmenau, Germany  Tandon School of Engineering, New York University,
New York City, NY 11201, USA

###### Abstract

Recurrent neural networks are machine learning algorithms which are suited
well to predict time series. Echo state networks are one specific
implementation of such neural networks that can describe the evolution of
dynamical systems by supervised machine learning without solving the
underlying nonlinear mathematical equations. In this work, we apply an echo
state network to approximate the evolution of two-dimensional moist Rayleigh-
Bénard convection and the resulting low-order turbulence statistics. We
conduct long-term direct numerical simulations in order to obtain training and
test data for the algorithm. Both sets are pre-processed by a Proper
Orthogonal Decomposition (POD) using the snapshot method to reduce the amount
of data. Training data comprise long time series of the first 150 most
energetic POD coefficients. The reservoir is subsequently fed by these data
and predicts of future flow states. The predictions are thoroughly validated
by original simulations. Our results show good agreement of the low-order
statistics. This incorporates also derived statistical moments such as the
cloud cover close to the top of the convection layer and the flux of liquid
water across the domain. We conclude that our model is capable of learning
complex dynamics which is introduced here by the tight interaction of
turbulence with the nonlinear thermodynamics of phase changes between vapor
and liquid water. Our work opens new ways for the dynamic parametrization of
subgrid-scale transport in larger-scale circulation models.

##  I Introduction

Machine learning (ML) algorithms have changed our way to model and analyse
turbulent flows including thermally driven convection flows Brenner et al.
(2019); Brunton et al. (2020); Pandey et al. (2020). The applications reach
from subgrid-scale stress models Ling et al. (2016); Duraisamy et al. (2019),
via the detection of large-scale patterns in mesoscale convection Fonda et al.
(2019) to ML-based parametrizations of turbulent transport and clouds in
climate and global circulation models Brenowitz and Bretherton (2018);
O’Gorman and Dwyer (2018); Gentine et al. (2018). The implementations of such
algorithms help to process and reduce increasing amounts of data coming from
numerical simulations and laboratory experiments Moller et al. (2020) by
detecting patterns and statistical correlations Goodfellow et al. (2016).
Moreover, the computational cost that is in line with a direct numerical
simulation (DNS) of the underlying highly nonlinear Navier-Stokes equations
can often be reduced significantly by running a neural network instead that
generates synthetic fields with the same low-order statistics as in a full
simulation Schneider et al. (2017); Mohan et al. (2020).

Turbulent convection, as all other turbulent flows, is inherently highly
chaotic so that specific flow states in the future are hard to predict after
exceeding a short horizon. The additional possibility of the fluid to change
its thermodynamic phase, as it is the case in moist turbulent convection
Stevens (2005); Mellado (2017), adds further nonlinearities and feedbacks to
the turbulence dynamics. Learning low-order statistics of such a turbulent
flow provides a challenge to an ML algorithm. For such a task, an internal
memory is required since statistical correlations decay in a finite time. This
necessitates the implementation of a short-term memory or internal cyclic
feedbacks in the network architecture. That is why a particular class of
supervised machine learning algorithms – recurrent neural networks (RNNs) –
will be in the focus of the present work Hochreiter and Schmidhuber (1997);
Lukoševičius and Jaeger (2009); Tanaka et al. (2019).

In this paper we apply a specific implementation of an RNN, the echo state
network (ESN) Jaeger and Haas (2004); Yildiz et al. (2012), to two-dimensional
(2d) turbulent moist Rayleigh-Bénard convection (RBC) flow. We use this RNN to
predict the low-order statistics, such as the buoyancy and liquid water fluxes
or fluctuation profiles of these quantities across the layer. Our present work
extends a recent application of ESNs to two-dimensional dry Rayleigh-Bénard
convection Pandey and Schumacher (2020) in several points. (1) The physical
complexity of turbulent convection is enhanced in the present moist case. This
is due to the total water content, an additional active scalar field which
comprises of vapor and liquid water contents. The total water content couples
as an additional field to the original dry Boussinesq model for temperature
and velocity. (2) The size of the convection domain is increased by a factor
of 4 such that the necessary degree of data reduction is significantly higher.
(3) Moist convection requires also the satisfying reproduction of new physical
quantities which are derived from different turbulence fields, e.g. the cloud
cover in or the liquid water flux across the layer. This can be considered as
a firmer test of the capabilities of the ESN to model complex dynamics. Such
quantities are of particular interest for larger-scale models of atmospheric
circulation that include mesoscale convection processes in form of parameters
or minimal models Grabowski and Smolarkiewicz (1999); Grabowski (2001). (4)
Finally, the hyperparameters of the ESN, in particular the spectral radius
$\rho(W^{r})$ of the reservoir matrix $W^{r}$, has been tested in more detail
(exact definitions follow). The grid search in our work thus sums up to a
total of more than 1800 different hyperparameter sets.

Echo state networks have recently received renewed attention for their
capability of equation-free modeling of several chaotic systems, such as of
the Lorenz96 model Vlachas et al. (2020) or the one-dimensional partial
differential Kuramoto-Sivashinsky equation Lu et al. (2017); Pathak et al.
(2018). Furthermore, low-order flow statistics in 2d dry Rayleigh-Bénard
convection with ${\rm Ra}=10^{7}$ have been successfully reproduced using an
ESN that trains the most energetic time coefficients of a Karhunen-Loéve
expansion (or proper orthogonal decomposition) of the convection fields Holmes
et al. (2012). This latter step can be considered as an autoencoder that
provides training and test data for the ESN which cannot take the data
directly, even for the present 2d case.

A second popular implementation of RNNs, which we want to mention here for
completeness, are long short-term memory networks (LSTM) Hochreiter and
Schmidhuber (1997) which have been also applied to fluid mechanical problems,
such as the dynamics of Galerkin models of shear flows in Srinivasan et al.
(2019). These models also demonstrated to capture the longer-term time
dependency and low-order statistics of important turbulent transport
properties well (see also Pandey et al. (2020) for a direct comparison). An
acceleration of the training, which requires the backpropagation of the errors
through the whole network in contrast to ESNs, were obtained recently with
Momentum LSTMs that apply the momentum extension of gradient descent search of
the cost function minimum to determine the weights at the network nodes Nguyen
et al. (2020).

In this work, a moist Rayleigh-Bénard convection model with moist Rayleigh
number ${\rm Ra}_{\rm M}\simeq 10^{8}$ and Prandtl number $\text{Pr}=0.7$ is
considered as an example case. We choose a 2d domain $\Omega=L\times H$ with
aspect ratio $A=L/H=24$. Here $L$ is the horizontal length and $H$ the height
of the simulation domain. The number of grid points was chosen as $N_{x}\times
N_{y}=7200\times 384$. The data are obtained by direct numerical simulations
(DNS) which apply a spectral element method Fischer (1997); Scheel et al.
(2013); nek (2017). Comprehensive studies of further data sets with different
parameters, such as Rayleigh numbers, to study the generalization properties
will be presented elsewhere. Our intention is to demonstrate the capability of
the ESN to deliver reliable low-order statistics for a turbulent convection
flow with phase changes.

The manuscript is structured as follows. The next section introduces the moist
RBC model and the central ideas of ESNs. Then the generation and analysis of
the DNS data, including a brief description of the numerical scheme and the
proper orthogonal decomposition (POD), is specified. Finally the results of
our machine learning approach to moist turbulent convection will be discussed
in detail. In section V we summarize our results and provide a conclusion and
an outlook.

##  II Methods

###  II.1 Moist Rayleigh-Bénard Convection Model

We now briefly review our model for moist Rayleigh-Bénard convection in two
spatial dimensions. A detailed derivation can be found in Pauluis and
Schumacher (2010); Weidauer et al. (2010); Schumacher and Pauluis (2010). The
framework is based on the mathematically equivalent formulation by Bretherton
Bretherton (1987, 1988). Similar simplified models of moist convection with
precipitation have been developed by Smith and Stechmann Smith and Stechmann
(2017) and Vallis et al. Vallis et al. (2019). Evaporative cooling and
buoyancy reversal effects were discussed for example by Abma et al. Abma et
al. (2013).

The buoyancy $B$ in atmospheric convection is given by Emmanuel (1994)

$\displaystyle
B=-g\frac{\rho(S,q_{v},q_{l},q_{i},p)-\overline{\rho}}{\overline{\rho}}$ (1)

with the gravity acceleration $g$, a mean density $\overline{\rho}$, the
pressure $p$, the entropy $S$ and the contents of water vapor $q_{v}$, liquid
water $q_{l}$ and ice $q_{i}$. We consider warm clouds only, i.e. $q_{i}=0$
and assume local thermodynamic equilibrium. From the latter assumption, it
follows that no precipitation is possible and the number of independent
variables in eq. (1) reduces to three. By introducing the total water content
$q_{T}=q_{v}+q_{l}$ the buoyancy can be expressed as $B(S,q_{T},p)$. In the
Boussinesq approximation, pressure variations about a hydrostatic equilibrium
profile are small, such that the buoyancy becomes $B(S,q_{T},y)$ with the
vertical spatial coordinate $y$ for the present 2d case. Furthermore, the
convection layer is near the vapor-liquid phase boundary. The buoyancy can
then be expressed as a piecewise linear function of $S$ and $q_{T}$ on both
sides of the saturation line. This step preserves the discontinuity of the
first partial derivatives of $B$ and therefore the physics of a first-order
phase transition. The advantage of this formulation is that locally the
saturation state of the air can be determined. In the last step, we substitute
the linear combinations of $S$ and $q_{T}$ on both sides of the phase boundary
by a dry buoyancy $D$ and a moist buoyancy $M$. Consequently the buoyancy
field $B$ is given by Pauluis and Schumacher (2010)

$\displaystyle B(x,y,t)=\max\left(M(x,y,t),D(x,y,t)-N_{s}^{2}y\right)$ (2)

where the fixed Brunt-Väisälä frequency
$N_{s}=\sqrt{g(\Gamma_{u}-\Gamma_{s})/T_{\rm ref}}$ is determined by the lapse
rates of the saturated and unsaturated moist air, $\Gamma_{s}$ and
$\Gamma_{u}$, and a reference temperature $T_{\rm ref}$. An air parcel at
height $y$ and time $t$ is unsaturated if $M(x,y,t)<D(x,y,t)-N_{s}^{2}y$ and
saturated if $M(x,y,t)>D(x,y,t)-N_{s}^{2}y$. Note that the newly introduced
dry buoyancy field $D$ is proportional to the liquid water static energy and
the moist buoyancy field $M$ to the moist static energy. As in dry Rayleigh-
Bénard convection with Dirichlet boundary conditions for the temperature, the
static diffusive profiles $\overline{D}(y),\overline{M}(y)$ are vertically
linear

$\displaystyle\overline{D}(y)$ $\displaystyle=D_{0}+\frac{D_{H}-D_{0}}{H}y$
(3) $\displaystyle\overline{M}(y)$
$\displaystyle=M_{0}+\frac{M_{H}-M_{0}}{H}y$ (4)

where $D_{0}$, $M_{0}$ and $D_{H}$, $M_{H}$ are the imposed values of $D$, $M$
at the bottom ($y=0$) and top ($y=H$) of the computational domain. Here,
$D_{0}=M_{0}$. The governing equations of the moist Boussinesq system are
given by

$\displaystyle\frac{d\mathbf{v}}{dt}$
$\displaystyle=-\nabla\tilde{p}+\nu\nabla^{2}\mathbf{v}+B\left(D,M,y\right)\hat{\mathbf{e}}_{y}$
(5) $\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=0$ (6)
$\displaystyle\frac{dD}{dt}$ $\displaystyle=\kappa\nabla^{2}D$ (7)
$\displaystyle\frac{dM}{dt}$ $\displaystyle=\kappa\nabla^{2}M.$ (8)

Here $\mathbf{v}=(v_{x}(x,y,t),v_{y}(x,y,t))^{T}$ is the two-dimensional
velocity field, $\tilde{p}=p/\rho_{\rm ref}$ the kinematic pressure, $\nu$ the
kinematic viscosity and $\kappa$ the scalar diffusivity. The term
$d/dt=\partial/\partial t+(\mathbf{v}\cdot\nabla)$ is the material derivative.
This idealized model describes the formation of warm, non-precipitating low
clouds in a shallow layer up to a depth of $\sim 1$km. The assumptions, which
are made here, hold for example to a good approximation over the subtropical
oceans.

The problem is made dimensionless using length scale $[x,y]=H$, buoyancy scale
$[B]=M_{0}-M_{H}$, and (free-fall) time scale $[t]=\sqrt{H/(M_{0}-M_{H})}$.
The characteristic velocity scale follows by
$[v_{x},v_{y}]=\sqrt{(M_{0}-M_{H})H}$. Four dimensionless numbers can be
identified: the Prandtl number, dry Rayleigh number and moist Rayleigh number
are given by

Pr $\displaystyle=\frac{\nu}{\kappa}$ (9) $\displaystyle{\rm Ra}_{\rm D}$
$\displaystyle=\frac{\left(D_{0}-D_{H}\right)H^{3}}{\nu\kappa}$ (10)
$\displaystyle{\rm Ra}_{\rm M}$
$\displaystyle=\frac{\left(M_{0}-M_{H}\right)H^{3}}{\nu\kappa}\,.$ (11)

An additional parameter arises from the additional phase changes, the
dimensionless form of (2)

${\rm CSA}=\frac{N_{s}^{2}H}{M_{0}-M_{H}}\,.$ (12)

The condensation in saturated ascent (CSA) controls the amount of latent heat
an ascending saturated parcel can release on its way to the top. The
saturation condition (2) implies that liquid water is immediately formed at a
point in space and time when $M>D-N_{s}^{2}y$. There is no supersaturation
considered in this model and the liquid water content field $q_{l}$ and thus
the clouds are given by

$q_{l}(x,y,t)=M(x,y,t)-(D(x,y,t)-N_{s}^{2}y)\,.$ (13)

Note that in this formulation, $q_{l}$ can become negative, as it is a measure
of the degree of saturation. In a nutshell, $q_{l}<0$ stand thus for a liquid
water deficit. When the atmosphere is saturated, $q_{l}\geq 0$ and the
conventional liquid water content is retained.

Here we study the case of $D_{0}>D_{H}$ and $M_{0}>M_{H}$. Both fields are
linearly unstable. For the case of a conditionally unstable moist layer with
$D_{0}\leq D_{H}$ we refer to refs. Bretherton (1987, 1988) or Pauluis and
Schumacher (2011).

###  II.2 Reservoir Computing

Reservoir computing (RC) is a special type of RNN implementation. Contrary to
standard feed forward networks, neurons in the hidden layers of RNN are
recurrently connected to each other. In this way, RNNs have a similar
architecture to biological brains and are said to posses an internal memory as
cyclic connections allow information to stay inside the network for a certain
amount of time before fading out, known as the echo state property. Yet the
training of such RNNs is exceedingly difficult, as common training schemes
like the back propagation through time struggle with fading error gradients,
slow convergence Jaeger (2002) and bifurcations Doya (1992). An alternative
training method was proposed by Jaeger Jaeger (2001) and in an independent
work by Maass Maass et al. (2002). Their idea, which is now known as reservoir
computing Lukoševičius and Jaeger (2009), was to train the weights of the
output layer only, which connect the RNN, denoted to as the reservoir, to the
output neurons. The weights of the input layer as well as the internal
reservoir weights should be initialized at random and then kept constant. This
training procedure reduces the computational costs for training significantly
and shifts the focus to an adequate initialization of the input and reservoir
weight matrix $W^{r}$ (which is an adjacency matrix in network theory). While
Jaeger’s approach is known as ESN, Maass’ framework is called liquid state
machine. They differ in their field of origin, as the ESN stems from the field
of machine learning and the liquid state machine from computational
neuroscience. We will stick to the ESN formulation, but note that RC refers to
the concept mentioned above and summarizes both models.

Despite their simple training scheme, ESNs have been said to tackle many
tasks, from forecasting closing prices of stock markets Lin et al. (2008) to
estimating the life span of a fuel cell Morando et al. (2013). Especially its
application to dynamical systems shows great promise. For instance it has been
demonstrated that the dynamics of two of the three degrees of freedom of the
Rössler system can be inferred from the evolution of the third one Lu et al.
(2017). Further, the Lyapunov exponents of the dynamical system that a trained
ESN represents have been shown to match the exponents of the data generating
system Pathak et al. (2017).

Figure 1: Sketch of the echo state network in the training phase (a) for time
steps $n<0$ and the prediction phase (b) for time steps $n\geq 0$.

Figure 1 shows the concept and components of the ESN, for the training phase
in panel (a) and for the prediction phase in panel (b). The input
$\mathbf{u}(n)\in\mathbb{R}^{\rm N_{\rm in}}$ at a time instance $n$ as well
as a constant scalar bias $b=1$ are passed to the reservoir via the input
weight matrix $W^{\rm in}\in\mathbb{R}^{\rm N_{\rm r}\times(1+\rm N_{\rm
in})}$. The weighted input contributes to the dynamics of the reservoir state
$\mathbf{s}\in\left[-1,1\right]^{\rm N_{\rm r}}$ at time $n$ which is given by

$\displaystyle\mathbf{s}(n)$ $\displaystyle=(1-\gamma)\mathbf{s}(n-1)$
$\displaystyle+\gamma\tanh\left[W^{\rm in}\left[b;\mathbf{u}(n)\right]+W^{\rm
r}\mathbf{s}(n-1)\right].$ (14)

Here $\left[b;\mathbf{u}(n)\right]$ stands for the vertical concatenation of
the scalar bias and the input 111In some cases the bias $b=0$ and
$\mathbf{u}(n-1)$ are used in eq. (II.2).. This update rule comprises external
forcing by the inputs $\mathbf{u}(n)$ as well as a self-interaction with the
past instance $\mathbf{s}(n-1)$. The reservoir weight matrix $W^{\rm
r}\in\mathbb{R}^{\rm N_{\rm r}\times\rm N_{\rm r}}$ blends the state
dimensions, while the nonlinear hyperbolic tangent $\tanh\left(\cdot\right)$,
applied to each component of its argument vector, is the nonlinear activation
function of the neurons in this model. The leaking rate
$\gamma\in\left(0,1\right]$ moderates the linear and nonlinear contributions
and assures that the state is confined to $\left[-1,1\right]^{N_{\rm r}}$. As
mentioned above, the existence of echo states, i.e., states that are purely
defined by the input history, is crucial. An ESN is said to include such echo
states when two different states $\mathbf{s}(n-1)$, $\mathbf{s}^{\prime}(n-1)$
converge to the same state $\mathbf{s}(n)$, provided the same input
$\mathbf{u}(n)$ is given and the system has been running for many iterations
$n$ Jaeger (2001). Therefore, the first few state iterations are considered as
a reservoir washout, even for a reservoir with echo state property. After this
transient phase, the updated state is concatenated with the bias and the
current input to form the extended reservoir state
$\tilde{\mathbf{s}}(n)=\left[b;\mathbf{u}(n);\mathbf{s}(n)\right]$. Finally,
$\tilde{\mathbf{s}}$ is mapped via the output matrix $W^{\rm
out}\in\mathbb{R}^{\rm N_{\rm in}\times\left(1+\rm N_{\rm in}+\rm N_{\rm
r}\right)}$ to form the reservoir output $\mathbf{y}(n)\in\mathbb{R}^{\rm
N_{\rm in}}$

$\mathbf{y}(n)=\mathrm{W}^{\rm out}\tilde{\mathbf{s}}(n).$ (15)

For our application the output dimension matches the input dimension $\rm
N_{\rm in}$, which generally does not need to be the case.

Before the ESN can be used in the prediction phase, as sketched in Fig. 1(b),
the elements of $W^{\rm out}$ have to be computed first. This process is known
as training phase of this supervised machine learning algorithm. Only when the
reservoir is properly trained it will produce reasonable output. A set of
$n_{\rm train}$ training data instances $\\{\mathbf{u}(n),\mathbf{y}^{\rm
target}(n)\\}$ (where $n=-n_{\rm train},-(n_{\rm train}-1),...,-1$) needs to
be prepared. The target output $\mathbf{y}^{\rm target}(n)$ represents the
desired output the ESN should produce for the given input $\mathbf{u}(n)$. The
reservoir state $\mathbf{s}$ is computed for all inputs in the training data
set and assembled into a mean square cost function with an additional $L_{2}$
regularization $C\left(W^{\rm out}\right)$ which is given by

$\displaystyle C\left(W^{\rm out}\right)$ $\displaystyle=\frac{1}{n_{\rm
train}}\sum_{n=-n_{\rm train}}^{-1}\left(W^{\rm
out}\tilde{\mathbf{s}}(n)-\mathbf{y}^{\rm target}(n)\right)^{2}$
$\displaystyle+\beta\sum\limits_{i=1}^{N_{\rm in}}\|w_{i}^{\rm
out}\|_{2}^{2},$ (16)

and has to be minimized corresponding to

$\displaystyle W_{\ast}^{\rm out}$ $\displaystyle=\arg\min C\left(W^{\rm
out}\right)\,.$ (17)

Here $w_{i}^{\rm out}$ denotes the $i$-th row of $W^{\rm out}$ and
$\|\cdot\|_{2}$ the $L_{2}$ norm. Equations (16) and (17) are known as ridge
regression with the ridge regression parameter $\beta$. The last term in (16)
suppresses large values of the rows of $W^{\rm out}$, which could
inadvertently amplify small differences of the state dimensions in (15). This
well known regression problem is solved by the fitted output matrix

$\displaystyle W_{\ast}^{\rm out}=Y^{\rm target}S^{\rm T}\left(SS^{\rm
T}+\beta{\rm Id}\right)^{-1}$ (18)

where $\left(\cdot\right)^{\rm T}$ denotes the transposed and ${\rm Id}$ the
identity matrix. $Y^{\rm target}$ and $S$ are matrices where the $n$-th column
is the target output $\mathbf{y}^{\rm target}(n)$ and the extended reservoir
state $\tilde{\mathbf{s}}(n)$, respectively.

As the output weights are the only parameters that are trained, RC is
computationally inexpensive. However, the algebraic properties of the
initially randomly generated matrices $W^{\rm in}$ and $W^{\rm r}$ are
hyperparameters which have to be tuned beforehand. In our approach we draw the
elements of $W^{\rm in}$ from a uniform distribution in
$\left[-0.5,0.5\right]$ and impose no further restrictions on this matrix. For
the generation of the reservoir weights in $W^{\rm r}$ it has been reported
that the proportion of non-zero elements, i.e., the reservoir density D and
the spectral radius $\varrho\left(W^{\rm r}\right)$, i.e., the largest
absolute eigenvalue of $W^{\rm r}$ are crucial parameters that determine
whether the desired echo state property holds Lukoševičius (2012). We choose a
sparse reservoir ($\text{D}<1$) with few internal node connections and draw
the elements from a uniform distribution in $\left[-1,1\right]$. We then
normalize $W^{\rm r}$ by its largest absolute eigenvalue and multiply it with
the desired spectral radius $\varrho\left(W^{\rm r}\right)$. This scaling
approach, initially proposed by Jaeger Jaeger (2002), has established itself
as one of the standard ways Morando et al. (2013); Lu et al. (2017); Pathak et
al. (2018); Pandey and Schumacher (2020) to control the spectral radius.
Nevertheless, other procedures have been proposed Yildiz et al. (2012);
Strauss et al. (2012). In addition, the size of the reservoir $\rm N_{\rm r}$
is a hyperparameter. Usually $\mathbf{s}$ should be a high-dimensional
extension of the inputs $\mathbf{u}$, so that $\rm N_{\rm r}\gg\rm N_{\rm in}$
is satisfied. Moreover, we consider the leaking rate $\gamma$ and ridge
regression parameter $\beta$ as further quantities that have to be adjusted to
our data.

##  III Echo state network for 2d Moist Convection

###  III.1 Direct Numerical Simulations

DNS using the spectral element solver Nek5000 Fischer (1997); Scheel et al.
(2013); nek (2017) were conducted to solve the two-dimensional moist Rayleigh-
Bénard system (5)-(8) in a domain $\Omega=L\times H$ with aspect ratio
$A=L/H=24$. The Rayleigh numbers are ${\rm Ra}_{\rm D}=2\cdot 10^{8}$, ${\rm
Ra}_{\rm M}=4\cdot 10^{8}$. The Prandtl number is ${\rm Pr}=0.7$ representing
moist air. The additional parameter is set to ${\rm CSA=0.3}$. In the vertical
direction $y$, Dirichlet boundary conditions are imposed for both buoyancy
fields at top and bottom in combination with free-slip boundary conditions for
the velocity field. Periodic boundaries are set for all fields in the
horizontal direction. We chose a spatial resolution of $N_{x}\times
N_{y}=7200\times 384$ grid points and a time step size of $5.0\cdot 10^{-4}$.
This setup corresponds to an absolutely unstable atmosphere, i.e. where both
unsaturated and saturated air are unstable w.r.t. vertical displacements. The
initial conditions are small perturbations around the diffusive equilibrium
state $\overline{M}(y)$ and $\overline{D}(y)$, which result in turbulent
convection. The flow statistics relaxes into a statistically stationary state
(see Fig. 2) which provides training and test data for further processing.

Figure 2: Turbulent kinetic energy $E_{\rm kin}(t)=\langle
v_{x}^{2}+v_{y}^{2}\rangle_{x,y}$ of the moist Rayleigh-Bénard flow versus
time $t$. After an initial transient the values such as those of $E_{\rm
kin}(t)$ become statistically stationary. In the statistically stationary
regime, $2000$ snapshots, each separated by $\Delta t=0.25$, are analyzed by a
POD (see Sec. III.2).

###  III.2 Data reduction by POD

We sample a total of $n_{s}=2000$ snapshots of $v_{x}$, $v_{y}$, $M$ and $D$
at a time interval $\Delta t=0.25$ in the statistically stationary regime. The
original DNS data snapshots have been sampled at a constant time interval such
that we stick to constant time steps throughout this work, including the
subsequent RC model. Furthermore, the data are spectrally interpolated on a
uniform grid with a resolution of $N_{x}^{\prime}\times
N_{y}^{\prime}=640\times 80$ points from the originally unstructured element
mesh. They are decomposed into temporal mean and fluctuations subsequently,

$\displaystyle v_{x}(x,y,t)$ $\displaystyle=$ $\displaystyle\langle
v_{x}\rangle_{t}(x,y)+v_{x}^{\prime}(x,y,t)$ (19) $\displaystyle v_{y}(x,y,t)$
$\displaystyle=$ $\displaystyle\langle
v_{y}\rangle_{t}(x,y)+v_{y}^{\prime}(x,y,t)$ (20) $\displaystyle D(x,y,t)$
$\displaystyle=$ $\displaystyle\langle D\rangle_{t}(x,y)+D^{\prime}(x,y,t)$
(21) $\displaystyle M(x,y,t)$ $\displaystyle=$ $\displaystyle\langle
M\rangle_{t}(x,y)+M^{\prime}(x,y,t)$ (22) $\displaystyle q_{l}(x,y,t)$
$\displaystyle=$ $\displaystyle\langle
q_{l}\rangle_{t}(x,y)+q_{l}^{\prime}(x,y,t)\,.$ (23)

Figure 3: Eigenvalue spectrum of the POD mode obtained from the analysis of
2000 snapshots. (a) Individual and cumulative contribution of each mode. The
shaded region indicates the first $\rm N_{\rm POD}=150$ modes which capture
$81\%$ of the total energy of the original snapshot data. (b) Time
coefficients $a_{p}(n)$ for the 1st, 2nd, 10th, 50th, and 100th mode are
shown. The first coefficients show a slow variation compared to higher
coefficients. Time series are shifted with respect to each other for better
visibility.

A grid dimension of $640\times 80$ in terms of ESN input dimensions ${\rm
N_{\rm in}}$ is still too big. We therefore follow the approach in Pandey and
Schumacher (2020) and introduce an intermediate step of data reduction before
handing the data to the ESN. We make use of the periodicity and expand the
data in a Fourier series in the horizontal $x$-direction and take the
Karhunen-Loéve expansion, also known as POD of our data. In particular we
choose the snapshot method Sirovich (1987) which decomposes the $k$-th
component of a vector field $\mathbf{g}(x,y,t)$ into

$\displaystyle g_{k}(x,y,t)$ $\displaystyle=\sum\limits_{p=1}^{\rm n_{\rm
s}}\sum_{n_{x}=-N_{x}^{\prime}/2}^{N_{x}^{\prime}/2}a_{p,n_{x}}(t)\Phi_{k,n_{x}}^{(p)}(y)\exp{\left(i\frac{2\pi
n_{x}x}{L}\right)}$ $\displaystyle=\sum\limits_{p=1}^{\rm n_{\rm
s}}a_{p}(t)\Phi_{k}^{(p)}(x,y).$ (24)

Here $a_{p}(t)$ and $\Phi_{k}^{(p)}(x,y)$ are the $p$-th time coefficient and
the corresponding spatial POD mode respectively. In our approach we take the
POD of $\mathbf{g}=(v_{x}^{\prime},v_{y}^{\prime},D^{\prime},M^{\prime})^{T}$.
The POD spectrum of the turbulent convection data can be seen in Fig. 3(a).
The eigenvalues of the covariance matrix fall off quickly and we therefore
truncate the POD expansion at $\rm N_{\rm POD}=150\ll\rm n_{\rm s}$ and
include only the most energetic modes (green shaded region). These capture
$81\%$ of the total energy of the original data. In Fig. 3(b) the time
coefficients $a_{p}(n)$ are shown for $p=1,2,10,50$, and 100 for all POD time
steps. The first time coefficients ($a_{1}$ to $a_{10}$) posses only few
temporal features opposed to higher coefficients. This range of active scales
is inherent to turbulence as kinetic energy of large-scale motion is
transferred to small eddies down to the Kolmogorov scale. Moreover, the
influence of the additional nonlinearity due to the phase changes impacts the
dynamics, as the first coefficients varied more in the dry RBC case with
aspect ratio $6$ at ${\rm Ra=10^{7}}$ Pandey and Schumacher (2020). This is
one order of magnitude below our Rayleigh number values. Nevertheless our RC
model will receive values for all $\rm N_{\rm POD}$ coefficients and therefore
for a wide range of temporal frequencies. We note here that the design of the
RC model could be adapted to the different frequencies in future approaches.

Figure 4 shows the spatial modes $\Phi_{2}^{(1)}$, $\Phi_{2}^{(50)}$ and
$\Phi_{4}^{(1)}$, $\Phi_{4}^{(50)}$.

Figure 4: Spatial structure of two POD modes for $v_{y}$: (a)
$\Phi_{2}^{(1)}(x,y)$, (b) $\Phi_{2}^{(50)}(x,y)$ and the moist buoyancy field
$M$: (c) $\Phi_{4}^{(1)}(x,y)$, (d) $\Phi_{4}^{(50)}(x,y)$. For visualization
purposes the aspect proportions do not match the actual aspect ratio of
$A=24$.

We limit ourselves to the last $1400={\rm n}_{\rm train}+{\rm n}_{\rm test}$
time instances of our data and use the first 150 time coefficients
$\mathbf{a}(n)=\left(a_{1}(n),a_{2}(n),...,a_{150}(n)\right)^{T}$ as the input
for the ESN. The first ${\rm n}_{\rm train}=700$ snapshots are assembled into
a training data set

$\displaystyle\\{\mathbf{u}(n),\mathbf{y}^{\rm
target}(n)\\}=\\{\mathbf{a}(n),\mathbf{a}(n+1)\\}$ (25)

with $-{\rm n}_{\rm train}\leq n\leq-1$. The training data span 175 free-fall
time units $T_{f}$ that correspond to 61 eddy turnover times. This time scale
is given by $\tau_{\rm eddy}=H/u_{\rm rms}\approx 2.9T_{f}$ with $u_{\rm
rms}=\langle u_{x}^{2}+u_{y}^{2}\rangle^{1/2}_{x,y,t}$. The ESN is trained to
predict the time instance $\mathbf{a}(n+1)$ when given the POD time
coefficients $\mathbf{a}(n)$ at the last time step as input. We use the first
$46$ time steps to initialize the reservoir. Using eq. (18) we compute $W^{\rm
out}$, which can then be used for prediction. In the prediction phase, we give
the initial input $\mathbf{u}(0)=\mathbf{a}(0)$ to the reservoir and redirect
the output to the input layer (see Fig. 1(b)) such that

$\displaystyle\mathbf{u}(n)$ $\displaystyle=\mathbf{y}(n-1)\qquad n=1,...,{\rm
n}_{\rm test}-1$ (26)

with ${\rm n}_{\rm test}=700$. This coupling creates an autonomous system that
generates new outputs without providing external inputs. Contrary to the
teacher forcing approach, where at each time step the input is given by the
actual time coefficients, this method is more suited for real world
application. Finally, the outputs at each time step are gathered and validated
by the last ${\rm n}_{\rm test}$ snapshots $\\{\mathbf{y}^{\rm
val}(n)\\}=\\{\mathbf{a}(n+1)\\}$.

###  III.3 Training of ESN and choice of hyperparameters

We quantify the quality of ESN predictions, with the set of hyperparameters
$\rm h=\\{\gamma,\beta,\rm N_{\rm r},{\rm D},\varrho(W^{\rm r})\\}$, by two
types of measures. The mean square error ${\rm MSE}_{\rm h}$ of ESN output to
the validation data

$\displaystyle{\rm MSE}_{\rm h}$ $\displaystyle=\frac{1}{{\rm n}_{\rm
test}}\sum\limits_{n=0}^{{\rm n}_{\rm test}}{\rm mse}(n)$ (27) where
$\displaystyle{\rm mse}(n)$ $\displaystyle=\frac{1}{\rm N_{\rm
POD}}\sum\limits_{i=1}^{{\rm N}_{\rm POD}}\left(y_{i}(n)-y^{\rm
val}_{i}(n)\right)^{2}$ (28)

is the mean square error at time step $n$ averaged over all ${\rm N}_{\rm
POD}$ modes. Additionally, we take the physically more relevant normalized
average relative error (NARE) as defined in Srinivasan et al. (2019) into
account. For the moist buoyancy field $M$, it is for example given by

$\displaystyle E_{\rm h}\left[\langle M\rangle_{x,t}\right]$
$\displaystyle=\frac{1}{C_{\max}}\int\limits_{0}^{1}\Big{|}\langle
M\rangle_{x,t}^{\rm ESN}(y)-\langle M\rangle_{x,t}^{\rm POD}(y)\Big{|}dy$ (29)
with $\displaystyle C_{\max}$
$\displaystyle=\frac{1}{2\max_{y\in[0,1]}(|\langle M\rangle_{x,t}^{\rm
POD}|)}.$ (30)

The superscript defines whether the field was reconstructed with $a_{i}(n)$
(POD) or $y_{i}(n)$ (ESN). It measures the integral deviation of the
reconstructed line-time average profile $\langle\cdot\rangle_{x,t}$ of a
specific field. We consider the three NAREs: $E_{\rm h}\left[\langle
M\rangle_{x,t}\right]$, $E_{\rm h}\left[\langle q_{l}^{\prime
2}\rangle_{x,t}\right]$ and $E_{\rm h}\left[\langle
v_{y}^{\prime}M^{\prime}\rangle_{x,t}\right]$, that is the NARE of the total
moist buoyancy field $M$, the liquid water content fluctuations
$q_{l}^{\prime}$, and the moist buoyancy flux fluctuations
$v_{y}^{\prime}M^{\prime}$.

$\gamma$ | $\beta$ | $D$ | $\varrho(W^{\rm r})$
---|---|---|---
$0.50$ | $5\cdot 10^{-4}$ | $0.1$ | $0.00$
$0.60$ | $5\cdot 10^{-3}$ | $0.2$ | $0.90$
$0.70$ | $5\cdot 10^{-2}$ | $0.3$ | $0.91$
$0.80$ | $5\cdot 10^{-1}$ | $0.4$ | $0.92$
$0.90$ |  | $0.5$ | $0.93$
$0.95$ |  | $0.6$ | $0.94$
|  | $0.7$ | $0.95$
|  |  | $0.96$
|  |  | $0.97$
|  |  | $0.98$
|  |  | $0.99$
|  |  | $1.00$
Table 1: Range of the four hyperparameters upon which a grid search was
conducted. For each of the $1848$ combinations, an ESN was trained and
validated with the training and validation data set. The reservoir dimension
$N_{\rm r}$ was set to $4000$ for all studies. The MSE and NARE measures were
computed and evaluated to find the optimal parameter set $h^{*}$.

$\gamma^{*}$ | $\beta^{*}$ | $N_{\rm r}^{*}$ | $D^{*}$ | $\varrho(W^{\rm r})^{*}$
---|---|---|---|---
$0.9$ | $5\cdot 10^{-1}$ | $4000$ | $0.1$ | $1.0$

MSE($\rm h^{*}$) | $E_{\rm h}\left[\langle M\rangle_{x,t}\right]$ | $E_{\rm h}\left[\langle q_{l}^{\prime 2}\rangle_{x,t}\right]$ | $E_{\rm h}\left[\langle v_{y}^{\prime}M^{\prime}\rangle_{x,t}\right]$
---|---|---|---
$8.18\cdot 10^{-4}$ |  $0.032$% |  $0.033$% |  $4.5$%

Table 2: Hyperparameter set $h^{*}$ that was chosen for the ESN setup and the
associated errors, which this ESN run has produced. The results of the
reservoir with these listed hyperparameters are presented in section IV.
Figure 5: Representative profiles taken from the error landscape for the
leaking rate $\gamma$ (a,b,c) and the spectral radius $\varrho(W^{\rm r})$
(d,e,f). The data are obtained by a grid search study (see Table 1). We find a
systematic dependence for the two quantities in this parameter domain. Note
the different magnitudes between single quantity-NARE in (b,e) and multiple
quantity-NARE in panels (c,f). The legends in (c) and (f) hold also for panels
(a,b) and (d,e), respectively.

A grid search for the four quantities $\gamma$, $\beta$, $D$, $\varrho(W^{\rm
r})$ was conducted in a suitable range (see Table 1), based on the results in
Pandey and Schumacher (2020). The reservoir size was fixed to $N_{\rm r}=4000$
for all runs. The resulting ${\rm MSE}_{h}$ and NAREs were computed to find an
adequate parameter setting. Figure 5 shows ${\rm MSE}_{h}$, $E_{\rm
h}\left[\langle q_{l}^{\prime 2}\rangle_{x,t}\right]$, $E_{\rm h}\left[\langle
v_{y}^{\prime}M^{\prime}\rangle_{x,t}\right]$ and their dependence on
$\varrho\left(W^{\rm r}\right)$ and $\gamma$. We detect a systematic
dependence of both, spectral radius and leaking rate, even when slightly
changing a third parameter (see legends). Interestingly, as both parameters
are increased, the mean square error increases as well, while the NARE either
decrease or pass a local minimum. We emphasize that our grid search is only an
excerpt of the much bigger error landscape. Further, we did not average over
multiple random realizations which would be the basis for a more rigorous
discussion of parameter dependencies. Nevertheless, a starting point of the
discussion of the hyperparameter dependencies is as follows: as the spectral
radius grows, the magnitude of the argument of the hyperbolic tangent builds
up too. This will saturate the activation function, which in turn will act in
an increasingly binary way since
$\lim\limits_{x\rightarrow\pm\infty}\tanh(x)=\pm 1$. In this limit of a fully
saturated activation function, eq. (II.2) simplifies to

$\displaystyle\mathbf{s}(n)\simeq(1-\gamma)\mathbf{s}(n-1)\pm\gamma\mathbf{1},$
(31)

where $\mathbf{1}=(1,1,...,1)^{T}\in\mathbb{R}^{N_{\rm r}}$. This corresponds
to a linear dependence of each reservoir state on its last instance plus the
constant leaking rate $\gamma$ with stochastically changing sign, depending on
the randomly generated weight matrices. As the leaking rate is increased
towards unity, the memory effect is lost as well and the reservoir state is
basically updated by the constant last term in (31). The resulting reservoir
output will lead to increased mean square deviations from the varying ground
truth signal. We thus speculate that such activation saturation is already
satisfied for several reservoir state components at $\rho(W^{r})\lesssim 1$
which in turn contribute to the increasing ${\rm MSE}_{\rm h}$ in panel (d) of
Fig. 5.

Finally, we chose the hyperparameter set $h^{*}$, listed in Table 2 as it
results in a minimum of $E_{\rm h}\left[\langle
v_{y}^{\prime}M^{\prime}\rangle_{x,t}\right]$. The reason for settling with
this measure is that it is susceptible to two ESN estimates. Quantities like
$E_{\rm h}\left[\langle q_{l}^{\prime 2}\rangle_{x,t}\right]$, which depend on
only one ESN estimate, exhibit small values for many parameter settings (see
Fig. 5 (b),(e)).

##  IV Results for the moist RBC case

Figure 6: Time evolution of the POD time coefficients $a_{i}(t)$. The gray
shaded area marks the training phase (reservoir output not shown). At the end
of the training phase the prediction phase starts. The curves labeled POD
stand for the ground truth of the evolution of the coefficient, while those
labeled ESN are the network predictions. Panels (f)–(j) show the initial part
of the forecast and correspond to (a)–(e). Figure 7: Fourier spectrum
$\|\mathcal{F}(a_{i})\|^{2}$ of the POD time coefficients and of the
corresponding reservoir prediction. The first 200 frequencies are shown only.

After the ESN receives the initial input $\mathbf{u}(0)=\mathbf{a}(0)$, the
autonomous predictor (see eq. (26)) produces estimates for the POD time
coefficients which can be seen in Fig. 6. From the predicted coefficients and
the known POD modes we can reconstruct all fields and compare these with the
ground truth which is the POD expansion of the test data.

Deviations for the first ten coefficients are detected while predictions for
subsequent POD coefficients associated with less energetic modes agree with
the values of the validation data for the first few time steps. Nevertheless,
the ESN accomplishes to produce a time series with matching temporal frequency
as the actual data, but shows bigger deviations to compute the trend of the
slowly varying first coefficients.

The frequency spectra of the ESN predictions for the coefficients $a_{i}(t)$
in comparison to those of the test data are shown in Figure 7 for 5 different
cases. While the spectral values of the first 100 frequencies are captured
well by the ESN, the higher frequency part starts to deviate in most cases. As
discussed already above, this might be due to a simple RC model architecture,
which does not differentiate between significantly different time scales that
are always present in a turbulent flow. Nevertheless, the result underlines
that the ESN is able to learn the time scales of the most significant POD
coefficients.

Figure 8(a) shows the mean square error ${\rm mse}(n)$ over all modes, as
defined in (28), as a function of time steps after training. The mean error
initially rises and then saturates. The fact that errors increase stems from
the coupling scheme of output to input; small errors will inadvertently be
amplified by the nonlinearity of the activation function in (II.2). Figure
8(b) shows the deviations $(y_{i}^{\rm val}(n)-y_{i}(n))$ for $i=1,10,50,100$.

Figure 8: ESN prediction error: (a) Mean square error over all modes ${\rm
mse}(n)$, see. eq. (28), a function of time steps $n$ after the training
phase. (b) The difference $(y_{i}^{\rm val}(n)-y_{i}(n))$, i.e. the deviation
of the ESN prediction $y_{i}(n)$ of $a_{i}$ at time steps $n$ after training.

Figures 9(a)–(c) show the three weight matrices of the RC model. As described
in section II.2, the input and reservoir weights are initialized randomly and
left unchanged. With a reservoir density of $D=0.2$ the reservoir matrix
$W^{\rm r}$ is a sparse matrix containing many weights equal to zero. Note
further, that the fitted output weights $W^{\rm out}$ have low magnitudes in
comparison to the entries of the input matrix $W^{\rm in}$. This is adjusted
according to the number of reservoir nodes $N_{\rm r}$ and the magnitude of
the training data. Moreover, the magnitude of the first $1+N_{\rm in}$ columns
are close to zero. This indicates that the contributions of the output bias
$b$ and the current input $\mathbf{u}$ to the reservoir output (see eq. (15))
are small.

Figure 10 (a,b) shows the training and prediction phase dynamics of three
exemplary hidden reservoir state components $s_{1},s_{1000}$, and $s_{4000}$.
As the first $46$ time steps of the training data were used to initialize the
reservoir state, the last $n_{\rm Train}-46$ time steps are shown in panel (a)
only. During both phases the individual time series $s_{i}$ are confined to a
certain subrange of the whole range $\left[-1,1\right]$. They have comparable
amplitudes. Nevertheless, the prediction phase time series differ from the
their training phase counterparts by slightly smoother variations with respect
to time, see also the corresponding Fourier frequency spectra in panels (c,d)
of the same figure. This might be explained by the fact that the states for
$n\geq 0$ experience the learned output matrix $W^{\rm out}_{\ast}$ via the
closed feedback loop. The states in the training phase, $n<0$, on the other
hand, do neither experience the fitted output matrix, nor is the last output
fed back to the reservoir. We suspect that this has a significant impact on
the evolution of the $s_{i}$.

Figure 9: Reservoir weight matrices: (a) Input weight matrix $W^{\rm in}$
which is a $4000\times 151$ matrix in our case. (b) Reservoir weight matrix
$W^{\rm r}$ which is a $4000\times 4000$ matrix in the present case. All
weights that are unequal to zero are marked as black dots. (c) Optimized
output weight matrix $W_{\ast}^{\rm out}$ which is a $150\times 4151$ matrix.
The aspect ratios of $W^{\rm in}$ and $W_{\ast}^{\rm out}$ have been adjusted
for illustration purposes. Figure 10: Reservoir state components
$s_{1},s_{1000},s_{4000}$ versus time step $n$ during (a) training and (b)
prediction phases. Panels (c) and (d) show their corresponding Fourier spectra
$\|\mathcal{F}(s_{i})\|^{2}$. Note that in (a), the first $46$ time steps are
not shown, as they were used for the initialization of the reservoir.

We now take a look at the reconstruction of the three fields $v_{x}(x,y)$,
$v_{y}(x,y)$ and $M(x,y)$ with the reservoir outputs as temporal coefficients
to see whether large-scale features are captured correctly. An instantaneous
snapshot at the half-time of the prediction phase at time step $n=350$ is
depicted in Fig. 11. We apply a POD-mode composition (24) to obtain the
fluctuation fields $v_{x}^{\prime}$, $v_{y}^{\prime}$ and $M^{\prime}$ from
the reservoir outputs $a_{p}(t)$. The mean profiles $\langle v_{x}\rangle_{t}$
and $\langle v_{y}\rangle_{t}$ and $\langle M\rangle_{t}$ are subsequently
added to obtain the full fields, see eqns. (20) and (23). The resulting ESN
predictions are displayed in panels (b), (d), and (f). For reference, the
validation (POD) fields are shown in panels (a), (c), and (e).

The horizontal velocity field $v_{x}(x,y)$ in (a) and (b) shows some
differences in the structure of the right- and left-going fluid patches, but
the large-scale structure as a whole is in surprisingly good qualitative
agreement. The structure of vertical velocity field $v_{y}(x,y)$ in panel (c)
and (d) does not show a systematic distinction, even though slight differences
in shape and maximum values of up- and downdrafts are detectable. Finally the
moist buoyancy field $M(x,y)$ in (f) does not fully reproduce all moist plumes
that detach from the bottom plate, see validation field in (e). Nevertheless
the predicted time coefficients lead to reconstructed fields that contain the
same features as the original fields.

Figure 11: Instantaneous snapshot of the fields (a, b) $v_{x}(x,y)$, (c, d)
$v_{y}(x,y)$ and (e, f) $M(x,y)$ at time step $n=350$ after the training
phase. Panels (a, c, e) are validation data from the POD and (b, d, f) the ESN
output data. The fields were reconstructed using the first 150 $a_{p}(n)$
(POD) and the predictions $y_{p}(n)$ (ESN). Here, $n$ is a discrete time step.
The ESN predictions deviate locally from the POD fields, but capture large-
scale features of the flow. The aspect ratio has been adjusted again for
illustration purposes. The corresponding colorbars can be seen on the right.

To get a better grasp on the time evolution of the error in the field
reconstruction, we compute a normalized field deviation at constant height $y$
of the vertical velocity field component in Fig. 12 which is given by

$\displaystyle{\rm Err}(x,n)=\frac{|v_{y}^{\prime\rm
POD}(x,n)-v_{y}^{\prime\rm ESN}(x,n)|}{\max_{x,y,n}\left(v_{y}^{\prime\rm
POD}\right)}\Bigg{|}_{y={\rm const}}$ (32)

where the superscript defines whether the field was reconstructed with
$a_{i}(n)$ (POD) or $y_{i}(n))$ (ESN). We find that the fast growing errors in
the time coefficients lead to fast amplifications of local field errors.
Furthermore, different horizontal and vertical positions in the domain show
different error magnitudes.

Figure 12: Time evolution of the prediction error ${\rm Err}(x,n)$, which is
given by eq. (32), at specified height $y$ (see top). As time progresses, the
errors start to grow in magnitude. Different positions $(x,y)$ in the domain
give rise to different magnitudes of the deviation.

We now discuss the ability of the ESN to reproduce the low-order statistical
properties of the turbulent flow. This is done by comparison of vertical line-
time average profiles $\langle\cdot\rangle_{x,t}(y)$. The averages are taken
along $x$-direction in combination with respect to time $t$. Such profiles are
for example of interest in larger-scale atmospheric models for the
parameterization of sub-grid scale transport Grabowski (2001); Khairoutdinov
and Randall (2001). Figure 13 depicts the profiles as a function of the domain
height $y$. The actual profiles obtained by the original DNS are plotted as a
dash-dotted, the POD reconstruction as a solid and the reconstruction from the
ESN outputs as a dashed line. Figure 13(a) shows the moist buoyancy $M$. Here
the time mean $\langle M\rangle_{t}$ was added to see whether the
reconstructed POD and ESN fields would deviate from the full DNS data. We
observe an excellent agreement and find that the ESN produces correct
fluctuations which preserve this profile. The fluctuations of the vertical
moist buoyancy flux $\langle v_{y}^{\prime}M^{\prime}\rangle$ are shown in
Fig. 13(b). Again, an excellent agreement between the curves in the bulk of
the domain and small deviations in the boundary layers only are observable,
despite the fact that this quantity is much more susceptible to errors since
it consists of two ESN estimates. Finally the fluctuations from the liquid
water content and the liquid water flux, further derived fields, are shown in
13(c) and (d), respectively. Here we see that POD and ESN curves match
throughout the whole of the domain. We thus conclude that the ESN is able to
reproduce essential low-order statistics very well.

For comparison, we add a comparison of test data with the output of an LSTM
network. The network parameters are as follows: the number of hidden states is
300 and the number of stacked LSTM layers 3. The loss function is again the
mean-squared error, the optimization applies the method of adaptive moments,
and the learning rate is $10^{-3}$ Goodfellow et al. (2016). A training of the
network proceeds over 1000 epochs. We can conclude that the LSTM performs
similarly well as the ESN even though the reproduced profiles deviate a bit
for both fluxes in the center of the convection layer.

Figure 13: Line-time averaged vertical profiles $\langle\cdot\rangle_{x,t}$ of (a) the full moist buoyancy $M$, (b) the vertical moist buoyancy flux $v_{y}^{\prime}M^{\prime}$, (c) the fluctuations of the liquid water content $q_{l}$ and (d) the vertical liquid water flux $v_{y}^{\prime}q_{l}^{\prime}$. The time average for DNS (dash-dotted) and POD (solid) were computed over the whole range of 1400 time steps, while for the ESN (dashed) and LSTM (dotted) only the range of the prediction phase ($700$ time steps) was taken into account. $\langle{\rm CC}^{\rm DNS}\rangle_{t}$ | $\langle q_{l}^{\rm DNS}\geq 0\rangle_{x,y,t}\cdot 10^{3}$ | $\langle{\rm CC}^{\rm POD}\rangle_{t}$ | $\langle q_{l}^{\rm POD}\geq 0\rangle_{x,y,t}\cdot 10^{3}$ | $\langle{\rm CC}^{\rm ESN}\rangle_{t}$ | $\langle q_{l}^{\rm ESN}\geq 0\rangle_{x,y,t}\cdot 10^{3}$
---|---|---|---|---|---
$89.43$% | $3.24$ |  $83.25$% | $2.72$ |  $82.49$% | $2.72$
Table 3: Time mean average $\langle\rm CC\rangle_{t}$ and $\langle q_{l}\geq
0\rangle_{x,y,t}$ of the cloud cover CC and the volume average of liquid water
for the DNS, POD and ESN case. See also Fig. 14. Figure 14: Clouds, i.e.
$q_{l}(x,y,n)\geq 0$ at time step $n=350$ of (a) the POD and the (b) ESN
prediction. The liquid water content $q_{l}$ differs in the magnitude and
shape of its envelope, the isosurface $q_{l}=0$. (c) shows the cloud cover as
defined in eq. (33) computed for the DNS data (brown), POD data (green) and
the ESN predictions (blue). The POD approximation does not capture all of the
original cloud cover; the value of the DNS exceeds the one of the POD by a few
per cent. The cloud cover prediction of the ESN itself deviates by a few per
cent in comparison to that of the POD. (d) shows the change of the volume
average of positive liquid water content $\langle q_{l}\geq 0\rangle_{x,y}$
with time.

Motivated by these results, we now investigate whether quantities such as the
cloud cover can be also modeled by the ESN. We define the cloud cover CC of
the two-dimensional domain $N_{x}^{\prime}\times N_{y}^{\prime}$ as the ratio
of the number of vertical grid lines $N_{q_{l}>0}$ that contain at least one
mesh point with $q_{l}>0$ along their vertical line of sight and the total
number of vertical grid lines, $N_{x}^{\prime}$. Thus follows

${\rm CC}=\frac{N_{q_{l}>0}}{N_{x}^{\prime}}\times 100\%.$ (33)

The time average $\langle{\rm CC}\rangle_{t}$ and volume-time average of
positive liquid water content $\langle q_{l}\geq 0\rangle_{x,y,t}$ are given
in Table 3. The truncation to the first $\rm N_{\rm POD}$ POD modes leads to a
loss of about $6.9\%$ of the original DNS CC and a $16.1\%$ loss of $\langle
q_{l}\geq 0\rangle_{x,y,t}$. We find good agreement between ESN estimate and
POD results. In Fig. 14, the POD and ESN results of the cloud distribution at
time step $n=350$ are shown. Despite small discrepancies in the local
distribution of the liquid water content and the shape of the cloud
boundaries, i.e. the isosurfaces $q_{l}=0$, the overall distribution is
comparable. In panels (c) and (d) of the same figure, the time evolution of
cloud cover and volume-time average positive liquid water content are
displayed. While the predicted cloud cover does not deviate too much from the
reference case, the variations in the amount of liquid water are less well
reproduced.

##  V Summary and Conclusion

In the present work, we have applied a machine learning algorithm to two-
dimensional turbulent moist Rayleigh-Bénard convection, in order to infer the
large-scale evolution and low-order statistics of convective flow undergoing
phase changes. We apply a specific learning scheme for recurrent neural
networks, called echo state network, which has been applied for learning the
dynamics of nonlinear systems, such as turbulent shear flows. Here, we test
its capabilities successfully by fitting a reservoir to complex convection
flow dynamics which results by the interaction of turbulence with the
nonlinear thermodynamics originating from the first order phase changes
between vapor and liquid water as present for example in atmospheric clouds.
We therefore generate comprehensive data by means of a simple moist convection
model in the Boussinesq framework, the moist Rayleigh-Bénard convection model.

We obtain moist convection data from direct numerical simulations. As the 2d
set of data has still a large amount of degrees of freedom and therefore
cannot be passed directly to the echo state network, we introduce the POD as
an additional dimensionality reduction step. We therefore decompose the data
into a data-driven, spatially dependent basis with temporal coefficients,
where the latter are fed to the reservoir. We truncate the POD to its most
energetic modes and coefficients, reducing the degrees of freedom of the
dynamical system at hand considerably. This reduced data set serves as the
ground truth for the training of the echo state network as well as validation
of its outputs. The network setup is tuned by conducting a comprehensive grid
search of important hyperparameters. By coupling the output of the trained
network back to its input, the autonomous system estimates the evolution of
the temporal coefficients. Reconstructing the velocity and thermodynamic
fields from these estimates, allow us to check whether the dynamics have been
learned. We find an excellent agreement of the vertical profiles of moist
buoyancy, vertical moist buoyancy transport as well as liquid water content.
Furthermore, we report a good agreement of essential parameters in moist
convection such as the fraction of clouds covering the two-dimensional
atmosphere, as well its content of liquid water.

This first approach of our reservoir computing model to moist Rayleigh-Bénard
convection shows, its potential to infer low-order statistics from a set of
training data. Though the reservoir output quickly diverges from the actual
system trajectory, time averaged quantities are robustly reproduced. This
result might seem trivial at first glance, yet the reservoir produces velocity
and thermodynamic fluctuation fields which do not deviate too strongly from
those of the original flow, even for combined quantities such as the liquid
water flux across the layer. This indicates that the present echo state
network did not just learn the statistics, but the dynamical system itself.
Our additional comparison with an LSTM network gives a similar outcome. A more
detailed comparison of both RNN implementations has to be left however as a
future work.

Our approach can be considered as a first step of applying reservoir computing
as a machine learning-based parameterization. General circulation models
already use multi-scale modeling methods where small-scale resolving models
interact with the large-scale motion by their low-order statistics,
essentially relaxing one to each other, e.g. in superparametrizations
Grabowski and Smolarkiewicz (1999); Grabowski (2001); Khairoutdinov and
Randall (2001). An echo state network can serve as a simple dynamical
substitute for the unresolved subgrid scale transport.

Even though the present results are promising, the development is still in its
infancy. We state that for example the mathematical foundations of reservoir
computing, which could provide deeper insights on the role of the
hyperparameters on the prediction quality, are still mostly unexplored.
Moreover, we reckon that for an extension of the ESN approach to a three
dimensional flow, the data reduction step via POD will not suffice to cope
with the large amount of simulation data. For this scenario one might propose
the usage of a convolutional autoencoder/-decoder neural network in
combination with the RC model. Furthermore, we mention that the machine
learning algorithm is supposed here to learn dynamics of a nonlinear system
which incorporates processes on different spatial and temporal scales. This
circumstance is so far not fully captured by the network architecture.
Particularly for turbulence, this might imply that the different spatial
scales which interact with each other and exchange their energy, could be
trained separately allowing for a subsequent coupling. The exploration of such
ideas is currently under way and will be reported elsewhere.

## Acknowledgments

This work is supported by the project ”DeepTurb – Deep Learning in and of
Turbulence” which is funded by the Carl Zeiss Foundation. The authors
gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-
centre.eu) for funding this project by providing computing time through the
John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS
at Jülich Supercomputing Centre (JSC). We thank Martina Hentschel, Erich
Runge, and Sandeep Pandey for helpful comments.

## References

  * Brenner et al. (2019) M. P. Brenner, J. D. Eldredge, and J. B. Freund, Phys. Rev. Fluids 4, 100501 (2019).
  * Brunton et al. (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Annu. Rev. Fluid Mech. 52, 477 (2020).
  * Pandey et al. (2020) S. Pandey, J. Schumacher, and K. R. Sreenivasan, J. Turbul. 21, 567 (2020).
  * Ling et al. (2016) J. Ling, A. Kurzawski, and J. Templeton, J. Fluid Mech. 807, 155 (2016).
  * Duraisamy et al. (2019) K. Duraisamy, G. Iaccarino, and H. Xiao, Annu. Rev. Fluid Mech. 51, 357 (2019).
  * Fonda et al. (2019) E. Fonda, A. Pandey, J. Schumacher, and K. R. Sreenivasan, Proc. Natl. Acad. Sci. 116, 8667 (2019).
  * Brenowitz and Bretherton (2018) N. D. Brenowitz and C. S. Bretherton, Geophys. Res. Lett. 45, 6289 (2018).
  * O’Gorman and Dwyer (2018) P. A. O’Gorman and J. G. Dwyer, J. Adv. Model Earth Sy. 10, 2548 (2018).
  * Gentine et al. (2018) P. Gentine, M. Pritchard, S. Rasp, G. Reinaudi, and G. Yacalis, Geophys. Res. Lett. 45, 5742 (2018).
  * Moller et al. (2020) S. Moller, C. Resagk, and C. Cierpka, Exp. Fluids 61, 111 (2020).
  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, 2016).
  * Schneider et al. (2017) S. Schneider, T. ans Lan, A. Stuart, and J. Teixeira, Geophys. Res. Lett. 44, 12396 (2017).
  * Mohan et al. (2020) A. Mohan, D. Tretiak, M. Chertkov, and D. Livescu, J. Turbul. 21, 525 (2020).
  * Stevens (2005) B. Stevens, Annu. Rev. Earth Planet. Sci. 33, 605 (2005).
  * Mellado (2017) J. P. Mellado, Annu. Rev. Fluid Mech. 49, 145 (2017).
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber, Neural Comput. 9, 1735 (1997).
  * Lukoševičius and Jaeger (2009) M. Lukoševičius and H. Jaeger, Comp. Sci. Rev. 3, 127 (2009).
  * Tanaka et al. (2019) G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Neural Netw. 115, 100 (2019).
  * Jaeger and Haas (2004) H. Jaeger and H. Haas, Science 304, 78 (2004).
  * Yildiz et al. (2012) I. B. Yildiz, H. Jaeger, and S. J. Kiebel, Neural Netw. 35, 1 (2012).
  * Pandey and Schumacher (2020) S. Pandey and J. Schumacher, Phys. Rev. Fluids 5, 113506 (2020).
  * Grabowski and Smolarkiewicz (1999) W. W. Grabowski and P. K. Smolarkiewicz, 133, 171 (1999).
  * Grabowski (2001) W. W. Grabowski, J. Atm. Sci. 58, 978 (2001).
  * Vlachas et al. (2020) P. R. Vlachas, J. Pathak, B. R. Hunt, T. P. Sapsis, G. M., E. Ott, and P. Koumoutsakos, Neural Netw. 126, 191 (2020).
  * Lu et al. (2017) Z. Lu, J. Pathak, B. Hunt, M. Girvan, R. Brockett, and E. Ott, Chaos 27 (2017).
  * Pathak et al. (2018) J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Phys. Rev. Lett. 120, 024102 (2018).
  * Holmes et al. (2012) P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley, _Turbulence, Coherent Structures, Dynamical Systems and Symmetry_ , Cambridge Monographs on Mechanics (Cambridge University Press, Cambridge, UK, 2012), 2nd ed.
  * Srinivasan et al. (2019) P. A. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, and R. Vinuesa, Phys. Rev. Fluids 4, 054603 (2019).
  * Nguyen et al. (2020) T. Nguyen, R. Baraniuk, A. Bertozzi, S. Osher, and B. Wang, in _Advances in Neural Information Processing Systems_ , edited by H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (Curran Associates, Inc., 2020), vol. 33, pp. 1924–1936.
  * Fischer (1997) P. F. Fischer, J. Comput. Phys. 133, 84 (1997).
  * Scheel et al. (2013) J. D. Scheel, M. S. Emran, and J. Schumacher, New J. Phys. 15, 113063 (2013).
  * nek (2017) _nek5000 version 17.0_ (2017), URL https://nek5000.mcs.anl.gov.
  * Pauluis and Schumacher (2010) O. Pauluis and J. Schumacher, Commun. Math. Sci. 8, 295 (2010).
  * Weidauer et al. (2010) T. Weidauer, O. Pauluis, and J. Schumacher, New J. of Phys. 12, 105002 (2010).
  * Schumacher and Pauluis (2010) J. Schumacher and O. Pauluis, J. Fluid Mech. 648, 509–519 (2010).
  * Bretherton (1987) C. S. Bretherton, J. Atmos. Sci. 44, 1809 (1987).
  * Bretherton (1988) C. S. Bretherton, J. Atmos. Sci. 45, 2391 (1988).
  * Smith and Stechmann (2017) L. M. Smith and S. N. Stechmann, J. Atmos. Sci. 74, 3285 (2017).
  * Vallis et al. (2019) G. K. Vallis, D. J. Parker, and S. M. Tobias, J. Fluid Mech. 862, 162–199 (2019).
  * Abma et al. (2013) D. Abma, T. Heus, and J. P. Mellado, J. Atmos. Sci. 70, 2088 (2013).
  * Emmanuel (1994) K. A. Emmanuel, _Atmospheric Convection_ (Oxford University Press, 1994).
  * Pauluis and Schumacher (2011) O. Pauluis and J. Schumacher, Proc. Natl. Acad. Sci. USA 108, 12623 (2011).
  * Jaeger (2002) H. Jaeger, GMD-Forschungszentrum Informationstechnik (2002).
  * Doya (1992) K. Doya, [Proceedings] 1992 IEEE International Symposium on Circuits and Systems 6, 2777 (1992).
  * Jaeger (2001) H. Jaeger, GMD-Forschungszentrum Informationstechnik Technical Report 148 (2001).
  * Maass et al. (2002) W. Maass, T. Natschläger, and H. Markram, Neural Computation 14, 2531 (2002).
  * Lin et al. (2008) X. Lin, Z. Yang, and Y. Song, in _Advances in Knowledge Discovery and Data Mining_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008), pp. 932–937.
  * Morando et al. (2013) S. Morando, S. Jemei, R. Gouriveau, N. Zerhouni, and D. Hissel, in _IECON 2013 - 39th Annual Conference of the IEEE Industrial Electronics Society_ (2013), pp. 1632–1637.
  * Pathak et al. (2017) J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, Chaos 27, 121102 (2017).
  * Lukoševičius (2012) M. Lukoševičius, LNCS 7700, 659 (2012).
  * Strauss et al. (2012) T. Strauss, W. Wustlich, and R. Labahn, Neural Comput. 24, 3246 (2012).
  * Sirovich (1987) L. Sirovich, Q. Appl. Math. XLV, 561 (1987).
  * Khairoutdinov and Randall (2001) M. F. Khairoutdinov and D. A. Randall, Geophys. Res. Lett. 28, 3617 (2001).