Low-dimensional representation of intermittent geophysical turbulence
  with High-Order Statistics-informed Neural Networks (H-SiNN)
We present a novel machine learning approach to reduce the dimensionality of
state variables in stratified turbulent flows governed by the Navier-Stokes
equations in the Boussinesq approximation. The aim of the new method is to
perform an accurate reconstruction of the temperature and the three-dimensional
velocity of geophysical turbulent flows developing non-homogeneities, starting
from a low-dimensional representation in latent space, yet conserving important
information about non-Gaussian structures captured by high-order moments of
distributions. To achieve this goal we modify the standard Convolutional
Autoencoder (CAE) by implementing a customized loss function that enforces the
accuracy of the reconstructed high order statistical moments. We present
results for compression coefficients up to 16 demonstrating how the proposed
method is more efficient than a standard CAE in performing dimensionality
reduction of simulations of stratified geophysical flows characterized by
intermittent phenomena, as observed in the atmosphere and the oceans.
Raffaello Foldes, Enrico Camporeale, Raffaele Marino
# Low-dimensional representation of intermittent geophysical turbulence with H-SiNN

###### Abstract

We present a novel machine learning approach to reduce the dimensionality of state variables in stratified turbulent flows governed by the Navier-Stokes equations in the Boussinesq approximation. The aim of the new method is to perform an accurate reconstruction of the temperature and the three-dimensional velocity of geophysical turbulent flows developing in homogeneities, starting from a low dimensional representation in latent space, yet conserving important information about non-Gaussian structures, captured by high-order moments of distributions, needed to reproduce embedded dynamics To achieve this goal we modify the standard Convolutional Autoencoder (CAE) by implementing a customized loss function that enforces accuracy of the reconstructed high order statistical moments. We present results for compression coefficients up to 16 demonstrating how the proposed method is more efficient than a standard CAE in performing dimensionality reduction of simulations of stratified geophysical flows characterized by intermittent phenomena, as observed in the atmosphere and the oceans.

+
Footnote †: preprint: Low-dimensional representation of intermittent geophysical turbulence with H-SiNN

## I Introduction

