Introduction

Turbulent flows are ubiquitous in many geophysical systems, including atmospheric and oceanic circulations, and play an important role, for example, greatly enhancing mixing and transport. Direct numerical simulation (DNS) of high-dimensional turbulent flows often becomes computationally intractable. Therefore, numerical simulations of most geophysical turbulent flows cannot resolve all the relevant scales (Fox-Kemper et al., 2019;Palmer, 2001;Schneider, Teixeira, et al., 2017). Large-eddy simulation (LES) is a practical approach to balance computational cost and accuracy: the large scales of the flow are explicitly resolved, while the effects of the small-scale features which cannot be resolved by the given grid resolution, called subgrid-scale (SGS) features, are parameterized as a function of the resolved flow (Pope, 2000;Sagaut, 2006;Smagorinsky, 1963). However, the performance of the LES models strongly depends on the accuracy of the employed SGS closure. Over years, there have been extensive efforts focused on formulating physics-based and semi-empirical SGS closures using various techniques in many turbulent flows (Meneveau & Katz, 2000;Moser et al., 2021;Pope, 2000;Sagaut, 2006), including geophysical flows (Alexander & Dunkerton, 1999;Anstey & Zanna, 2017;Berner et al., 2017;Cessi, 2008;Gallet & Ferrari, 2020;Herman & Kuang, 2013;Jansen & Held, 2014;Khodkar et al., 2019;O'Kane & Frederiksen, 2008;Sadourny & Basdevant, 1985;Schneider, Teixeira, et al., 2017;Sridhar et al., 2022;Sullivan et al., 1994;Tan et al., 2018;Zanna et al., 2017).

The challenge of modeling SGS closures lies in faithfully representing the two-way interactions between the SGS processes and the resolved large-scale dynamics. There are two general approaches to SGS modeling: (a) functional and (b) structural (Sagaut, 2006). The functional SGS closures are developed by considering the interscale interactions (e.g., energy transfers). This is often achieved by introducing a dissipative term. Hence, functional SGS closures generally take an eddy-viscosity form to mimic the average function of the SGS eddies. Among the first and most-used functional closures is the Smagorinsky model (Smagorinsky, 1963). Later, dynamic formulations of this model were proposed, in which the key coefficient is dynamically adjusted to the local structures of the flow (Chai & Mahesh, 2012;Germano, 1992;Ghosal et al., 1993;Lilly, 1992). Existing functional closures, most of which are the eddy-viscosity type, can be excessively dissipative (Guan et al., 2022;Vreman et al., 1996). Furthermore, they cannot capture the structure of the SGS terms, leading to a low correlation coefficient (CC < 0.5) with the true SGS terms, that is, those diagnosed from the DNS data (Carati et al., 2001;Guan et al., 2022;Moser et al., 2021).

On the contrary, structural closures tend to have much higher CC with the true SGS terms. Structural closures approximate the SGS terms by constructing it from an evaluation of large-scale motions or a formal series expansion. One of the most common structural closures is the nonlinear gradient model (Clark et al., 1979;Leonard, 1975), referred to as NGM hereafter (it is also known as the tensor diffusivity model (Leonard & Winckelmans, 1999)). The NGM can be derived analytically: the SGS term is approximated using a first-order truncated Taylor-series expansion of the SGS stress' convolution integral (details discussed later). However, despite CC > 0.9, LES with NGM closure has been found to be unstable in many studies of two-dimensional (2D) and three-dimensional (3D) turbulence. These instabilities are often attributed to insufficient dissipation and more importantly, to the presence of too-strong backscattering in NGM (S. Chen et al., 2003Chen et al., , 2006;;Fabre & Balarac, 2011;Leonard, 1997Leonard, , 2016;;Liu et al., 1994;Lu & Porté-Agel, 2010;Meneveau & Katz, 2000;Moser et al., 2021;Prakash et al., 2021;Vollant et al., 2016). As a result, while backscattering (basically anti-diffusion or up-gradient flux) is an important process to represent in closure models (Grooms et al., 2015;Guan et al., 2022;Hewitt et al., 2020;Nadiga, 2010;Shutts, 2005), it is ignored in most practical SGS closures in favor of stability (though there has been some new exciting progress in this area; see, e.g., Jansen et al. (2015) and Juricke et al. (2020)). In fact, currently operational climate models do not account for backscattering in their ocean parameterizations (Hewitt et al., 2020). Consequently, a framework for developing SGS closures with the right amount of diffusion and backscattering, that can capture both the structure and function of the SGS terms, has remained elusive (Moser et al., 2021;Pope, 2000;Sagaut, 2006).

Before moving forward, it should be pointed out that while the discussion so far has been focused on closure for geophysical turbulence, many other critical processes in the Earth system (in atmosphere, ocean, land, cryosphere, biosphere and at their interfaces) require parameterizations in Earth system models (Schneider, Jeevanjee, & Socolow, 2021;Stensrud, 2009). Thus, the discussion below and as clarified later, the findings of this paper, are broadly relevant to parameterization efforts in Earth science.

Recently, machine learning (ML) has brought new tools into SGS closure modeling (Balaji, 2021;Brunton et al., 2020;Duraisamy, 2021;Gentine et al., 2021;Schneider, Lan, et al., 2017;Zanna & Bolton, 2021). The strength of ML techniques is their ability to handle high-dimensional data and learn strongly nonlinear relationships. Therefore, ML techniques are attractive tools that might be able to extract more hidden knowledge from data, potentially providing better SGS closures and even new insights into SGS physics. Data-driven SGS closures, for example, based on deep neural networks trained on high-fidelity simulation data such as DNS data, have been developed for canonical geophysical flows such as 2D and quasi-geostrophic turbulence (Bolton & Zanna, 2019;Frezat et al., 2022;Guan et al., 2022Guan et al., , 2023;;Maulik et al., 2018;Pawar et al., 2020;Srinivasan et al., 2023) and oceanic and atmospheric circulations (Beucler et al., 2021;Brenowitz & Bretherton, 2018;Cheng et al., 2022;Guillaumin & Zanna, 2021;Rasp et al., 2018;Yuval & O'Gorman, 2020;X. Zhang et al., 2022). While some of these studies found the learned data-driven SGS closures to lead to stable and accurate LES (Frezat et al., 2022;Guan et al., 2022Guan et al., , 2023;;Yuval & O'Gorman, 2020), a number of major challenges remain (Balaji, 2021;Schneider, Jeevanjee, & Socolow, 2021). Perhaps the most important one is interpretability, which is difficult for neural networks, despite some recent advances in explainable ML for climate-related applications (Clare et al., 2022;Mamalakis et al., 2022), including for SGS modeling (Pahlavan et al., 2024;Subel et al., 2023). The black-box nature of neural network-based closures aside, there are also challenges related to generalizability, computational cost, and even implementation (Balaji, 2021;Chattopadhyay et al., 2020;Guan et al., 2022;Kurz & Beck, 2020;Maulik et al., 2019;Subel et al., 2021;Xie et al., 2019;Zhou et al., 2019), limiting the broad application of such closures in operational climate and weather models, at least for now.

An alternative approach that is rapidly growing in popularity involves using ML techniques that provide interpretable, closed-form equations, for example, using sparse linear regression. The underlying idea of this equationdiscovery approach is that given spatial, temporal, or spatio-temporal data from a system, one can discover the governing (algebraic or differential) equations of that system (Brunton et al., 2016;Y. Chen et al., 2022;Goyal & Benner, 2022;Mojgani et al., 2022;Rudy et al., 2017;Schaeffer, 2017;Schmidt & Lipson, 2009;Schneider, Stuart, & Wu, 2021;Schneider et al., 2020Schneider et al., , 2022;;Udrescu & Tegmark, 2020;S. Zhang & Lin, 2018). Most of the aforementioned studies are focused on discovering the entire governing equations from data, though few recent studies have used this approach to discover SGS closures (see below). This approach has the following advantages over more complex methods such as neural networks in the context of SGS modeling: (a) the learned closure is significantly easier to interpret based on physics (Zanna & Bolton, 2020), (b) the number of required training samples and the training costs are often considerably lower (Brunton et al., 2020;Mojgani et al., 2022Mojgani et al., , 2023)), and (c) the computational cost of implementation in conventional solvers is lower, as the discovered closures often involve traditional operations, for example, gradients and Laplacians (Ross et al., 2023;Udrescu & Tegmark, 2020).

