Subgrid-scale (SGS) closures, or parameterizations, are essential for simulating real-world systems involving turbulent flows [1], [2]. A prominent example is climate modeling, where closures are required for simulating the geophysical turbulence in the atmosphere and ocean, along with various other nonlinear, multi-scale processes [3]–[5]. Despite over a century of research, a general framework for deriving turbulence closures from first principles has remained elusive, resulting in the widespread reliance on semi-empirical and ad-hoc models [2], [6], [7]. In fact, the shortcomings of these closures are the main source of epistemic uncertainty in climate change projections and weather forecasts, especially for extreme events [8]–[11]. This has led to extensive efforts in recent decades to pursue innovative physical, mathematical, and, recently, artificial intelligence (AI) tools to develop SGS closures for geophysical flows [11]–[14].

An ideal closure 1) should accurately capture the detailed structure of the SGS fluxes, 2) should accurately reproduce the interactions between the SGS processes and the resolved, large-scale dynamics, e.g., the inter-scale energy and enstrophy transfers, and 3) when coupled to the numerical solver of the resolved scales, e.g., in large-eddy simulation (LES), the simulated flow should have key characteristics such as energy and enstrophy spectra matching those of the direct numerical simulation (DNS). In the context of weather and climate prediction, of particular interest is also capturing the statistics of rare (extreme) events. However, currently, a framework for developing closures satisfying (1)-(3) is lacking [1], [2], [6].

The “structural modeling” approach [2], e.g., based on a second-order truncated Taylor-series expansion of the SGS flux, yields closures such as the nonlinear gradient model (NGM2), which satisfies (1) with above 0.9 pattern correlation with filtered DNS [15]–[17]. However, these closures fail (2)-(3), e.g., producing no inter-scale energy transfer in 2D turbulence and leading to unstable LES [18]–[20]. The “functional modeling” approach [21]–[24] leads to eddy-viscosity closures such as Smagorinsky [25] and Leith [24], which produce stable LES and some aspects of (2) though they are often overly dissipative and miss backscattering (the transfer of energy from the SGS processes to the resolved scales), an important process in atmospheric and oceanic flows [4], [26]–[35]. Excessive dissipation also degrades the representation of extreme events in LES with these closures [36], [37]. Furthermore, these closures fail (1), e.g., lower than 0.5 pattern correlation [36], [38], [39].

Recently, AI approaches for developing data-driven closures have received significant attention. Most efforts have focused on supervised learning with filtered DNS data, which is an approach akin to structural modeling [7], [11], [40], [41]. While some studies have demonstrated that these data-driven closures can satisfy conditions (1)–(3), they suffer from three major limitations: (i) the need for large amounts of DNS snapshots for training, (ii) lack of interpretability, and (iii) poor generalization to out-of-distribution regimes, e.g., extrapolation to much higher Reynolds numbers [36], [42]–[45]. Self-supervised learning approaches, akin to functional modeling, only use DNS statistics and have also shown success; however, they still suffer from shortcomings (ii) and (iii) [46]–[48].

A number of studies have pursued an alternative approach: “equation discovery”, a class of AI techniques that aim to find a closed-form mathematical representation of the data [49]–[56]. The major appeals of this approach are data efficiency, interpretability, and generalizability, thus addressing the shortcomings of deep learning. In their pioneering work, Zanna and Bolton [19] used Bayesian sparse regression to find a closed-form equation for SGS momentum stress in geophysical turbulence. They robustly discovered a closure that resembled NGM2 but found that while this closure accurately captured the detailed structure of the SGS fluxes, the LES was unstable (unless the predicted SGS flux was attenuated). Using the same methodology, [20] discovered closures of the same structure for SGS momentum and heat fluxes across various setups of 2D turbulence and Rayleigh-Bénard convection and confirmed that they all exactly match the outcome of second-order truncated Taylor-series expansion; i.e., the discovered closure is NGM2. They also demonstrated that the LES instability results from NGM2’s inability to capture any interacscale energy transfer (neither diffusion nor backscattering). However, it remained unclear how a better equation discovery approach can be devised, though the need for a physics-informed loss function, evolving library (e.g., through genetic programming [32], [49]) were speculated.

In this Letter, we inform the criterion for “discovery” of a sparse representation with the fundamental understanding of turbulence physics. We robustly find, for diverse setups of 2D turbulent flows and with an expansive library, a closure that is NGM2 plus an additional term. We further show that this additional term is, in fact, the 4th-order term in the Taylor-series expansion of the SGS flux. We demonstrate that this new closure (NGM4, hereafter), accurately captures the structure of the SGS flux (pattern correlation \(\approx 0.99\)), accurately reproduce the inter-scale energy and enstrophy transfers, and leads to stable LES, that among other things, accurately reproduces the statistics of the rare, extreme events in diverse test cases that cover a broad range of dynamics mimicking atmospheric and oceanic turbulence.

