Machine-Learned Closure of URANS for Stably Stratified Turbulence: Connecting Physical Timescales & Data Hyperparameters of Deep Time-Series Models
Muralikrishnan Gopalakrishnan Meena Demetri Liousas Andrew D. Simin Aditya Kashi Wesley H. Brewer James J. Riley Stephen M. de Bruyn Kops




Spectral Density and Sum Rules for Second-Order Response Functions
    Peter Abbamonte
    April 29, 2024
==================================================================




We develop time-series machine learning (ML) methods for closure modeling of the Unsteady Reynolds Averaged Navier Stokes (URANS) equations applied to stably stratified turbulence (SST). SST is strongly affected by fine balances between forces and becomes more anisotropic in time for decaying cases. Moreover, there is a limited understanding of the physical phenomena described by some of the terms in the URANS equations. Rather than attempting to model each term separately, it is attractive to explore the capability of machine learning to model groups of terms, i.e., to directly model the force balances. We consider decaying SST which are homogeneous and stably stratified by a uniform density gradient, enabling dimensionality reduction. We consider two time-series ML models: Long Short-Term Memory (LSTM) and Neural Ordinary Differential Equation (NODE). Both models perform accurately and are numerically stable in a posteriori tests. Furthermore, we explore the data requirements of the ML models by extracting physically relevant timescales of the complex system. We find that the ratio of the timescales of the minimum information required by the ML models to accurately capture the dynamics of the SST corresponds to the Reynolds number of the flow. The current framework provides the backbone to explore the capability of such models to capture the dynamics of higher-dimensional complex SST flows.

Notice: This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-planhttp://energy.gov/downloads/doe-public-access-plan).









§ INTRODUCTION


Stably stratified turbulence (SST) is a model for understanding fluid flows that involve an ambient stable density stratification in the oceans due to temperature or salt gradients.  Studying SST has important implications for climate modeling, pollution mitigation, and deep sea mining.  In general the dynamics of SST reflect a competition between many forces including inertial, viscous, buoyancy, Coriolis, and shearing.  Even in the simplest case, though, SST can be parameterized by a minimum of two dimensionless groups defined by ratios of the inertial, buoyancy, and viscous forces.  These ratios can be chosen in a variety of ways, and we use Reynolds and Froude numbers defined later.  Consistent with the parameter space of SST being, at a minimum, two dimensional, the flows are highly intermittent and anisotropic at large scales.  <cit.> provide an overview of SST, <cit.> review many theoretical, laboratory, and numerical studies in SST, and <cit.> present a framework for parameterizing a three parameter variant of SST that includes mean shear in order to describe the flow dynamics in terms of where the flow resides in parameter space.

The extreme scales required to capture the evolution of such flows in nature hinder computational experiments even with the latest high-performance computing platforms <cit.> and thus, the atmospheric and oceanographic communities have relied on parametrizations or surrogate models to capture the effect of the small scale flow features in SST. Turbulence closure models in Unsteady Reynolds Averaged Navier Stokes (URANS) simulations and large-eddy simulations (LES) of the atmosphere, ocean, and their sub-systems are broad classes among many such turbulence models. Particularly, URANS second-moment closure (SMC) based modeling has been prevalent in simulations of stratified shear layers and wakes <cit.>. Such models are limited by the appropriate choice of the parametrization, interpretation of the parametrizations, and generalizability with different flow regimes and types. Data-driven techniques are attractive to overcome such challenges. Over the past decade, there is a history of using machine learning (ML) methods to tackle fluid flow problems <cit.>, and specifically, to assist in the augmentation of turbulence models <cit.>. The atmospheric and oceanographic communities have been using ML to formulate surrogate models to capture the effects of sub-grid scale structures. Such efforts have been focused on turbulence closures for LES and global climate models <cit.>. There have been various efforts to use data-driven models in modeling other applications of stratified turbulent flows as well, such as shear layers <cit.>, stratified wakes <cit.>, engineering applications <cit.>, optical turbulence <cit.>, and mixing in SST <cit.>.

Developing predictive capability for SST is the first motivation for this study.  A second is to apply machine learning (ML) to the difficult problem of closing the URANS equations or the related approach of LES.  In URANS, differential equations for averaged quantities are solved numerically with models for terms involving fluctuations about the averages.  In LES, equations are solved for quantities resolved on the numerical grid and subgrid-scale quantities are modeled.  With URANS, even in configurations involving just two types of forces, e.g., inertial and viscous, all but the simplest closure models typically involve differences between model terms.  SST flows involve more than two forces, so there are more terms to be modeled in either URANS or LES, and there is limited understanding of the physical phenomena described by some of these terms. Rather than attempting to model each term separately, as is done in SMC modelling, it is attractive to see how well ML can be used to model groups of terms, that is, to directly model the force balances.

Furthermore, the inherent nonlinear dynamics of the SST, even in a reduced-order sense, pose the challenge of finding interpretable and generalized surrogates. Thus, the third and final motivation of this study is to formulate interpretable ML models that accurately predict force balances. Specifically, we will focus on interpreting the data requirements of time-series-based ML models for modeling a reduced dimensional representation of canonical SST.