A number of equation-discovery techniques and test cases have been recently employed for structural modeling of the SGS stress. In the first study of its kind, Zanna and Bolton (2020) used relevance vector machine (RVM), a sparsity-promoting Bayesian linear regression technique, with a library of second-order velocity derivatives and their nonlinear combinations, to learn a closed-form closure model for the SGS momentum and buoyancy fluxes from filtered high-resolution simulations of ocean mesoscale turbulence. They found a closure that resembled the NGM, with close connections to earlier physics-based modeling work by Anstey and Zanna (2017). Although, the discovered closure performed well in a priori (offline) tests, it was unstable a posteriori (online), that is, when it was coupled to a low-resolution ocean solver. Following the same general approach, more recently, Ross et al. (2023) proposed a novel equation-discovery approach combining linear regression and genetic programming. This hybrid approach uses genetic programming to discover the structure of the equation followed by linear regression to fine-tune the coefficients. In contrast to methods such as RVM, genetic programming does not require an explicit library of features, instead, it uses a simple set of features and operations, and constructs expressions by successively applying operators and combining expressions. Similarly, in other disciplines, Reissmann et al. (2021) and Li et al. (2021) recently used are the number of grid points on the DNS grid in the horizontal and vertical directions, respectively. L = 6π is the length of the domain in the horizontal direction. Note that the lowest N LES is chosen such that the LES resolution resolves at least 80% of the DNS kinetic energy (Pope, 2000). Filters are only applied along the horizontal direction. (Pope, 2000). Filters are applied in both spatial dimensions for 2D-FHIT.

gene expression programming to discover SGS stress for the Taylor-Green vortex and the 3D isotropic turbulence, respectively. They developed a nonlinear closure consisting of the local strain rate and rotation rate tensors, based on what is known as Pope tensors (Pope, 1975), which will be discussed later. Overall, these more recent studies found that gene expression programming-and genetic programming-based closures often outperform common baselines such as the Smagorinsky and the mixed models when turbulence statistics and flow structures are considered (Li et al., 2021;Reissmann et al., 2021;Ross et al., 2023). Note that there also have been a number of studies focused on equation-discovery for functional modeling, for example, using techniques such as Ensemble Kalman inversion (Schneider, Stuart, & Wu, 2021;Schneider et al., 2020); see the Section 4.

In this study, we build on the work by Zanna and Bolton (2020) and use 2D forced homogeneous isotropic turbulence (2D-FHIT) and Rayleigh-Bénard convection (RBC) to extend and expand their analysis in several directions:

1. We use RVM with an expansive high-order library to discover closures from DNS data for the SGS momentum flux tensor (2D-FHIT and RBC) and the SGS heat flux vector (RBC). 2. We conduct extensive robustness analysis of the discovered closures across a variety of flow configurations, filter types, and filter sizes, and examine the potential effects of numerical errors. 3. Further clarify the connections between the robustly discovered SGS momentum and heat flux closures, and the SGS closures obtained analytically from the truncated Taylor-series expansion of the filter's convolution integral, the NGM (Leonard, 1975). 4. Explain the physical reason for the unstable a posteriori LES with the discovered SGS closures, despite their high a priori accuracy in some metrics (such as CC). 5. Present a decomposition of the SGS tensor to the Leonard, cross, and Reynolds components, showing their relative importance and dependence on the filter type/size. 6. Based on these findings, we present a number of ideas for discovering stable and accurate SGS closures from the data in future work. Note that while we focus on the use of RVM here, our findings and conclusions in Equations 1-6 are applicable to any equation-discovery effort, and not just for SGS momentum and heat fluxes in geophysical turbulence, but for SGS modeling in any nonlinear dynamical system.

This paper is organized as follows. In Section 2, we provide an introduction to the methodology, including the governing equations of test cases (2D-FHIT and RBC), filtering procedure for data and equations, RVM algorithm, and the employed library of the basis functions. Section 3 includes the discussion on the discovered closures, a priori and a posteriori tests, connection with the physics-based closures, and contribution of the Leonard, cross, and Reynolds components. Summary and Discussion are in Section 4.

Models, Methods, and Data

Filtering Procedure

In DNS, the velocity field, u(x, t), is resolved using high spatio-temporal resolutions down to all relevant scales. In LES, a low-pass filtering operation, denoted by (.), is performed on the equations and flow fields. The resulting filtered fields, for example, filtered velocity, u(x,t) , can be adequately resolved using relatively coarse spatiotemporal resolutions: the required grid spacing is proportional to the specified filter width, Δ, which is analogous to the size of the smallest eddies resolved in the LES (Pope, 2000;Sagaut, 2006). Using u(x, t) as an example, the general spatial filtering operation is defined by Sagaut ( 2006)

where * is the convolution operator, and the integration is performed over the entire domain. The specified filter kernel, G, satisfies the normalization condition

∫

Subsequently, any flow field such as velocity can be decomposed into a filtered (resolved) part and SGS (residual) part:

where u′ is the SGS field. While this appears to be analogous to the Reynolds decomposition, an important distinction should be noted: the filtered residual field may not be strictly zero (u′ ≠ 0, thus u ≠ u), depending on the choice of the filter function (Sagaut, 2006). Further details about the filters used in this work (Gaussian, box, Gaussian + box, and sharp-spectral) are given in Appendix A.

Two-Dimensional Forced Homogeneous Isotropic Turbulence (2D-FHIT)

We consider 2D-FHIT as the first test case. This canonical flow has been extensively used for testing novel physics-based and ML-based SGS closures for geophysical turbulence in the past decades (Boffetta & Ecke, 2012;Chandler & Kerswell, 2013;Guan et al., 2022;Tabeling, 2002;Thuburn et al., 2014;Vallis, 2017;Verkley et al., 2019). The dimensionless continuity and momentum equations for 2D-FHIT in (x, y) spatial dimensions are:

where u = (u,v) is the velocity, p is the pressure, F represents a time-constant external forcing, R is the Rayleigh drag, and Re is the Reynolds number. The domain is doubly periodic with length L = 2π.

The equations for LES are obtained by applying a homogeneous 2D filter (Equation 1) to Equations 4 and 5. The filtered continuity and momentum equations are:

where τ is the SGS stress tensor:

A closure model is needed to represent τ xx , τ xy (=τ yx ), and τ yy , in terms of the resolved flow (u,v,p) . However, currently, this is not possible just using the first principles due to the presence of the u 2 ,uv, and v 2 terms.

We study three cases of 2D-FHIT (Table 1), creating a variety of flows that differ in dominant length scales and energy/enstrophy cascade regimes. For DNS, as discussed in Appendix B, Equations 4 and 5 are numerically solved at high spatio-temporal resolutions using a Fourier-Fourier pseudo-spectral solver. For the LES, the same solver at lower spatio-temporal resolution is used (Appendix B).

Note. For Cases R1-R3 and different filter types, the mean and standard deviation of the discovered coefficients over different N LES are reported (see Table 2). The average CC of the closure for each element of J is shown in parentheses. The last column shows the analytically derived coefficients for the NGM (see Section 3.2).

Turbulent Rayleigh-Bénard Convection (RBC)

As our second test case, we use 2D turbulent RBC, a widely used canonical flow for buoyancy-driven turbulence (Chillà & Schumacher, 2012;Dabbagh et al., 2017;Hassanzadeh et al., 2014;Kooloth et al., 2021;Lappa, 2009;Sondak et al., 2015), which in addition to the SGS (momentum) stress, requires closure modeling of the SGS heat flux (Pandey et al., 2022;Peng & Davidson, 2002;Wang et al., 2008). Under the Oberbeck-Boussinesq approximation, the dimensionless governing equations for the flow between horizontal walls at fixed temperatures (the bottom wall being warmer than the top) in (x, z) spatial dimensions are:

where v = (u,w) is the velocity, θ is the temperature (T ) departure from the conduction state, ẑ is the unit vector in the vertical direction, and Ra and Pr are the Rayleigh and Prandtl numbers, respectively. The domain is periodic in the horizontal direction with horizontal domain size L = 6π and vertical length of 1. No-slip boundary conditions are applied at the walls. We use three cases of turbulent RBC (Table 2) in which the Ra and Pr are varied.

To properly resolve the thin boundary layers in turbulent RBC, a pseudo-spectral solver with (non-uniform) Chebyshev collocation points in the vertical direction is used. However, filtering variables on a non-uniform grid can cause major errors in the diagnosed SGS terms, because the filters will not commute with spatial derivatives (Yalla et al., 2021). As a result, following the common practice for LES, we only filter the equations in the horizontal direction, where (uniform) Fourier collocation points are used. The LES equations obtained by applying a 1D filter along the horizontal direction, x, to Equations 9-11 are:

where τ is the SGS (momentum) stress tensor Note. The CC of P τ for both forward transfer (>0) and backscatter (<0) of SGS energy is "undefined" since P NGM τ = 0 everywhere for 2D-FHIT (in general, P FDNS τ ≠ 0). On the contrary, the forward transfer and backscatter of SGS enstrophy are captured well by the NGM. The Gaussian filter is used in FDNS.

and J is the SGS heat flux vector

Here, in addition to τ, J needs a closure model too.

For DNS, as discussed in Appendix C, Equations 9-11 are numerically solved at high spatio-temporal resolutions using a Fourier-Chebyshev pseudospectral solver. For LES, the same solver with lower spatial resolution is used (Appendix C).

Filtered Direct Numerical Simulation (FDNS) Data

It should be highlighted that in this study with two canonical test cases, we consider DNS data as the "truth", and use filtered DNS (FDNS) data to discover the closures. However, in reality, performing DNS for many geophysical flows is computationally prohibitive. In such cases, highresolution LES that adequately resolves the process of interest (e.g., ocean eddies, gravity waves, etc.) is often used as the truth to train the ML algorithms for SGS modeling (Shen et al., 2022;Sun et al., 2023;Yuval & O'Gorman, 2020;Zanna & Bolton, 2021).

Here, we compute FDNS variables on the LES grids, which are 4-64 times coarser than the DNS grid in both spatial dimensions for 2D-FHIT and one spatial dimension for RBC (see Tables 1 and 2). More specifically, we first apply the respective filter's transfer function (Tables A1 and A2) to the DNS data, and then coarse-grain the results onto the LES grid. Note that following some of the recent papers (Grooms et al., 2021;Guan et al., 2022), we define "filtering" as an operation that removes the small scales but keeps the grid resolution (e.g., DNS), and "coarse-graining" as an operation that changes the grid size, for example, from the DNS resolution to LES resolution. Note that τ and J in Equations 7, 13, and 14 need to be on the LES grid.

The filtering and coarse-graining are performed following Sagaut (2006) and Guan et al. (2022). Briefly, using the velocity u(x DNS ,t) as an example, and denoting the DNS grid and wavenumber as x DNS and k DNS , we first transform the DNS velocity into the spectral space û(k DNS ,t), where (.) means Fourier transformed. This is followed by applying the filter in the spectral space:

Here, Ĝ(k DNS ) can be any of the transfer functions listed in Tables A1 and A2, and ⊙ is the Hadamard (elementwise) multiplication. After the filtering operation, coarse-graining is performed to transform the filtered variable from the DNS to the LES grid. In this study, we perform coarse-graining in spectral space with cutoff, k c = π/ Δ LES , which for example, in 2D, yields

It is important to highlight that the filtering operation can be reversible under certain conditions: The DNS data can be recovered (via deconvolution) from the filtered data (still on the DNS grid) if the filter is a Reynolds operation, not a projection (Sagaut, 2006). This is the case for Gaussian, box, and Gaussian + box filters, whose transfer function only attenuates the signal. However, for the sharp-spectral cutoff filter (which is a projection

Table 7 The Average Correlation Coefficient (CC) Between Inter-Scale Kinetic Energy Transfer (P τ ) or Enstrophy Transfer (P Z ) or Potential Energy Transfer (P J ) Patterns of the SGS Fluxes From FDNS and From NGM Closure (Equations 25 and 26) for Cases R1-R3 and Different N LES Cases N LES = 128 N LES = 256 N LES = 512 CC for P τ (P τ > 0, P τ < 0) R1 0.94(0.96, 0.85) 0.99(0.99, 0.98) -R2 0.97(0.81, 0.98) 0.98(0.91, 0.98) 0.99(0.97, 1.00) R3 0.79(0.81, 0.74) 0.88(0.92, 0.81) 0.96(0.97, 0.93) CC for P Z (P Z > 0, P Z < 0) R1 1.00(1.00, 0.99) 1.00(1.00, 1.00) -R2 0.99(0.96, 1.00) 0.99(0.98, 1.00) 1.00(1.00, 1.00) R3 0.96(0.95, 0.96) 0.99(0.99, 0.98) 1.00(1.00, 0.99) CC for P J (P J > 0, P J < 0) R1 0.89(0.89, 0.15) 0.97(0.97, 0.46) -R2 0.76(0.75, 0.65) 0.91(0.91, 0.63) 0.98(0.98, 0.76) R3 0.77(0.75, 0.40) 0.87(0.86, 0.39) 0.94(0.94, 0.44)

Note. Note that for RBC, filtering is conducted in only one direction (x), therefore, P τ is not identically zero. Here, the forward transfer and backscatter of SGS kinetic energy and enstrophy are overall captured well, specially as N LES increases. However, the backscatter of SGS potential energy is not well captured, specially at low LES resolutions. The Gaussian filter is used in FDNS. See Appendix F for the definition of P J .

operator), and for coarse-graining, information is lost beyond the cutoff wavenumber k c and therefore data recovery is limited only up to k c . This is further discussed in Section 3.2.

Hereafter, for brevity, we use the term "filtered" (still denoted by .

-) to mean "filtered" and then "coarse-grained."

The spectrum of a filtered and coarse-grained variable can be found in Figure A1b in Appendix A.

Figure 1 shows the effects of filtering on the vorticity and temperature fields for 2D-FHIT and RBC, illustrating that the small-scale structures of ω and T are removed due to filtering and the fields are smoothed out.

The Equation-Discovery Method

In this study, we employ the RVM (Tipping, 2001) to discover closed-form closures for each element of the τ tensor and J vector from the FDNS data. RVM is a sparsity-promoting Bayesian (linear) regression technique that has shown promise in applications involving dynamical systems (Mojgani et al., 2022;Zanna & Bolton, 2020; S. Zhang & Lin, 2018). RVM relies on a pre-specified library of basis functions Φ; each column of this matrix is a basis, for example, a linear or nonlinear combination of relevant variables such as velocity and temperature and/or their derivatives. The library should be expressive enough so that s, a vectorized snapshot of a element of any τ or J, could be completely represented as

The vector of regression weights, c, is computed by minimizing the mean-squared error (MSE)

where vector S consists of n samples of s stacked together. RVM assumes Gaussian prior distributions for each weight, and the width of the Gaussian posterior provides a measure of the weight's uncertainty. Sparsity is enforced via an iterative process: basis functions whose weights' uncertainties exceed a pre-specified hyperparameter (threshold), α, are removed (pruned), and Equation 20 is minimized again. The iterations stop when all the remaining basis functions have uncertainties smaller than α. Larger α results in a lower MSE but more terms in the discovered model (see below).

A critical step in using RVM (and most equation-discovery methods) is the choice of the library. Here, we have chosen the following libraries. For momentum stress, we use

where A,B = u or v (2D-FHIT) and C, D = u or w (RBC). Note that experiments with including θ in D yield the same results. For heat flux, we use

where A = u, w, or θ (RBC). These libraries are expansive, with integers 0 ≤ q ≤ 8 and 0 ≤ p ≤ 2, though the total derivative order is limited to 8th (there are a total of 546 and 614 terms in the libraries used for momentum and heat fluxes, respectively). The form of these libraries is motivated by the Galilean-invariant property of the SGS terms, and by past studies. For example, these libraries include Pope's tensors (Pope, 1975), which have been used in physics-based (Anstey & Zanna, 2017;Gatski & Speziale, 1993;Jongen & Gatski, 1998;Lund & Novikov, 1993) and equation-discovery (Li et al., 2021;Reissmann et al., 2021;Ross et al., 2023) approaches in the past (and include the structure of the Smagorinsky model; see below). Our library also includes the basis functions used by Zanna and Bolton (2020).

Note that all calculations for the libraries (and any computation in this work) is performed using the same spectral methods used for DNS and LES.

We have found it useful for interpretability of the outcome and improving the robustness of the algorithm to remove redundant terms using the continuity equation (e.g., using ∂v/∂y = -∂u/∂x, ∂ 2 v/∂y∂x = -∂ 2 u/∂x 2 , etc.). Also, we have found it essential to normalize each basis in Φ to have a zero mean and a unit variance, because the amplitude of higher-order derivatives can be much larger than that of the lower-order ones.

Like any method, equation discovery using RVM has a number of strengths and weaknesses:

1. It is data efficient (Mojgani et al., 2022;Zanna & Bolton, 2020). For example, here, we report the results with n = 100 FDNS samples, but even with n = 1, the results remain practically the same. 2. It is more robust, in terms of convergence, compared to similar sparsity-promoting techniques (Zanna & Bolton, 2020;S. Zhang & Lin, 2018). 3. A pre-specified library is needed and it is assumed that the true answer (e.g., the SGS stress) can be represented as a linear combination of the chosen basis functions. 4. The pre-specified hyper-parameter α determines how parsimonious the discovered model is. Decreasing α leads to a smaller (likely, more interpretable) model at the expense of increasing the MSE. Here, we follow the model-selection literature (Mangan et al., 2017;Mojgani et al., 2022) and objectively choose α using the Lcurve, as shown later. 5. The answer can depend on the choice of the loss function. The RVM's MSE loss (Equation 20) is strictly following the principle of structural modeling, matching the flux between the FDNS and the discovered model.

Note that the above strengths ( 1) and ( 2) are highly desirable while these weaknesses (3)-( 5) are common among many equation-discovery methods, although techniques such as genetic programming and gene expression programming can address (3) and ( 5), for example, using an evolving library. We will further discuss (3)-( 5) in Section 4.

Results

In this section, we present and discuss the discovered closures, and analyze them a priori (offline) and a posteriori (online, coupled with LES). We then uncover the connections between the discovered closure and the NGM. For all results presented here, we use n = 100 FDNS samples from a training set and 20 FDNS samples from an independent testing set.

The Discovered Closures for SGS Momentum and Heat Fluxes

For each of the six cases in Tables 1 and 2, we separately discover closures for three elements of the SGS stress tensor, that is, τ xx , τ xy = τ yx , and τ yy for 2D-FHIT, and τ xx , τ xz = τ zx , and τ zz for RBC. Additionally, we discover two elements of the SGS heat flux vector, that is, J x and J z for RBC. We discover individual closures for four filter types: Gaussian, box, sharp-spectral, and Gaussian + box. The first three are common filter types, while the last one is motivated by a few recent studies (Guillaumin & Zanna, 2021;Zanna & Bolton, 2020). We also examine several filter sizes, Δ (see Tables 1 and 2), and the effect of varying α, which as mentioned earlier, is a key hyperparameter in RVM.

We analyze the a priori performance of the discovered closures using the most commonly used metric: the average of CCs for testing samples (Guan et al., 2023;Maulik et al., 2019;Sagaut, 2006). For each element of τ or J, denoted below by τ for convenience, the CC for each testing sample is calculated between 2D patterns of τ from FDNS and τ predicted by the RVM-discovered closure for the corresponding filtered flow variables (e.g., u, v etc.):

Journal of Advances in Modeling Earth Systems

where 〈⋅〉 is domain averaging. The same equation is also used for computing CC values of 2D patterns of interscale energy or enstrophy transfer, P (defined later).

As a representative example of the findings, Figures 2a and 2b shows the averaged CC for τ yy (K1-K3) and J x (R1-R3) as α is increased. Figures 2c and 2d presents the number of terms in the discovered closures. With small α, the discovery is unsuccessful (CC = 0; zero terms). However, as α is further increased, for all cases, CC abruptly jumps to above 0.8 -0.9 with 1 and 2 discovered terms, and then gradually converges to 1 but with exponentially growing number of terms in the discovered closure. The CC-α relationship forms an "L-curve." The elbow of this curve indicates the α that balances accuracy and model size, and is extensively used in the model-selection and 2). The Gaussian filter is applied in both cases.

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874 equation-discovery literature to objectively choose α (Calvetti et al., 2000;Goyal & Benner, 2022;Lawson & Hanson, 1995;Mangan et al., 2017;Mojgani et al., 2022). Examining all cases with other filter sizes and filter types reveals the same behavior as shown in Figure 2, with the exception of the sharp-spectral filter as shown in Figure D1. For this filter, the discovery is unsuccessful, leading to low CC and non-robust results; we will explain the reason for this failure later in this section. ) is used, but the same behavior is observed with any other N LES and filter type (except for the sharp-spectral, see the text). In general, for small α (<1), no closure is discovered (CC = 0, zero term). With increasing α, the CC converges to ∼1 (a more accurate a priori closure) but at the expense of a larger closure with many more terms (note the logarithmic scale of the y axes in panels (c) and (d)). However, the CC-α relationship forms an "L-curve," whose elbow indicates the α that balances accuracy and model size (see the text).

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874

