12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

SUBGRID-SCALE SURROGATE MODELING OF IDEALIZEDATMOSPHERIC FLOWS: A DEEP LEARNED APPROACH USINGHIGH-RESOLUTION SIMULATION DATA

Muralikrishnan Gopalakrishnan Meena
Matthew R. Norman
National Center for Computational Sciences National Center for Computational Sciences
Oak Ridge National Laboratory
Oak Ridge National Laboratory
Oak Ridge, TN, 37831-6008, USA
Oak Ridge, TN, 37831-6008, USA
gopalakrishm@ornl.gov
normanmr@ornl.gov

David M. HallNVIDIA
Lafayette, CO, USA
dhall@nvidia.com

ABSTRACTlocation in the host model.

We introduce a deep learned subgrid-scale surrogate
model for dry, stratified idealized atmospheric flows fromhigh-resolution simulation data. Deep neural networks (NNs)
are used to model the full state differences between a coarse
resolution simulation and a high-resolution simulation, run simultaneously with the coarse resolution simulation forced by
the high-resolution simulation, hence capturing both dissipative and anti-dissipative effects. The NN model is able to accurately capture the state differences in a priori tests outside thetraining regime. In a posteriori tests intended for productionuse, the NN coupled coarse simulation is accurate comparedto the high-resolution simulation over a finite period in time.
With the accumulation of the errors, the NN-coupled simulation becomes computationally unstable after a while. Thesesurrogate models further pave the way for formulating stable,
complex, physics-based NN models which are driven by traditional subgrid-scale turbulence closure models.

Such models are limited by the appropriate choice of theparametrization. Moreover, simulating such models with theglobal climate model at large climate time scales is still acomputational challenge, even with the current generation ofaccelerator-based high-performance computing systems (Norman et al., 2021). Recently, the use of artificial intelligence has
been at the forefront to form surrogate models parametrizing
the subgrid-scale processes (Duraisamy et al., 2019; Brunton
et al., 2020; Kochkov et al., 2021), taking advantage of theavailability of high-resolution simulation data. These models
generally rely on components in the eddy-viscosity term of thefiltered Navier-Stokes equation, which concentrate on capturing the dissipative effects (Gamahara & Hattori, 2017; Raspet al., 2018; Maulik et al., 2019; Yuval & O’ Gorman, 2020).
Moreover, deep learning techniques have been used to capture and parametrize the subgrid-scale processes in global climate models (Krasnopolsky et al., 2013; Gentine et al., 2018;
Brenowitz & Bretherton, 2019).

Introduction

As stratified turbulent flows encountered in cloud physics
also are comprised of significant anti-dissipative effects, it isbeneficial to capture both the dissipative and anti-dissipative
effects through the SGS closures. We approach this problem
by directly modeling the full state difference between a coarse
simulation and a simultaneously running very fine resolution
simulation, enabling us to capture both the dissipative andanti-dissipative effects from the subgrid-scale processes. An
overview of the modeling framework is shown in Fig. 1. We
use the two-dimensional (2D) inviscid, non-hydrostatic, compressible Euler equations with density stratification for solving
canonical LES and mesoscale atmospheric flows and to generate subgrid-scale data as the difference between high- andlow-resolution flow field data.

The challenge for numerically analyzing atmospheric turbulence is tied to the broad range of multi-scale physics constituting such flows. More specifically, the added complexityof density stratification involving cloud microphysics and turbulence leads to greater challenges (Wyngaard, 1992). Capturing the complex subgrid-scale processes at the scale of clouds
is important for constructing accurate global climate models(Schneider et al., 2017). There have been various efforts in theatmospheric science community to capture the effects of thesubgrid-scale processes using some form of parametrization.
Such parametrizations, like large-eddy simulations, rely on filtering the high-wavenumber (subgrid-scale) structures and resolving only a coarser domain with embedded subgrid-scalemodeling (SGS modeling) (Pressel et al., 2017). The SGS
models parametrize the viscous dissipation from the filteredsubgrid-scale structures. Cloud resolving models (CRM) areanother example for such a modeling framework used in simulating global climate models (Randall, 2013; Hannah et al.,2020). CRMs are comprised of a high-resolution cloud model
which parametrizes deep convection for the fluid state at that

