{'fraction_non_alphanumeric': 0.036143367013522346, 'fraction_numerical': 0.010782687377469461, 'mean_word_length': 4.462996018124399, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Echo State Network for two-dimensional turbulent moist Rayleigh-B\'enard convection

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\'enard 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. The training data comprise long time series of the first 150 most energetic POD coefficients. The reservoir is subsequently fed by the data and results in predictions of future flow states. The predictions are thoroughly validated by the data of the original simulation. 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 [1][2][3]. The applications reach from subgrid-scale stress models [4,5], via the detection of large-scale patterns in mesoscale convection [6] to MLbased parametrizations of turbulent transport and clouds in climate and global circulation models [7][8][9]. The implementations of such algorithms help to process and reduce increasing amounts of data coming from numerical simulations and laboratory experiments [10] by detecting patterns and statistical correlations [11]. 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 [12,13].
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 [14,15], 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 shortterm 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 [16][17][18].
In this paper we apply a specific implementation of an RNN, the echo state network (ESN) [19,20], to twodimensional (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 [21] 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 [22,23]. (4) Finally, the hyperparameters of the ESN, in particular the spectral radius ρ(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.
low-order flow statistics in 2d dry Rayleigh-Bénard convection with 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 [27]. 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 shortterm memory networks (LSTM) [16] which have been also applied to fluid mechanical problems, such as the dynamics of Galerkin models of shear flows in [28]. These models also demonstrated to capture the longer-term time dependency and low-order statistics of important turbulent transport properties well (see also [3] 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 [29].
In this work, a moist Rayleigh-Bénard convection model with moist Rayleigh number Ra M 10 8 and Prandtl number Pr = 0.7 is considered as an example case. We choose a 2d domain Ω = L × 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 × N y = 7200 × 384. The data are obtained by direct numerical simulations (DNS) which apply a spectral element method [30][31][32]. 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.

A. 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 [33][34][35]. The framework is based on the mathematically equivalent formulation by Bretherton [36,37]. Similar simplified models of moist convection with precipitation have been developed by Smith and Stechmann [38] and Vallis et al. [39]. Evaporative cooling and buoyancy reversal effects were discussed for example by Abma et al. [40].
The buoyancy B in atmospheric convection is given by [41] with the gravity acceleration g, a mean density ρ, 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 [33] B(x, y, t) = max M (x, y, t), D(x, y, t) − N 2 s y where the fixed Brunt-Väisälä frequency N s = g(Γ u − Γ s )/T ref is determined by the lapse rates of the saturated and unsaturated moist air, Γ s and Γ u , and a reference temperature T ref . An air parcel at height y and time t is unsaturated if M (x, y, t) < D(x, y, t) − N 2 s y and saturated if M (x, y, t) > D(x, y, t) − N 2 s 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 D(y), M (y) are vertically linear where Here v = (v x (x, y, t), v y (x, y, t)) T is the two-dimensional velocity field,p = p/ρ ref the kinematic pressure, ν the kinematic viscosity and κ the scalar diffusivity. The term d/dt = ∂/∂t + (v · ∇) is the material derivative. This idealized model describes the formation of warm, nonprecipitating low clouds in a shallow layer up to a depth of ∼ 1km. 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 An additional parameter arises from the additional phase changes, the dimensionless form of (2) 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 2 s y. There is no supersaturation considered in this model and the liquid water content field q l and thus the clouds are given by 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 ≥ 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 ≤ D H we refer to refs. [36,37] or [42].

B. 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 [43] and bifurcations [44]. An alternative training method was proposed by Jaeger [45] and in an independent work by Maass [46]. Their idea, which is now known as reservoir computing [17], 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 [47] to estimating the life span of a fuel cell [48]. 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 [25]. 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 [49]. 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 u(n) ∈ R Nin 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 in ∈ R Nr×(1+Nin) . The weighted input contributes to the dynamics of the reservoir state s ∈ [−1, 1] Nr at time n which is given by Here [b; u(n)] stands for the vertical concatenation of the scalar bias and the input [54]. This update rule comprises external forcing by the inputs u(n) as well as a self- interaction with the past instance s(n − 1). The reservoir weight matrix W r ∈ R Nr×Nr blends the state dimensions, while the nonlinear hyperbolic tangent tanh (·), applied to each component of its argument vector, is the nonlinear activation function of the neurons in this model. The leaking rate γ ∈ (0, 1] moderates the linear and nonlinear contributions and assures that the state is confined to [−1, 1] Nr . 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 s(n − 1), s (n − 1) converge to the same state s(n), provided the same input u(n) is given and the system has been running for many iterations n [45]. 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 statẽ s(n) = [b; u(n); s(n)]. Finally,s is mapped via the output matrix W out ∈ R Nin×(1+Nin+Nr) to form the reservoir output y(n) ∈ R Nin y(n) = W outs (n).
For our application the output dimension matches the input dimension N 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 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 train training data instances {u(n), y target (n)} (where n = −n train , −(n train − 1), ..., −1) needs to be prepared. The target output y target (n) represents the desired output the ESN should produce for the given input u(n). The reservoir state 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 (W out ) which is given by and has to be minimized corresponding to Here w out i denotes the i-th row of W out and · 2 the L 2 norm. Equations (16) and (17) are known as ridge regression with the ridge regression parameter β. The last term in (16) suppresses large values of the rows of W 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 T denotes the transposed and Id the identity matrix. Y target and S are matrices where the n-th column is the target output y target (n) and the extended reservoir states(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 in and W r are hyperparameters which have to be tuned beforehand. In our approach we draw the elements of W in from a uniform distribution in [−0.5, 0.5] and impose no further restrictions on this matrix. For the generation of the reservoir weights in W r it has been reported that the proportion of non-zero elements, i.e., the reservoir density D and the spectral radius (W r ), i.e., the largest absolute eigenvalue of W r are crucial parameters that determine whether the desired echo state property holds [50]. We choose a sparse reservoir (D < 1) with few internal node connections and draw the elements from a uniform distribution in [−1, 1]. We then normalize W r by its largest absolute eigenvalue and multiply it with the desired spectral radius (W r ). This scaling approach, initially proposed by Jaeger [43], has established itself as one of the standard ways [21,25,26,48] to control the spectral radius. Nevertheless, other procedures have been proposed [20,51]. In addition, the size of the reservoir N r is a hyperparameter. Usually s should be a high-dimensional extension of the inputs u, so that N r N in is satisfied. Moreover, we consider the leaking rate γ and ridge regression parameter β as further quantities that have to be adjusted to our data.

III. ECHO STATE NETWORK FOR 2D MOIST CONVECTION
A. Direct Numerical Simulations DNS using the spectral element solver Nek5000 [30][31][32] were conducted to solve the two-dimensional moist Rayleigh-Bénard system (5)

B. 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 ∆t = 0.25 in the statistically sta-tionary 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 × N y = 640 × 80 points from the originally unstructured element mesh. They are decomposed into temporal mean and fluctuations subsequently, A grid dimension of 640 × 80 in terms of ESN input dimensions N in is still too big. We therefore follow the approach in [21] 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 [52] which decomposes the k-th component of a vector field g(x, y, t) into Here a p (t) and Φ (p) k (x, y) are the p-th time coefficient and the corresponding spatial POD mode respectively. In our approach we take the POD of g = (v x , v y , D , M ) 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 N POD = 150 n 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 Ra = 10 7 [21]. This is one order of magnitude below our Rayleigh number values. Nevertheless our RC model will receive values for all N 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 Φ . We limit ourselves to the last 1400 = n train + n test time instances of our data and use the first 150 time coefficients a(n) = (a 1 (n), a 2 (n), ..., a 150 (n)) T as the input for the ESN. The first n train = 700 snapshots are assembled into a training data set {u(n), y target (n)} = {a(n), a(n + 1)} (25) with −n train ≤ n ≤ −1. The training data span 175 freefall time units T f that correspond to 61 eddy turnover times. This time scale is given by τ eddy = H/u rms ≈ 2.9T f with u rms = u 2 x + u 2 y 1/2 x,y,t . The ESN is trained to predict the time instance a(n + 1) when given the POD time coefficients 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 out , which can then be used for prediction. In the prediction phase, we give the initial input u(0) = a(0) to the reservoir and redirect the output to the input layer (see Fig. 1(b)) such that u(n) = y(n − 1) n = 1, ..., n test − 1 with n 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 n test snapshots {y val (n)} = {a(n + 1)}.

C. Training of ESN and choice of hyperparameters
We quantify the quality of ESN predictions, with the set of hyperparameters h = {γ, β, N r , D, (W r )}, by two types of measures. The mean square error MSE h of ESN output to the validation data where is the mean square error at time step n averaged over all N POD modes. Additionally, we take the physically more relevant normalized average relative error (NARE) as defined in [28] into account. For the moist buoyancy field M , it is for example given by with .
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 · x,t of a specific field. We consider the three NAREs: that is the NARE of the total moist buoyancy field M , the liquid water content fluctuations q l , and the moist buoyancy flux fluctuations v y M .
A grid search for the four quantities γ, β, D, (W r ) was conducted in a suitable range (see Table I), based on the results in [21]. The reservoir size was fixed to N r = 4000 for all runs. The resulting MSE h and NAREs were computed to find an adequate parameter setting. Figure 5 shows MSE h , E h q 2 l x,t , E h v y M x,t and their dependence on (W r ) and γ. 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,   FIG. 5: Representative profiles taken from the error landscape for the leaking rate γ (a,b,c) and the spectral radius (W r ) (d,e,f). The data are obtained by a grid search study (see Table I). 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. 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 in-creasingly binary way since lim x→±∞ tanh(x) = ±1. In this limit of a fully saturated activation function, eq. (14) simplifies to where 1 = (1, 1, ..., 1) T ∈ R Nr . This corresponds to a linear dependence of each reservoir state on its last instance plus the constant leaking rate γ 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 ρ(W r ) 1 which in turn contribute to the increasing MSE h in panel (d) of Fig. 5. Finally, we chose the hyperparameter set h * , listed in Table II as it results in a minimum of E h v y M x,t . The reason for settling with this measure is that it is susceptible to two ESN estimates. Quantities like E h q 2 l x,t , 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
After the ESN receives the initial input u(0) = 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 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 (14). Figure 8(b) shows the deviations (y val i (n) − y i (n)) for i = 1, 10, 50, 100. Figures 9(a)-(c) show the three weight matrices of the RC model. As described in section II B, the input and reservoir weights are initialized randomly and left unchanged. With a reservoir density of D = 0.2 the reservoir matrix W r is a sparse matrix containing many weights equal to zero. Note further, that the fitted output weights W out have low magnitudes in comparison to the entries of the input matrix W in . This is adjusted according to the number of reservoir nodes N r and the magnitude of the training data. Moreover, the magnitude of the first 1 + N in columns are close to zero. This indicates that the contributions of the output bias b and the current input 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 Train 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 .
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 PODmode composition (24) to obtain the fluctuation fields v x , v y and M from the reservoir outputs a p (t). The mean profiles v x t and v y t and M 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.
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 where the superscript defines whether the field was reconstructed with a i (n) (POD) or y i (n)) (ESN). We find 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 · x,t (y). The averages are taken along xdirection 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 [23,53]. 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 M 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 v y M 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  (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. 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 [11]. 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.
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 twodimensional domain N x × N y 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 . Thus follows The time average CC t and volume-time average of positive liquid water content q l ≥ 0 x,y,t are given in Table  III. The truncation to the first N POD POD modes leads to a loss of about 6.9% of the original DNS CC and a 16.1% loss of q l ≥ 0 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 largescale 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 [22,23,53]. 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.