Geophysical fluids are characterized by the interplay of non-linear vortices and waves, developing very complex turbulent dynamics due to the density stratification and the Earth's solid body rotation. Large scales motions at planetary scales (of the order of 1000 km) emergence from the so-called quasi-geostrophic balance [1]. At mesoscale \(O(100\,\mathrm{km})\) and even more at the sub-mesoscale \(O(10\,\mathrm{km})\), when the turbulent eddy turnover time becomes comparable to the characteristic time scales associated to internal waves, motions feel less the constrains of the force balance, a regime that is generically identified as _stratified turbulence_[2]. This is characterized by the presence of shears leading to instabilities, the vertical extension of strata being controlled by the Brunt-Vaisala frequency \(N\). In stratified geophysical flows the probability density function (PDF) of fields fluctuations at the scales comparable with that of the mean flow can be characterized by fat tails, with a departure from Gaussianity. In particular, this behavior has been observed in the vertical velocity and/or (potential) temperature fluctuations in both the atmosphere [3; 4; 5] and the ocean [6; 7], and extensively in a variety of numerical investigations [8; 9; 10; 11; 12; 13] through the forth-order moment (i.e., kurtosis) of the vertical velocity. Marino _et al._[13] have shown how the extreme drafts responsible for the non-Gaussian behavior of the vertical component of the velocity do generate local turbulence and enhance the classical small-scale intermittency [11]. They provided as well evidence that these structures are associated with patches of enhanced kinetic and potential energy dissipation, making stably stratified flows within the regime analyzed in [9] more efficient at dissipating energy but also non-homogeneous due to the presence of strong vertical drafts. The community resorted to a variety of different techniques to model the scales at which turbulence dominates in geophysical flows, from the implementation of codes based on Reynolds averaged Navier-Stokes (RANS) equations and Large Eddy Simulations schemes (LES) [14; 15; 16], to the simulation of geophysical turbulence using a reduced models [17]. On the other hand, progresses made in the field of High Performance Computing have made it possible to perform direct numerical simulations (DNS) of stratified turbulent flows in a parameter space of geophysical interest [18; 19; 20]. Disregarding the physical model implemented, common trait of high-resolution three-dimensional simulations of geophysical flows is the large amount of data produced. For example, a DNS run with resolution of \(512^{3}\) grid points in single precision requires more than 0.5Gb of storage space per field and per time-step. Retaining the output with high temporal cadence is therefore in many case not doable, which requires to find a compromise between spatial and temporal resolutions, though investigating the evolution of the fields in space and time is often needed in order to assess the complex dynamics developing in from the need to perform run-time analysis, as well as for the implementation of machine learning techniques able to learn features of the dynamical systems from the latent space [21; 22; 23]. The rationale behind the dimensionality reduction of the output of turbulent flows simulations is therefore twofold: on the one hand this is done in order to store data for post-processing at higher temporal frequency, while it also enables the efficient training of data-mining and machine learning techniques. The study of low dimensional manifolds has a long history in dynamical systems and in turbulence [24; 25], laying at the foundation of so-called reduced order methods (ROMs) [26]. In turn, ROMs have also recently benefited by progresses in machine learning based approaches [27; 28]. Dimensionality reduction techniques such as proper orthogonal decomposition [29; 30] and dynamic mode decomposition [31; 32] have been extensively used to address is
sues concerning a variety of flows [33; 34; 35; 36; 37]. Other examples of model order reduction techniques applied to fluid frameworks are Galerkin-projection based nonlinear methods [38] and system identification based auto-regressive models [39; 40]. Techniques based on machine learning principles that have been employed in modeling the equations of motion are Gaussian process regression [41; 42], symbolic regression [43; 44], field inversion [45; 46], and artificial neural networks (ANNs). In data-driven methods, ANNs are the most powerful tools in terms of their complexity, ability to generalize and to capture highly nonlinear phenomena, typical of turbulent processes. Other developments have focused on improving reliability and accuracy of low fidelity models (i.e., RANS and LES), by using data from high fidelity simulations, e.g. DNS, to learn proper closures [47; 48]. Convolutional neural networks (CNNs) have been successfully applied to the identification of flow structures [49] and to perform nonlinear modal decomposition in turbulent flows [50; 51]; more recently, CNN-based architectures have been extended to fully 3D DNS of turbulent channel flow [52] and to numerical simulations of flows characterized by the emergence of extreme events, in terms of boundary coherent structures [53]. These implementations emphasized the advantages of using CNNs over traditional methods based on principal component analysis (PCA), the former being able to capture the intrinsic nonlinear dynamics of turbulent flows. More generically, deep learning proved to be a powerful tool for the analysis and generation of reduced-order models of turbulent systems. In terms of machine learning applications to stratified flows, CNN-based deep learning has recently been applied by [54] to get parametrization of diapycnal mixing using data from DNS. Convolutional AutoEncoders (CAE) are a family of autoencoders that proved successful in extracting information from two- and three-dimensional data [55] and their use has recently been proposed in the context of fluid-dynamics [56; 57; 58]. However, most of these approaches present severe limitations if implemented without proper tuning, and their use is often constrained to idealized and/or simplified cases due to the major constraints that come from capturing transient behavior [59]. Conversely, the presence of large-scale transient phenomena developing in geophysical flows, such as hurricanes, tornadoes, oceanic fronts, etc., can be viewed as an obstacle to the creation of reduced-order manifolds sufficiently informative to reliably recover properties and dynamical features of geophysical fluid systems. Multi-scale models, such as global circulation models (GCM), are powerful tools for the investigation of the Earth's atmosphere and the oceans which would greatly benefit from the possibility to simulate complex fluid dynamics on a dimensional phase space lower than what can be achieved numerically. The efficiency of machine learning approaches is reduced when applied to systems with strong transients, such as the non-stationary states resulting from the presence of intermittent phenomena in turbulent flows. The purpose of the study presented here is precisely to develop a machine learning tool based on CAE able to produce reliable low-order representations of stratified geophysical turbulent flows in a range of parameters characterized by developing intermittent (in space and time) vertical velocity drafts, leading to the departure of the vertical velocity statistical moments from the Gaussian reference. We show that standard CAE architectures are well suited to learn a latent space containing enough information to reconstruct the general features of the original stratified turbulent flow. However, strong vertical drafts developing in a certain regime of the governing parameter, and in general high order statistics of velocity and temperature in presence of large- and small-scale intermittency, will typically not be well recovered. To overcome this issue, we introduce here a novel implementation of CAE: statistics-informed convolutional autoencoder (or SiCAE). The idea is based on a more general approach, applicable to any neural network, which we refer to as high-order statistics-informed neural network (H-SiNN). The concept of statistics-informed neural network (SiNN) has been first introduced by [60], where it has been proposed the addition of two terms to the loss function whose effect is to constrain both the PDF and the auto-correlation function of the reconstructed fields. On the other hand, in our implementation, the loss function explicitly enforces the preservation of high-order moments and this is done in order to ensure consistency between the statistics of the original and the reconstructed flow fields. Indeed, it is worth recalling that statistical moments of the velocity distribution function play a crucial role in the characterization of turbulent frameworks, being directly related to invariants such as energy and enstrophy, to other fundamental quantities, such as the dissipation occurring at the small scales. We will demonstrate here the enhanced capability of the H-SiNN in recovering finer statistical features as well as flow inhomogeneitis and intermittent dynamics with respect to a standard CAE implementation.

The manuscript is organized as follows: in Sec. II the numerical simulations performed and analyzed are briefly introduced along with a description of the equations and parameters governing the dynamics of the stratified flows under study; Sec. III explains the classical CAE architecture, and in Sec. IV our implementation of statistical-informed CAE is thoroughly described; finally, Sec. V shows in detail the outcome of the application of the model reduction tool we developed to DNS of stratified turbulent flows developing large-scale intermittent events in the vertical component of the velocity and in the temperature field [9].

## II DNS solver for stratified turbulence

The CAEs are trained using the output of a direct numerical simulation of a stably stratified turbulent flow obtained by integrating the Navier-Stokes equations in the Boussinesq approximation in which the velocity field \(\mathbf{u}=(\mathbf{u}_{\perp},w)\) remains incompressible, \(\nabla\cdot\mathbf{u}=0\), while small density variations are taken into account only in the buoyancy term. Such a model can be written as:

\[\partial_{t}\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p-N\theta\mathbf{\hat{z}}+\nu\nabla^{2}\mathbf{u}+ \mathbf{F} \tag{1}\] \[\partial_{t}\theta+(\mathbf{u}\cdot\nabla)\theta = Nw+\kappa\nabla^{2}\theta. \tag{2}\]

with \(\theta\) being the temperature fluctuations evaluated relative to a mean temperature profile \(\theta_{0}\), \(N=[-g\partial_{z}\bar{\theta}/\theta_{0}]^{1/2}\) is the
Brunt-Vaisala frequency. Such a simulation has been initialized with zero temperature fluctuations and random velocity modes applied at large-scale, in a Fourier shell centered at \(k_{0}=[2,3]\); the size of the computational box is associated with \(k_{min}=1\), with \(k_{min}=2\pi/L_{0}\) and the resolution is \(512^{3}\) grid-points. A random forcing \(F\) is imposed to the momentum equation at \(k_{F}=2\pi/L_{f}\in[2,3]\), continuously injecting kinetic energy into the system and allowing to reach a turbulent stationary state. In the above equations, \(\nu\) and \(\kappa\) are the kinematic viscosity and the thermal diffusivity, respectively and we take the Prandtl number \(\nu/\kappa=1\), with \(\nu=10^{-3}\); finally, \(p\) is the pressure. We adopt the following definitions for the dimensionless Reynolds and Froude numbers,

\[\mathrm{Re}=U_{rms}L_{int}/\nu,\ \mathrm{Fr}=U_{rms}/[L_{int}N]\, \tag{3}\]

with \(U_{rms}\) and \(L_{int}\) the characteristic (rms) velocity and the integral scale of the fluid, respectively. These parameters, and in particular the buoyancy Reynolds \(R_{B}=\mathrm{ReFr}^{2}\), measure the relative strength of buoyancy to dissipation and are commonly used to distinguish between wave-dominated and turbulence-dominated regimes [61]. The equations (1-2) are integrated numerically using the Geophysical High-Order Suite for Turbulence (GHOST), a pseudo-spectral code which employs a hybrid parallelization combining MPI, OPENMP and CUDA [62; 63]. It allows for a variety of physical solvers and supports non-cubic geometry [64], and non-periodic boundary conditions [65]. In the following, we will make use of the statistical moments of a distribution function up to the forth-order. The definition of the third- and forth-order moments, namely the skewness \(Sk_{w}\) and kurtosis \(K_{w}\) respectively, of the distribution of vertical velocity \(w\) is:

\[Sk_{w}=\frac{\left\langle(w-\bar{w})^{3}\right\rangle}{\left\langle(w-\bar{w} )^{2}\right\rangle^{3/2}}, \tag{4}\]

\[K_{w}=\frac{\left\langle(w-\bar{w})^{4}\right\rangle}{\left\langle(w-\bar{w} )^{2}\right\rangle^{2}}, \tag{5}\]

where averages can be either computed over the simulation volume or others ensembles, such as horizontal planes. It is worth recalling that, \(Sk=0\) and \(K=3\) are the reference values for a Gaussian distribution. On the contrary, for values of \(K_{w}\) higher than 3, PDFs will display fat (non-Gaussian) tails, indicating that the variable presents extreme values which deviate consistently from the mean value. For this reason, the kurtosis \(K\) has been largely employed to characterize the emergence of large-scale intermittency in the vertical velocity field \(w\) (and temperature \(\theta\)) in numerical simulations [9]. As observed by Feraco _et al._[9], \(K_{w}\) presents strong non-Gaussian values in a narrow range of the Froude number, with an evident peak at \(Fr\approx 0.076\), being compatible with values observed in the upper atmosphere and oceans. The following analysis is performed using the same run which has been comprehensively studied in Feraco _et al._[9], and subsequent works [11; 13], presenting the highest level of intermittency (see Fig.2 of Feraco _et al._[9]).

## III Convolutional autoencoder theory

An autoencoder is an unsupervised feed forward artificial neural network (ANN) that aims to reconstruct a given data set through a process involving data compression and recovery [66]. Indeed, its main objective is to create a reliable and reduced representation of the input data set, being particularly helpful for either creating reduced model or simpler but still informative representation of high-dimensional data.

The network consists of an encoder part \(P:\mathbb{R}^{D}\rightarrow\mathbb{R}^{d}\) (with \(d<D\)), which compresses the input data into a smaller space representation (or latent space), and a decoder part \(Q:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}\), which reconstructs the encoded data back to the original input dimension. When the mapping kernel is linear, the autoencoder can be considered as a singular value decomposition analysis. An autoencoder that employs convolutional layers to perform the encoding and decoding operations on images, or more in general on 2D (or higher) data, is generally termed Convolutional AutoEncoder (CAE). During the convolution operation, the kernel swipes the data domain to extract features and learn spatial and/or temporal dependencies. The process is then carried out across multiple layers, obtaining representative features in a hierarchical sense [67]; CAEs with many hidden layers can be considered as deep neural networks (DNN). The initial layers learn a low-level representation or particular details about the data, while subsequent layers focus on larger and larger-scale information. In our work we developed four CAEs with different compression factors (CFs), defined as the ratio between the input and the latent dimensions, from 2 to 16; an example of a general CAE architecture used in our implementations is shown in Fig. 1. In Tab. 1 we report the main characteristics of the four CAEs implemented here; starting from CF2, we obtained higher CFs by adding more hidden layer to the initial architecture, but trying to keep the number of parameters as fixed as possible in order to make a more fair comparison of the four networks performance. Indeed, the number of weights and biases, as highlighted by the values in Tab. 1, varies only by \(\approx 10\%\). On the other hand, more hidden layers allow networks with higher CF to learn more features of the original fields during the compression step. The number of hidden layers is intended for the entire CAE architecture, comprising both the layers of the encoding \(n_{e}\) and decoding \(n_{d}\) part; they usually have an amount of layers which follow the relation \(n_{d}=n_{e}+1\), while the CAE is symmetric if the total number of layers is even.