Advances in ML have shown the capability of ML models to implicitly learn the underlying dynamics of a time-dependent physical system <cit.>. Recurrent neural networks, such as Long short-term memory (LSTM) <cit.> and gated recurrent units (GRUs) <cit.>, have been shown to model and predict time-series data effectively.
Thus, time-series ML models could provide a computationally inexpensive way to simulate complex time evolving force balance of the SST URANS equations. The task of interpreting such models still remains a challenge <cit.>. State-of-the-art approaches involve using techniques such as model visualization, layer-wise relevance propagation, and attention mechanisms <cit.>. While such techniques provide valuable insights into the working of the model, extraction of physical knowledge about the system being modeled is limited. One step in this direction is to use ML models for extracting the timescales of the physical system, which are fundamental to understanding complex dynamical systems.

In this paper, we report on the use of time-series ML to close the URANS equations for one of the simplest configurations of SST, namely, homogeneous and decaying turbulence.  The SST configuration and URANS equations are presented in <ref>. In <ref>, we review time-series ML approaches, the training data sets, and techniques for interpreting the data requirements of the ML models. The results from our URANS models are discussed in <ref>. Some conclusions are drawn in <ref>.








§ THEORETICAL BACKGROUND FOR URANS OF SST



For this work, we consider temporally evolving SST that is axisymmetric in the sense that the statistics are independent of horizontal direction.  The flow is stably stratified by a uniform density gradient so that it is homogeneous in the vertical direction.  While the configuration is formally homogeneous in all directions, we show in <ref> with contour plots from direct numerical simulations (DNS) that the flow is intermittent at length scales that are a significant fraction of the domain size. We take advantage of the homogeneity to reduce the RANS grid down to one point in each direction in <ref>.  Here, let us write the governing equations and their Reynolds-averaged versions for the general case.

The modeled flow evolves with density stratification and is governed by the dimensionless Navier–Stokes equations satisfying the Boussinesq assumption, given by




    ∇·u   = 0, 
    ∂u/∂ t + u·∇u   = -( 2π/Fr)^2 ρe_z -∇ p + 1/Re∇^2u, 
    ∂ρ/∂ t + u·∇ρ - w    = 1/Re Pr∇^2 ρ.

Here, u = (u,v,w) represents the velocity vector in the Cartesian coordinate system x = (x,y,z) evolving in time t, and ρ and p are the deviations of density and pressure from their ambient conditions, respectively. The velocity and length scales for non-dimensionalzing the variables are the r.m.s. velocity Û and integral length scale L̂, respectively. Note that the dimensional quantities are represented with ·̂. Time t is non-dimensionalized as t = (t̂ - t̂_0)/t̂_LE, where t̂_0 is the initial time at which the turbulence decay process starts and t̂_LE is the large-eddy advective time scale. The non-dimensional quantities, Froude number Fr, Reynolds number Re, and Prandtl number Pr, are defined respectively as

    Fr = 2πÛ/N̂L̂, Re = ÛL̂/ν̂,  and  Pr=ν̂/α̂.

Here, N̂ = [-(ĝ/ρ̂_0) (dρ̂/dẑ)]^1/2 is the Brunt–Väisälä frequency or the buoyancy frequency, ĝ is the acceleration due to gravity in the -z direction, ρ̂_0 is the density at a reference height of ẑ_0, ν̂ is the kinematic viscosity, and α̂ is the thermal diffusivity. The density gradient dρ̂/dẑ is a uniform, time-invariant, stable ambient density gradient applied to impose stratification.

For the homogeneous case, we can approximate the Reynolds average with a spatial average over the entire domain, denoted by < ·>. Then the time evolution of the average flow can be described in terms of energy equations by decomposing kinetic energy into its horizontal and vertical contributions, E_H and E_V, which are coupled via pressure. The potential energy E_P is a scalar multiple of the variance of the buoyancy. An additional equation is required for the buoyancy flux, which couples E_V and E_P. The resulting set of evolution equations is




    d/dt<E_H>    = < p ∇_H·u_H >_needs model - ν̂< ( ∂u_H/∂ x_j)^2 >_needs model
    d/dt<E_V>    = < p ∂ w/∂ z>_needs model - < bw > - ν̂< ( ∂ w/∂ x_j)^2 >_needs model
    d/dt<E_P>    = < bw > - α̂/N̂^2< ( ∂ b/∂ x_j)^2 >_needs model
    d/dt<bw>    = < p ∂ b/∂ z>_needs model + < w^2 >N̂^2 - < b^2 > - ( ν̂ + α̂) < ( ∂ b/∂ x_j) ( ∂ w/∂ x_j) >_needs model

