Espen Gaarder
email: espenhaug@mac.com

Gianfranco Spavieri
email: gspavieri@gmail.com

The Haug-Spavieri metric introduces in the Schwarzschild metric an observable term verifiable even in weak gravitational fields

Keywords: Mass-charge metric, time dilation, red-shift, light bending, black holes

Haug and Spavieri have recently presented a new exact solution to Einstein's field equation. In this paper, we will list multiple implications of this model concerning predicted: gravitational time dilation, redshift, light bending, escape velocity, and more. All the new results contain a relativistic correction term added to the Schwarzschild metric. The Schwarzschild metric can almost be seen as a weak gravitational field approximation of our mass-charge metric. The di↵erence in predictions of the two metrics, can be verified even in the weak gravitational field of the Sun and also in the gravitational field of the Earth by utlizing atomic or optical clocks on board of satellites. Yet, in the gravitational field of the Earth the e↵ect is so small, that to be detected, we must have access to the most advanced state of the art in modern optical clocks. However in the Sun's gravitational field it seems that even current cesium atomic clocks are good enough to test the di↵erences in predictions between the two metrics.

Also, photons, sent out in a strong gravitational field close to the event horizon of black holes, can likely be used to verify what metric derived from general relativity theory is the best to explain observed phenomena in black holes, such as the lack of velocity time dilation in high-z quasars.

1 Predictions from the Haug-Spavieri metric versus the Schwarzschild metric Einstein [START_REF] Einstein | Die grundlage der allgemeinen relativitätstheorie[END_REF] field equation is given by

R µ⌫ Rg µ⌫ = 8⇡G c 4 T µ⌫ (1)
The field equation must be solved based on certain boundary conditions to yield equations capable of predicting gravitational e↵ects, which can then be verified through observations. The most well-known and thoroughly tested solution is, by far, the Schwarzschild metric. It has undergone extensive testing in weak gravitational fields, confirming its accuracy in predicting gravitational time dilation, redshift, light bending, and more.

Recently, Haug and Spavieri [START_REF] Haug | Mass-charge metric in curved spacetime[END_REF] have derived a new exact solution to Einstein's field equation. This solution takes into account a stress-energy tensor that considers both the electric field energy of the charge and the gravitational field energy of the mass. To determine if this metric better describes the real world compared to the Schwarzschild metric, one must calculate predictable di↵erences.

The new metric derived by Haug and Spavieri is given by in S.I. units:

ds 2 = ✓ 1 2GM c 2 r + G 2 M 2 c 4 r 2 + Q 2 r 2 + P 2 r 2 ◆ c 2 dt 2 + ✓ 1 2GM c 2 r + G 2 M 2 c 4 r 2 + Q 2 c 2 r 2 + P 2 c 2 r 2 ◆ 1 dr 2 +r 2 ⌦ 2
(2) We will focus here on the case when the charge Q and magnetic moment P are set to 0 (Q = P = 0). Then, the metric simplifies to

ds 2 = ✓ 1 2GM c 2 r + G 2 M 2 c 4 r 2 ◆ c 2 dt 2 + ✓ 1 2GM c 2 r + G 2 M 2 c 4 r 2 ◆ 1 dr 2 + r 2 ⌦ 2 (3)
From this, we can easily derive predictions for gravitational time dilation, redshift, and the gravitational bending of light. The results, when compared to the Schwarzschild metric, are presented in Table 1.

Table 1:

The table displays well-known predictions from the Schwarzschild metric and also from the metric proposed by Haug and Spavieri. In all the predictions, there is a second-order term in the Haug and Spavieri metric that is absent in the Schwarzschild metric. This additional term primarily leads to predictable di↵erences in strong gravitational fields, especially for observations from photons emitted close to a black hole horizon. Alternatively, extremely precise measurement devices are required to detect the predicted di↵erence in a weak gravitational field, such as experiments conducted in the gravitational field of the Earth or the Sun.

Prediction

Schwarzschild:

Haug-Spavieri:

Gravitational red shift z = q 1 2GM R 1 c 2 q 1 2GM R 2 c 2 1 z = r 1 2GM r 1 c 2 + G 2 M 2 r 2 1 c 4 r 1 2GM r 2 c 2 + G 2 M 2 r 2 2 c 4 1 Time dilation T 1 = T 2 q 1 2GM r 2 c 2 q 1 2GM r 1 c 2 T 1 = T 2 r 1 2GM r 2 c 2 + G 2 M 2 r 2 2 c 4 r 1 2GM r 1 c 2 + G 2 M 2 r 2 1 c 4 Gravitational deflection ✓ = 4GM c 2 r ✓ = 4GM c 2 r 2G 2 M 2 c 4 r 2 Escape velocity v e = q 2GM r v e = q 2GM r G 2 M 2 r 2 c 2 Event horizon black-hole R s = 2GM c 2 R h = GM c 2 Mass from horizon radius M = c 2 Rs 2G M = c 2 R h G 2
The Haug-Spavieri metric is likely verifiable relative to the Schwarzschild metric even in weak gravitational field

Theoretically, there is a clear distinction in the formulas between the Schwarzschild and Haug-Spavieri metric. However, it is important to inquire about the magnitude of the pre-dicted di↵erences. Since it is a second-order term that creates the di↵erence in predictions between the two metrics, it is evidently a relativistic e↵ect that will predominantly have significant implications in strong gravitational fields. Nevertheless, it will also exhibit some minor observable e↵ects in weak gravitational fields, the magnitude of which we can easily calculate. In Earth's gravitational field, if we examine the gravitational time dilation between a clock placed at sea level and another one at a 1000-meter tower (or mountain top), we would only expect to observe approximately 7.6 ⇥ 10 23 seconds of di↵erence in time dilation over a second. This is considerably outside the precision of the current best atomic and optical clocks, which have a precision of approximately 10 18 per second. See, for example, [START_REF] Hinkley | An atomic clock with 10 18 instability[END_REF][START_REF] Bloom | An optical lattice clock with accuracy and stability at the 10 18 level[END_REF][START_REF] Nicholson | Systematic evaluation of an atomic clock at 2 ⇥ 10 18 total uncertainty[END_REF][START_REF]Frequency ratio measurements at 18-digit accuracy using an optical clock network[END_REF]. To increase precision, one can increase the distance between the clocks, for example, by synchronizing an atomic clock on the surface of the Earth with an atomic clocks placed on board of a geostationary satellite. Geostationary satellites are about 35,786 km out in space, so this is the distance between the clocks. Then the predicted di↵erence is about 2.367 ⇥ 10 19 seconds per second. Over a month (30 days), it amounts to 6.13 ⇥ 10 13 seconds. This could be possible with clocks better than the o↵-the-shelf atomic clocks, as it is well within the limit of precision in today's best atomic and optical clocks (10 18 ). The predicted di↵erence arises from the di↵erence between the gravitational time-dilation foreseen by the Schwarzschild metric and that of the Haug-Spavieri metric. This predicted di↵erence is after properly adjusting for velocity time dilation as well as the Sagnac e↵ect [START_REF] Sagnac | The demonstration of the luminiferous aether by an interfermoeter in uniform rotation[END_REF]. Table 2 gives examples of predicted di↵erences for gravitational time dilation from the Schwarzschild metric and the Haug and Spavieri metric in clocks placed distance d apart from each other. In the Sun's gravitational field, at a distance of 6.9 million kilometers from the center of the sun-this is the distance at which observational solar satellites orbit the sun-we find a predicted time-dilation di↵erence of about 3.1 ⇥ 10 14 seconds per second when using two satellites with clocks 1 km apart from each other (one further out from the sun). If the satellites are placed 100 km apart, one at 6.9 million kilometers from the sun and one 100 km further out, then the predicted time dilation di↵erence between the Schwarzschild metric and the Haug-Spavieri metric is about 3.1 ⇥ 10 12 seconds per second.