Instead of handling the problem on the global spatial domain all at once, we choose to handle it locally using stencils ofdata (not unlike the nature of convolutions), to reduce the deeplearning model size, improve generalization, and increase thenumber of samples available for training. To avoid chaotic divergence between the two simulations, the states of coarse sim-

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
scalars, ρ 0 and (ρθ )0 ), removing the dominant underlying hydrostatic balance described by Eq. (4). Although inviscid, theEuler equations form the basis of atmospheric flow simulations, and inherent numerical dissipation maintains stability.
Moreover, the 2D setting is an apt choice for idealized simulation as the majority of idealized test cases in literature are
2D.
The flow is perturbed from a neutrally stratified dry atmosphere. The domain is of size (x, z) ∈ [0, 20] × [0, 10] km.
Slip, solid wall boundary conditions are prescribed to the top
and bottom walls. The left and right walls are prescribed withperiodic boundary conditions. The time steps are adjusted sothat the Courant–Friedrichs–Lewy number is 0.8. We choose
the grid resolution, nx × nz , of the coarse simulation to be
200 × 100 (grid spacing of 100 m in all directions). The fineresolution domain has a resolution of 1000 × 500 (grid spacing of 20 m in all directions), resulting in a grid mapping ratioof 5× between the coarse and fine domains. We simulate theflows till 1000 secs. After the initial laminar regime dominated with anti-dissipative effects, the flow evolves to a highly
turbulent regime from about 500 secs onward. The temporal regimes consisting of both anti-dissipative and dissipative
effects makes the colliding thermal test case a useful model
problem to test the capability of the current framework to capture these effects. Further details of the test case, solver, andnumerical schemes can be found in Norman (2021).

ulation are updated by the states of fine simulation after eachtime step. The dissipative and anti-dissipative components areembedded into the coarse simulation as it is integrated withthe high-resolution simulation as a parallel driver. The current
modeling approach for 2D flows is useful as a pilot project todevelop the surrogate modeling techniques, particularly since
2D CRMs have been found to be accurate and computationally cheaper parametrizations for global climate models. Incorporating cloud microphysics to the present supervised approach could serve as a novel alternative to the SGS models of
CRMs, capturing both dissipative and anti-dissipative effectsof the subgrid-scale processes.

General procedure
We use NNs to augment the flow states of the coarse resolution simulation, qcoarse , with states of a simultaneously running fine resolution simulation, qfine . The majority of subgridscale effects missing from the coarse resolution simulation are
captured by the full state difference∆q = qfine − qcoarse ,(1)

where qfine denotes qfine interpolated to the coarse resolution
grid stencil by summing over the stencil around a given cellin the fine grid to reduce the dimension to that of the coarse
grid. A sample portrayal of the procedure is shown in Fig. 1.
Given a coarse resolution data, the objective is to obtain a NN
learned state difference, ∆qNN , and correct the coarse states toget
qNN = qcoarse + ∆qNN .

Neural network setup(2)

The procedure resembles a mapping from fine grid to coarse
grid. Modeling the full state difference enables us to capture
both the dissipative and anti-dissipative effects of the subgridscale processes. In the following section, we will introduce theflow problem we use to demonstrate the modeling technique.
We will discuss the details of the NN architecture and training
procedures in Sections and .

Fluid flow problem
We consider a 2D thermal collision of two hot and cold
thermals as a sample problem to demonstrate the modelingframework (Norman, 2021). The colliding thermals createstrong discontinuities, strong winds, and significant turbulent
regimes as the flow evolves in time. The flow evolution is described by the 2D, dry, compressible, non-hydrostatic Euler
equations, given by
 



 

ρu
ρw
ρ
0

 2


 

∂ 
ρuw
0

 ρu  + ∂ ρu + p + ∂ 
=
2







−(ρ − ρH )g
ρuw
∂t ρw
∂x
∂ z ρw + p − pH
ρθ
0
ρuθ
ρwθ
(3)
ρH = −