where E_H = (1/2)|u|^2_H, E_V = (1/2)w^2, and E_P = b^2 / (2N̂^2). Here, b=ĝ/ρ̂_0ρ̂ is the buoyancy, and |u_H| = ( u^2 + v^2 )^1/2 is the horizontal velocity. Several terms on the right-hand-sides (RHSs) of these equations are marked to indicate that they must be modeled in order to close the system of ordinary differential equations. Rather than model each term, or even the RHS of each equation, our objective is to model the entire set of RHSs using ML techniques, as broadly depicted in <ref>.









§ APPROACH






 §.§ Procedure


As mentioned in <ref>, the objective is to model the temporal evolution of the energies (horizontal and vertical kinetic energies, and potential energy) and buoyancy flux for a single point in space. Approximating the derivatives of the energy variables are of interest for modeling the RHS of <ref>-<ref>, needed to close the system. The terms to be modeled can be represented as

    f = [d/dt<E_H>   d/dt<E_V>   d/dt<E_P>   d/dt<bw>]^T.

The terms are non-dimensionalized by the total energy, E_T = <E_H>+<E_V>+<E_P>, and buoyancy frequency, N̂, at the initial condition, giving <E_H> = (1/2)|<u_H>|^2/E_T_0, <E_V> = (1/2)<w>^2/E_T_0, <E_P> = <b>^2 / (2N̂^2E_T_0), and <bw>/(2N̂E_T_0) as the terms to be modeled. Here, we can represent the independent variables as q = [ |<u_H>|^2   <w>^2   <b>^2   <bw> ]^T. Thus, the problem can be generically represented in discrete time as

    dq(t_i)/dt = f(q(t_i), t_i)

where t_i represents a discrete time instant. We use supervised deep learning techniques <cit.> to model the relationship between q and f, and thus, model the time evolution of q, as shown in <ref>.

Neural networks have the capacity to approximate any Borel measurable function if it has the proper architecture <cit.>. We are interested in neural networks capable of modeling the RHS of the system of ODEs in <ref>, representing the reduced form of <ref>-<ref>, so that any numerical ODE solver can solve the time evolution of the states, q(t). 
Moreover, as the problem is fundamentally dependent on time and the time history of the states, modeling the memory effect present in the time series using time-series ML models is attractive. Thus, our objective revolves around using sequences of time-series data to predict the time evolution of q. Here, sequence refers to the ordered snapshots of state time history (or memory) used to model the dynamics of the system.

Generally, sequences of time-series data have two parameters analogous to timescales in physical time-evolving systems, as illustrated in <ref>. One timescale is the sequence length, which is the time span of the data in a single sequence. The second timescale is the sampling period of the data within one sequence. We will refer to latter timescale by its inverse, the sampling frequency. With respect to physical timescales, the sequence length can be considered as the largest timescale of the system being modeled whereas the sampling period can be treated as the shortest timescale of the system. We will elaborate more on what these timescale mean for SST in <ref>.



We use these underlying timescales of the temporal data as hyperparameters of time-series ML models. Furthermore, we determine the minimum information, based on these timescales, required by the  ML models to effectively model the system dynamics. These techniques are used specifically for interpreting the data requirements for the time-series ML models.

We choose to model the evolution of the flow states using two routes: (1) model the time derivatives (or the entire RHS of <ref>-<ref>) using an ML model and pass this information into an ODE solver to obtain the time evolution of q and (2) directly model the time evolution of q using an ML framework.
We achieve these two routes using the Long Short-Term Memory (LSTM) <cit.> and Neural Ordinary Differential Equation (NODE) <cit.> models, respectively.
In the following sections, we describe the architecture of the LSTM and NODE models, and the dataset used to train these models.





 §.§ Neural network setup


Without a loss of generality, the system of differential equations governing the SST, <ref>, are modeled as a neural network (NN), 𝒩, given by

    dq(t_i)/dt=f̃(q(t_i),t_i,θ)_𝒩

where the NN is parameterized by a tensor, θ, containing learnable weights and biases defining the map between the input and output layers of the network. The symbol ·̃ represents the quantities obtained from the learned ML model.





  §.§.§ Long Short-Term Memory (LSTM) model


In the first route for modeling the system of ODEs, given the components of q at an instant in time t_i, we use NNs of the form 𝒩(q, θ) to model the RHS of <ref>, with θ representing the parameters of the NN (“weights”). Specifically, we use the Long Short-Term Memory (LSTM) <cit.> network as the NN model. An LSTM architecture is a sub-class of the recurrent neural network (RNN) architecture, which uses internal states of nodes (memory) to process arbitrary sequences of inputs and thus, allows for time-series to be modeled <cit.>. LSTM is commonly used for time-series prediction because of its improved stability and predictive capabilities over generic RNNs <cit.>. A sample depiction of the model is shown in <ref>.



The input to the LSTM model is the time history q(t_i-n), q(t_i-(n-1)), …, q(t_i) where n is the sequence length. The output of LSTM is f̃(q(t_i), t_i,θ), which is fed into an ODE solver (like a Runge–Kutta scheme) to obtain the next value in the series, q(t_i+1). The process is repeated to predict the entire time evolution of q. Thus, the entire workflow can be viewed as a LSTM+ODE numerical simulation, where the LSTM serves as the RHS of the ODE solver. The LSTM model is trained using the loss function comprised of the RHS of <ref>, given by

    ℒ[f(q(t_i),t_i), f̃(q(t_i),t_i,θ)].