The predicted di↵erences are well within the measurement capabilities of the best atomic and optical clocks with precision of approximately 10 18 s. If the satellites are placed at 100,000 km apart, then we only need clocks with nanosecond accuracy, 10 9 , which is achiev-Table 3: The table shows predicted time dilation di↵erences between clocks placed a distance d apart on board solar satellites. The inner satellite is assumed to be placed at a distance of 6.9 million kilometers from the sun, and the next one is located a distance d further out from that. The clocks on the two satellites are Einstein synchronized and then compared after 1 second or after 1 day or after 30 days, the longer one wait the bigger the predicted di↵erence between the two metrics Then, we go down to the need of a reduced accuracy of 10 2 seconds per month. There are multiple methods to reduce errors in clocks. For example one would typically have at least three or even whole clusters of initially synchronized clocks on each satellite for control and error reduction. This should be no big task for NASA and could be one of the most important breakthroughs in gravity for decades if the Haug-Spavieri metric is verified to be a more accurate model of reality.

3 The Haug-Spavieri metric is likely verifiable relative to the Schwarzschild metric in strong gravitational fields When close to a black hole horizon, the term G 2 M 2 r 2 c 4 that appears in most formulas as an additional relativistic term relative to the Schwarzschild predictions, suddenly becomes very significant. Another important way to test the Haug-Spavieri metric against the Schwarzschild metrics is therefore from photons sent out close to the horizon of black holes.

There are multiple interesting types of black holes that could be worth investigating. For example, in high-z quasars (which all likely have black holes at their center), surprisingly, velocity time dilation [START_REF] Hawkins | On time dilation in quasar light curves[END_REF][START_REF] Hawkins | Time dilation and quasar variability[END_REF] has not been observed. Based on the results foreseen by our metric, photons can be sent out closer to the center of the mass of the quasar black holes than what possible in the Schwarzschild metric. The reason is that in our metric, the event horizon is only half of that of the Schwarzschild horizon, and photons can then be sent out closer to the center of the black hole without being inside the horizon. So, this is one possibility that can be more carefully investigated by researchers.

It has also been pointed out by many researchers [START_REF] Pathria | The universe as a black hole[END_REF][START_REF] Stuckey | The observable universe inside a black hole[END_REF][START_REF] Easson | Universe generation from black hole interiors[END_REF][START_REF] Christillin | The machian origin of linear inertial forces from our gravitationally radiating black hole universe[END_REF][START_REF] Pop | The universe in a black hole in einstein-cartan gravity[END_REF][START_REF] Akhavan | The universe creation by electron quantum black holes[END_REF][START_REF] Lineweaver | All objects and some questions[END_REF] that the Hubble sphere has many mathematical properties similar to a black hole, but the Schwarzschild black hole does not seem to be able to fit observations that well. One possibility is therefore to investigate if the Hubble sphere resembles more closely to a Haug-Spavieri black hole; we have outlined this possibility in a tentative theory for it in [START_REF] Haug | New exact solution to einsteins field equation gives a new cosmological model[END_REF], which should also be investigated further.

When it comes to gravitational light bending (deflection), our preliminary result (which we will have to carefully check later on) shows that the deflection is

✓ = 4GM c 2 r 2G 2 M 2 c 4 r 2 .
In a weak gravitational field, the extra term 2G 2 M 2 c 4 r 2 has too little impact to be tested, for example, from the bending of light around the sun. Here, the two models are currently indistinguishable. In a very strong gravitational field, it is interesting that the Schwarzschild metric predicts light bending at the event horizon itself (or just outside) as ✓ = 4GM c 2 rs = 2. A light bending of 2 basically means a photon could orbit the black hole. We get the identical result at the horizon from our new metric, namely

✓ = 4GM c 2 r h 2G 2 M 2 c 4 r 2 h = 2.
However, for a Planck mass black hole, in our theory, we can go all the way down to the Planck length radius, while the Schwarzschild metric can only go down to twice the Planck length. This could possibly indicate that our model is better suited to fit the Planck scale properties than the Schwarzschild metric; see also the table about micro black holes in [START_REF] Haug | Generalized rotation metric in curved spacetime[END_REF].