We use the L-curve to determine the optimal α. In 2D-FHIT, there are two kinks in the curve around the elbow, corresponding to the discovery of closures with 1 and 2 terms, respectively (Figures 2a and 2c).

Given the robust and asymptotic behavior in α after the second kink, we use the corresponding α to identify the discovered closure (see the black circles). We find that consistently, across Cases K1-K3, filter types, and filter sizes, this closure is of the form τ = [ τ xx τ xy τ yx τ yy ] = Δ 2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ a xx ( ∂u ∂x ) 2 + b xx ( ∂u ∂y ) 2 a xy ∂u ∂x ∂v ∂x + b xy ∂u ∂y ∂v ∂y a xy ∂u ∂x ∂v ∂x + b xy ∂u ∂y ∂v ∂y a yy ( ∂v ∂x ) 2 + b yy ( ∂v ∂y

where a xx , b xx , a xy , b xy , a yy , and b yy are the discovered coefficients (Δ 2 is factored out to further highlight the independence of these coefficients from the filter size). Table 3 shows that these six coefficients are the same, and the same for Cases K1-K3, although they can depend on the filter type. This table also shows the average CC values of the discovered closure, which are around 0.99, demonstrating the accurate prediction of each element of the stress tensor and the excellent a priori (offline) performance of the discovered closure for a broad range of LES resolutions.

Following the same approach, we discover basically the same closure for τ in RBC Journal of Advances in Modeling Earth Systems

where, as before, d xx , d xz , and d zz are the coefficients with Δ 2 factored out. Note that Equation 25 is the same as Equation 24, except that here, there is one term rather than two in each element of the tensor, which is a result of filtering (in RBC) performed only in the horizontal, x, direction. As before, Table 4 shows that these d coefficients are the same, and the same for Cases R1-R3, though varying with filter type. Like before, the discovered closure has fairly high CC values.

Again, following the same approach, we determine the optimal α for discovering the closure of J. In Figure 2b, Case R1 has a clear elbow while Cases R2-R3 have two kinks around the elbow. Examining all cases and the number of discovered terms (Figure 2d), we find that the single-term closures discovered at the first kink (circled) provide consistent and robust results. The first kink is observed for all the Cases R1-R3. In contrast, the second kink is observed only in Cases R2-R3. The closure at the first kink is

where d x and d z are the discovered coefficients with Δ 2 factored out. Table 5 shows that these d coefficients are the same, and the same for Cases R1-R3, but varying with filter type. As before, the discovered closure has a good a priori performance.

To summarize the findings, Equations 24-26 and Tables 3-5 show that 1. Closures of the same form are robustly discovered for τ in two vastly different systems, 2D-FHIT and RBC.

Even the closure for J overall has the same form, consisting of the products of the first-order derivatives of the variables involved in the nonlinearity of the SGS term.

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874

2. Not just the form, but even the coefficients of the terms in the closures, are consistently the same as parameters such as Re, forcing wavenumber, Ra, or Pr are changed in Cases K1-K3 and R1-R3, leading to different dynamics. The coefficients are independent of the fluid and even the flow properties. 3. The form of the closures is independent of the filter type unless the sharp-spectral filter is used. The coefficients, once normalized by Δ 2 , are independent of filter size, but depend on filter type. 4. The discovered closures have outstanding a priori performance, often with CC > 0.95 and even as high as 0.99.

It should be noted that the CCs reported in these tables are averaged over a broad range of N LES . The values of CC are higher for larger N LES , that is, smaller Δ.

The Nonlinear Gradient Model (NGM): Taylor-Series Expansion of the SGS Term

A closer examination of Equation 24 reveals that this closure is indeed the NGM (this includes both the form and the coefficients, within the uncertainty range). This connection was already pointed out by Zanna and Bolton (2020), although the implications and findings such as 1-4 mentioned in the previous subsection were not further discussed in their short letter.

First, let's briefly review the NGM (Clark et al., 1979;Leonard, 1975;Sagaut, 2006). As a simple illustration of the idea behind this model, we have reproduced the derivation of NGM in Appendix E using a 1D arbitrary field, a (x), for the reader's convenience. Taylor-series expansion of a(xr x ) around a(x) (Equation E2) simplifies the convolution integral of the filtering operation (Equation E1) such that a(x) can be written in terms of a(x) and its derivatives, with coefficients that depend only on the moments of the filter's kernel, G (Equation E4). Using u 2 and u 2 as a(x), we eventually arrive at an analytically derived closure for τ xx with error O(Δ 4 ) (Equation E12). In 2D with filtering applied in both directions (like our 2D-FHIT), the NGM is (Sagaut, 2006)

∂y ) 2 ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y ( ∂v ∂x ) 2 + ( ∂v ∂y ) 2