Hyperparameter selection for the LSTM model is done through experimentation (broad grid search). We found that for modeling the SST problem, the best predictions are obtained with four LSTM layers with a hidden size of 10. We also use a multi-layer perceptron (MLP) consisting of one hidden layer with 15 neurons before the output layer. A Leaky Rectified Linear Unit (LeakyReLU) activation function with a slope of 0.1 is applied for both the LSTM and MLP layers. We measure the mean squared error as the loss function, ℒ, for optimizing the network weights during the training process. We select a learning rate of 0.001 and train the model for 4000 epochs, after which the model training converges.




  §.§.§ Neural Ordinary Differential Equation (NODE) model


In the second route for modeling the system of ODEs, given the components of q at an instant in time t_i, we use an ML framework devised by <cit.>, Neural Ordinary Differential Equations (NODE), to directly predict the time evolution of the system. The NODE framework models an arbitrary system of ODEs as an NN given discrete observations. <cit.> refer to this framework as a time-stepper model compared to direct solution models which map a given time, t_i, to its observation, q_i, such as deep neural networks (DNN), recurrent neural networks (RNN), and residual neural networks (ResNet) <cit.>. Motivated by the similarity between ResNets and Euler's method for numerical integration, wrapping a black-box ODE solver around an NN offers many modeling advantages. The main advantage being that, unlike ResNets or traditional sequential learning architectures, transformations between hidden states evolve continuously. Many systems described by differential equations evolve as such, so the combination of well known mathematical integration tools and neural networks can be a powerful paradigm in scientific applications. In recent years, NODE has been used to model complex systems in fluid dynamics, such as to model and predict thermoacoustic instability <cit.> and learn subgrid-scale models <cit.>. Related work also includes the implementation of ODE-based ML solutions for modeling dynamical systems <cit.>, reduced-order modeling <cit.>, and turbulent flow control <cit.>.

As depicted in <ref>, NODE takes as input the current state of the system q(t_i) and the desired time instants t_i+1, t_i+2, …, t_i+n at which to evaluate the state. The output is the state variables q(t_i+1), q(t_i+2), …, q(t_i+n) at those discrete time instants. What makes a time-stepper model unique compared to direct solutions, is that the loss is computed after integrating the learned approximation of the differential equation, <ref>. As a result, the input to the scalar loss function ℒ is the output of the ODE solver, given by

    ℒ[q(t_i+n), q̃(t_i+n)] = ℒ[q(t_i+n), q(t_i)+  ∫_t_i^t_i+nf̃(q(t),t,θ)dt] = ℒ[q(t_i+n), ODESolve(q(t_i), f̃,t_i,t_i+n,θ)]

where t_i and t_i+n represents the initial and final time instants for a given sample of the time-series data. Gradients can be computed using the adjoint sensitivity method, solving an augmented ODE backwards in time. The loss function can then be minimized by comparing the integrated prediction, q̃(t), and observed solution, q(t), in relation to the tensor, θ, which parameterizes the differential equation. Using this framework, recovering both the approximation of the state variables and their derivatives are convenient. Derivatives are recovered by inputting the predictions back into the neural network defining the temporal derivative of the system, f̃(q(t),t,θ) (<ref>).



In literature, NODE has been shown to be less sensitive to hyperparameter selection compared to traditional NN models, like MLPs and ResNets <cit.>. This pertains to the hyperparameters of the model architecture and not the two data hyperparameters of interest in the current study – sequence length and sampling frequency. For the model hyperparameters of the NODE architecture, we chose an absolute tolerance of 10^-3 and a relative tolerance of 10^-4 for the integrated ODE solver <cit.>. Furthermore, the MLP within the NODE model has 10 layers with 40 neurons per layer. We use a Sigmoid Linear Unit (SiLU) activation function for these layers. Mean squared error is used as the loss function, ℒ, for optimizing the network parameters, θ, during the training process. For this configuration, we used a higher learning rate of 0.05 compared to LSTM and lower number of epochs of 1000 to train the NODE model. This is because the NODE model converged faster with similar performance.




  §.§.§ Comparing LSTM and NODE models


As described above, the time-series prediction that we formulate the LSTM and NODE models to perform are fundamentally different. In general, the inputs and outputs of the LSTM and NODE models can be summarized as shown in <ref> and their time-dependence can be summarized as shown in <ref>. While the LSTM model inherently relies on the time history of the data, the NODE model utilizes the future time evolution of the data to model the behaviour of the system. Nonetheless, there are two timescales for the input data to both the models - the sequence length and sampling frequency. Using these two timescales as the main hyperparameter of the two models, we investigate the minimum information (with respect to these two timescales) that is required to effectively model a system of ODEs - assessed in terms of the accuracy and numerical stability of the resulting ODE solver.








 §.§ Training the models