The dimensionless continuity and momentum equations for LES of 2D turbulence in spatial dimensions (\(x,y\)) are given by: \[\begin{align} \frac{\partial \overline{u}_i}{\partial x_i} & = & 0, \tag{1}\\ \frac{\partial \overline{u}_i}{\partial t} +\frac{\partial\overline{u}_i\overline{u}_j}{\partial x_j} & = & -\frac{\partial \overline{p}}{\partial x_i} + \frac{1}{Re}\frac{\partial^2 \overline{u}_i}{\partial x_j x_j} -\frac{\partial \tau_{ij}}{\partial x_j} + \overline{\mathcal{Q}}_i. \tag{2} \end{align}\] Here, \(\overline{(.)}\) represents a low-pass filtering operation, \(u_i\) denotes velocity, \(p\) is pressure, \(Re\) is the Reynolds number, \(\overline{\mathcal{Q}}_i\) represents a time-invariant external forcing, Rayleigh drag, and the Coriolis force. The SGS stress tensor, \(\tau_{ij}=\overline{u_iu}_j-\overline{u}_i\overline{u}_j\), requires a closure model that expresses \(\tau_{11}, \tau_{12}=\tau_{21},\) and \(\tau_{22},\) as a function of the resolved flow variables (\(\overline{u}_i, \overline{p}\)).

We investigate four cases of 2D turbulence (see Figure 4), generating a diverse flow dynamics that vary in dominant length scales, energy (\(E=\frac{1}{2}u_iu_i\)), and enstrophy (\(Z=\frac{1}{2}\omega^2\)) cascade regimes (\(\omega\) is the vorticity). Here, we regard direct numerical simulation (DNS) data as the ground truth and utilize filtered DNS (FDNS) data to learn closure for \(\tau_{ij}\). The FDNS data are obtained by first filtering the DNS data using a Gaussian filter and then coarse-graining the results to the \(4\) to \(64\times\) coarser LES grid following [20].

The equation discovery method used in this study is a sparsity-promoting Bayesian linear regression approach [19], [52], [57] based on the relevance vector machine (RVM) [58], employed to derive closed-form closures for each element of \(\tau_{ij}\) using FDNS data. This method operates on a library of basis functions, \(\boldsymbol{\Phi}\), comprising linear or nonlinear combinations of relevant variables, such as velocity or its derivatives: \[\begin{align} \left[ \frac{\partial^{\left(q_1+q_2\right)} A}{\partial x^{q_1} \partial y^{q_2}}\right]^{p_1} \left[ \frac{\partial^{\left(q_4 + q_5\right)} B}{\partial x^{q_4} \partial y^{q_5}}\right]^{p_2} \left[ \frac{D C}{D t}\right]^{p_3},\label{eq10} \end{align}\tag{3}\] where \(A, B = \overline{u}\) or \(\overline{v}\), and \(C = \overline{\omega}\) or an element of \(\overline{S}_{ij} = \frac{1}{2}\left(\frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial y_i} \right)\), the strain rate. Material derivative \(\frac{D}{Dt} = \frac{\partial}{\partial t} + u_j\frac{\partial}{\partial x_j}\) is included to account for non-Markovian (memory) effects, motivated by Mori-Zwanzig formalism [59], [60].

This library is extensive, with integers \(0 \le q \le 8\) and \(0 \le p \le 2\), although the total derivative order is constrained to 8th (yielding a total of 930 terms in the library). To accurately compute higher-order spatial derivatives required for this study, we used arbitrary-precision methods [61]. The construction of this library is motivated by the Galilean-invariant and symmetry properties of the SGS terms; e.g., this library includes the Pope’s tensors [32], [62]–[68].

The library must be comprehensive enough to express \(\boldsymbol{{s}}\), a vectorized element of \(\tau_{ij}\), as \(\boldsymbol{s}^{\text{RVM}} = \boldsymbol{\Phi}\boldsymbol{c}\). The regression weights, \(\boldsymbol{c}\), are optimized by minimizing the error \(\text{MSE} = \lVert \mathit{S}^{\text{RVM}} - \mathit{S}^{\text{FDNS}} \rVert^2_2\), where \(\mathit{S}\) consists of \(n\) stacked samples of \(\boldsymbol{s}\). The RVM enforces sparsity by iteratively pruning basis functions with weights exhibiting uncertainties above a predefined threshold, \(\alpha\), and recalculating the MSE until all remaining functions have uncertainties below \(\alpha\). A higher \(\alpha\) increases model complexity but reduces the MSE.