where c τ depends on the filter's kernel. Similarly, for the 2D RBC with filtering only in the x direction, the NGM is

As emphasized in Appendix E, there is nothing specific to momentum flux or even turbulence (or even physical systems) in the derivation of NGM. In fact, for the filtered quadratic nonlinearity of any two arbitrary variables, one arrives at the same expression with O(Δ 4 ) accuracy. For example, following this derivation, for the SGS heat flux, we obtain

where like c τ and d τ , d J only depends on the filter's kernel. Note that we are referring to NGM as a general class of closure models, representing the leading term of the Taylor-series expansion of the SGS term (whether it is for momentum or heat or any other flux).

Computing c τ , d τ , d J for each of the filter types used in this study, we confirm that the discovered closures for the SGS stress are basically the NGM (Equations 27 and 28), and in the case of the SGS heat flux, an NGM (Equation 29) closure (see Tables 345).

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874

Based on the above analyses, we can now explain the findings (1)-( 4) in Section 3.1. Closures of the same structure are robustly discovered for both SGS momentum and heat fluxes in two vastly different turbulent flows (and independent of parameters such as Re, Ra, Pr, and forcing) because the first term in the Taylor-series expansion dominates the SGS flux. As a result, in equation-discovery using common loss functions such as MSE and evaluation metrics such as CC, which aim at closely matching τ or J, NGM or NGM-like closures are discovered (if the library is expansive enough to include all the relevant terms). We emphasize that this would be the case with discovering the representation of the filtered nonlinearity of any two arbitrary variables. As already observed, the coefficients of the discovered closure become even closer to those of NGM as Δ is decreased (thus reducing potential contributions from the truncated O(Δ 4 ) terms).

The connection to the analytical derivation also explains why the coefficients in the discovered models are independent of the fluid or even the flow properties (Ra, Re, Pr) and only depend on the filter size (Δ) and filter type.

For the Gaussian and box filters we obtain c τ = d τ = d J = 1/12: this is because the parameters of the Gaussian filter are chosen such that Gaussian and box filters' kernels have the same second moment (see Equation E5) for the definition of moment) (Pope, 2000). For the Gaussian + box filter, the coefficients are 1/6 because the kernel of this filter is convolution of the Gaussian and box filter kernels. For the sharp-spectral filter, the moments are indefinite, this is why there is no NGM discovery with this filter (and we will discuss later why the equation discovery fails altogether). Note that coarse-graining done here via cutoff in the spectral space does not change c τ , d τ , and d J ; however, if coarse-graining is done by other techniques such as box averaging, then the coefficients might change (note that the NGM coefficient for Gaussian + box filter is twice the coefficient of either filter, as the coefficients are additive; see Tables 345).

In short, one can explain the effects of different filter kernels and coarse-graining strategies on the discovered closures following the analytically derivable NGM (see Appendix E and Sagaut ( 2006)).

Note that as mentioned earlier, reconstructing SGS fluxes from the filtered solution is a solvable problem as long as the filtering operation is reversible. This is the case for filters with compact support, which include Gaussian, box, and Gaussian + box filters, but not sharp-spectral (cutoff) filters. However, once coarse-graining is applied in addition to any of these filters, the full operation is no longer reversible. Still, the part of the spectrum above the Nyquist frequency can be recovered. It should be noted that NGM is an approximation of the deconvolution of SGS fluxes (Sagaut, 2006). The deconvolution procedure can be done only for invertible filters, that is, nonprojective filters. Projective filters, such as sharp-spectral cutoffs, induce an irreversible loss of information, rendering the original data unrecoverable. This is another way of understanding why the NGM exists for Gaussian and box filters (and their combinations) but not for sharp-spectral cutoff filters.

An important implication of the above findings and discussions is that the discovered closure may not be unique and can depend on the filtering and coarse-graining procedure: it depends on the filter type (and up to a factor, on the filter size). This is not a problem of equation-discovery; in fact, the SGS fluxes diagnosed from FDNS are not unique and depend on the filtering and coarse-graining procedure (this is further shown in Figure 4 and discussed at the end of this section). This has implications not just for equation-discovery, but more broadly, for the ongoing efforts on learning SGS parameterizations for various processes from high-fidelity data using ML. See Sun et al. (2023) for extensive discussions about this issue focused on the data-driven SGS modeling of atmospheric gravity waves.

The next key question is about the accuracy and stability of LES of the 2D-FHIT and RBC with the NGM closures, τ NGM and J NGM . However, before discussing the a posteriori (online) performance of NGM closures, we address one more issue, and that is about any potential influence from numerical calculations in our equation discovery.

Effects of Numerical Discretization

The appearance of gradients of velocity (or temperature) in Equations 27-29 might suggest to some that the discovered equations represent the truncated terms of finite difference/volume discretization schemes (the methods used in Zanna and Bolton (2020)). The discussions in their paper and the comprehensive analyses here should leave no ambiguity that Equations 27-29 represent the physics of the SGS fluxes, rather than any numerical error. Still, we wish to discuss a few more points here, as numerical errors from finite difference/volume discretizations or from aliasing (in spectral calculations) can certainly contaminate equation discovery.

All numerical calculations in this study are performed using Fourier and Chebyshev spectral methods. Moreover, we have repeated our calculations of the SGS fluxes and of the basis functions in the library after de-aliasing based on the 2/3 rule (Orszag, 1971). Furthermore, we have repeated the discovery on fluxes that are only filtered but not coarse-grained (thus they remain on the high-resolution DNS grid). The outcomes of all these experiments are Equations 27-29, demonstrating that the discovered closures do not contain any contributions from numerical errors.

A Posteriori (Online) Tests and Inter-Scale Energy/Enstrophy Transfer

For all 6 cases and all tested N LES , the LES runs with NGM closures are unstable: High-wave number features appear and the simulations eventually blow up (not shown). This is consistent with the findings of Zanna and Bolton (2020), who only found stable LES once the SGS momentum fluxes predicted by the discovered closure were attenuated (also see recent work by Perezhogin et al. (2023) who showed that further filtering the output of NGM can lead to promising online results in an ocean model). More generally, this is also consistent with extensive studies in the 1990s (though mainly focused on 3D turbulence), which found that LES with the NGM closure is unstable (Borue & Orszag, 1998;S. Chen et al., 2003S. Chen et al., , 2006;;Leonard, 1997;Liu et al., 1994;Meneveau & Katz, 2000;Pope, 2000;Vreman et al., 1997). The exact reason(s) for the instabilities remain unclear but these studies found that in general, in 3D turbulence, NGM has insufficient dissipation and/or too much backscattering; see, for example, the discussions in Leonard (1997Leonard ( , 2016) ) and Sagaut (2006). Over the years, several methods, including positive clipping (removing backscattering), regularization, modulation, and the dynamic procedure, have been employed to overcome this instability problem and enhance the NGM's performance in a posteriori tests (Balarac et al., 2013;Khani & Dawson, 2023;Khani & Porté-Agel, 2017, 2022;Liu et al., 1994;Prakash et al., 2022). Later studies focused more on eddy-viscosity closures, or on NGM combined with eddy-viscosity, the so-called mixed models (Balarac et al., 2013;Cottet, 1996;Winckelmans et al., 1998). Such versions of NGM have been used in some geophysical flows, for example, for atmospheric boundary layer (Khani & Porté-Agel, 2017, 2022;Khani & Waite, 2020;Lu & Porté-Agel, 2010, 2014) and oceanography (Khani & Dawson, 2023).

In 2D turbulence with filtering done in both directions, such as our 2D-FHIT cases, the NGM has a clear major shortcoming: it cannot capture any energy transfer between the subgrid and resolved scales, despite capturing the enstrophy transfer well (S. Chen et al., 2003Chen et al., , 2006;;Nadiga, 2008). To further explore this issue, first note that the rate of kinetic energy transfer between the resolved and subgrid scales, P τ , is (Pope, 2000)

