Classifying Turbulent Environments via Machine Learning

The problem of classifying turbulent environments from partial observation is key for some theoretical and applied fields, from engineering to earth observation and astrophysics, e.g. to precondition searching of optimal control policies in different turbulent backgrounds, to predict the probability of rare events and/or to infer physical parameters labelling different turbulent set-ups. To achieve such goal one can use different tools depending on the system's knowledge and on the quality and quantity of the accessible data. In this context, we assume to work in a model-free setup completely blind to all dynamical laws, but with a large quantity of (good quality) data for training. As a prototype of complex flows with different attractors, and different multi-scale statistical properties we selected 10 turbulent 'ensembles' by changing the rotation frequency of the frame of reference of the 3d domain and we suppose to have access to a set of partial observations limited to the instantaneous kinetic energy distribution in a 2d plane, as it is often the case in geophysics and astrophysics. We compare results obtained by a Machine Learning (ML) approach consisting of a state-of-the-art Deep Convolutional Neural Network (DCNN) against Bayesian inference which exploits the information on velocity and enstrophy moments. First, we discuss the supremacy of the ML approach, presenting also results at changing the number of training data and of the hyper-parameters. Second, we present an ablation study on the input data aimed to perform a ranking on the importance of the flow features used by the DCNN, helping to identify the main physical contents used by the classifier. Finally, we discuss the main limitations of such data-driven methods and potential interesting applications.


Introduction
Extracting information about the statistical context from partial measurements of a turbulent system is a fundamental problem in fluid dynamics. Unless the data available are large and of very high quality, any linear regression trying to statistically distinguish different turbulent set up from limited data is destined to fail. This is evident if you need, e.g., to estimate the turbulent intensity, given by a combination of the velocity root mean square, viscosity and flow correlation length, out of sparse and temporal measurements, and without any prior on the dissipative properties [1]. Similarly, many turbulent set-ups are subjected to transitions at changing control parameters or boundary conditions, developing macroscopically similar behaviour but different subtle multi-scale statistical properties and/or transient phases [2,3,4,5,6]. As a result, classifying turbulent scenarios from spot measurements is a key problem to reduce the complexity of predicting the flow evolution, to select optimal policies, pre-trained in different set-ups, for optimal navigation [7,8,9,10,11,12,13] and/or estimating the probability of extreme events [3,14,15,16,17,18]. Similarly in the context of data-driven turbulence modelling of non-stationary and transient flows, an instantaneous classification of the environment would be beneficial to perform an ad-hoc finetuning of the closure model depending on the specific local turbulence condition, e.g. high/low shear, stable/unstable conditions or strong/weak updrafts, such as to cure the 'curse-of-dimensionality' and improve the training performances [19,20,21,22]. As a prototype of complex flows with different dynamical attractors, and different multiscale statistical properties we selected 10 turbulent 'ensembles' by changing the rotation frequency of the frame of reference of a 3d domain and we supposed to have access to a set of partial observations limited to the instantaneous kinetic energy distribution in a 2d plane, as it is often the case in geophysics and astrophysics [23,24,25,26,27,28,29,30]. Increasing rotation rate, turbulence undergoes global and local changes, concerning the development of quasi-2d columnar structures, extreme non-Gaussian small-scale anisotropic fluctuations, long-range wave-wave correlations [31,32,33,34], as a result, it can be considered a paradigmatic test-bed where spatial and Fourier features have different fingerprints in the statistical distributions. Classification methods can be split in two different classes. The equation-driven tools, where the data analysis is enhanced by the physical knowledge of the underlying dynamical laws [35,36,37], and the purely data-driven methods, based on the statistical analysis of the available data. Some attempts to classify turbulence using equation-driven tools have been developed in the recent years, using Fourier spatio-temporal decomposition [38] or Nudging [37,39,40]. and Ω = 10 in the right panel. In the first case the rotation frequency is small enough that the flow develops a dynamics similar to the standard 3D non-rotating case, with forward cascade of energy and small scales fluctuations. In the second case instead rotation is large enough to produce an inverse energy transfer that give rise to large-scales structures almost constant in the rotation direction.
However, all these techniques on top of the knowledge of the equation-of-motion require access to important numerical resources to run digital twins of the observed system. In this work we discuss the model-free approaches, where we are completely agnostic of the system dynamical laws and we aim to implement purely data-driven tools trained on a large quantity of good quality data. In this paper, we attack this general inverse problem following recent progresses in Machine Learning (ML) augmented fluid dynamics, [41,42,43,44,45,46,7,47,48,49,1], by training a state-of-the-art Deep Convolutional Neural Network (DCNN) in a supervised way, [50,51,52,53], thanks to a large dataset build up with the use of high performance Direct Numerical Simulations (DNS) and accessible at the webpage http://smart-turb.roma2.infn. it. We compare the non-linear ML regression against a Bayesian baseline, based on the estimation of velocity and vorticity flow moments, in order to have comparison with both large-and small-scale physics inputs. Furthermore, we critically analyze the performance of the DCNN at changing the amount of training data and the number of optimized hyper-parameters. Finally, we try to open-theblack-box by performing an ablation study on the input data aimed to rank the flow features used by the DCNN, and to help to identify the main physical contents used by the classifier. Such 'inverse-engineering' approach, identifying key degrees-of-freedom for the classification can also have important applications concerning active control to favour or contrast flow transitions in other general setups. The paper is organized as follow. In Sect. 2 we give details on the physics of turbulent flows and on the regression problem we are focusing on. In Sect. 3 we describe the datasets used. In Sect. 4 we describe the methodologies implemented, the DCNN and the BI. Results are presented in Sect. 5. In Sect. 6 we present the ablation study on the input data, and in Sect. 7 we draw our concluding remarks.