\begin{table}
\begin{tabular}{c c c c c} \hline Id & Compression factor & No. Layers & No. params. & Latent space \\ \hline CF2 & 2 & 7 & 321,156 & \(64\times 64\times 128\) \\ CF4 & 4 & 8 & 339,892 & \(32\times 32\times 256\) \\ CF8 & 8 & 9 & 324,772 & \(32\times 32\times 128\) \\ CF16 & 16 & 11 & 349,804 & \(16\times 16\times 256\) \\ \hline \end{tabular}
\end{table}
Table 1: Description the four CAEs architecture: the compression factor (CF), the number of layers referring to both the encoding and decoding part of the network. The number of parameters account for all the weights and biases which have to be optimized during the training phase, and the latent space is the reduced space after the encoding network is applied.
## IV High-order statistics informed neural network: H-Sinn

The aim of this work is to devise a method that attains a compression of the data volume of DNS for stratified turbulence flows, yet preserving their statistics (up to the fourth-order moment, kurtosis). Mathematically, this goal can be cast as a multi-objective optimization (MOO) problem, and one of the greatest advantage of using a neural network to solve MOO problems is that, in principle, it is sufficient to include any term that one wish to minimize in the cost function and, as long as such terms are differentiable, the neural network will find a solution that approximately minimizes the objective (possibly a local minimum) through the back-propagation procedure. This is known as _scalarization method_, where the different terms are s simply added together, each one weighted by its own scalar factor [68]. Determining those weights is crucial, especially when employing ANNs, to make sure that the cost function is not dominated by one or a few terms only and that the solution actually minimizes each term. Intuitively one can think of the weights as normalization factors when different terms in the cost function have different order of magnitudes. Our implementation of H-SiNN involves the use of a statistics-aware loss function into a CAE, thus creating what we call statistics-informed convolutional autoencoder (SiCAE). As mentioned above, we focus on up-to-fourth order moments of the vertical velocity distribution, since we already know that it presents the highest variability in our simulations. In particular, our aim is to preserve standard deviation, skewness, and kurtosis of the vertical component of the velocity field \(w\), which we indicate with the symbols \(\sigma_{w}\), \(Sk_{w}\), and \(K_{w}\), respectively. As a result, the cost function of the CAE computed between the original data set \(y\) and its reconstruction \(\tilde{y}\) is defined as:

\[\begin{split}\mathcal{L}(y,\tilde{y})&=\frac{1}{D} \left(W_{1}\sum_{y}^{[u,v,w,\theta]}\sum_{i=1}^{N}\frac{(y_{i}-\tilde{y}_{i}) ^{2}}{N}+W_{2}\sum_{i=1}^{N}\frac{(\sigma_{w_{i}}-\sigma_{\tilde{y}_{i}})^{2} }{N}\right.\\ &+W_{3}\sum_{i=1}^{N}\frac{(Sk_{w_{i}}-Sk_{\tilde{w}_{i}})^{2}} {N}+W_{4}\sum_{i=1}^{N}\frac{(K_{w_{i}}-K_{\tilde{w}_{i}})^{2}}{N}\Bigg{)}\\ &=\frac{1}{D}\sum_{j=1}^{4}W_{j}R_{j}\end{split} \tag{6}\]

with \(D=3(W_{1}+W_{2}+W_{3}+W_{4})\), being the total sum of the weights. The second, third and forth term in (6) are the mean square errors (MSEs) between the moments of the original and the reconstructed vertical velocity \(W\) computed over horizontal planes (\(x\),\(y\)), namely the standard deviation \(\sigma_{w}\), the skewness \(Sk_{w}\) and kurtosis \(K_{w}\) respectively. We recall that for a normal distribution the reference values are: \(\sigma=1\), \(Sk=0\) and \(K=3\). In addition, here, we propose a method for dynamically weight the four terms in the cost function; each term is properly weighted with values \(W_{i}\) varying at each epoch during the training phase. In order to obtain a loss function with terms having the same order of magnitude during the whole training procedure at a given epoch \(m\), the weight of the \(i\)th term \(W_{i}\) is proportional to the sum of the other three terms at the previous epoch \(m-1\), e.g. \(W_{i}^{m}=\sum_{j\neq i}W_{j}^{(m-1)}R_{j}^{(m-1)}\) (and normalized such that their sum is equal to 1). Even though the various terms of the statistics-informed loss function are well balanced, the training procedure with additional constraints either coming from the knowledge of the physical system, as it happens for physics-informed neural networks, or from statistical information can be difficult, slower and sometimes very noisy. Finally, we have heuristically verified that the strategy which is computationally faster to train a H-SiNN is in two stages. In the first stage we treat the autoencoder as a classical CAE, with only the first term in Eq. 6. After we train for 50 epochs, we add the remaining terms in the cost function, i.e. starting from a network that achieves a low reconstruction error, in terms of mean error. During the second stage we train the CAEs for 100 more epochs with the addition of statistical constraints on the first four moments of the vertical velocity PDF; this is shown in panels (b)-(f) of Fig. 2. The behavior of the different terms of the loss function during one hundred epochs of training is shown in Fig. 3. We can appreciate from the four panels how there are essentially three weights (\(W_{MSE}\), \(W_{\sigma_{w}}\) and \(W_{Sk_{w}}\)) which are nearly equal for the entire training phase, whereas the coefficient weighting the kurtosis term is significantly smaller indicating a larger error on the forth-order moment, as expected.