where summation over repeated indices is implied. S and τ r are the 2D filtered rate of strain tensor and the anisotropic part of the SGS stress tensor (see Appendix F for details). In homogeneous 2D turbulence, globally (domain averaged), there is no net forward cascade of kinetic energy, which has been used to develop successful parameterizations that are globally energetically consistent (Boffetta & Ecke, 2012;Jansen et al., 2015). However, locally (at a given grid point), there can be both forward transfer and backscatter of energy, which can be important to capture by the SGS closures (Thuburn et al., 2014). In 2D turbulence with filtering done in both directions, using τ NGM in the above equation shows that P NGM τ (x,y,t) is identically zero at every grid point (see Appendix F). This is demonstrated numerically in Table 6, which also shows that NGM captures both forward transfer and backscatter of SGS enstrophy fairly well (CC > 0.95). Therefore, despite the high CC of τ NGM with τ FDNS , and even a fairly accurate inter-scale enstrophy transfer, NGM cannot capture any inter-scale energy transfer, indicating a major failure from a functional modeling perspective (note that in this context, "inter-scale" means between the resolved and subgrid scales). A physical/mathematical interpretation of this failure is that while NGM reproduces the structure of τ remarkably well, it does not at all capture the correlations between the τ and S tensors, for example, the angles between their principle directions (Leonard, 2016). This inability to represent any inter-scale energy transfer is likely the reason for the instabilities of LES with NGM closure in Cases K1-K3 (and generally, in 2D turbulence). But how about in RBC? In Cases R1-R3, filtering is conducted only in the horizontal direction, and as a result, P NGM τ is not identically zero. In fact, in these cases, the forward transfer and backscatter of both kinetic energy and enstrophy are captured fairly well by NGM, with CC often above 0.95 (Table 7). However, a deeper examination shows that the backscatter (anti-diffusion) of interscale SGS potential energy, measured as P J (see Appendix F), may not be captured well, specially at low N LES Journal of Advances in Modeling Earth Systems 10.1029/2023MS003874 (Table 7). Poor representation of backscattering can certainly lead to instabilities, as for example, shown by Guan et al. (2022) for 2D turbulence.

To further explore other potential shortcomings of NGM, we have also examined the spectra of elements of τ NGM and J NGM in comparison to those from FDNS (Figure 3). This analysis shows that the spectra of SGS momentum and heat fluxes are captured well across wavenumbers, even at high wavenumbers, indicating that NGM performs well in this a priori (offline) metric.

To summarize the above analyses: we find all LES with NGM closures for 2D-FHIT and RBC cases to become unstable even at high LES resolutions. Understanding the reason(s) for this poor a posteriori (online) performance is essential to make further progress. Examining a few functional and structural metrics beyond CC of SGS fluxes (e.g., inter-scale energy/enstrophy transfers, spectra) point to only one major shortcoming that is relevant to 2D-FHIT (and any 2D turbulent flow): NGM cannot capture any inter-scale kinetic energy transfer, which is likely the reason for the instabilities. This is not an issue in RBC, for which we only identify one shortcoming, and that is the poor representation of backscatter (anti-diffusion) of potential energy, specially at low LES resolution. These findings indicate that the poor a posteriori (online) performance of NGM might have different causes in different flows and requires more extensive investigations.

Before discussing ideas for addressing these challenges in future work, we present more analyses in two areas: a closer examination of the physics included in the library (Section 3.4) and the decomposition of the SGS fluxes and the sensitivity of the diagnosed fluxes to the filter type/size (Section 3.5).

A Physics-Guided Library: Pope Tensors

In Section 3.1, we consider an expansive library of basis functions combining the low-and high-order derivatives and polynomials of velocity and temperature. Under certain assumptions, smaller but physics-informed libraries can be devised. For example, Boussinesq (1877) hypothesized that for a nearly homogeneous, incompressible, high-Re flow, the anisotropic SGS stress τ r (Equation F2) is only a function of the filtered rates of strain S (Equation F1) and rotation Ω (Equation 32) tensors:

Note that in Equations 7 and 13 and in general, only ∇ ⋅ τ r has to be parameterized as the rest of ∇ ⋅ τ can be absorbed into ∇p (Sagaut, 2006). Owing to Cayley-Hamilton theorem (Gantmakher, 2000), τ r can be represented as a linear combination of a finite number of tensors, the so-called Pope tensors (Pope, 1975). In 2D, there are only three Pope tensors Z, thus

The three Pope's tensors are Z (0) = I, Z (1) = S, and

∂x ) 2 2( ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y ) 2( ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y ) -( ∂u ∂y ) 2

Journal of Advances in Modeling Earth Systems 10.1029/2023MS003874

which is related to the anisotropic part of the NGM stress. In fact, τ NGM-r = -Δ 2 Z (2) ⁄12 (see Equation F3). Note that this is also the physics-based closure derived in Anstey and Zanna (2017). Coefficients ζ (n) are functions of invariants I 1 = tr(S 2 ) and I 2 = tr(Ω 2 ). The standard Smagorinsky model is ζ (1) (I 1 )Z 1 .

Our expansive library, described in Equations 21 and 22, includes the individual terms to discover Z (n) (n = 0, 1, 2); however, we have always found the NGM stress, τ NGM . To see whether the results will change with a discovery only done on the anisotropic SGS stress, τ r , and with a smaller library that only has the terms relevant to the Pope tensors, we have conducted more experiments with three libraries for Cases K1-K3. The first library only includes the three Pope tensors {Z (0) , Z (1) , Z (2) }, the second library only includes the six non-zero elements of these tensors, and the third library only includes the eight terms that compromise these six:

( ∂v ∂x ) 2 , ∂u ∂x ∂v ∂x , ∂u ∂y ∂v ∂y }.

The RVM with any of these libraries robustly discovers τ NGM-r = -Δ 2 Z (2) ⁄12, without Z (1) (or Z (0) ) showing up (thus, no Smagorinsky/eddy viscosity-like term). Needless to say, LES with this closure is unstable.

The above analyses show the prevalence of NGM: it emerges whether the full or just the anisotropic part of the SGS stress tensor is discovered, and whether an expansive or a small physics-guided library is used.

Decomposition of SGS Fluxes: Leonard, Cross, and Reynolds Stresses

As discussed earlier, whether a closure could be successfully discovered from FDNS data or if the NGM could be derived depend on the choice of the filter. The latter was explained based on the dependence of the derived closure on the moments' of the filter kernel. Furthermore, the coefficients of the discovered closure and the analytically derived coefficients of the NGM depend on the choice of the filters (Tables 345). Here, we further demonstrate the sensitivity of the diagnosed FDNS SGS flux (which is treated as truth in offline/supervised learning data-driven modeling approaches) to the choice of the filter, and then decompose the flux into its three components to gain further insight.

The top row of Figure 4 shows examples of SGS τ in 2D-FHIT diagnosed from the FDNS data using different filter types. It is clear that the diagnosed fluxes are not unique and particularly different between Gaussian/box filters and sharp-spectral filter (similar differences can be seen in SGS momentum and heat fluxes in RBC). This sensitivity, which has important implications for data-driven SGS modeling efforts (Sun et al., 2023), has been known for a long time in the LES community (Leonard, 1975;Sagaut, 2006). The Gaussian and box filters extract fairly similar features, even of almost the same amplitude (which is due to their matched kernels' second moments). The Gaussian + box filter captures similar features but with a factor of ∼2 difference in amplitude (related to the factor of two difference in NGM coefficients). However, the sharp-spectral filter extracts very different features that have much smaller length scales and amplitudes. We speculate that the failure of the discovery with the cutoff filter is a result of the inability of the current library in representing these features (at least in sparsely representing these features, if it exists). We also point out that in Guan et al. (2022Guan et al. ( , 2023)), deep convolutional neural networks (CNNs) could not be successfully trained on FDNS data obtained using sharp-spectral cutoff filters, while high CC and stable/accurate LES runs in different systems were achieved using CNNs trained on FDNS data obtained through the Gaussian filter. Note that Ross et al. (2023) successfully trained CNNs (and performed equation-discovery) using a "smoothed" sharp-spectral filter that had exponential decay at high wavenumbers (rather than a cutoff). These findings further show the importance of how the "true" SGS fluxes are diagnosed for offline/supervised learning.

To see the reason for this difference, we decompose the SGS tensor using u = u + u′. Leonard (1975) introduced a decomposition of the SGS tensor into three components. However, since two of these components were not Galilean-invariant (Speziale, 1985), a Galilean-invariant decomposition was later proposed by Germano (1986):

τ and J of RBC can be decomposed in the same fashion. The most familiar component, the Reynolds stress, represents interactions in the unresolved scales that project onto the resolved scale. The cross stress represents the direct interactions between the unresolved and resolved scales that project onto the resolved scale. The Leonard stress includes interactions between the resolved scales not captured by the low-resolution LES grid. See Leonard (1975), McDonough (2007), and Perezhogin and Glazunov (2023) for more discussions.

The relative importance of these three components in τ and J depends on the filter type and size (and even the flow characteristics). Rows 2-4 of Figure 4 show examples of the Leonard, cross, and Reynolds stress components of τ xy . For Gaussian, box, and Gaussian + box filters, the Leonard stress dominates, followed by cross and then Reynolds stress. However, for sharp-spectral, only the Reynolds stress has coherent structures that more or less resemble the Reynolds stress from Gaussian/box filters. The strong dependence on filter type comes from the fact that for Gaussian box filters, u′ ≠ 0 and u ≠ u, leading to non-zero Leonard and cross stresses.

However, for the sharp-spectral filters, u′ = 0 and u = u. As a result, with de-aliasing applied, the (coarsegrained) Leonard stress may or may not be zero, depending on details of the calculations as discussed next (the cross stress is non-zero). This is because of the nuances in how the cutoff wavenumber for the sharp-spectral kernel is defined for both filtering and coarse-graining; for example, we find the coarse-grained Leonard stress to be non-zero in Figure 4. Here, following the common approach in the literature (Pope, 2000), for filtering (Table A2), we compare the total wavenumber

) with the cutoff wavenumber (k c ). However, for coarse-graining, we have to compare k x and k y with k c (Equation 18). This subtle difference results in the Leonard stress not being zero here. Note that our conclusion about the failure of discovering a robust closure model with the sharp-spectral filter (at least with the current library and algorithm) is not affected by these choices. We have found that we still cannot discover a closure when we repeat the equation-discovery process on the total stress or the Reynolds stress calculated with filtering that also uses Equation 18.