Physics Background
Before moving towards the details of the regression problem let us discuss the physical background of turbulence under rotations. Rotating turbulence is described by the Navier-Stokes equations (NSE) equipped with the Coriolis force, here ν is the kinematic viscosity, 2Ωẑ × u is the Coriolis force produced by rotation, and Ω is the frequency of the angular velocity around the rotation axisẑ. The fluid density is constant and absorbed into the definition of pressure p. f , represents an external forcing mechanism, see Sect. 3 for details on the type of forcing used here. The relative importance of non-linear (advection) with respect to linear terms (viscous dissipation) is given by the Reynolds number, Re = ε , where ε f is the rate of energy input at the forcing scale L f ∼ 1/k f . Rotation introduces a second non-dimensional control parameter, the Rossby number, which represents the ratio between rotation frequency and the flow characteristic time scale at the forcing wavenumber:    While it is well known that increasing Re leads to transitions to different turbulent regime, the effects of varying Ro are more subtle and the results of contrasting flow features, including wave-wave interactions, large-scale energy condensation in quasi-2d structures, small-scale anisotropic intermittency [16]. In the limit of large Ro or small Ω, standard 3d homogeneous and isotropic turbulent dynamics is observed, on the other hand, increasing Ω, Ro ≤ 1, the flow becomes almost bidimensional [54,55,56,57,58,34], moving with uniform velocity along the direction parallel to the rotation axis, i.e.ẑ. Simultaneously, at high rotation rates, the flow develops coherent/largescales vortex like structures that live on the horizontal plane, (x, y), as can be seen in the renderings of Fig. 1. These structures are fed by an inverse energy cascade that transports energy from the injection scales towards the largest scales available in the physical domain [31,59,60].
In Fig. 2 we present a visualization of our problem: classifying the turbulent environments from one realization of the velocity module on 2d planes extracted from 3d simu-lations at varying the reference rotation rate Ω. From the visual comparison of the 2d planes for different Ω, (from Ω = 2 to Ω = 30), and the results shown for the spectra and for the temporal evolution of the total kinetic energy (panels 'a' and 'b') one can immediately understand how it is not so difficult to classify the systems in two rough classes, strong/weak rotations. On the other hand, it looks much harder to distinguish different scenarios inside one of the two rough classes. The aim of the paper is to find a tool capable to disentangle also the subtle correlations that distinguish among all turbulent ensembles. We aim to solve the inference problem with the new data analyses paradigms proposed in computer vision such as a non-linear classifier based on DCNN [61,62,63,64,65,66,67,68]. We imagine to have only a partial knowledge of the system, and in particular, for the sake of a possible realistic implementation, we make the four following assumptions: -Time evolution is not accessible -We have access only to one scalar measurement, e.g. the velocity amplitude, |u| -We can observe the system on a 2d plane and not on the full 3d domain As a baseline, we have repeated the same investigation using a Bayesian Inference (BI) analysis, in order to compare the 'supremacy' of the black-box DCNN with a tool based on correlations with physically meaningful observables, i.e. moments of the velocity field and of its derivatives. Let us notice that we have also attempted to perform the same classification, using the unsupervised Gaussian Mixture Models (GMM) classification [69,70,71]. However we did not find any relation between the unsupervised classification and the rotation value, so it has not been included in the discussion with DCNN and BI.