### Plane-by-plane approach

We mentioned that the main objective is to obtain a tool which is able to reproduce most of the features of the velocity \(\mathbf{u}\) and temperature \(\theta\) fields of a fluid presenting high anisotropy (stratification) and non-homogeneity, due to the presence of large-scale intermittent structures. For these last

Figure 1: Schematic representation of a general Convolutional Autoencoder (CAE). The dimension of the _latent space_ and the number of encoding/decoding layers depend on the compression factor (CF). For each time step, the input is represented by horizontal slices (\(512^{2}\)) of the four physical variables produced by the simulations (\(u\), \(v\), \(w\) and \(\theta\)) stacked together forming a volume of \(512^{2}\times 4\). The number of filters and layers composing the encoder (and decoder) part also depend on the CF, all the CAE implemented in this work are symmetric.
reasons, we further implemented the standard CAE algorithm to make it more suitable and better performing in our particular application. A key ingredient of stratified flows is anisotropy, introduced by gravity, which tend to suppresses vertical motions with respect to the horizontal. Therefore, even though we deal with fully three-dimensional simulations the original cubic volume is split in horizontal planes along the vertical direction (\(z\)) before being used as input for the CAEs. This means that, for each time and height value \(z\), planes for every velocity component (\(u\), \(v\) and \(w\)) and for \(\theta\) are stacked together creating three-dimensional input data with dimensions \(512\times 512\times 4\), as indicated in Fig. 1. Before being divided, the data are normalized \(\hat{x}=(x-\mu)/\sigma\) using the average \(\mu\) and standard deviation \(\sigma\) computed on each simulation cube. In this approach, we believe that passing the information from the three components of the velocity, together with the temperature, is essential for the CAE to better reconstruct the flow dynamics at a given time \(t\) and altitude \(z\), since their variation is strictly correlated within the primitive equations. This approach has the shortcoming that it is not possible to add an important information we know about the velocity field, that is the incompressibility \(\mathbf{\nabla}\cdot\mathbf{u}=0\), since working plane by plane we cannot add constraints on the vertical derivative of the velocity \(\partial_{z}w\). This leads to values for the maximum and average divergence of \(\mathbf{u}\) reported in Tab. 2, where we can see that even if on average the condition is well satisfied in the reconstructed fields, locally this is not true while the condition is satisfied also point-wise for the original field.

## V Results

We analyze in details the output of a single simulation run for more than \(45\tau_{NL}\) corresponding to 100 time steps. Since in our approach we divide each simulation cube in planes along the vertical direction which are then treated as independent during training, the data set we use for training and testing our CAEs comprises \(N=51,200\) samples, each with dimension \(512^{2}\times 4\) (\(512^{2}\) grid-points for four variables \(u\), \(v\), \(w\) and \(\theta\)). We divide the data set in training and test set on a temporal basis, taking the first 70% snapshots for the training phase and the remaining 30% for testing the performance.

The strategy described above has been used to train four different statistics-preserving CAE, with compression factors (CF) of 2, 4, 8, and 16. Figure 2 shows the history of the different terms in the loss function, for increasing number of epochs. Red, green, blue, and cyan lines are respectively for CF 2, 4, 8, and 16. The first stage is made of 50 epochs, and in the second stage the networks are trained for an additional 100 epochs. The interesting feature, valid for all CFs, is that after the first stage (panel a) the MSE does not further decreases (panel c), while all other moments decrease by one order of magnitude or so (panels d-f). This shows that the second stage of training achieves the goal of minimizing the error of standard deviation, skewness, and kurtosis, without significantly degrade the point-wise accuracy of the reconstruction.

\begin{table}
\begin{tabular}{c|c|c c c c}  & Original & CF2 & CF4 & CF8 & CF16 \\ \hline \hline \(|(\mathbf{\nabla}\cdot\mathbf{u})|\) & \(|\times 10^{-11}|\) & 0.42 & 1.83 & 3.26 & 2.37 & 2.06 \\ \(max\left\{|\mathbf{\nabla}\cdot\mathbf{u}|\right\}\) & \(5.9\cdot 10^{-4}\) & 538 & 410 & 573 & 530 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Average and maximum value of the divergence of the velocity field after the reconstruction with the different CAE.

Figure 3: History of the weights multiplying the different terms of the loss functions: the mean squared error (MSE, black), the standard deviation (\(\sigma_{w}\), orange), the skewness (\(Sk_{w}\), blue) and the kurtosis (\(K_{w}\), gray).