For each of the four cases in Fig. 4, we separately discover closures for three elements of the SGS stress tensor. We also examine several LES grid sizes, and the effect of varying \(\alpha\), which as mentioned earlier, is a key hyper‐parameter in equation discovery. We use \(n\) = 100 FDNS samples from a training set and 20 FDNS samples from an independent testing set. As an illustrative example, Fig. 1 (a) displays the mean CC for all elements of \(\tau_{ij}\) as \(\alpha\) increases. Note that this structural modeling approach is what all past equation-discovery studies followed [19], [20]. For small \(\alpha\) values (\(\alpha < 1\)), no closure is discovered (CC = 0, zero terms). As \(\alpha\) increases, the number of discovered terms increases, and the CC value rises and eventually plateaus, forming an “L-Curve.” The elbow of this curve is commonly used to identify the \(\alpha\) that balances accuracy and model size in equation discovery [19], [52], [57]. In Fig. 1 (a), the CC-\(\alpha\) relationship for each element of \(\tau_{ij}\) converges to approximately 1 and robustly discovers the same model at the elbow. Recent studies [19], [20] found that the elbow of this curve (red circles) corresponds to the analytically derivable NGM2 [15], [16]. However, as these studies and earlier work found, LES with NGM2 leads to unstable a posteriori (online) simulations. [20] suggested that this instability is due to the inability of NGM2 in capturing any inter-scale energy transfers; i.e., \(P_{\tau}^{\text{NGM2}}(x,y) = 0\) (NGM2 captures neither diffusion nor backscattering).

Here, we introduce a physics-informed discovery criterion that requires not only accurate reconstruction of \(\tau_{ij}\) but also capturing the inter-scale energy transfer (both diffusion and backscattering), thus combining structural and functional modeling. Figure 1 (b) demonstrates that the L-curve in this new approach does form an elbow arond NGM2 (red circules), but rather, discovers a new closure (black circles) that has CC of \(\tau_{ij}>0.99\) (Fig. 1 (a)) and non-zero inter-scale energy transfer (Fig. 1 (b)). We have found that the closure discovered at the new elbow is NGM2 plus the second term from the Taylor-series expansion involving the higher orders (\(\mathcal{O}(\Delta^4\))), where \(\Delta\) is the LES grid spacing (see Eq. 4 ). We refer to this new closure as NGM4. As shown below, NGM4 outperforms existing physics-based closures in a priori and a posteriori tests.

It is noteworthy that the L-curves in Fig. 1 are consistently observed for all cases with any other \(N_{\text{LES}}\). Note that increasing \(\alpha\) leads to another closure with more terms but negligible improvements in CC; we have found this closure to be NGM6; see Eq. 4 (a generalized Taylor-series expansion of \(\tau_{ij}\) is provided in the supplementary information). The equation discovery with the expansive 930-term library confirms that Taylor-series expansion offers the best representation of both \(\tau_{ij}\) and the inter-scale energy transfer. Note that the memory terms never emerged during the discovery.

We first compare the performance of NGMs (NGM2-6) with common baseline closures, the eddy-viscosity Smagorinsky (Smag) [25], its dynamic version (DSmag) [22], and the dynamic Leith (DLeith) [24] models (see Appendix B). Figure 2 shows that NGM4 and 6 perfectly capture the structure of the SGS stress and inter-scale enstrophy transfers, and very well capture the interscale energy transfers (both diffusion and backscattering), significantly outperforming NGM2 and the physics-based closures. An important feature of NGM4 and 6 is that they capture the energy backscatter that is seen in the FDNS. While the overall representations of inter-scale energy transfer of NGM4 and 6 are not perfect (around \(62\%\) and \(78\%\) of FDNS’, respectively), they are substantially better than the energy transfers of NGM2 (\(=0\)) and the eddy-viscosity closures, which are purely and excessively diffusive (as indicated by the negative CC values). The advantages of NGM4 and 6, especially over the eddy-viscosity closures, can also be seen in the inter-scale energy and enstrophy transfers’ spectra (Fig. S1). The a priori tests show that NGM4 and 6 are failry comprable, though NGM6 has a noticeably better \(P_\tau\) (see Table S1).

NGM4 and 6’s strong a priori performance is reflected in their superior a posteriori performance. LES [69] with NGM 4 and 6 remain long-term stable when the LES resolution is high enough to capture over \(80\%\) of DNS enstrophy (see Fig. S2). This is reminiscent of findings in 3D turbulence, where LES typically needs a resolution high enough to cover at least \(80\%\) of DNS energy [1].