Numerical Simulations
To generate the database of turbulent flows on a rotating frame, we have performed a set of Direct Numerical Simulations (DNS) of the NSE for incompressible fluid in a triply periodic domain of size L = 2π with a resolution of N = 256 gird-points per spatial direction. We used a fully dealiased parallel pseudospectral code, with time integration implemented as a second-order Adams-Bashforth scheme and the viscous term implicitly integrated. A linear friction term, 0.1∆ −1 u, is added at wavenumbers with module |k| ≤ 2 to reach a stationary state. At the same time, viscous dissipation, ν∇ 2 u, is replaced by a hyperviscous term, ν∇ 4 u, with ν = 1.6 × 10 −6 , to increase the inertial range of scales dominated by the non-linear energy transfer. The forcing mechanism, f , is a Gaussian process delta-correlated in time, with support in wavenumber space centered at k f = 4. The total energy time evolution for the 10 different simulations varying Ω is presented in Fig. 2(a), while in panel (b) of the same figure we show the energy spectra, E(k) = 1 2 k≤|k|<k+1 |û(k)| 2 , averaged over time in the stationary regime for the same set of 10 simulations. From the energy spectra we can observe that while for Ω = 2 and 4 (in the direct cascade regime) there is a depletion of energy at wavenumbers smaller than |k| < k f , for Ω values above the transition the split cascade regime, also the small wavenumbers are filled in with energy driven by the inverse cascade, namely from the forcing to the largest system's scales.

Dataset extraction
The ML training and validation datasets, accessible at the webpage http://smart-turb.roma2.infn.it, are extracted from the DNS described above, following [72], such as to match all assumption discussed in the introduction (Sect. 1), namely; -In order to vary the Rossby number, we performed 10 numerical simulations each with a different rotation frequency. Two simulations well below the transition in the direct cascade regime, Ω = 2, 4, three of them in the transition region, Ω = 5, 6, 7, and the re- The network starts with four convolutional layers that map the input plane to 32 different 8x8 images of aggregated features. The second part of the network is composed by two fully connected layers that give in output the rotation prediction, Ω pred . All layers are followed by a ReLu activation function, a dropout of 30% and 20% follows respectively the two last fully connected layers.
maining five well inside the split cascade regime with Ω = 8, 10, 15, 20, 30. -For each of the 10 simulations we have dumped a number of roughly 600 snapshots of the full 3d velocity field with a temporal separation large enough to decrease correlations between two successive data-points (see Fig. 2(a) for the total energy evolution). -For each configuration, 64 horizontal planes (x, y), (perpendicular to the rotation axis), are selected at different z and for each of them we computed the module of the velocity field as; |u| = u 2 x + u 2 y + u 2 z . -To increase the dataset, we shift each of the 64 planes by choosing randomly a new center of the plane and using periodic boundary conditions, such as to obtain a total number of 50k planes for each Ω value.
The full dataset, in this way, is composed of half a million different velocity module planes, out of which the 80% of them are used as training set while 20% are used as validation set. To be sure that validation data were never seen during training, none of the planes used in the validation set is obtained from a random shift of a plane contained in the training set. Instead, the two sets are built using the same protocol but from different temporal realization of the flow splitting its time evolution in two separate parts as also displayed by the vertical dashed lines in panel (a) of Fig. 2.