Figure 2: Panel (a) shows the mean squared error (MSE) during the first 50 epochs of the training stage where no other terms are added to the loss function. In panel (b) the total loss during 100 epochs of the second training stage, in which the loss function is composed by four terms which are represented in panels (c)–(f) for the various CAE with different compression factor (CF).
### Global statistical properties of reconstructed fields

The trend of the statistical moments error seems to be inversely proportional to the order of moment, meaning that the improvement on the reconstructed standard deviation (2nd-order moment) is significantly larger than what is achieved for the kurtosis (4th-order moment) within the same amount of epochs. This is also pointed out by Fig. 3 where, for the four different SiCAE, the trend of the various weights balancing the loss function terms is reported as a function of the epoch; let us recall that each weight \(W_{i}\) is inversely proportional to the sum of the three contribution, given from the other moments, taken at the previous epoch \(~{}1/\sum_{(i\neq j)}W_{j}\). In addition, by looking at the behavior of the loss weights (Fig. 3), we observe that the error on the kurtosis is approximately between one and two orders of magnitude greater than the others. Such a large difference may be due to the high variability of \(K_{w}\) in this simulation, whereas standard deviation and skewness are expected to vary less and to be closer to reference values of a Gaussian distribution. Obviously a trend with the compression factor is observed, and, as it is expected, the higher CF the larger the reconstruction error both on the mean field and the statistical moments; this is not related to the newly introduced loss function since already during the first training stage (panel (a)) this trend can be clearly observed. The details of the results obtained on the test set for the various terms involved in the loss function during the two training stages are summarized in Tab. 3. We can see that including other terms into the loss function results in an average reconstruction error 2-3 times greater than what obtained after the first training stage, even though it remains on the order of MSE\(\sim 10^{-4}\), consistent with other results found in literature [69, 70]. In addition, when constraints on the statistical moments are added to the model, we observe that the reconstruction of high-order moments improves up to 10 times for any compression factor.

We selected an interval of the numerical simulation where isolated extreme bursts occur, producing evident peaks of the kurtosis \(K_{w}\) (see Fig. 4, panel (a)). In Fig. 4 the vertical velocity kurtosis computed on the original data (black dashed curve) is shown as a function of the turnover time for the entire test set, and compared with the kurtosis obtained from the CAEs reconstructions (different colors for different CF

Figure 4: Behavior of the kurtosis of the vertical velocity \(K_{w}\) over time for the test set. Panels (a) and (c) show \(K_{w}\) obtained for the different compression factors (CFs) in different colors (top) and the mean absolute relative error on the reconstructed kurtosis (bottom) during the first training stage when the loss function is composed by the MSE only. The same results are shown in panels (b) and (d) for the second training stage where the modified loss is adopted. Gray circles in panel (a) represent three times analyzed in the following figures with very different values of kurtosis \(K_{w}\) of the vertical velocity.

Figure 5: For each CAE, the behavior over time of the difference between the mean absolute percentage errors (MAPE) on the kurtosis \(K_{w}\) of the vertical velocity between the first \(\text{MAPE}^{(1)}_{K_{w}}\) and second training stage \(\text{MAPE}^{(2)}_{K_{w}}\). The MAPE is defined as \(\text{MAPE}^{(i)}_{K_{w}}=|(K^{\prime(i)}_{w}-K^{\prime\prime}_{w})/K^{\prime \prime}_{w}|\times 100\) with \(K^{\prime\prime}_{w}\) and \(K^{\prime\prime(i)}_{w}\) being the kurtosis of the original and reconstructed fields, respectively.

Figure 6: PDF of the vertical velocity \(P(w)\) computed after the first training stage (left panel) and at the end of the entire training stage (right panel), taken at \(t=47.9\). We can clearly notice that, despite the PDF core show a perfect agreement during the stage one yet, the fat tails of the distribution are consistently recovered only after the introduction of the custom loss function during the second training stage (right panel).

\begin{table}
\begin{tabular}{l c c c c c c}  & CAE & MSE & \(\sigma_{w}\) err. & \(S_{w}\) err. & \(K_{w}\) err. & Loss \\  & & \([\times 10^{-4}]\) & \([\times 10^{-5}]\) & \([\times 10^{-4}]\) & \([\times 10^{-2}]\) & \([\times 10^{-4}]\) \\ \hline \multicolumn{6}{l}{1st stage} \\ “standard loss” & & & & & \\ \hline \multirow{3}{*}{} & CF2 & 0.39 & 0.13 & 10.9 & 18.4 & 0.39 \\  & CF4 & 0.55 & 0.45 & 13.2 & 22.4 & 0.55 \\  & CF8 & 0.54 & 0.82 & 19.1 & 19.0 & 0.54 \\  & CF16 & 1.1 & 5.3 & 44.9 & 41.7 & 1.1 \\ \hline \multicolumn{6}{l}{2nd stage} \\ “statistical-informed loss” & & & & & \\ \hline \multicolumn{6}{l}{CF2} & 0.67 & 1.03 & 1.9 & 0.04 & 0.44 \\  & CF4 & 1.0 & 1.0 & 3.2 & 1.4 & 4.0 \\  & CF8 & 1.3 & 0.074 & 3.6 & 1.7 & 2.3 \\  & CF16 & 2.3 & 0.32 & 12.0 & 5.9 & 10.0 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Results of the four CAE on the test set for both the first (top) and second (bottom) training phase.
values). Panels (a) and (b) show a comparison obtained respectively after the first and second training phase; the relative error between the original \(K_{w}^{a}\) and recovered \(K_{w}^{r}\) kurtosis is shown in panels (c, first phase) and (d, second phase). This figure clearly highlights how the second stage of training (panels (b) and (d)) significantly improves the reconstruction of the fields up to the forth-order moment of the vertical velocity. All the CAEs show a significant reduction of the relative error, with a particular improvement in times where volume averaged \(K_{w}\) exhibits highly non-Gaussian values (\(t\approx 503\tau_{NL}\) and \(t\approx 520\tau_{NL}\)). In addition, with the application of the modified loss we are able to obtain percentage errors on the kurtosis smaller than \(10\%\) on average. However, the improvement in the reconstruction of the forth-order moment is obtained also for values close to the Gaussian reference (\(K_{w}\approx 3\)). Indeed, it is evident that our approach can be successfully extended to other situations where the PDFs of the dynamical fields are far from Gaussian-like shape, as it is the case for the concentration of a passive scalar [71] or for the velocity in the boundary layer [72; 3]. The statistical moments up to the forth-order embedded into the loss function guarantee an optimal reconstruction of the statistical properties of the physical fields without any prior knowledge about the shape or variability of their distribution function. The improvement we obtain with the custom loss is quantified in Fig. 5 where we show the difference between the mean absolute percentage error (MAPE) computed on \(K_{w}\) after the first and second training stage, as a function of time (once again, note that the training is performed by treating each plane and each time independently, i.e. no correlation along the \(z\) direction or in time is taken into account). The improvement reaches up to \(\approx 40\%\) at times with high values of \(K_{w}\) but in general an increasing of nearly \(10\%\) is observed at all times. In Fig. 6 we show the PDFs computed on the entire volume at time \(t\approx 520\tau_{NL}\) (third gray circle in Fig. 4), when \(K_{w}\simeq 6.8\), using both the original data (black) and data obtained from the four CAEs (colored curves). We notice how after the first training stage (left panel) the PDF core is already reliably recovered by the decoder, even though the tails significantly differ. In particular, the difference between original and reconstructed seems increasing for higher \(|w|\). The PDF in the right panel of Fig. 6, obtained after the second training stage, completely resembles the one computed from the original data, confirming that enforcing the statistical moments up to the forth-order is sufficient to enforce compatible PDFs. In