1 ∂ pHg ∂z(4)where ρ is density, u and w are wind velocities in x- andz-directions (horizontal and vertical directions, respectively),θ is potential temperature, p = C0 (ρθ )γ is pressure (C0 is
a constant, γ = c p /cv , c p is specific heat of dry air at constant pressure, and cv is specific heat of dry air at constantvolume), pH and ρH are the hydrostatic pressure and density,and g is acceleration due to gravity. Equation 3 can be represented in vector form as ∂t q + ∂x f + ∂z g = s with the state vector q = [ρ ρu ρw ρθ ]T . We only model the perturbed

2

We deploy a supervised deep learning technique to model
the relation between the inputs, a local stencil of coarse resolution states (qcoarse ), and the outputs, state difference in thecenter cell of the stencil with the fine resolution states (∆q).
We model this mapping between qcoarse and ∆q using deepNNs of the form N (qcoarse , w) with w representing the parameters (“weights”) of the NN. The learning process involves
optimizing the parameters w of the NN by minimizing the loss
computed by the function L[∆q, N (qcoarse , w)]. In the current analysis, we use a squared L2 norm as the loss function,L[∆q, N (qcoarse , w)] = k∆q − N (qcoarse , w)k22 .
A stencil-based approach is used to build the architectureof N (qcoarse , w). We use the stencil around each grid cellin the coarse simulation to model the ∆q at the respective locations as shown in Fig. 2. The framework of N (qcoarse , w)
is composed of 36 inputs, the 4 flow states of the 3 × 3 gridstencil around a given coarse cell. The outputs (4 values) arethe ∆q at the given grid cell (the center of the stencil of inputs). The hidden layer(s) of the NN are comprised of neurons of various configurations, enabling the nonlinear activation function to capture the relation between the inputs andoutputs (Goodfellow et al., 2016). We choose 3 different typesof NN architectures for the hidden layer(s) with a Leaky ReLU
activation function having a slope of 0.1 (Glorot et al., 2011):(a) a single layer of 45 neurons (single-layer model), (b) 10
layers with 45 neurons/layer in a ResNet-based configuration(ResNet model) (He et al., 2016), and (c) 10 layers with 45
neurons/layer in a DenseNet-based configuration (DenseNet
model) (Huang et al., 2017). We tune various hyperparametersof the NN architecture to arrive at these architectures. Thesehyperparameters for the current approach involve the numberof neurons per layer, number of layers, activation function, values for regularization at the input, hidden, and output layers,
optimization routine, and batch size of training samples.
The single-layer model is used as a shallow, simple NNmodel to test the capability of NNs to capture the grid mapping. Drop-out regularization is applied to the input layer tostabilize the single-layer model. Adding more layers to this12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Coarseresolution
ﬂow

Machinelearnedcorrection

Corrected
highresolution
ﬂow

ML
Model

Figure 1.

An overview of the NN-based approach to capture the subgrid-scale effects missing in a coarse resolution simulation.

Input Layer

inputs

NN modelOutput Layer

MLP

outputs