Deep Convolutional Neural Network (DCNN)
Fig. 3 schematizes the ML problem set-up. A plane of velocity module randomly extracted from the dataset of 10 simulations is given in input to a DCNN with the goal to infer the Ω value of the simulation from which the plane has been generated. The training is performed in a supervised setting, having pre-labelled all planes with their corresponding value, Ω true . As solution of the regression problem, the DCNN gives in output a real number corresponding to the predicted rotation value, Ω pred . In this way the inferring task is not limited to a predefined set of classes and can be generalized to estimate rotation values never observed during training. The network is trained to minimize the following loss function; which is the mean-squared-error between the target value Ω true and the predicted value Ω pred . E is the mean over a mini-batch of 256 different planes. To implement the DCNN we used TensorFlow [73] libraries with a Python interface, optimized to run on GPUs. Because of total dataset size, 245GB for the 500k different planes at a resolution of 256 2 grid points, the training had to be parallelized over eight A100 Nvidia GPUs. As described in Fig. 4, the overall architecture is composed of four convolutional layers that encode the input plane of velocity module on 32 different planes of size 8 × 8. This first   part of the network is asked to aggregate the relevant features for the regression problem. After, the results for the aggregated features are reshaped into a single vector of 2048 elements. The second half of the network, is composed by two fully connected layers that transform the features vector into one real number representing the network prediction, Ω pred . To reduce overfitting the training of these two final fully connected layers is performed using a dropout respectively of 30% and 20% [74]. All layers are followed by a non-linear ReLu activation functions. Using a mini-batch of 256 planes, and a training dataset of N tra = 400k planes, each epoch consists of 1563 backpropagation steps. The backpropagation has been performed with an Adam optimizer equipped with a Nesterov momentum [75] using a decreasing learning rate starting from 10 −4 with a scheduler that consists in a linear decrease of the learning rate up to reach 10 −5 in 20 epochs and then to remain constant. To achieve a good training with the learning protocol just introduced, we have normalized the input data to be in the [0; 1] range, namely we rescaled the module intensity of each gird-point by normalizing it to the maximum value calculated over the whole training set. The same normalization constant has been used to rescale the validation set. In the same way the training labels Ω true have been rescaled by their maximum value in the training. In Here Ω true = 20 is used, but similar behaviours are observed for all the other Ω. At epoch 200 the PDF becomes close to a delta-like distribution, suggesting very good accuracy of the network, with errors only every few cases over onethousand, and within 30% the correct value. In the central column, panels (b,e) of Fig. 5, we compare the same results obtained varying the size of the DCNN, precisely the number of its trainable weights keeping fixed the architecture. Here the label DCNN-f8 and DCNN-f64 represent two neural networks that have respectively 8 and 64 different filters in each convolutional layer, instead of the 32 as in the case of DCNN. The two networks have a total number of parameters of roughly 145k for the smaller one and 1.4M for the larger network. We can see that when the network is too small the accuracy decreases on both the validation and the training sets while if the network is too large there is some overfitting error on the validation set. Also the PDFs is optimal for the intermediate network labeled as DCNN. In the last column of Fig. 5, panels (c,f), we compare the DCNN results trained on smaller datasets. Here we can appreciate the quick deterioration of the prediction when reducing the quantity of the data used in the training.

Bayesian Inference
The second statistical inference approach considered is the Bayesian Inference (BI) [76,77,78]. The idea is to infer properties of the flow assuming the knowledge of some underlying probability distributions. More specifically, the Bayes' Theorem tells us that the posterior probability, P (Ω i |O), of observing Ω i given a measure of another generic observable, O, can be estimated as follow; P (Ω i ) is the prior probability, in our case for all Ω i equal to 0.1 because the dataset is composed by the same number of planes for each of the 10 Ω i . P (O|Ω i ), instead, is the likelihood to get the measure O of a generic observable for a fixed Ω i . Once known the posterior probability in Eq. (4), a first Bayes estimation can be obtained as follows; which correspond to the mean value on the posterior probability. A second possible estimation can be done as the most likely Ω i given O, namely; where in the last step the denominator of (4) is neglected because it does not affect the maximization on Ω i . For our problem, the two more relevant physical observable are the total energy and the total squared velocity gradients, representing large-scale (Gaussian) and small-scale (non-Gaussian) features, respectively. Adding the information of flow time evolution would be beneficial to have a measure of the waves propagation inside the flow, but as discussed this is in general a much more complicated measure and we are not assuming to rely on it. For this reason, in the sequel we will proceed consider u 2 and (∂ x u) 2 , where u = |u| is the module of the horizontal velocity and · · · stands the average over the plane. For each of these two observable we have calculated the likelihood P (O|Ω i ) from the 40k planes for each Ω i of the training dataset. Namely we have measured P ( u 2 |Ω i ), P ( (∂ x u) 2 |Ω i ) and the joint-likelihood, P ( u 2 , (∂ x u) 2 |Ω i ) on the training set. Measuring the same observable on each plane of the validation dataset, and knowing the likelihoods for each Ω i , we could estimate the posterior probability P (Ω i |O) (Eq. (4)) and than infer Ω Bay (O) (Eq. (5)). Fig. 6 shows the joint-likelihood of the mean energy and mean gradients. From this plot we can see that also combining large and small scales observables does not allow to have a clean statistical separation of the different Ω i for all the 10 'ensembles' considered here. In the next section we present and discuss results of Ω inference via DCNN and standard Bayesian approach.