As shown in Figs. 3 and S3, when stable, LES with NGM4 and 6 capture both the bulk and tails of the probability distribution function (PDF) of vorticity (the tails of this PDF represent “extreme weather” in this prototype of atmospheric circulation). LES with eddy-viscosity closures fail to capture the extremes, especially for Cases 3 and 4, which have high-wavenumber forcings. Note that NGM2 for any case and any LES resolution leads to instabilities.

In addition to better capturing extreme events, LES with NGM4 and 6 outperforms eddy-viscosity models in short-term forecasting of vorticity (see insets of panels (a)-(b) in Figs. 3 and S3). further demonstrate the correlation coefficient (CC) between the vorticity of FDNS, \(\overline{\omega}_{\text{FDNS}}\), and the SGS closures, \(\overline{\omega}_{\text{LES}}\). All models start from the same initial condition, with NGM4/6 maintaining a trajectory closer to FDNS than the eddy-viscosity models.

LES with NGM4 and 6 also outperform LES with other closures in capturing the domain-averaged kinetic energy and enstrophy (Table S2). Figures 3 and S3 (c)-(d) show the enstrophy spectra of FDNS and LES. The insets magnifying the tails of the spectra do not clearly show the superiority of LES with any closure. However, the insets magnifying the spectra at lower wavenumbers reveal that LES with eddy-viscosity closures struggle to capture the bulk of the enstrophy. This observation is also consistent with the energy spectra (Fig. S4). Note that in a posteriori tests, LES with NGM4 and NGM4 show fairly similar performances, suggesting that in practice, using NGM4 might be enough.

In conclusion, a physics-informed data-driven discovery identified an analytically derivable, yet previously ignored, representation of SGS closure that satisfies both structural and functional modeling requirements. This closure (NGM4) is “interpretable” as it is derivable from the Taylor-series expansion and it represents the Leonard and cross stresses (see supporting information). This closure is “generalizable”, as its coefficients are entirely determined by LES resolution (\(\Delta\)). NGM4 is accurate: it is the first closure to simultaneously excel in both a priori and a posteriori tests. LES with NGM4 is stable only when the numerical resolution is high enough to capture \(80\%\) of DNS enstrophy; however, this is consistent with general limitations of LES [1].

A major difference between NGM4 and NGM6 is that the latter represents the Reynolds stress too (see supporting information). However, while a priori tests show that NGM6 better represents the inter-scale energy transfer, a posteriori tests do not demonstrate any advantage over NGM4 in terms of stability or accuracy. This is encouraging as implementing models with high-order derivatives in numerical solvers can be challenging, but these findings suggest that NGM4 might suffice in practice, e.g., in ocean modeling.

We note that NGM4 and 6 violate the Boussinesq hypothesis as they depend on second-order velocity gradients, highlighting the need to go beyond classical turbulence tensors to represent backscatter and anisotropic effects [62].

While here we only focused on the SGS stress tensor in 2D turbulence, these findings are relevant to multi-scale modeling in other dynamical systems too, as coarse-graining any quadratic nonlinearity can lead to an NGM-like representation. This work also presents the importance of combining physical and mathematical insight to better harness he power of AI methods in accelerating scientific discovery.

Acknowledgment↩︎

This work is supported by NSF grant AGS-2046309 and Schmidt Sciences. Computational resources were provided by NSF ACCESS (allocation ATM170020) and NCAR’s CISL (allocation URIC0009).

Data availability↩︎

The code for the numerical solver “py2d” is available at [69]. The codes for equation discovery and data analysis in this work can be found at [70].

Appendix↩︎

Appendix A: Figure 4 shows the cases considered in this study and their physical and numerical parameters.

Appendix B: The baseline physics-based closures used here are Smag, DSmag, and DLeith. These eddy viscosity closures introduce dissipation and do not account for backscattering: \(\tau_{ij} = -2\nu_{e}\overline{S}_{ij}\), with \(\overline{S}_{ij}\) representing the filtered rate of strain and \(\nu_{e}\) denoting the eddy viscosity [1]. \(\nu_e = \left(C_s\Delta\right)^2\sqrt{2\overline{S}_{ij}\overline{S}_{ij}}\) for Smag and \(\nu_{e} = \left(C_l \Delta\right)^3 |\nabla \overline{\omega}|\) for Leith. We use \(C_s = 0.17\) for Smag, as proposed in a previous work [21]; \(C_s\) and \(C_l\) are estimated dynamically based on the local flow structure for DSmag and DLeith, respectively. This procedure can yield \(\nu_{e} \leq 0\), potentially resulting in numerical instabilities; therefore, the commonly used “positive clipping" is applied to enforce diffusion (\(\nu_e \geq 0\)) [71].

References↩︎