The time-series models are trained using data from a decaying homogeneous SST flow simulation described in <cit.>. A sample portrayal of the time evolution of the SST flow is shown in <ref>. The flow is initialized using Taylor–Green vortices with low-level noise (10% of initial energy) to break the symmetries of the vortices and trigger instabilities for the flow to evolve into turbulence. The ambient stratification is constant and the initial density perturbation field is zero. This configuration is a significant test for URANS because it starts out laminar, then transitions to turbulence, and finally becomes increasingly more strongly stratified as it decays such that the inertial force decreases relative to the buoyancy force.   We consider the case with initial Reynolds and Froude numbers Re=3200 and Fr=4, which is discussed in detail in <cit.>.



As mentioned in <ref>, we take advantage of the homogeneity of the flow and spatially average the spatio-temporal data to a single grid point and obtain time series of the state variables. The training data consist of sequence of these time-series data, defined by two time scales: (1) the length of the sequence, sequence length, and (2) the time interval between samples within the sequence, the sampling frequency. We use these two variables as hyperparameters to interpret the data requirements of the models, where each model is trained for a single type of sequence defined by the hyperparameters. Note that these data hyperparameters are different from the model hyperparameters, which are chosen as mentioned in <ref> and kept constant throughout the analysis. For each hyperparameter set, the time-series is interpolated to have the same time-step as that of the sampling frequency. Thus, the number of samples used for training the models is determined by the sequence length as well as the sampling frequency. 

For a given set of samples defined by the data hyperparameters, we split the time-series to training and testing regimes based on the temporal evolution of the system. We chose to train the models untill t=20, as the data captures both the transient as well as the beginning of the turbulent decay regimes. The testing is done both for the unseen data as well as from the beginning (in the training regime) to test both the accuracy and numerical stability of the model. The full training data for a given set of hyperparameters is split with a ratio of 90:10 for training and validation at each training epoch. As the terms modeled are already non-dimensionalized, we did not have to perform additional scaling of the data, which is typical for NN training due to the limited range in which activation functions can operate.








§ RESULTS AND DISCUSSION



Here, we evaluate the ability of time-series ML to learn the dynamics of the reduced-dimensional SST by modeling the system of equations governing the force balance of SST URANS equations. We demonstrate that data requirements of these models are interpretable by extracting the physical timescales corresponding to the information learned by the model. We study the relationship between the physical timescales and the ML hyperparameters analogous to timescales of the time-series data – the data hyperparameters sequence length and sampling frequency. 

To demonstrate the interpretability framework, we first model the time evolution of a damped pendulum <cit.> as a sample model problem. The damped pendulum is a very simple and well-studied nonlinear dynamical system. It has two timescales that can be derived analytically, which makes it a suitable choice for this study. Although the damped pendulum problem is not as complex as the nonlinear SST that motivates this study, we use the damped pendulum problem as a canonical model problem to demonstrate our framework of extracting physical timescales from time-series ML models <cit.>. Herein, we first show the results for the data hyperparameters (sequence length and sampling frequency) for which we obtained the best results and then we discuss the effects of these hyperparameters on the model performance and how they can be used to interpret the data requirements of the models.



 §.§ Sample test problem: Simple pendulum


Machine learning models for both the LSTM and NODE are tested and validated with a simpler, well studied, synthetic dataset describing the motion of a damped pendulum. A damped pendulum exhibits some of the fundamental behaviours observed in the results for homogeneous SST, making it a relevant dataset for evaluating model performance. The governing equations of the pendulum are typically written in terms of its angular position θ and angular velocity ω. However, many dynamical systems in engineering and nature can be described by energy equations. In particular, the motivation for this study, URANS modeling of SST, is formulated in terms of energies as described in <ref>-<ref>. The system of ODEs governing the time evolution of the pendulum in terms of its kinetic energy (E_K=1/2mℓ^2ω^2), potential energy (E_P=mgℓ(1-cosθ)), and the flux (F=mgℓωsinθ) between them is given by




    dE_K/dt   = -2 b/m E_K - F 
    dE_P/dt   = F 
    dF/dt   = 2g/ℓ(E_K-E_P) + E_P/m ℓ^2(2E_K - E_P) -b/mF ,

where the parameters of the pendulum are the damping coefficient b, mass m, acceleration due to gravity g, and length ℓ. Thus, the state of the system is defined as q = [E_K    E_P   F]^T. Solving the system of ODEs requires the initial conditions E_K_0, E_P_0, and F_0. The similarities of the states of the pendulum and the SST problem are notable.

The pendulum has two fundamental timescales that govern its motion. They are the oscillating timescale τ_o and the damping timescale τ_d. These can be derived analytically using the small-angle approximation <cit.> and are given by




    τ_o    = √(ℓ/g)
    τ_d    = m/b.


We can use these timescales to determine the theoretical minimum information that a data-driven model would need to effectively learn the dynamics of a pendulum. An accurate model must capture the information corresponding to both the damping as well as oscillating timescales.



  §.§.§ LSTM and NODE results