Figure 7: Left column: 3D renderings of the vertical velocity \(w/\sigma_{w}\), where the most extreme structures have been highlighted in solid colors. Right column: Vertical profile of the vertical velocity kurtosis \(K_{w}(z)\) of the original field (black dashed) and of the fields reconstructed with the SiCAE (colored, see legend). The vertical gray dash-dotted line (mid panel) represent the time and height taken for Fig. 9. Let us notice that the kurtosis is in logarithmic scale and the three panels show a different y-range for a clearer visualization.
Fig. 7 we show an example of large-scale intermittent structures observed at three different times, indicated by the gray circles in Fig. 4. The three-dimensional renderings represent the vertical velocity field \(w\) for low (top, \(t\simeq 506\tau_{NL}\)), medium (center, \(t\simeq 520\tau_{NL}\)) and high (bottom, \(t\simeq 502.5\tau_{NL}\)) values of the kurtosis \(K_{w}\), where values \(|w/\sigma_{w}|>4\) are highlighted in solid color while smaller values are depicted as a transparent blue. Aside each rendering the vertical profile of the kurtosis \(K_{w}\) computed plane-by-plane with the original data (black), as well as with the reconstructed (colored for the different CFs) is shown. Also for the top panel, when the kurtosis is small \(K_{w}\simeq 3.8\) there are planes reaching values of \(K_{w}\approx 6\), thus indicating a strong spatial variability even when the global kurtosis is close to the Gaussian reference value. These effects are enhanced when more vertical drafts develop within the flow (center and bottom) with kurtosis more than one order greater than the Gaussian reference values, as it shown in logarithmic scale in Fig. 7.

### CAE reconstruction of intermittent structures

We observed that the statistical properties, and in particular the kurtosis, of the vertical velocity are well recovered by the proposed implementation of statistics-informed CAEs. However, the presence of large-scale extreme events represents a challenge for their intermittency in time but also in space, creating patches of enhanced turbulence and regions populated by transient coherent structures. In addition, the low-dimensional representation of turbulent flows would represent a valuable result if the reconstruction produces truthful physical fields with local features which realistically resemble the dynamical quantities produced by the evolution of the governing equations. Therefore, in this section we analyze how the introduction of global statistical terms to the loss function locally affects the reconstruction of the entire velocity field \(\mathbf{u}\), as well as the potential temperature \(\theta\). In Fig. 8 we report a three-dimensional visualization of the vertical velocity \(w\) obtained with the four CAEs, as well as the original field for comparison. Data are taken at time \(t\simeq 520\tau_{NL}\), with relatively high kurtosis, that is \(K_{w}\approx 6.8\). Strong extreme vertical drafts, identified as values with \(|w/\sigma_{w}|>4\), have developed in the mid-bottom part of the domain. The same kind of structures can be noted in the reconstructed renderings, up to the maximum level of compression (i.e., CF16); in fact, at \(CF=16\) the shape and location of drafts starts being more noisy and deviating from the original field. The gray transparent shading represents velocity values below the threshold of \(4\sigma_{w}\), and therefore the majority of the volume. The reconstruction made by the implemented CAEs involves also the other components of the velocity field (\(u\) and \(v\)), as well as the (potential) temperature fluctuations \(\theta\). In order to have a general overview