data randomly sampled in space and time, particularly in regions with high gradients of the states. Thus, we use a total variation (TV) weighted random sampling technique. The
TV of a sample at the cell (i, j) in the 2D domain is given
by V (qi, j ) = ∑0k=−1 {|qi+(2k+1), j − qi, j | + |qi, j+(2k+1) − qi, j |}.
TV identifies regions with high variations in states. As the
majority of the domain is comprised of quiescent regions withsmall changes in the states during the initial transient regime,these regions are not captured by the TV weighting. This limitation makes the models biased against the background mean.
To alleviate this bias issue, we use a 50 : 50 split to curate samples with and without high TV weighting. Selecting the threshold to identify the 50 : 50 split poses further challenges as theflow evolves through both laminar and turbulent regimes, leading to a broad range of values for the states over time. Thus,
the threshold for identifying the 50 : 50 split is manually adjusted for various flow regimes.
The input and output data for the NN is scaled according
to min-max normalization. For a vector of input and output
variables of the network, Q = [q ∆q]T , the normalization
is given by Q̃ = (Q − mintrain Q)/(maxtrain Q − mintrain Q),
where Q̃ is the rescaled vector and the operators mintrain andmaxtrain refers to the minimum and maximum operators over
the training sample space of the variables, respectively.
We train the model in Python environment using PyTorchand then deploy the model in the C++ solver with GPU acceleration for production usage. The “NAdam” optimizer isused for optimizing the weights of the NNs (Dozat, 2016).
For training the 1M data points (with a 70 : 30 split between
training and validation loops, and random shuffling of training
samples), we use mini-batches of 1024 samples. Even withthe use of mini-batches, the model incurs high variance in accuracy over the training epochs. We have learned that using
a higher learning rate at initial training epochs till saturationof the model accuracy and then lowering the value helps toreduce this high variance in accuracy. Furthermore, for the
deep NNs, we also rely on an averaging technique to ensemble
the model weights across various training epochs to alleviate
the high variance in the training accuracy. After a regime ofstable accuracy is reached over epochs, we collect ensembles
of model weights from models with error below a particularthreshold. The weights of the final model are an ensemble ofthese collected weights. The procedure helps average the gradient basin in the optimization manifold, aiding to form a stable solution. This procedure has similarities with the stochastic weighted averaging (SWA) method introduced in Izmailov
et al. (2018).

ResNetDenseNet

Figure 2. Illustration of the network model for modeling thefull state difference using stencil data.

model reduced the stability of the model during testing without considerable increase in accuracy and ability to capture
the nonlinear physics. This lead us to the last two models. The
last two models are deep NNs, used to increase the complexity of the network and accurately capture the highly non-linear
physics. These deep NNs take advantage of the so-called skipconnections to enable sparse, complex, non-linear relations between the input and output (He et al., 2016).
We have found that the simple, small framework of using
a stencil-based input and output to the network is advantageousover convolutional neural networks for higher accuracy andbetter stability of the model during testing in unknown dataregions. The architecture intrinsically is similar to convolution
operations. Moreover, the stencil-based approach mimics ahigher-order finite difference scheme. The current formulation
does not assume a dissipation- or advection-based mapping forthe relationship between qcoarse and ∆q. Our approach is to let
the NN capture the complex mathematical relations between
the states of the center cell and neighboring cells with the state
difference at the center cell in the coarse domain.
It is important to note that the actual domain of constraint
for the fine-scale model is wider than the 3×3 coarse cells over
the temporal domain of a coarse-grid time step. Therefore, fullconstraint of the state difference is likely not possible. Still, we
wished to maintain a simpler NN model with minimal inputs,
and results will show that the 3 × 3 stencil is effective in providing a constraining set of inputs.

Training and validation of model
We have learned that sampling the data for training is a
crucial step for building an accurate and stable NN model,
which we elaborate here. The training dataset is comprised of
1M samples collected between the temporal regime of [0, 900]
secs (both laminar and turbulent regimes). The small sizeof the stencil-based framework allows us to use significantly
fewer samples for training. We have found that the accuracyof the model is significantly deteriorated when trained with