As for the dependency on filter size, as Δ increases, the relative importance of Reynolds stress increases: See Figure 5 for examples from Cases K3 (τ xx ) and R3 (J z ). Finally, note that the relative importance of these three components might depend on the flow itself. For example, in 3 km-resolution regional simulations of the tropics, Sun et al. (2023) found that the vertical (horizontal) flux of the SGS zonal momentum is dominated by the Reynolds (Leonard) stress, which was attributed to the substantial differences of the filtered vertical wind and the filtered zonal or meridional winds.

The above analyses further explain the strong dependency of the diagnosed "true" SGS flux and the discovered closures on the filter type and size. These analyses also show that depending on the filter type/size, the Reynolds stress may not be the only component of the SGS flux that needs to be parameterized. In fact, the Leonard and cross stresses might be even larger and have to be included in the calculation of the total SGS flux and in the closure. Thus, there is no unique way to quantify the SGS fluxes, and these sensitivities have major implications for the "true" SGS flux that is fed into the RVM or any equation-discovery algorithm (and more broadly, any ML algorithm). Needless to say, these sensitivities have major impacts on the choices to make for developing datadriven closures in real-world applications (see item (f) in the next section).

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874

Summary and Discussion

In this work, we have used relevance vector machine (RVM) to discover subgrid-scale (SGS) closures from filtered direct numerical simulation (DNS) data for both the SGS momentum flux tensor (in 2D forced homogeneous isotropic turbulence, 2D-FHIT, and Rayleigh-Bénard convection, RBC) and the SGS heat flux vector (in RBC). The expansive library includes derivatives of velocity (and temperature) up to 8th order (calculated using spectral methods) and their quadratic combinations. We have conducted extensive robustness analysis of the discovered closures across a variety of flow configurations (changing Re, Ra, Pr, and the forcing wavenumber), filter types (Gaussian, box, Gaussian + box, and sharp-spectral cutoff), and filter sizes.

Based on these analyses, except for when the sharp-spectral filter is used (see below), we have robustly discovered the same closure for the SGS stress in 2D-FHIT and RBC. We have further shown that this closure model is in fact the NGM, which can be derived analytically from the first term of the Taylor-series expansion of the convolution integral. The discovered SGS heat flux in RBC is also consistent with the truncated Taylor-series expansion. We have demonstrated a few important points about these discovered closures:

1. They all have high CC (often >0.9 -0.95) with the true SGS terms obtained from filtered DNS data, that is, excellent performance based on this commonly used a priori test metric. The same closure is discovered regardless of the system because the expansion's first term dominates the MSE loss function of RVM. 2. Despite this high CC, all a posteriori (online) tests result in unstable LES. This is consistent with the past findings about the NGM in the LES community (mainly for 3D turbulence) and in the climate community (for geophysical turbulence). Here, we argue that the inability of NGM to capture any inter-scale kinetic energy transfer in 2D-FHIT (or any 2D flow filtered in both directions) is likely the reason for the instability. For RBC, where filtering is done only in one direction, deeper investigations into the spectra of the SGS fluxes and inter-scale enstrophy and potential energy transfer, pointed to another likely reason for the instability: poor representation of the backscatter of SGS potential energy. This suggests that the poor a posteriori (online) performance of NGM might have different reasons in different flows. 3. The exact form of the discovered closure depends on the filter type and the filter size, Δ. For filters with compact support (i.e., all filters used here except for sharp-spectral cutoff), the structure of the closures is the same, but the coefficients are different (still, consistent with the Taylor-series expansion, as shown in the appendices). For the sharp-spectral cutoff filter, the equation discovery fails, again, consistent with the fact that the Taylor-series expansion cannot be conducted, a known issue in the LES literature (Sagaut, 2006). Finally, we would like to emphasize that what is reported here, for example, based on low CC, is the "failure" of the "equation discovery process," whose goal was to minimize the MSE loss to match the stress (i.e., structural modeling). As discussed earlier in the paper as well as further below, functional closures with low CC (such as Smagorinsky (1963)) can produce stable and relatively reasonable LES, thus low CC does not necessarily mean a failed closure. That said, we could not discover any parsimonious closure with the sharp-spectral filtered data that was consistent across several cases and resolutions, and thus, we do not perform any a posteriori LES with these low-CC closures. Note that as mentioned before, with a "smoothed" sharp-spectral filter, Ross et al. (2023) successfully performed equation-discovery.

As a side note, while the terms of the discovered closures might look like truncation error of finite difference/ volume discretization, in our work, all calculations (DNS solver, SGS terms, library) are done using spectral methods. This further shows, along with the Taylor-series expansion results, that the discovered closures are indeed representing the physics of the SGS terms, rather than any numerical error.

As an additional piece of analysis, we also present the decomposition of the SGS terms to the Leonard, cross, and Reynolds terms. We show that the Leonard and then cross terms often dominate the total SGS term, though the relative amplitude of these terms decreases as the filter size increases. However, this analysis shows that only computing the Reynolds momentum stress or heat flux can lead to discovering an inaccurate closure (and in general, in data-driven SGS modeling, in too-small SGS fluxes). That said, the relative importance of these three terms depends on the filter type and size, and likely, on the flow's spatial spectrum (Sun et al., 2023).

The analyses presented in this paper are aimed at highlighting the promises and challenges of the equationdiscovery approach to SGS modeling. On one hand, it is promising that this approach robustly discovers closures that could be closely connected with those mathematically derived, and could be easily interpreted and analyzed in terms of turbulence physics. On the other hand,

Journal of Advances in Modeling Earth Systems

10.1029/2023MS003874