Figure 8: Three-dimensional visualizations, taken at \(t/\tau\approx 520\), of the extreme events in the vertical velocity \(w\) from the original data (big cube), alongside the field reconstructed by the four SiCAE (small cubes). The vertical drafts (\(|w/\sigma_{w}|>4\)) are highlighted as solid colors (red/blue for positive/negative), whereas the other regions are represented as gray shaded areas
of how the neural networks recover all the physical quantities interested in the analyzed DNSs we represent them with several panels in Fig. 9. This figure shows an horizontal cut (\(x\),\(y\)) of the simulation domain taken at the same time of the previous figure (Fig. 8 at the height \(z^{*}\) indicated with a dash-dotted line in Fig. 7 (middle panel). The columns of Fig. 9 refer to the three components of the velocity field \(\mathbf{u}=(u,v,w)\) and the temperature \(\theta\). All the images are represented by the same color bar which is not shown since the main objective of this figure is the comparison between the first row, the original data, and the others, being the reconstructed physical fields for increasing values of the compression factor. Since we are looking at the domain from top (gravity is a vector entering the page), it is correct to have horizontal components strongly dominated by a nearly constant positive (red) for \(v\) and negative (blue) for \(u\) mean wind (see color bar); this is indeed the

Figure 9: Horizontal slices of the four physical fields, from left to right: horizontal component of the velocity \(u\) and \(v\) (left and mid-left), the vertical component \(w\) (mid-right) and the potential temperature \(\theta\) (right) taken at \(t/\tau\approx 520\) with a kurtosis \(K_{w}(z^{*})=10.0\). From the second line to bottom the reconstructed field obtained with different CAE having increasing compression factor (CF).
effect of Vertical Sheared Horizontal Winds (VSHWs) which are ubiquitous in stratified flows. Nevertheless, also the horizontal components of the velocity show small-scale perturbations where extreme vertical drafts develop, as already seen throughout this manuscript. The extreme events developed in this snapshot are clearly visible in the vertical velocity \(w\) and partially from the temperature renderings (third and forth columns in Fig. 9), and the same detail is captured also by the field reconstructions. Indeed, as already noticed for the three-dimensional visualization, the reconstruction is very reliable up to \(\mathrm{CF}=16\) where a significant checkerboard effect starts developing everywhere in the domain. This is probably due to the combined effect of the high compression factor and of the statistical-informed loss function presenting large-scale statistical constraints. In fact, by looking at the reconstructed fields after the first stage of training (not shown here), we can observe the same artifact at high compression but it is slightly reduced for the absence of other terms into the loss function.

## VI Conclusion

The work presented in this manuscript is an example of how the capability of machine learning tools to reproduce dynamical fields in fluid mechanics can be enhanced by incorporating physical and/or statistical knowledge about the system under study. In particular, we proposed a novel implementation of convolutional autoencoder (CAE) which include a modified loss function able to preserve the statistical moments of the vertical velocity field up to the forth order. We focused on the vertical component of the velocity field produced by the DNS of a stratified turbulent flow developing strong velocity drafts along the z direction, as is observed in the atmosphere and in oceans, which poses a challenge in finding an informative low-order representation of the physical fields. The comparison between a standard implementation of CAE and the SiCAE shows how the introduction of additional terms in the cost function, enforcing high-order moments of the PDF, do significantly improve the reconstruction of peculiar features of stratified geophysical flows, yet maintaining an acceptable level of the error in the reconstruction of the mean field. The various diagnostics we implemented all confirm that the novel implemented algorithm is able to reconstruct with good accuracy physical fields characterized by highly variable PDF (both in space and time). In the physical framework considered, the emergence of vertical drafts affect the local dynamics of stratified turbulent fluids, affecting their statistical description which not only can be non-stationary but is also non homogeneous. These results can potentially be generalized to other contexts and physical variables, to address for instance the reconstruction of intermittent passive scalar fields [73], in order to investigate their diffusion properties, or the velocity field obtained from observations of the planetary boundary layer [72; 3] and the upper atmosphere [74; 75]; strongly variables PDFs have also been detected in the oceans by means of GOMs [76] and using data from observations [6], further widening the application of the analysis proposed here, which can be very helpful in managing large volume of geophysical data from high-resolution numerical simulations and state-of-art observational campaigns. The addition of other terms to the loss function makes it necessary to have longer training phase and larger data set in order to converge to a well-trained model. Indeed, very often a compromise needs to be made between having detailed reconstructions of the mean fields dynamics and the reproducing peculiar structures of the turbulent flows, which is precisely where the algorithm presented improve the classical methodological approach based on CAE.

###### Acknowledgements.

 R.M. and R.F. acknowledge support from the project "EVENTFUL" (ANR-20-CE30-0011), funded by the French "Agence Nationale de la Recherche" - ANR through the program AAPG-2020. E.C. was partially supported by NASA grants 80NSSC20K1580 "Ensemble Learning for Accurate and Reliable Uncertainty Quantification" and 80NSSC20K1275 "Global Evolution and Local Dynamics of the Kinetic Solar Wind".

## Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.