Results
For an effective model predicting the subgrid-scale effects, the model should be accurate as well as computationally stable when tested in temporal regimes outside the training
regime. As with subgrid-scale turbulence models, we perform
both a priori and a posteriori testing (Piomelli et al., 1988) of3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
bly in the results of the L1 and L2 difference norms. The main
coherent structures are accurately captured and the larger L1
errors are from the small-scale fluctuations, isolated to a few
regions in the flow field. Note that these regions of high local errors are located at high-gradient regions. Note that we
also test the models in the initial laminar regime where antidissipative effects dominate, which are not discussed here. We
have found the models effectively capture high-gradient regions encountered in the laminar flow regimes as comparedto the coarse simulations.
With spatial accuracy comes the next challenge of temporal stability of the NN-coupled solver. The high spatial accuracy attained by the NN does not remain numerically stable indefinitely. As the flow evolves, the simulation becomes
unstable. The normalized error of the rms value of windspeeds and total kinetic energy (KE) over time are shown in
Fig. 6. The rms value (spatial) of wind speeds and the total kin ,nnetic energy are defined by urms = [1/(nx nz ) ∑i,xj z u2i, j ]0.5 andnx ,nz2
2
Total KE = 0.5[∑i, j (ρu)i, j + (ρw)i, j ], respectively, whereu = [ρu ρw]T . Overestimates of variability, as seen previously in Fig. 5 (d), accumulate over time until the solver becomes unstable at t0 + 17 secs, which is after ∼ 85 coarse-gridtime steps. The error in (ρu)rms by the coarse simulation withcorrection using the NN model is always lesser than that of thecoarse simulation without correction while in the stable regime(see Fig. 6 (a) solid and dashed blue lines, respectively). The
error in (ρw)rms by NN model is comparable to that of thecoarse simulation without correction. The total KE shows that
the NN-based correction indeed outperforms the coarse simulation without correction during most of the stable temporal
regime.
The results suggests that the NN-based correction ishighly accurate over a finite time interval, whereas further research needs to be done to make the solver stable from theaccumulation of the errors at high-gradient regions. Note the“spike” at the far right of Fig. 6. After this time, the model
becomes unstable. We note that the DenseNet model is able toperform similarly to that by the ResNet model, becoming unstable after t0 + 12 secs but with higher levels of error before it
becomes unstable. Whereas, the single-layer model is highly
unstable, as it becomes unstable after only t0 + 3 secs. Thisinability of the single-layer model emphasizes the need for acomplex formulation of the NN model with multiple layers
and skip-connections. We provide insights from our ongoingwork on improving the stability of the NN-coupled solver in
the Conclusion section.

the NN model in turbulent regimes outside the training regime.
The results reveal that the ResNet-based model gives themost accurate and computationally stable NN-coupled solver.
Herein, we only show the results of the ResNet-based model,
while we also discuss insights gained from the other models.

A priori testing For the a priori testing, the results fromthe NN model is directly compared with the expected solution,
evaluated independently without coupling the NN model with
the numerical solver. We observe the behavior of the model byusing data from a full snapshot in a turbulent regime outside
the training regime, as shown in Fig. 3 for θ 0 . The ResNetmodel predicts the state difference accurately with an overall
Euclidean error norm of L2 = k∆q − ∆qNN k2 /k∆qk2 = 0.195.
We also plot the L1 error, ε = |∆q − ∆qNN |/ max(|∆q|), tostudy the results locally in space, as shown in Fig. 3 (c). The
main structures are accurately captured whereas the smallerscale features have lower accuracy. Nonetheless, the maximum error is 10%, which is isolated to a small region in thelower-left side of the domain. The DenseNet model performed
close to that by the ResNet model for both the random samples (normalized errors have a mean and standard deviationof O(10−3 ) and O(10−2 ), respectively) and full flow field(L2 = 0.243) test datasets. Surprisingly, the single-layer model
did give accurate results (normalized errors have a mean andstandard deviation of O(10−2 ) and O(10−2 ), respectively, andL2 = 0.629) but did not perform as well as the ResNet model.

A posteriori testing The a posteriori testing is done byimplementing the NN model with the flow simulation. A sample depiction of the NN-coupled solver is shown in Fig. 4. The
fine simulation is used as correction until the testing time, t0 ,after which the NN model is invoked. Therefore, at time t0 ,the model state is essentially perfect with respect to the finegrid model. We use the PyTorch C++ frontend API1 to couple
the NN model with the C++ solver, performing the so-called
in-the-loop ML integration with the scientific solver. A GPU
kernel is used to update the qcoarse with the machine learned∆qNN after every time step, following Eq. 2. Even without special batching consideration of online inferencing and refactoring the code for GPU performance, we attain speed-ups up to
8× when using the NN-coupled solver compared to a fine resolution simulation. Since the current discussion is concentrated
more towards the scientific aspects of the model, we reserve
the full computational aspects of the framework for future investigation.
We invoke the NN model at time t0 = 900 secs, which
is outside the training regime and well within the turbulent
regime. Recall that the training data was sampled from t ∈[0, 900] seconds. The results we show here are those after thesimulation has run for t0 + 5 secs, which is approximately 25
coarse-grid time steps after t0 . The results for θ 0 of the coarse
simulation with correction using fine simulation, without correction (after t0 ), and with correction using the NN model are
shown in Figs. 5 (a-c), respectively. Without correction after t0(Fig. 5 (b)), diffusive effects start to dominate the coarse simulation and the subgrid-scale effects are not captured. The NN
model accurately predicts the flow evolution with the sharp
gradients, maintaining a difference of L2 = 0.118 and maximum L1 difference of 20%, as shown in Figs. 5 (c-d), respectively. Note that, the chaotic divergence of the modeled flow
and its difference from the ideal states are convolved insepara-