We select a damped pendulum with damping coefficient b=0.5, mass m=1, acceleration due to gravity g=25, length ℓ=l, and initial conditions E_K_0=0, E_P_0=1, and F_0=0. Note that for these parameters, τ_d >> τ_o. Thus, the pendulum is underdamped and oscillates considerably before damping, similar to SST. The data for training and testing the models were generated numerically using these settings of the pendulum. Training and testing of the models were performed at different instants in time.


We show the results for these best LSTM and NODE models in <ref>. Both models perform accurately and are numerically stable. We observe that the NODE model slightly outperforms the LSTM model. It is interesting to note that the total energy of the pendulum E_T, <ref> (bottom-right), is accurately replicated by both the LSTM and NODE models. E_T was not explicitly modeled by the two time-series ML models, nor were any constraints imposed on the model training to follow the decay of E_T. The values for E_T are computed post-simulation, by summing the kinetic and potential energies from the ML models.





  §.§.§ Interpreting the ML model data requirements: minimum information required by model


The interpretability of the data requirements of the time-series ML models is studied by attempting to extract physically realizable timescales from the ML models. For this, the data hyperparameters, sequence length and sampling frequency of the input data, are varied to see the effects on the model performance. The objective is to identify the minimum values of the hyperparamters which results in good model performance, thus finding the minimum information that is required by the ML models to effectively learn the system dynamics. 

The results of the parametric study of the two ML timescales are shown in <ref>. The data hyperparameters are varied and the corresponding performance of the LSTM and NODE models are observed using 3-D contour plots (<ref> (a-b)) and their corresponding 2-D projections (<ref> (c-d)). The accuracy of the models (vertical axis in <ref> (a-b) and the contour levels in <ref> (c-d)) is measured by the maximum normalized root mean square error of all the states, max_q(NRMS).



For the LSTM model, small values of sampling frequency and sequence length have extremely high errors as noted by the lower-left corner in <ref> (a) and (c). This intuitively makes sense because when enough information is not provided to the model, the model is not able to learn the dynamics of the system. The 2-D projection of max_q(NRMS) shown in <ref> (c) reveals that this region lives below the minimum information line computed analytically (black dashed curve). The minimum information line corresponds to when the ratio of the sequence length and the sampling frequency is equal to the ratio of the damping timescale, τ_d, and oscillating timescale, τ_o, of the pendulum. 

Although there are local regions of high relative error, the general trend is that the LSTM model converges to accurate models beyond the minimum information required to represent the dynamics of the system. Thus, the inherent timescales of the minimum data input for accurate LSTM models corresponds to the physical timescale of the system.

The model error for the NODE models has a spike below the minimum information line similar to LSTM (lower-left corner of <ref> (b) and (d)), whereas closely after the minimum information line, the NODE model performs accurately. Again, these results demonstrate that analyzing the inherent timescales of the input data to the time-series ML model can enable one to extract the timescales of the physical system. This ML-learned timescale correspond to the minimum information required by the time-series ML models to accurately predict the system dynamics.

The NODE model performs notably better than LSTM in the parametric study, as seen from <ref> (a-d). Firstly, the highest error observed for the NODE model is an order of magnitude smaller than that of the LSTM model. Furthermore, NODE is much more numerically stable above the minimum information line since the surface plot is much flatter (and smoother) than that of LSTM, which has many local minima and maxima.




 §.§ SST initialized with Taylor–Green vortices


We next demonstrate the capability of the time-series ML models to capture the nonlinear dynamics of the energy equations of the SST flow initialized with Taylor–Green vortices. We perform both a priori (offline test) and a posteriori (online test) testing of the LSTM and NODE models.



  §.§.§ a priori testing


The a priori or offline test involves testing the accuracy of the models for a given set of inputs. For the LSTM model, we provide values of q at random time instants as inputs to the model and compare the outputs of the model with true or expected values – the time derivative of q (i.e.,  f) at the respective instants. The testing is performed at conditions within and outside the training data regime. The results of the a priori test on the LSTM model are shown in <ref>, sorted in the order of the time evolution. We only show the results from the best model identified from the timescales analysis, which will be discussed more in <ref>. The overall accuracy of the model for each variable is measured using the coefficient of determination, R^2, with respect to the true time derivatives of q. The results demonstrates the high accuracy of the LSTM model to capture the complex nonlinear relationship between q and its time derivatives. For the NODE model, the outputs are the time evolution of q and thus corresponds to a posteriori or online testing.





  §.§.§ a posteriori testing


The a posteriori or online test involves testing the accuracy and stability of the models over the time evolution of the system. For the LSTM models, this involves coupling the output of the model with the simulations, which in this case is an ODE solver run using a fourth order explicit Runge–Kutta scheme. The outputs from the LSTM model, time derivatives of q at an instant (or the RHS of <ref>, f), serves as the RHS function for the ODE solver. We validate the numerical accuracy and stability of the ODE solver by comparing the time evolution of q obtained from DNS with that obtained from the ODE solver when using the true values of time derivative of q as the RHS of the solver.