We hope that other researchers will closely investigate our metric and compare carefully its predictions with observations from black holes and the Schwarzschild metric. The comparison can, hopefully, be extended to other metrics such as the Reissner-Nordström [START_REF] Reissner | Über die eigengravitation des elektrischen feldes nach der einsteinschen theorie[END_REF][START_REF] Nordström | On the energy of the gravitation field in einstein's theory[END_REF], Kerr [START_REF] Kerr | Gravitational field of a spinning mass as an example of algebraically special metrics[END_REF] metric, the Kerr-Newman [START_REF] Newman | Note on the kerr spinning-particle metric[END_REF][START_REF] Newman | Metric of a rotating, charged mass[END_REF] metric, and even our new generalized rotation metric [START_REF] Haug | Generalized rotation metric in curved spacetime[END_REF].

Conclusions

We have shown that, with present technology, the metric derived by Haug-Spavieri appears to be verifiable relative to the Schwarzschild metric, even in the weak gravitational field of the Sun or the Earth. This can be achieved using just two Einstein-synchronized standard atomic clocks at normal orbital altitudes around the Sun. If one wishes to test the metrics in the Earth's gravitational field, it is necessary to maximize the distance and employ very modern equipment-likely requiring the use of state-of-the-art optical or atomic clocks, as achieving an observable minimum precession of 10 13 seconds per month is essential. Solar satellite experiments are possibly preferable, as such experiments meet accuracy requirements well within the reach of standard atomic clocks.

Additionally, observations from photons sent out close to the horizon of a black hole can likely help us distinguish which metric is best suited to describe reality. It seems that the Haug-Spavieri metric possibly interprets more e↵ectively the lack of observed time dilation in High z quasars. Furthermore, if the Hubble sphere is taken to represent a black hole, then photons sent out close to the Hubble horizon could possibly be better described through the Haug-Spavieri metric. Our tentative considerations about some cosmological and gravitational phenomena, related to the Haug-Spavieri metric, need to be scrutinized by multiple researchers over time before drawing any conclusions.

Table 2 :

 2 The table shows predicted gravitational time dilation di↵erences between initial Einstein synchronized clocks placed distance (altitude) d apart in the Earth gravitational field. The di↵erence represents the prediction of the Schwarzscihld metric versus the Haug-Spavieri metric.

	Altitude distance d	Predicted di↵erence	Predicted di↵erence	Predicted di↵erence
	between clocks		per second :	per day:	per month (30 days):
	one ground based
	Altitude of tower or mountain clock:
	1 km Altitude of airplane clock:	7.6 ⇥ 10 23 s	6 .57 ⇥ 10 18 s	1 .97 ⇥ 10 16 s
	10 km Altitude of air-baloon clock:	7.59 ⇥ 10 22 s	6 .56 ⇥ 10 18 s	1 .97 ⇥ 10 15 s
	30 km Altitude satelite clock:		2.27 ⇥ 10 21 s	1 .96 ⇥ 10 16 s	5 .88 ⇥ 10 15 s
	20,000 km Altitude geostationary satelite clock: 2.28 ⇥ 10 19 s	1 .97 ⇥ 10 14 s	5 .91 ⇥ 10 13 s
	35 786 km		2.37 ⇥ 10 19 s	2 .05 ⇥ 10 14 s	6 .14 ⇥ 10 13 s

Availability of data and materials

The study contains no new data.

Declarations Competing interests

The authors have no competing interests to report.


R Foldes

E Camporeale

R Marino