a) The commonly used MSE loss function, or similar loss functions, will be always dominated by the leading term(s) of the Taylor-series expansion. Thus, sparsity-promoting equation-discovery techniques, at least with the common derivative/polynomial-based libraries, will always find the NGM (if all the relevant terms exists in the library). Note that this is true for the closure of any SGS process, as the Taylor-series expansion applies to any compact filter. Thus, this point and many of the main points of this paper are relevant beyond just SGS modeling for turbulence, but also SGS modeling of other nonlinear, multi-scale processes in the Earth system. b) Given that our diagnoses show shortcomings of the NGM with functional modeling metrics (e.g., inter-scale energy transfer), one idea is to include such physics constraints in the loss function. For example, Guan et al. (2023) demonstrated that a loss function that combines structural and functional modeling constraints can enhance the a priori and a posteriori performance of the data-driven closure model in the small-data regime. More functional-modeling physics constraints (as domain averaged or wavenumber-dependent quantities) can be included in the loss function, which can potentially close the gap between a priori (offline) and a posteriori (online) performance. While the loss function of some techniques such as RVM may not be flexible to change beyond MSE, other methods such as genetic programming, gene expression programming or symbolic regression provide such flexibility (Cranmer, 2023;Ross et al., 2023). Also, equation-discovery using neural network-based algorithms has gained popularity recently, as for example, their loss functions can be very flexible given the use of backpropagation for training (Z. Chen et al., 2021). That said, "spectral bias" (Xu et al., 2019), the fundamental challenge of neural networks in learning high wavenumbers, can become an issue when equation-discovery is the goal; see Mojgani et al. (2023) for an example and a solution in a quasigeostrophic turbulence test case. c) The fault may not entirely (or at all) lie with the MSE loss function. Guan et al. (2022) showed that a deep CNN with basically the same MSE loss function as the one used here (which only accounts for structural modeling) can learn a closure for 2D turbulence that has CC > 0.95 and leads to stable and accurate LES (and accurate inter-scale transfers; see Guan et al. (2023)). But a major difference between the CNN and RVM is that the former does not use a pre-defined set of basis functions, but rather, learns them. Recent work by Subel et al. (2023) has shown that the CNN of Guan et al. (2022) learned a set of low-pass, high-pass, and band-pass Gabor filters. As another major difference, the CNN's sparsity is not user-defined, but rather, comes from overparameterization. d) Related to (c), the discovered closures can depend on the choice of the library. This issue can be addressed by trying more expansive libraries (though this can lead to non-robust discoveries) or as mentioned earlier, by using methods such as genetic programming or gene expression programming, which allow the library to evolve (see Schmidt and Lipson (2009), Udrescu andTegmark (2020), andRoss et al. (2023)). Libraries inspired by the CNNs' basis functions or distilled from other deep neural networks could be explored as well (Cranmer et al., 2020;Subel et al., 2023). Furthermore, there are studies, for example, based on the Mori-Zwanzig formalism, suggesting that memory has to be included in closures (Parish & Duraisamy, 2017;Wouters & Lucarini, 2013). Hence, basis functions that include temporal information (as already used in Ross et al. (2023)) should be further explored in future work. e) Choosing the hyper-parameter(s) that determine the level of sparsity might require more thought too. While the L-curve criterion has shown success in many problems, the metrics for which the curve is calculated for should be further investigated. The common a priori metrics such as CC of SGS fluxes are completely incapable of identifying shortcomings from a functional modeling perspective, such as lack of inter-scale energy transfer or poor representation of backscattering, which can be diagnosed using additional metrics. Note that a high CC of SGS fluxes has been found as the necessary but not sufficient condition for a successful closure (Meneveau, 1994). f) Aside from all of the above issues related to the discovery algorithm, what needs to be discovered (the "truth") should be further examined. The discovered closures can depend on the filter type/size and the methodology (e.g., calculating Reynolds stress or the full SGS stress), because what is diagnosed as the "true" SGS flux from DNS has such dependencies. This has important implications for any data-driven SGS modeling approach, including those using deep neural networks or any other statistical learning method (Fatkullin & Vanden-Eijnden, 2004;Grooms et al., 2021;Guan et al., 2022;Sun et al., 2023;Zanna & Bolton, 2021). In fact, this non-uniqueness and uncertainty of the true SGS term is a major shortcoming of the supervised/offline learning approach to data-driven closure modeling in real-world applications, and is one of the main motivations to pursue online or at least offline-online learning approaches (Pahlavan et al., 2024;Schneider, Stuart, & Wu, 2021;Schneider et al., 2023).

∂ω ∂t + N(ω, ψ) = 1 Re ∇ 2 ωχωf , (B3) where N(ω, ψ) is N(ω, ψ) = ∂ψ ∂y ∂ω ∂x -∂ψ ∂x ∂ω ∂y . (B4)

f is a deterministic forcing (Chandler & Kerswell, 2013;Kochkov et al., 2021):

where f k x and f k y are the forcing wavenumbers and χ = 0.1 represents the Rayleigh drag coefficient. Comparing Equation 5with Equation B3, it is evident that ∇ × R =χω and ∇ × F =f .

In DNS, Equations B2 and B3 are solved in a doubly periodic domain using a Fourier-Fourier pseudo-spectral solver with second-order Adams-Bashforth and Crank Nicholson time-integration schemes for the advection and viscous terms, respectively (time step Δt DNS ). For LES, we use the same solver with lower spatio-temporal resolution: We use N LES that is 8-64 times smaller than N DNS , and Δt LES = 10Δt DNS .

Appendix C: The RBC Numerical Solver

We solve Equations 9-11 using a Fourier-Chebyshev pseudo-spectral solver (Khodkar & Hassanzadeh, 2018;Khodkar et al., 2019). Briefly, using the ωψ formulation, the dimensionless governing equations become For DNS, we solve Equations C1-C3 in domain (6π,1) . Periodic boundary conditions are imposed in the horizontal direction and no-slip and fixed temperature boundary conditions are imposed on the horizontal walls. Second-order Adams-Bashforth and Crank Nicholson time integration schemes are used for the advection and viscous terms, respectively. Table 2 presents the N DNS and N LES for each case. For LES, we use the same solver but with lower resolution in the horizontal direction.

Appendix E: Taylor-Series Expansion of the SGS Flux for a 1D Field

Let's focus on a spatially 1D field a(x) (dependence on t is not explicitly written for brevity). The filtering operation's convolution integral (Equation 1) becomes

G(r x ) a(xr x ) dr x , (E1)

The Taylor-series expansion of a(xr x ) around a(x) gives a(xr x ) = a(x) -

Substituting this into Equation E1 and using a = a(x), a = a(x) for brevity yields

The second line follows the first line considering that a and its derivatives do not depend on the variable of integration, r x . In Equation E4, a depends on a and its derivatives, with coefficients that only depend on the filter type and size through moments of the kernel, G. For example, for a Gaussian filter (Table A1)

Note that all the odd moments are 0, resulting in O(Δ 4 ) as the order of the truncated terms once the moments in Equation E5 are substituted in Equation E4: Journal of Advances in Modeling Earth Systems 10.1029/2023MS003874 a = a + 1 2! Δ 2 12 ∂ 2 a ∂x 2 + O(Δ 4 ). (E6) To calculate a term like τ xx = u 2u 2 , we first use a = u in Equation E6 and then square it to arrive at u 2 = u 2 + 2u( 1 2! Δ 2 12 ∂ 2 u ∂x 2 ) + O(Δ 4 ). (E7) Next, we use a = u 2 in Equation E6 to obtain u 2 = u 2 + 1 2! Δ 2 12 ∂ 2 u 2 ∂x 2 + O(Δ 4 ), (E8) = u 2 + 2 2! Δ 2 12 (( ∂u ∂x ) 2 + u ∂ 2 u ∂x 2 ) + O(Δ 4 ). (E9) Using Equations E7 and E9 we find τ xx = u 2u 2 = Δ 2 12 ( ∂u ∂x ) 2 + O(Δ 4 ). (E10) Note that this expression depends on u rather than u, which is what we desire. Next, we use a = ∂u/∂x in Equation E6 to obtain ∂u ∂x = ∂u ∂x + 1 2! Δ 2 12 ∂ 3 u ∂x 3 + O(Δ 4 ). (E11) Using this expression in Equation E10 yields an analytically derived closure for τ xx with error O(Δ 4 ) τ NGM xx = u 2u 2 = Δ 2 12 ( ∂u ∂x ) 2 .

(E12)

This is the NGM (Leonard, 1975;Sagaut, 2006). Four issues should be emphasized here i. This procedure can be followed for any filter type. However, the Taylor series does not exist for filters such as sharp-spectral (cutoff), whose kernel's second-order moment is indefinite. Thus, for such filters, NGM does not exist (Meneveau & Katz, 2000;Sagaut, 2006). ii. The same procedure can be followed to derive NGM for higher dimensions, for example, τ NGM xx , τ NGM xy , and τ NGM yy in 2D; see Sagaut (2006). iii. The coefficients in NGM depend on the filter's kernel and its moments (Equation E5). For Gaussian and tophat, the parameters of the kernels are chosen to match their first moment, leading to Δ 2 /12 coefficient for both. However, the coefficients differ for higher-order terms (Sagaut, 2006). iv. The procedure presented above is not specific to turbulence or even dynamical systems. The procedure and its outcome are valid for the filtered quadratic nonlinearity of any two variables, even random variables.

S = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂u ∂x 1 2 ( ∂u ∂y + ∂v ∂x ) 1 2 ( ∂u ∂y + ∂v ∂x ) ∂v ∂y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (F1) τ r = τ -1 2 tr(τ)I, (F2) where I is the identity matrix. In 2D with filtering in both directions, the anisotropic part of the SGS stress tensor from the NGM is τ NGM-r = τ NGM -1 2 tr(τ NGM ) I. (F3) τ NGM-r 2D = Δ 2 12 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2 (( ∂u ∂y ) 2 -( ∂v ∂x ) 2 ) ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂y -1 2 (( ∂u ∂y ) 2 -( ∂v ∂x )

2

Inserting this and S (Equation F1) into Equation 30shows zero point-wise inter-scale (kinetic) energy transfer in NGM: P NGM τ (x, y, t) = 0.

In buoyancy-driven turbulence such as RBC, the total inter-scale energy transfer rate P E is the sum of the rate of transfer of kinetic energy (P τ ) due to SGS momentum fluxes and potential energy (P J ) due to SGS heat fluxes (Eidson, 1985;Peng & Davidson, 2002):