We show the results from the best LSTM and NODE models (based on the timescales analysis) in <ref>. Both the LSTM and NODE models are able to accurately predict the time evolution of q. In terms of the overall accuracy measured by the R^2 values, the NODE model offers better accuracy. Moreover, the predictions are numerically stable over time. Note that although the results show the model initialized from t=0, which is within the training regime, we also tested the model at various initial conditions within and outside the training regime. All the tests resulted in similarly accurate and numerically stable predictions. We wanted to demonstrate that the error acquired by the models when tested in the training regime (see the predictions for E_P in the interval t = [0,12]) does not influence the predictions in the testing regime drastically. The time-series models are stable even with accrual of small errors over time, which has been shown to be detrimental for the numerical stability of models focused only on spatial information <cit.>.

The results demonstrate the capability of the time-series ML models to accurately predict the time evolution of complex nonlinear dynamical systems like the reduced-order SST. Similar to the pendulum results, the capability of the ML time-series models to accurately predict the total energy of the turbulent flow is notable (<ref> bottom row), especially given that no information regarding the total energy was provided to the models. 
Note that, as a separate training campaign, we also trained the models with a physics-informed approach by imposing constraints of time rate of change of total energy, dE_T/dt, on the loss functions, thereby explicitly preventing the system from gaining energy. However, these results were not notably different from the current analysis. It is interesting to see the time-series models are able to predict the total energy of the system when the correct timescales are used for the inputs to the ML models.





  §.§.§ Interpreting the ML model data hyperparameter for SST


We apply the physical timescale extraction framework discussed in <ref> on the highly nonlinear SST problem. The current study is motivated by this turbulence modeling problem because there are no analytical solutions for the timescales in turbulent flow. In SST, there are, at a minimum, three important timescales: the buoyancy time scale, the inertial time scale, and the viscous time scale. Currently, the timescales can only be deduced from a combination of analytical work, laboratory or field measurements, and DNS – extremely computationally expensive campaigns, e.g. <cit.>. Even with these tools, the dominant timescales are hypothesized and then tested against data. This process becomes increasingly difficult when the flow becomes more complex due to, for example, mean shear. Of interest here is whether the relevant timescales could emerge from the ML modeling framework without any hypothesis. The pendulum results demonstrated this claim for a simple dynamical system.

We conduct a similar ML timescales parametric study on the SST data using both the LSTM and NODE models, as shown in <ref>. The parametric plot of timescales identifies unique regions of minimum error: between (1) sampling frequency 30-50 and sequence length 0.4-0.6 for LSTM, and (2) sampling frequency 20-40 and sequence length 0.5-1 for NODE. Interestingly, while the sampling frequency region is similar for both the models, the NODE model is accurate and numerically stable above a sequence length of 0.5. This may pertain to the increased numerical stability and accuracy observed for NODE. These local regions correspond to the minimum information required by the inputs of the ML models to effectively predict the SST dynamics. Based on our inferred knowledge gained from studying the pendulum problem, we can potentially claim that these timescales are related to the physical timescales of the system. 



To validate the above claim, we identify the points of lowest error in these regions, denoted by the blue and red dots in <ref> (a-b), respectively for the LSTM and NODE models. The ratio of the ML timescales, sequence length to sampling period (1/sampling frequency), at these points of lowest error are approximately 16 for LSTM and 32 for NODE. We note that these ML timescale ratios are integer multiples of the Reynolds number for the SST data, Re=3200. Physically, the Reynolds number of the SST flow can be defined by the ratio of the fluctuating timescale due to the buoyancy flux and decay timescale due to viscous dissipation <cit.>. Broadly, they represent the smaller and larger timescales of the system, which are what the time-series ML models extract as well. 
Thus, the ratio of the ML timescales are integer multiples of the physical timescale ratio known for SST flows. 
These observations validate our claim that our ML timescales analysis is able to interpret the behavior of the ML models and extract physical timescales of even complex turbulent fluid flow systems.

We generalize the above timescales extraction framework for various parameter ranges of SST flows initialized by Taylor–Green vortices <cit.>. These flow are parameterized by the Fr and Re. We chose to perform the timescale analysis on two flows: (1) Fr=2, Re=3200 and (2) Fr=4, Re=6400, thereby testing the timescale extraction framework for change in Fr and Re. The time series of the flow states for the three SST flows are shown in <ref> (a-e). Our objective is to compare the timescales extracted using the original/baseline analysis with the dataset Fr=4 and Re=3200 (blue lines in <ref> (a-d)) with that when (1) the Re is increased (orange lines in <ref> (a-d)) and (2) the Fr is decreased (green lines in <ref> (a-d)).