Results
A first quantitative comparison of the two techniques is presented in Fig. 7, where in the four panels we report the scatterplot between the reference rotation value, Ω true and the corresponding prediction obtained via the DCNN (panel a), and BI with the three different likelihoods, namely the one based on the mean energy (panel b), on the mean square gradients (panel c) and on the joint likelihood combining energy and gradients (panel d). The four panels are all obtained with the same validation dataset of 10k planes for each of the 10 Ω i values. The first important observation is that the DCNN prediction is surprisingly close to the correct values not only on average but also for each single sample. Indeed, from the top left panel of Fig. 7 we can see that all points in the scatter plot fall well inside the hexagon symbol, just few exceptions are observed for the cases of Ω true = 8 and 10, and for Ω true = 2 and 4. Let us stress that this result is far from being trivial and provides a first example where a ML data-driven technique succeed in the inferring of an hidden parameter from the analysis of real fully developed turbulent configuration. The reason of this success is double, the intrinsic potential of DCNN and the good quality of the dataset produced and implemented in the training of the network. Indeed as already shown a dataset reduction quickly leads to a deterioration of the training results, using a dataset up to eight times smaller we have seen that the DCNN predictions does not change on average but the points distribution The second important consideration can be drawn by the observation of the BI results. Here, looking at panels (bc-d) we can assess the ranking of different features for the classification task. When using only the energy, BI is able to distinguish among the two rough classes (below/above the energy transition) only, while it is completely blind to differences among the intermediate Ω values. In panel (c) we see that measuring small-scale features is more useful for the classification, at least on average. Panel (d) shows that combining small and large scales input improves the results, keeping the BI close to the ML approach. Keeping adding different observable to condition the BI approach is clearly doomed of fail, because of the difficulties to have a large statistical sample to probe the conditional probability on a larger and larger dimensional space. The detailed comparison between DCNN and the best BI with joint likelihood on velocity module and gradients is presented in Fig. 8, where 10 different panels report the PDF of the DCNN and BI predictions, for all planes and all Ω of the validation set. Each PDF is normalized to its reference rotation value Ω true such as Ω pred /Ω true = 1 indicates the perfect prediction. Predictions in the small rotation regime (Ω = 2, 4) are the more difficult because no large-scale coherent structures are expected to appear in the flow and the resulting field is very little affected by the presence of rotation, so the predictions are easily confused among these two values. Comparing the two PDFs panels (a) and (b) we notice that while Bayes peaks always at the center between the two values, the DCNN is correctly peaked at the right value. In general for all the analyzed Ω (panels a-l) the DCNN accuracy is always very high and the larger errors happen only for a few cases on the 10k planes of the validation set. The last comparison between DCNN and BI is presented in Fig. 9. Here in panel (a) there is a summary with the statistics of the different DCNN and BI predictions. To do so we measure the Mean Squared Error (MSE) prediction normalized to the target value, From this figure we can see that DCNN is superior to the BI in a consistent way for all rotation rates. On the other hand, in panel (b), we focus on the sensitivity of   Fig. 8. PDF of rotation value predicted by the neural network, Ω pred (DCNN, solid lines full triangles), compared with the mean, ΩBay (red histograms), and most likely, Ω * Bay (black histograms), Bayes inference using joint-likelihood P ( u 2 , (∂xu) 2 |Ωi), on the validation set. For the sake of data presentation each of the PDF measured on the 10 different Ω values are reported in 10 different panels (from 'a' to 'l'). As indicated in the axis, the PDF are all normalized by their corresponding Ωtrue indicate in the labels of each panel.
the two inferring tools to distinguish parameters used in different simulations. In particular we aim to identify from the analysis of all planes of each simulation the ones that belong to a specific rotation class, that in this case was chosen to be Ω = 8. In Fig. 9(b) we measure the mean  normalized distance, between the prediction Ω pred and the value of interest, Ω = 8, averaging over all planes from taken from the same simulation as a function of the different Ω values used in the simulations. We can see that the distance is sharply minimized when the planes come from the DNS with the chosen rotation Ω = 8. In other words, Fig. 9(b) highlights the sensitivity of the data-driven tools in classifying complex turbulent conditions from 'ensembles' of different environments.