1 https://pytorch.org/cppdocs/. A sample implementation of the
API is provided in https://github.com/muralikrishnangm/pytorch-cppexample.

Concluding remarks and future work
We use deep learning to generate surrogate models fromhigh-resolution data for capturing subgrid-scale effects in dry,stratified turbulence in idealized atmospheric flows. Startingfrom the inviscid Euler equations with density stratification,
we generate subgrid-scale data as a difference between highand low-resolution flow field data. The coarse resolution simulation is kept in sync with the fine resolution simulation via
coarsening interpolation after each time step to avoid chaotic
non-linear divergence. The results show that a deep ResNet
type model using the stencil information around a given gridcell is able to accurately perform a posteriori tests and is stable over a significant time period. Such accuracy is consistent over both laminar as well as turbulent regimes, which are
dominated by anti-dissipative and disspative effects, respectively. The present supervised approach serves as a novel alternative to the subgrid-scale models in large eddy simulations
and CRMs, capturing both dissipative and anti-dissipative effects from the subgrid-scale phenomena.

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
(a)

(b)

(c)

Figure 3. Full flow field prediction by the neural network model. Potential temperature corrections are shown along with the L1 error
norm, ε = |∆q − ∆qNN |/ max(|∆q|) and L2 norm of the full flow field, L2 = k∆q − ∆qNN k2 /k∆qk2 .

REFERENCES
Coarse simulation

Brenowitz, N. D. & Bretherton, C. S. 2019 Spatially extended
tests of a neural network parametrization trained by coarsegraining. J. Adv. Modeling Earth Sys. 11 (8), 2728–2744.
Brunton, S. L. et al. 2020 Machine learning for fluid mechanics. Ann. Rev. Fluid Mechanics 52, 477–508.
Dozat, T. 2016 Incorporating Nesterov momentum into Adam.
Duraisamy, K. et al. 2019 Turbulence modeling in the age ofdata. Ann. Rev. Fluid Mechanics 51, 357–377.
Gamahara, M. & Hattori, Y. 2017 Searching for turbulencemodels by artificial neural network. Phy. Rev. Fluids 2 (5),054604.
Gentine, P. et al. 2018 Could machine learning break the convection parameterization deadlock? Geophysical Research
Letters 45 (11), 5742–5751.
Glorot, X. et al. 2011 Deep sparse rectifier neural networks.
Goodfellow, I. et al. 2016 Deep learning. MIT press.
Hannah, W. M. et al. 2020 Initial results from the superparameterized E3SM. J. Adv. Modeling Earth Sys. 12 (1).
He, K. et al. 2016 Deep residual learning for image recognition.
Huang, G. et al. 2017 Densely connected convolutional networks.
Izmailov, P. et al. 2018 Averaging weights leads to wider optima and better generalization.
Kochkov, D. et al. 2021 Machine learning–accelerated computational fluid dynamics. Proc. National Academy of Sciences 118 (21).
Krasnopolsky, V. M. et al. 2013 Using ensemble of neural networks to learn stochastic convection parameterizations forclimate and numerical weather prediction models from data
simulated by a cloud resolving model. Adv. Artificial Neural
Systems 2013.
Ling, J. et al. 2016 Reynolds averaged turbulence modellingusing deep neural networks with embedded invariance. J.
Fluid Mechanics 807, 155–166.
Maulik, R. et al. 2019 Subgrid modelling for two-dimensional
turbulence using neural networks. J. Fluid Mechanics 858,