We do not change the model hyperparameters of the LSTM and NODE models, and obtain two different models by training the two datasets independently. The focus of generalizability here is not to have a single model to solve parameter ranges of the flow, but to investigate the generalizability of our timescale extraction framework. The results of the timescale analyses are shown in <ref> (f-h). When the Re is doubled and Fr is kept constant, there is marginal change in the decay rate and buoyancy flux as seen from the time-series in <ref> (a-d). The timescales extracted by the ML models also reflect this, depicted by the slight changes in the sampling frequency and sequence length as seen from <ref> (h). When the Fr is halved and Re is unchanged, the flow is characterized by longer times to decay and thus, lower decay rate, as shown in <ref> (a-c) and (e). Thus, one could expect the time-series ML models to require more time lengths of data to learn the dynamics – i.e., longer sequence length. This is what is reflected in the timescales extracted by the ML models as shown in <ref> (f). While the sampling frequency remains same, the sequence length increases, denoting the capability of the framework to capture the change in the decay dynamics of the system, thus, demonstrating the generalizability of the timescales extraction framework for various parameter regimes of SST flows.








§ CONCLUDING REMARKS



We present a time-series based machine learning (ML) modeling framework for second moment closure modeling of Unsteady Reynolds Averaged Navier Stokes (URANS) simulations of stably stratified turbulence (SST). Instead of individually formulating higher-order closures for the nonlinear terms in the URANS equations, we train time-series ML models using high fidelity direct numerical simulation (DNS) data to model the entire right-hand-side (RHS) of the URANS equations. The current approach is applied to one-dimensional SST data, which is spatially reduced by taking advantage of the homogeneity of the flow in three-dimensions. We choose two ML architectures for modeling: (1) Long short-term memory (LSTM) and (2) Neural Ordinary Differential Equation (NODE). The interpretability of the data requirements for the models are assessed by testing the accuracy of the models against two data hyperparameters, the sampling frequency and sequence length of the input data to the models, revealing two ML timescales of the data. The accuracy of the models based on these hyperparameters corresponds to the minimal information required by the ML models to effectively capture the dynamics of the complex system. We find that variations/trends in the ML timescales of the best model corresponds to variations/trends in physical timescales known for such flows.

The framework is first demonstrated for a simple dynamical system – a damped pendulum. The time-series ML models accurately capture the dynamics of the pendulum. Moreover, the ratio of the ML timescales extracted based on the sampling frequency and sequence length of the best models correspond to the ratio of the oscillating frequency to the damping frequency of the pendulum. Physically, this finding denotes that the minimum information required by the time-series ML models to accurately capture the dynamics of a system also captures fundamental timescales of the system. 

We extend the demonstration to a decaying SST flow initialized using Taylor–Green vortices. The time series data of the SST involves laminar, transition, and turbulent decay, posing a much more complex nonlinear problem for the time-series ML modeling framework. The LSTM and NODE models are not only accurate for predicting the time evolution of the SST, but also are numerically stable. Moreover, the ratio of the timescales extracted by the ML models is related to the Reynolds number of the flow. We also generalize the timescale extraction framework for various SST, with different Froude and Reynolds numbers. The ML timescales extracted are validated with the relative change in dynamics of the SST according to Froude and Reynolds numbers.

There are a couple of avenues to extend the current framework. First, while the LSTM model has shown great capability to model the dynamics of the SST, there are better time-series models in literature, e.g. GRU <cit.>, which have been shown to be more accurate than LSTM. The current work is a preliminary effort to demonstrate that such time-series models can indeed be used to formulate interpretable models for SST dynamics. There is a clear avenue to extend the analysis for more complex models. Second, the current analysis is not generalized for predicting various parameter ranges of flow using a single model. This requires further investigation and potentially, a framework to incorporate the variation in the nature of the dynamics. Finally, while the prediction of the models are deterministic, the behavior of SST in nature warrants uncertainty quantification of the results of the models. The unique opportunity provided by GPU enabled high-performance computing platforms to perform ML inference at scale can further enable probabilistic approaches for modeling. Such complimentary advancements to the current framework can enable modeling of flows with higher complexity and spatial dimensions.












§ ACKNOWLEDGEMENTS

This research was supported by the Office of Naval Research via grant N00014-19-1-2152. High performance computing resources were provided via the U.S. Department of Energy(DOE) INCITE program by the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. Computing resources were also provided through the U.S. Department of Defense High Performance Computing Modernization Program. 





§ DATA AND CODE AVAILABILITY

Correspondence and requests for material should be addressed to M.G.M. 
Our algorithms will be made freely available in https://github.com/muralikrishnangm/sst-rans-ml.gitLink after peer-review.




§ AUTHOR CONTRIBUTIONS

M.G.M, S.M.deB.K, and J.J.R initiated and designed the project. S.M.deB.K curated the stably stratified turbulence direct numerical simulation data. M.G.M., D.L., and A.D.S. designed and performed the machine learning analysis with feedback from A.K. and W.H.B.. M.G.M., D.L., A.D.S., and S.M.deB.K wrote the manuscript with feedback from all the coauthors.



§ AUTHOR DECLARATION

The authors declare no competing interests.












unsrtnat