Occlusion
Here we present the ablation study on the input data aimed to perform a ranking on the importance of the flow features used by the DCNN to get a correct prediction, see [79] where a similar approach has been used on plane Couette flow to characterize the relevance of different regions for the prediction of relaminarisation events. The idea consists on keeping fixed the best DCNN trained model and to provide in input data where a fraction of the information were occluded such as to evaluate the deterioration of the inferring performance as a function of the different features occluded. As illustrated in Fig. 10 the first type of ablation performed is done in real space and consists in applying a filter (zero-mask) over all high or low energy regions up to a given percentage of points in the plane. In this way we aim to occlude a particular spatial input region, characterized by its energy, in order to evaluate its relative importance for the network predic-tion. The real space occlusion fields can be defined as; where u < o is occluded field on the high energy regions, i.e. the mask removes all points where the energy is larger than e h . The same idea holds for u > o , but now this field is masked on the low energy regions, namely on all points with energy lower than e l . The two thresholds (e h , e l ), are chosen such that the percentage of points occluded is equal in the two cases. In Fig. 10 there is a visualization for both the high energy and low energy ablation for three different percentages of occlusion and three different rotation values.
To evaluate the importance of the occluded regions in terms of the network prediction, in Fig. 11 we measure the MSE, as defined in Eq. (7), obtained for the predictions on the occluded fields, varying the percentage of occlusion and comparing the ablation of high or low energy regions. Also in this case it is useful to distinguish between the weak and the strong rotation rates. The panels 'a' and 'b' show that for Ω = 2 the high energy regions are more important, and their ablation deteriorates the prediction. On the other hand for Ω = 4 the network focuses on the low energy regions which are found to be more important by this ablation study. All the remaining panels from 'c' to 'l' correspond to the strong rotation rates, and support in a consistent way that in the inverse cascade regime the intense vortical regions are more important for the correct classification. It is interesting to observe that for few cases, Ω = 5, 6 and 20, the occlusion of a small percentage of the less energetic regions improved the DCNN prediction compared to the one obtained on full input. The second type of occlusion, illustrated in Fig. 12, is a Fourier-space ablation. Here we remove all the degrees-offreedom (DOF) contained in a specific circumference of radius |k| = k f ilt from the original plane in Fourier-space. Let us stress that the Fourier-space ablation is non-local in real space. Hence, all points of the occluded inputs are modified by the Fourier-space ablation. It is also important to notice that when increasing the filter wavenumber, k f ilt , the percentage DOF removed from the system becomes larger. To quantify the relevance of a particular Fourier frequency k f ilt we have evaluated the relative error between the DCNN predictions evaluated on the occluded and on the full input, defined as; where we have indicated as DCNN[u(x)] the network prediction and as u o (x, k f ilt ) the occluded field depending on the filter wavenumber. R(k f ilt ) > 1 indicates highrelevance for the occluded features because their removal deteriorates the prediction. On the contrary R(k f ilt ) < 1 suggest low relevance for the DOF at k f ilt whose removal improved the DCNN prediction. In Fig. 13 we report the results for the relevance R(k f ilt ) averaged over all planes of the validation set. Panel (a) shows results for a representative subset of the high rotation rates, Ω = 5, 8 and 30, while panel (b) shows results for the two lower rotation rates. The main plot in both panels shows the mean relevance R(k f ilt ) , as a function of the filter wavenumbers, averaged over all planes in the validation set of each Ω. As already observed from the real space occlusion, the role of the large scales depends on the rotation rate. When the rotation is strong, Ω ≥ 5, the relevance of the large and more energetic scales is preeminent. On the other hand, as expected, they are not relevant for the weak rotation rates. Furthermore, with the Fourier space occlusion we have access to the role of the smaller scales and we found that the DCNN is sensitive to the occlusion on a wide range of scales up to k f ilt ∼ 40. This supports the idea that the neural network has learned to evaluate the complex multi-scale statistical properties to improve its pre-  Mask on high |u| 2 Mask on low |u| 2 Fig. 11. MSE, as defined in Eq. (7), obtained for the predictions on the occluded fields, varying the occlusion percentage for both masks applied to the higher energy (green triangles) and to the lower energy regions (red circles). Panels (a) and (b) present results for the weak rotation rates (Ω = 2 and 4), while the remaining eight panels form (c) to (l) report results for the strong rotation rates where the fields are dominated by the intense vortical structures. diction. Only in the particular case of Ω = 2 we have not observed any particular relevance at any scale. In the inset of both panels (a) and (b) of Fig. 13 we report the R(k f ilt ) variance as a function of the filter scale. Again the variance is evaluated over all planes of the validation set of each Ω. From the variance we can highlight the presence of an exponential cutoff in the relevance of all wavenumbers larger than k f ilt ∼ 40. This observation tells us that such scales are dominated by viscous dissipation which is independent of the effect of rotation. In conclusion we have shown an attempt to use ablation to perform a features ranking on the flow properties more important for prediction of the DCNN. In our case we have found that the large scales swirling structures generated by rotation are generally more relevant for the correct classification of the high rotation rates. On the other hand we have seen that there is a full range of scales up to wavenumbers around k f ilt ∼ 40 that must be considered to improve the inferring prediction.