Low-dimensional representation of intermittent geophysical turbulence with H-SiNN Low-dimensional representation of intermittent geophysical turbulence with High-Order Statistics-informed Neural Networks

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.

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 socalled quasi-geostrophic balance [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF] . At mesoscale O(100 km) and even more at the sub-mesoscale O(10 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 [START_REF] Lilly | Stratified Turbulence and the Mesoscale Variability of the Atmosphere[END_REF] . This is characterized by the presence of shears leading to instabilities, the vertical extension of strata being controlled by the Brunt-Väisälä 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 [START_REF] Mahrt | Observations of turbulence in stratified flow[END_REF][START_REF] Mahrt | Intermittency of Atmospheric Turbulence[END_REF][START_REF] Lenschow | A comparison of higher-order vertical velocity moments in the convective boundary layer from lidar with in situ measurements and large-eddy simulation[END_REF] and the ocean [START_REF] D'asaro | High-Frequency Internal Waves on the Oregon Continental Shelf[END_REF][START_REF] Capet | Mesoscale to Submesoscale Transition in the California Current System. Part I: Flow Structure, Eddy Flux, and Observational Tests[END_REF] , and extensively in a variety of numerical investigations [START_REF] Rorai | Turbulence comes in bursts in stably stratified flows[END_REF][START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF][START_REF] Sujovolsky | Invariant manifolds in stratified turbulence[END_REF][START_REF] Feraco | Connecting large-scale velocity and temperature bursts with small-scale intermittency in stratified turbulence[END_REF][START_REF] Rosenberg | Correlation between buoyancy flux, dissipation and potential vorticity in rotating stratified turbulence[END_REF][START_REF] Marino | Turbulence generation by large-scale extreme vertical drafts and the modulation of local energy dissipation in stably stratified geophysical flows[END_REF] through the forth-order moment (i.e., kurtosis) of the vertical velocity. Marino et al. [START_REF] Marino | Turbulence generation by large-scale extreme vertical drafts and the modulation of local energy dissipation in stably stratified geophysical flows[END_REF] 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 smallscale intermittency [START_REF] Feraco | Connecting large-scale velocity and temperature bursts with small-scale intermittency in stratified turbulence[END_REF] . 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 [START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF] 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) [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF][START_REF] Girimaji | Partially-Averaged Navier-Stokes Model for Turbulence: A Reynolds-Averaged Navier-Stokes to Direct Numerical Simulation Bridging Method[END_REF][START_REF] Spalart | Detached-eddy simulation[END_REF] , to the simulation of geophysical turbulence using a reduced models [START_REF] Julien | Rapidly rotating turbulent rayleigh-bénard convection[END_REF] . 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 [START_REF] Pouquet | Geophysical turbulence and the duality of the energy flow across scales[END_REF][START_REF] Marino | Resolving the paradox of oceanic large-scale balance and small-scale mixing[END_REF][START_REF] Rosenberg | Evidence for Bolgiano-Obukhov scaling in rotating stratified turbulence using high-resolution direct numerical simulations[END_REF] . 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 [START_REF] Mahrt | Observations of turbulence in stratified flow[END_REF] 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 runtime analysis, as well as for the implementation of machine learning techniques able to learn features of the dynamical systems from the latent space [START_REF] Parente | Identification of lowdimensional manifolds in turbulent flames[END_REF][START_REF] Wan | Data-assisted reduced-order modeling of extreme events in complex dynamical systems[END_REF][START_REF] De Jesús | Data-driven low-dimensional dynamic model of kolmogorov flow[END_REF] . 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 datamining and machine learning techniques. The study of low dimensional manifolds has a long history in dynamical systems and in turbulence [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. I. Coherent structures[END_REF][START_REF] Zdybał | Cost function for low-dimensional manifold topology assessment[END_REF] , laying at the foundation of so-called reduced order methods (ROMs) [START_REF] Lassila | Model order reduction in fluid dynamics: Challenges and perspectives[END_REF] . In turn, ROMs have also recently benefited by progresses in machine learning based approaches [START_REF] Mohan | A deep learning based approach to reduced order modeling for turbulent flow control using lstm neural networks[END_REF][START_REF] Chen | Physics-informed machine learning for reduced-order modeling of nonlinear problems[END_REF] . Dimensionality reduction techniques such as proper orthogonal decomposition [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF][START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF] and dynamic mode decomposition [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Kutz | Dynamic Mode Decomposition[END_REF] have been extensively used to address is-sues concerning a variety of flows [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF][START_REF] Garicano-Mena | A composite dynamic mode decomposition analysis of turbulent channel flows[END_REF][START_REF] Williams | Hybrid reduced-order integration with proper orthogonal decomposition and dynamic mode decomposition[END_REF][START_REF] Marensi | Symmetry-reduced dynamic mode decomposition of near-wall turbulence[END_REF][START_REF] Cheng | Large eddy simulation and dynamic mode decomposition of turbulent mixing layers[END_REF] . Other examples of model order reduction techniques applied to fluid frameworks are Galerkin-projection based nonlinear methods [START_REF] Carlberg | The gnat method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows[END_REF] and system identification based auto-regressive models [START_REF] Raveh | Identification of computational-fluid-dynamics based unsteady aerodynamic models for aeroelastic analysis[END_REF][START_REF] Sarma | Uncertainty reduction in aeroelastic systems with time-domain reduced-order models[END_REF] . Techniques based on machine learning principles that have been employed in modeling the equations of motion are Gaussian process regression [START_REF] Raissi | Deep multi-fidelity gaussian processes[END_REF][START_REF] Zhang | Machine learning methods for data-driven turbulence modeling[END_REF] , symbolic regression [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Schmelzer | Discovery of algebraic reynolds-stress models using sparse symbolic regression[END_REF] , field inversion [START_REF] Duraisamy | New approaches in turbulence and transition modeling using data-driven techniques[END_REF][START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF] , and artificial neural networks (ANNs). In datadriven 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 [START_REF] Angriman | Generation of Turbulent States using Physics-Informed Neural Networks[END_REF][START_REF] Angriman | Assimilation of statistical data into turbulent flows using physics-informed neural networks[END_REF] . Convolutional neural networks (CNNs) have been successfully applied to the identification of flow structures [START_REF] Morimoto | Convolutional neural networks for fluid flow analysis: toward effective metamodeling and low dimensionalization[END_REF] and to perform nonlinear modal decomposition in turbulent flows [START_REF] Fukami | Convolutional neural network based hierarchical autoencoder for nonlinear mode decomposition of fluid field data[END_REF][START_REF] Murata | Nonlinear mode decomposition with convolutional neural networks for fluid dynamics[END_REF] ; more recently, CNN-based architectures have been extended to fully 3D DNS of turbulent channel flow [START_REF] Nakamura | Convolutional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow[END_REF] and to numerical simulations of flows characterized by the emergence of extreme events, in terms of boundary coherent structures [START_REF] Khoa Doan | Convolutional autoencoder for the spatiotemporal latent representation of turbulence[END_REF] . 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 [START_REF] Salehipour | Deep learning of mixing by two 'atoms' of stratified turbulence[END_REF] 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 [START_REF] Wang | An efficient and effective convolutional auto-encoder extreme learning machine network for 3d feature learning[END_REF] and their use has recently been proposed in the context of fluiddynamics [START_REF] Gonzalez | Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems[END_REF][START_REF] Xu | Multi-level Convolutional Autoencoder Networks for Parametric Prediction of Spatio-temporal Dynamics[END_REF][START_REF] King | From deep to physicsinformed learning of turbulence: Diagnostics[END_REF] . 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 [START_REF] Kutz | Deep learning in fluid dynamics[END_REF] . 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 statis-tical 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 smallscale 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 highorder statistics-informed neural network (H-SiNN). The concept of statistics-informed neural network (SiNN) has been first introduced by [START_REF] Zhu | Learning stochastic dynamics with statistics-informed neural network[END_REF] , 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 largescale intermittent events in the vertical component of the velocity and in the temperature field [START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF] .

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 u = (u ⊥ , w) remains incompressible, ∇ • u = 0, while small density variations are taken into account only in the buoyancy term. Such a model can be written as:

∂ t u + (u • ∇)u = -∇p -Nθ ẑ + ν∇ 2 u + F (1)
∂ t θ + (u • ∇)θ = Nw + κ∇ 2 θ . (2)
with θ being the temperature fluctuations evaluated relative to a mean temperature profile θ 0 , N = [-g∂ z θ /θ 0 ] 1/2 is the Brunt-Väisälä 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π/L 0 and the resolution is 512 3 grid-points. A random forcing F is imposed to the momentum equation at k F = 2π/L f ∈ [2, 3], continuously injecting kinetic energy into the system and allowing to reach a turbulent stationary state. In the above equations, ν and κ are the kinematic viscosity and the thermal diffusivity, respectively and we take the Prandtl number ν/κ = 1, with ν = 10 -3 ; finally, p is the pressure. We adopt the following definitions for the dimensionless Reynolds and Froude numbers,

Re = U rms L int /ν, Fr = U rms /[L int N] , (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 = ReFr 2 , measure the relative strength of buoyancy to dissipation and are commonly used to distinguish between wave-dominated and turbulence-dominated regimes [START_REF] Ivey | Density stratification, turbulence, but how much mixing?[END_REF] . 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 [START_REF] Mininni | A hybrid mpi-openmp scheme for scalable parallel pseudospectral computations for fluid turbulence[END_REF][START_REF] Rosenberg | Gpu parallelization of a hybrid pseudospectral geophysical turbulence framework using cuda[END_REF] . It allows for a variety of physical solvers and supports non-cubic geometry [START_REF] Sujovolsky | Generation of turbulence through frontogenesis in sheared stratified flows[END_REF] , and non-periodic boundary conditions [START_REF] Fontana | Fourier continuation method for incompressible fluids with boundaries[END_REF] . 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 = (w -w) 3 ⟨(w -w) 2 ⟩ 3/2 , (4)
K w = (w -w) 4 ⟨(w -w) 2 ⟩ 2 , (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 θ ) in numerical simulations [START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF] . As observed by Feraco et al. [START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF] , K w presents strong non-Gaussian values in a narrow range of the Froude number, with an evident peak at Fr ≈ 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. [START_REF] Feraco | Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number[END_REF] , and subsequent works [START_REF] Feraco | Connecting large-scale velocity and temperature bursts with small-scale intermittency in stratified turbulence[END_REF][START_REF] Marino | Turbulence generation by large-scale extreme vertical drafts and the modulation of local energy dissipation in stably stratified geophysical flows[END_REF] , presenting the highest level of intermittency (see Fig. 2 I. 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.

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 [START_REF] Demers | Non-linear dimensionality reduction[END_REF] . 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 : R D -→ R d (with d < D), which compresses the input data into a smaller space representation (or latent space), and a decoder part Q : R d -→ 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 [START_REF] Masci | Stacked convolutional auto-encoders for hierarchical feature extraction[END_REF] ; 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. I 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. I, varies only by ≈ 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.

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 fourthorder 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 ares simply added together, each one weighted by its own scalar factor [START_REF] Gunantara | A review of multi-objective optimization: Methods and its applications[END_REF] . 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 σ 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 ỹ is defined as: (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 σ 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.

L (y, ỹ) = 1 D W 1 [u,v,w,θ ] ∑ y N ∑ i=1 (y i -ỹi ) 2 N +W 2 N ∑ i=1 (σ w i -σ wi ) 2 N +W 3 N ∑ i=1 (Sk w i -Sk wi ) 2 N +W 4 N ∑ i=1 (K w i -K wi ) 2 N = 1 D 4 ∑ j=1 W j R j (

A. 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 u and temperature θ fields of a fluid presenting high anisotropy (stratification) and non-homogeneity, due to the presence of large-scale intermittent structures. For these last 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 θ are stacked together creating three-dimensional input data with dimensions 512 × 512 × 4, as indicated in Fig. 1. Before being divided, the data are normalized x = (xµ)/σ using the average µ and standard deviation σ 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 ∇ • u = 0, since working plane by plane we cannot add constraints on the vertical derivative of the velocity ∂ z w. This leads to values for the maximum and average divergence of u reported in Tab. II, 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τ 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 × 4 (512 2 grid-points for four variables u, v, w and θ ). 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.

A. 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 (2ndorder 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/ ∑ (i̸ = 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. III. 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∼ 10 -4 , consistent with other results found in literature [START_REF] Glaws | Deep learning for in situ data compression of large turbulent flow simulations[END_REF][START_REF] Momenifar | Dimension reduced turbulent flow data from deep vector quantisers[END_REF] . 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

K w = |(K r,(i) w -K o w )/K o w | × 100 with K o w and K r, (i)
w being the kurtosis of the original and reconstructed fields, respectively. (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 ≈ 503τ NL and t ≈ 520τ 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 ≈ 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 [START_REF] Pumir | Exponential tails and random advection[END_REF] or for the velocity in the boundary layer [START_REF] Mahrt | Observations of turbulence in stratified flow[END_REF][START_REF] Brun | Large-eddy simulation of a katabatic jet along a convexly curved slope. part i: Statistical results[END_REF] . 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 ≈ 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 ≈ 520τ NL (third gray circle in Fig. 4), when K w ≃ 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

B. 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 repre-sents 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 u, as well as the potential temperature θ . 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 ≃ 520τ NL , with relatively high kurtosis, that is K w ≈ 6. 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 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 [START_REF] Orsi | Scalar mixing in homogeneous isotropic turbulence: A numerical study[END_REF] , in order to investigate their diffusion properties, or the velocity field obtained from observations of the planetary boundary layer [START_REF] Mahrt | Observations of turbulence in stratified flow[END_REF][START_REF] Brun | Large-eddy simulation of a katabatic jet along a convexly curved slope. part i: Statistical results[END_REF] and the upper atmosphere [START_REF] Vierinen | Observing mesospheric turbulence with specular meteor radars: A novel method for estimating second-order statistics of wind velocity[END_REF][START_REF] Chau | Radar observation of extreme vertical drafts in the polar summer mesosphere[END_REF] ; strongly variables PDFs have also been detected in the oceans by means of GOMs 76 and using data from observations [START_REF] D'asaro | High-Frequency Internal Waves on the Oregon Continental Shelf[END_REF] , further widening the application of the analysis proposed here, which can be very helpful in managing large volume of geo-physical 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.

Figure 1 .

 1 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 θ ) stacked together forming a volume of 512 2 × 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.

  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 σ w , the skewness Sk w and kurtosis K w respectively. We recall that for a normal distribution the reference values are: σ = 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 ith term W i is proportional to the sum of the other three terms at the previous epoch m -1, e.g. W m i = ∑ j̸ =i W

4 4 × 10 - 5 6 × 10 - 5 2 × 10 - 4 3 × 10 - 4 MSE

 4105105104104

Figure 2 .

 2 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).

Figure 3 .

 3 Figure3. History of weights multiplying the different terms of the loss functions: the mean squared error (MSE, black), the standard deviation (σ w , orange), the skewness (Sk w , blue) and the kurtosis (K w , gray).

Figure 4 .Figure 5 . 1 )

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

  Stage 2: moments-informed loss original: Kw = 6.8 CF 2 : Kw = 6.6 CF 4 : Kw = 6.2 CF 8 : Kw = 6.1 CF 16 : Kw = 5.9

Figure 6 .

 6 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).

Figure 7 .

 7 Figure 7. Left column: 3D renderings of the vertical velocity w/σ 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.

Figure 8 .

 8 Figure 8. Three-dimensional visualizations, taken at t/τ ≈ 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/σ w | > 4) are highlighted as solid colors (red/blue for positive/negative), whereas the other regions are represented as gray shaded areas

8 .Figure 9 .

 89 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 θ (right) taken at t/τ ≈ 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).

Table II .

 II Average and maximum value of the divergence of the velocity field after the reconstruction with the different CAE. CAE MSE σ w err. Sk w err. K w err. Loss [×10 -4 ] [×10 -5 ] [×10 -4 ] [×10 -2 ] [×10 -4 ]

	1st stage
	"standard loss"
	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
	2nd stage
	"statistical-informed loss"
	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
							Original CF2 CF4 CF8 CF16
						|⟨∇ • u⟩| [×10 -11 ] max {|∇ • u|}	0.42 5.9 • 10 -4 538 410 573 530 1.83 3.26 2.37 2.06

Table III .

 III Results of the four CAE on the test set for both the first (top) and second (bottom) training phase.

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