Coarse simulation
with correction
Correction method after
each time stepFine simulation

NN modelFigure 4. Sample depiction of a posteriori testing using theneural network model.

We are currently investigating alternate formulations toenhance the stability of the model. Complex diffusion-based
approaches inspired from traditional SGS models are an alternative approach in the physics-informed realm. Such models
along with enforcing conservation of various conserved quantities can help stabilize the NN coupled solver (Ling et al.,2016). Moreover, recurrent neural networks and generative
adversarial networks can be used for correcting the emulation.
One aspect we have not explored in this manuscript is the explainability and interpretability of the NN models. Recently,
various techniques have been used to understand the physicalimplications of machine learning, particularly in the meteorological domains (McGovern et al., 2019). Finally, performant
and scalable integration of the model to the solver in hybrid architectures (Partee et al., 2021) is of utmost urgency for practical production use of such models in climate models.

Acknowledgments This research was supported by the
Exascale Computing Project (17-SC-20-SC), a collaborative
effort of the U.S. Department of Energy Office of Science andthe National Nuclear Security Administration. This researchused resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supportedby the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
(a)

(b)

(c)

(d)

Figure 5. Full flow state prediction using neural network based emulation in turbulent regime outside of training regime. Snapshotsof potential temperature perturbation (θ 0 ) after 5 secs from the initial state of testing regime are shown for: (a) coarse simulation with0 ), (b) coarse simulation without correction (θ 0correction using fine simulation (θfinecoarse-no correction ), and (c) coarse simulation with0
correction using NN model (θcoarse-NN ), along with the L2 = kqfine − qcoarse-NN k2 /kqfine k2 difference norm of the flow field. (d) L100 , ε = |qdifference of θcoarse-NN
with respect to θfine
fine − qcoarse-NN |/ max(|qfine |).
(a)(b)
Total Kinetic Energy
No correction fromResNet from

No correction fromResNet from

Figure 6. Comparison of L1 errors in the (a) rms values of flow velocity and (b) total kinetic energy over time. The dashed lines
are the variables from the coarse simulation without correction and the solid lines are those of the coarse simulation with correction
using NN model. The L1 error norm for the time series data is defined by ε = |fcoarse-no correction/coarse-NN − ffine |/|ffine |, where f =[urms Total KE]T .

ulation of turbulent channel flows. Physics of Fluids 31 (7),
1884–1891.
Pressel, K. G. et al. 2017 Numerics and subgrid-scale modeling in large eddy simulations of stratocumulus clouds. J.
Adv. Modeling Earth Sys. 9 (2), 1342–1365.
Randall, D. 2013 Beyond deadlock. Geophysical Research
Letters 40 (22), 5970–5976.
Rasp, S. et al. 2018 Deep learning to represent subgrid processes in climate models. Proc. National Academy of Sciences 115 (39), 9684–9689.
Schneider, Tapio et al. 2017 Climate goals and computing thefuture of clouds. Nature Climate Change 7 (1), 3–5.
Wyngaard, J. C 1992 Atmospheric turbulence. Ann. Rev. Fluid
Mechanics 24 (1), 205–234.
Yuval, J. & O’ Gorman, P. A. 2020 Stable machine-learningparameterization of subgrid processes for climate modelingat a range of resolutions. Nature Comm. 11 (1), 1–10.

122–144.
McGovern, A. et al. 2019 Making the black box more transparent: Understanding the physical implications of machinelearning. Bulletin of the American Meteorological Society
100 (11), 2175–2199.
Norman, M. R. 2021 A high-order WENO-limited finitevolume algorithm for atmospheric flow using the ADERdifferential transform time discretization. Quarterly J.
Royal Meteorological Society 147 (736), 1661–1690.
Norman, M. R. et al. 2021 Unprecedented cloud resolution ina GPU-enabled full-physics atmospheric climate simulationon OLCFs summit supercomputer. Int. J. HPC Applications
36(1), 93–105.
Partee, S. et al. 2021 Using machine learning at scale in HPC
simulations with SmartSim: An application to ocean climate modeling.
Piomelli, U. et al. 1988 Model consistency in large eddy sim-

6