Conclusions
We have studied the problem of inferring physical parameters of turbulent flows with purely data-driven tools without assuming any knowledge on the dynamical laws governing the systems. The case of study is 3d turbulence on a rotating frame with the goal of inferring the rotation frequency from the analysis of the flow velocity module on a 2d plane observed at a fixed time from the domain top. This regression problem, without considering the temporal dynamics of the flow, has not a clear solution even considering equation-driven approaches, and it represents an open question for theoretical physicists working on rotating flows. In this work we take advantage of the new ML tools and of a big pre-labeled dataset from high quality DNS [72] to tackle such problem. In particular we have trained a Deep Convolutional Neural Network in a supervised way to obtain a ML model able to solve this complex task. To compare the 'supremacy' of DCNN results and to guide the reader towards a better understanding of the DCNN solution we have repeated the same game using standard Bayesian Inference (BI). For the latter estimation we have used two physically meaningful observables, moments of the velocity and of its gradients. In this way we have shown how the combination of more independent observables can allow to reach a better prediction, however simply increasing the number of observables would lead to a failure because of the difficulties to have a enough statistics to probe on the larger and larger dimensional space of the conditional probability. Performing an ablation study on the input data and fixing the best trained network we have shown that the DCNN is sensitive to the removal of small fraction of the input data, and it can be used to rank different flow features in terms of their importance for the prediction. Performing an ablation in real space we have shown that intense energy regions are relevant for a correct classification of high rotation rates. At the same time with a Fourier space ablation we have shown that DCNN is sensitive to changes in the input applied over a large range of wavenumbers, indicating the capability of the network to identify not only regions in space more relevant than other but also that it is measuring non-trivial and multiscale observables to improve its prediction. As briefly discussed in the introduction, an unsupervised classification of the same dataset has been attempted using a GMM classification algorithm, however no relation between the classes identified and the rotation value has been found. Indeed GMM, as more in general kmeans algorithms [80], are mainly sensitive to the number of vortexes inside each plane or to their spatial location and not to observables that are related to the rotation rate. For this reason a detailed analysis using this tool has been omitted in the paper. The classification of complex turbulent environments can be extended in several ways and has many real world applications. One example, consists in inferring of the intensification of extreme cyclones in the atmosphere [81,82], where learning how to combine standard information with more complex data as satellites images can lead to better predictions of such phenomena. Another potential application is in self-navigation problem, where autonomous vehicles need to know the on-time flow conditions before to choose the best policy to optimize their navigation [7,9,8,12,13]. Also in the realm of turbulence modelling it would be beneficial to perform an ad-hoc fine-tuning of the closure model depending on the specific local turbulence condition to improve the training performances [19,20,21]. Our work show that DCNN can solve this kind of inferring problems with unprecedented accuracy still under the limitation of having a large quantity of good quality data.