Subgrid-Scale Parameterization in Burgers’ Equation Using Structure-Preserving Neural Networks and Entropy Variables Aijaz Nazir ∗,
Ilya Timofeyev†,
July 17, 2026
arXiv:2607.14855v1 [math.NA] 16 Jul 2026
Abstract We present a machine learning approach for developing subgrid-scale (SGS) parametrizations in coarse simulations of partial differential equations. We utilize structure-preserving neural networks and entropy variables to learn subgrid fluxes in coarse simulations of the Burgers’ equation. In particular, we employ a decoupled neural network architecture explicitly separating the subgrid corrections into two distinct components: a conservative Flux Potential network and an Eddy Viscosity network. We demonstrate that this reduced-order framework maintains high physical fidelity, accurately reproducing the energy spectrum, spatial and temporal correlation functions, and dynamical characteristics of the full-scale system. Furthermore, we show that our approach is robust and applicable to parameters outside the training regime.
Keywords: Burgers’ Equation, Entropy variables, Deep Learning, Sub-Grid Parametrization, Structure-preserving Neural Networks. MSC: 65M99, 68T07
1
Introduction
Developing models for unresolved degrees of freedom in coarse-grained simulations of fluid models has been an active research topic for many decades. The literature on this subject is extensive; some examples from computational fluid dynamics (CFD) include Large Eddy Simulations (e.g., [40, 25]), subgrid-scale viscosity models (e.g., [42, 17]), POD-ROM (Proper Orthogonal Decomposition Reduced-Order Models) approach (e.g., [35, 3, 38]), tensorial ROM for parametrized systems [33, 32, 31, 4, 5]. Various extensions of these techniques have also been proposed, which led to advances in computational methods for fluid dynamics. In addition, stochastic mode reduction [30] has been used to develop effective computational models for coarse variables for the 1D Burgers [15] and the 1D shallow-water [50] equation. Recently, neural networks have been used to develop parametrizations of subgrid terms (e.g. [28, 37, 6, 10, 8, 23]) with many examples demonstrating the validity and effectiveness of this approach. However, with a complete theory of neural networks still in the early stages of development, there is no guarantee of the mathematical properties ∗
Dept. of Mathematics, University of Houston, Houston, TX 77204, ORCID: 0000-0002-3494-4101, [email protected] † Dept. of Mathematics, University of Houston, Houston, TX 77204, ORCID: 0000-0002-3978-4047, [email protected]
1
of ML reduced models, such as stability, convergence, response to external perturbations, etc. Therefore, one of the most fruitful future research directions is to combine the advantages of rigorous mathematics and machine learning. Hyperbolic systems have a rich mathematical structure that has been extensively explored in the literature. Analytical advancements for hyperbolic systems have led to a better understanding of nonlinear wave phenomena and advances in computational methods. One important analytical concept in the theory of hyperbolic systems is the existence of the entropy function η(u), where u is the vector of dependent variables in the system (see e.g., [14]). This leads to the definition of entropy variables v = ∇u η(u) and entropy fluxes that provide additional analytical structure and an alternative symmetric formulation in the space of entropy variables. This symmetric structure plays an important role in both the theoretical analysis of hyperbolic systems and the design of stable numerical methods. In particular, the entropy function and entropy variables provide a natural framework for constructing entropy-consistent and entropy-stable discretizations. Recently, the entropy function and entropy variables have been used to develop neural network models to learn hyperbolic conservation laws from solution trajectories [11, 29]. In this work, we use a similar approach where we learn the entropy function and construct a neural network model to learn subgrid fluxes in coarse simulations of hyperbolic equations. In contrast to [11, 29], our focus is not on learning the underlying conservation law, but on adapting the methodology to learn the subgrid parametrization for unresolved physical processes in a coarse discretization of hyperbolic conservation laws. Moreover, we focus on long stationary simulations reproducing the energy cascade, the spatial and temporal correlation functions, and individual solutions of the hyperbolic system. Therefore, our work is related to CFD techniques mentioned at the beginning of this introduction, since our goal is to represent unresolved subgrid dynamics through effective neural network closures. The main contribution of this paper is the formulation and numerical validation of a structure-preserving neural-network closure framework for coarse simulations of hyperbolic conservation laws. In this work, we utilize a structure-preserving neural network to model the subgrid processes in coarse simulations of the Burgers’ equation. The main contribution of this paper is to introduce and systematically test a general structure-preserving framework for constructing data-driven subgrid parametrizations for coarse discretizations of hyperbolic conservation laws. The proposed framework combines machine learning with analytical properties of hyperbolic systems, particularly the entropy function, entropy variables, and flux-potential formulation. Although the present study focuses on the one-dimensional Burgers’ equation, the objective is not merely to construct a model for this particular equation, but to establish a methodology that can subsequently be extended to more general hyperbolic systems and multidimensional fluid-dynamical problems. In particular, we use three neural networks to model different components. First, we use an Input-Convex Neural Network (ICNN) [2] to learn the entropy function η(u); entropy variables can then be obtained by automatic differentiation. Second, we use a feed-forward neural network to learn the flux potential. Third, we use a feed-forward neural network with feature variables to develop a Smagorinsky-type [42] Eddy-Viscosity approximation. Details are presented in section 4.2. An important part of the contribution is the numerical evaluation of the complete framework rather than of its individual neural-network components in isolation. We test whether the resulting closure can simultaneously reproduce individual solution trajectories, long-time statistical properties, and the transfer of energy across resolved scales. We also assess the stability of the proposed
2
coarse-grid discretization and compare its performance with standard numerical and subgrid-scale models. We add large-scale forcing to perform long-time stationary simulations and show that our neural network reduced model accurately reproduces the energy spectra and the spatial and temporal correlation structure of the fully-resolved simulations. We also demonstrate that the resulting coarse mesh discretization is stable, accurately reproduces individual solutions, and does not require any additional stabilization (e.g., flux limiters). Our results demonstrate that incorporating the mathematical properties of hyperbolic systems directly into a machine learning approach yields an efficient, stable, and accurate subgrid modeling framework. Thus, the principal contribution of this work is both the formulation of the entropy-based neural closure framework and its validation across several complementary measures of coarse-model performance. The rest of the paper is organized as follows. In Section 2, we introduce the Burgers’ equation and numerical discretization. In Section 3, we discuss the problem formulation, including the definition of coarse variables, a brief overview of the entropy function and entropy variables, and a summary of the neural network modeling approach. Section 4 presents details about the data preparation, network architecture, and training. Numerical results are presented in Section 6. Conclusions are discussed in Section 7.
2
Problem formulation
2.1
Burgers’ Equation with stochastic forcing
The Burgers’ equation [9] is one of the well-studied nonlinear partial differential equations, often used as a prototype model for turbulence. We consider the stochastically forced Burgers’ equation over a periodic domain x ∈ [0, L] with L = 2π given by ∂ ∂ u+ f (u) = ρ(x, t) ∂t ∂x
(1)
with the nonlinear flux f (u) = u2 /2 and u ≡ u(x, t). Here, ρ denotes a large-scale stochastic forcing term, defined as X ρ(x, t) = A [αk (t) cos(kx) + βk (t) sin(kx)] , (2) k∈K
where A represents the forcing amplitude and L is the length of the spatial domain. The coefficients αk and βk are random variables evolving according to the AR(1) (autoregressive model of order 1) process αk (t + ∆t) = ψαk (t) + σϵk,1 (t), βk (t + ∆t) = ψβk (t) + σϵk,2 (t), (3) where 0 < ψ = 1 − γ∆t with γ > 0 and σ > 0 are AR(1) parameters, and ϵk,i (t), i = 1, 2 are i.i.d. Normal N (0, ∆t) random variables. In this paper we use γ = 1 and σ = 1.41. The forcing magnitude is A = 1, unless otherwise specified. Coefficients αk and βk are time-correlated Normal N (0, σ 2 /(1 − ψ 2 )) random variables. Equation (3) can be viewed as a temporal discretization of the Ornstein-Uhlenbeck process. In this paper, we use K = {1, 2, 3}, i.e., the first three wavenumbers are forced. The stochastic forcing ρ(x, t) is introduced to counteract the numerical dissipation inherent in finite-volume schemes.
3
2.2
Space–time discretization
Next, we consider a finite-volume discretization of the PDE in (1). The spatial domain is discretized using a uniform finite-volume mesh with grid spacing ∆x = L/Nf . The computational cells are defined as Ci = [xi−1/2 , xi+1/2 ], i = 0, . . . , Nf − 1, where the cell interfaces are located at xi−1/2 = i∆x and xi+1/2 = (i + 1)∆x. The cell centers are given by xi = (i + 1/2)∆x. This defines the fine-scale computational mesh Mf = {xi = i∆x | i = 0, 1, . . . , Nf − 1}. The semi-discrete system becomes fi+1/2 − fi−1/2 d ui = − + ρi , (4) dt ∆x where ∆x = L/Nf , ui is the average over the cell Ci , and fi+1/2 is a suitable discretization of the flux function. There are various options for the numerical flux fi+1/2 with particular details found in standard literature on finite volume and finite difference methods for conservation laws (e.g., [26, 21]). In this paper, we adopt the Local Lax-Friedrichs flux (LLF) [39] given by fi+1/2 =
f (ui+1 ) + f (ui ) λi+1/2 − (ui+1 − ui ), 2 2
(5)
where λi+1/2 is the local bound for the maximum wave speed. In practice, λi+1/2 is often approximated using the Rusanov formula λi+1/2 = max(|f ′ (ui )|, |f ′ (ui+1 )|) = max(|ui |, |ui+1 |). For the temporal discretization, we employ Heun’s method, a second-order explicit strong stability preserving (SSP) [20] Runge–Kutta scheme, to advance the semi-discrete system in time. Let u = {ui (t), i = 0, . . . , Nf − 1}. Then the semi-discrete system in (4) can be schematically written as u̇ = g(u) + ρ(t). The Heun’s method can be written as ũn+1 = un + g(un )∆t,
un+1 = un +
1 g(un ) + g(ũn+1 ) ∆t + ρn , 2
where un denotes the numerical solution at time tn . We also include stochastic forcing ρn = {ρni , i = 0, . . . , Nf − 1} with ρni ≡ ρ(xi , tn ) into the second step.
3
Methodology
Fully resolved simulations are often prohibitively expensive, making it desirable to obtain a coarse description of the physical phenomena. To this end, it is necessary to develop an appropriate subgrid model that represents the interactions between the resolved and unresolved (subgrid) degrees of freedom. In this work, we formulate such a model by adapting an entropy-based, structurepreserving neural-network methodology to the problem of subgrid closure. The proposed framework combines three neural networks that learn the entropy function, the flux potential, and the eddy-viscosity component, respectively, thus incorporating analytical structure directly into the subgrid parametrization.
3.1
Coarse Variables and Deviations
We introduce an averaging operator and define coarse variables as local averages over q fine grid points, i.e., q(I+1)−1 1 X UI (t) = ui (t) (6) q i=qI
4
for each cell I in a coarse mesh Mc = {0, 1, 2, . . . , Nc − 1} with Nc = Nf /q. The unresolved degrees of freedom are then defined as the residuals or deviations yi = ui − UI(i) , where the index i for yi and ui refers to the coarse cell I where the fine cell i is located. The dynamics of UI can be obtained by averaging the fine-scale equation (4): d 1 UI = − dt q∆x
q(I+1)−1
X
fi+1/2 − fi−1/2
i=qI
1 + q
q(I+1)−1
X
ρi (t).
(7)
i=qI
Using the telescoping property, the sum in the equation above can be rewritten as q(I+1)−1
X
fi+1/2 − fi−1/2 = fq(I+1)−1/2 − fqI−1/2 .
i=qI
Next, since our goal is to remain within the flux-difference framework, we define the “True” coarse fluxes as T T FI+ FI− (8) 1/2 = fq(I+1)−1/2 , 1/2 = fqI−1/2 , P q(I+1)−1 1 ρi (t) where the superscript T stands for the ”True” flux. In addition, we also define ρU I (t) ≡ q i=qI as the average of the forcing term on the coarse mesh. The equation (7) is exact, but not closed, because fluxes fi±1/2 depend on small-scale variables. Note that fq(I+1)−1/2 represents the exact physical flux evaluated at the right boundary of the coarse cell I. Because it acts on the fine grid interface, it is a function of both the coarse variables and deviations. Substituting ui = yi + UI(i) into (5) and computing explicitly the coarse index I(i), the ”True” flux can be decomposed into a macroscopic component and a subgrid flux: T FI+ 1/2 = FI+1/2 (U) + GI+1/2 (U, y),
(9)
where
F (UI ) + F (UI+1 ) 1 2 = UI+1 + UI2 . (10) 2 4 The subgrid flux GI+1/2 (U, y) depends on both the coarse variables and deviations. The goal of subgrid modeling is to develop a suitable approximation of the subgrid flux expressed only in terms of the coarse variables, thus eliminating explicit dependence on the deviations. In this paper, we use structure-preserving neural networks to approximate the subgrid flux GI+1/2 . We assume that the forcing is slowly varying in space, and the averaged forcing ρU I (t) can be U well-approximated by the forcing evaluated at the mid-point of cell I, i.e., ρI (t) ≈ ρ(XI , t). Thus, we denote ρI (t) := ρ(XI , t) and omit the subscript U for the rest of the paper. FI+1/2 =
3.2
Entropy Variables and Symmetrization
Data-driven subgrid models consistent with the underlying physical properties of the corresponding partial differential equations would have a better chance of reproducing the properties of fullyresolved dynamics. Therefore, we design our neural network for estimating subgrid fluxes to be consistent with the theory of entropy-stable conservation laws (see e.g., [14, 18, 44, 45]), which we summarize briefly in this section. Consider the conservation law ∂ ∂ u+ f (u) = 0. (11) ∂t ∂x 5
Many conservation laws admit a strictly convex entropy function η(u). The existence of such a function implies that the system is symmetrizable. Moreover, the strict convexity of η ensures a one-to-one mapping between the conservative variables u and the entropy variables, defined as: v = η ′ (u).
(12)
Since η ′ (u) is strictly convex, then v = η ′ (u) is invertible, and we can define u(v). Transforming the governing equations into the dual space of entropy variables symmetrizes the system, since it defines a one-to-one change of variables (at least locally). The derivative u′ (v) is positive definite since the entropy function η(u) is strictly convex. Next, we can write fluxes in terms of entropy variables g(v) := f (u(v)) and recast the conservation law in the symmetric form using entropy variables ∂v ∂v u′ (v) + g ′ (v) = 0. ∂t ∂x We can also define the flux potential ϕ(u) = vf (u) − r(u) with r′ (u) = η ′ (u)f ′ (u). If we substitute the entropy variables into the flux potential, then we can rewrite the conservation law (11) as ∂ ∂ ′ u+ ϕ (v) = 0. ∂t ∂x
3.3
Neural Network Subgrid Modeling
In this work, we propose a hybrid formulation for the discrete dynamics of the coarse variable UI . Rather than utilizing a black-box neural network to approximate the right-hand side of the equation for coarse variables, we keep the flux formulation and decompose fluxes into three parts: a macroscopic component, a learned structurally conservative part, and a learned dissipative correction. We utilize the following function form for the reduced equations for coarse variables d 1 UI = − FI+1/2 − FI−1/2 + ρI (t), dt q∆x
(13)
where the numerical flux at the interface I + 1/2 is constructed to preserve the underlying physical structure of the conservation law: NN NN FI+1/2 = FI+1/2 + FI+ 1/2 (UI , UI+1 ) − DI+1/2 (UI , UI+1 ).
(14)
N N and D N N is described below. Here, FI+1/2 is given by equation (10) and the construction of FI+ 1/2 I+1/2 N N and Entropy Encoding. The structurally conservative part is Neural Subgrid Flux FI+ 1/2 estimated by learning two separate neural networks: the Entropy Neural Network ηθ and the Flux Potential Neural Network ϕθ . To ensure consistency, this term is calculated in the dual space of N N is computed as entropy variables v = η ′ (U ) [11, 29]. In particular, the flux FI+ 1/2 NN FI+ 1/2 =
1 ′ ϕθ (vI ) + ϕ′θ (vI+1 ) , 2
vI = ηθ′ (UI ),
where derivatives are computed using automated differentiation and θ denotes the network’s parameters. To mathematically guarantee that ηθ (U ) is a strictly convex entropy function, it is parameterized using an Input-Convex Neural Network (ICNN) [2]. Strict convexity (i.e., η ′′ > 0) is architecturally enforced by restricting the hidden-layer weights to the non-negative domain and 6
employing convex, non-decreasing activation functions. Additionally, a quadratic skip-connection is incorporated to maintain a strictly positive lower bound on the second derivative. Particular details about the architecture of ηθ (u) and ϕθ (v) are provided in sections 4.2.1 and 4.2.2, respectively. Dissipation DN N . The idea is similar to the Smagorinsky-type approximation [42]. Here, a neural network is used to learn the state-dependent viscosity coefficient. In particular, we use the following form for the dissipation term NN DI+ 1/2 = Cθ (ξ)
ΛI+1/2 × (UI+1 − UI ). 2
(15)
Here, ΛI+1/2 denotes the local bound for the maximum wave speed, defined as ΛI+1/2 = max(|UI |, |UI+1 |). The coefficient Cθ is predicted by an Eddy Viscosity Neural Network, which interrogates the local flow state through a feature vector ξ comprising the local average velocity and the amplified interfacial jump. We ensure that Cθ (ξ) is positive by applying the sigmoid function in the final layer. The sigmoid activation allows us to impose hard bounds on the viscosity coefficient, so that Cθ (ξ) ∈ [Cmin , Cmax ]. Unlike traditional schemes where the dissipation coefficient is a fixed constant, Cθ is a statedependent parameter learned through end-to-end backpropagation. This allows the model to perN N and C ) collaborate to resolve the form Residual Discovery, where the neural branches (FI+ θ 1/2 high-frequency discrepancies between the coarse macroscopic discretization and the filtered highfidelity Direct Numerical Simulation (DNS) reference. Details about the architecture of Cθ are provided in section 4.2.3.
4
Data Preparation, Network Architecture, and Training
We develop a reduced model for the local averages of the stochastically forced Burgers’ equation by introducing a structure-preserving hybrid ”Gray-Box” framework. Rather than relying on standard feed-forward neural networks to approximate the entire non-linear flux, we decompose the numerical flux into an analytical macroscopic component and learned structural and dissipative subgrid corrections. We note that standard mean-squared error formulations struggled to accurately capture shock dynamics without over-dissipating smooth waves. To address this limitation, we introduce a gradient-weighted loss function that emphasizes high-frequency interfacial jumps and train the model end-to-end using the total physical flux. Furthermore, the spurious oscillations that traditionally require flux limiters are natively suppressed through a state-dependent Eddy Viscosity Network, while consistency for nonlinear flux corrections is architecturally guaranteed via an Input-Convex Neural Network (ICNN). The performance of the reduced model is assessed based on how well equilibrium statistical properties are recovered compared to those of the Direct Numerical Simulation (DNS). In particular, we illustrate that the reduced model is able to accurately reproduce the spectral energy decay and the spatial and temporal correlation functions of the large-scale variables. The response of the reduced model with respect to changes in the amplitude of the forcing is also analyzed. In this section, we describe how the training dataset is generated, the structure of the threebranch neural architecture, and the latent learning procedure used to approximate the subgrid physics and resolve the dynamics on the coarse mesh.
7
4.1
Data Preparation
We conduct a high-resolution numerical simulation of the full model (4), sampling the solution at intervals of 100∆t with ∆t = 0.001. Thus, each run generates 10000 snapshots. To ensure stability and remove transient effects, the initial 5000 model time units are discarded. The initial velocity field is defined as a superposition of the first two Fourier modes, 2πx 4πx u(x, 0) = 0.1 sin + ϕ1 + 0.1 sin + ϕ2 , L L where L = 2π and the phase shifts ϕ1 , ϕ2 ∼ N (0, 1) are independently sampled from a normal distribution. Time series of variables U0 and U1 , the corresponding ”true” flux F1T/2 , and the smoothness indicator β1/2 (computed for coarse variables) are recorded. We generate three datasets for three different initial conditions, filter them as described below, and then combine them into one balanced dataset. The final dataset contains approximately 152K samples. Details about training and validation samples are provided in Table 1. Dataset Filtering. Smoothness indicators, introduced in WENO schemes [22], provide a quantitative measure of local solution regularity. For Burgers’ dynamics, the indicator β identifies regions of nonlinear steepening and shock formation. Thus, it enables targeted selection of dynamically active zones for model training. In this work, local smoothness is quantified using a standard WENO-type smoothness indicator. For a discrete solution UI , the indicator βI is defined as βI =
13 1 (UI−1 − 2UI + UI+1 )2 + (UI−1 − UI+1 )2 . 12 4
(16)
The first term approximates the local curvature of the solution through a second-order finite difference, while the second term measures the squared first-order gradient. In several previous works [1, 34, 46], we determined that near-shock regions are most challenging for subgrid modeling. Therefore, to emphasize regions with significant subgrid activity, the dataset is filtered using the smoothness indicator β, which quantifies the local variability of the solution. A threshold of β ∗ = 0.25 is applied, separating the data into high-β (active subgrid regions) and low-β (smoother regions) subsets. Samples with β ≥ β ∗ , which correspond to shocks, are retained in full. We remove 70% of the low-β samples from the dataset. This selective downsampling prevents the over-representation of low-β values, preserving diversity in the dataset while ensuring sufficient representation of higher-β values. This approach ensures that the training dataset effectively captures different behaviors, allowing the neural network to generalize well across different conditions.
4.2
Network Architecture and Training
We employ a structure-preserving hybrid neural architecture to learn the subgrid structural correction and dynamic dissipation as functions of the local state, which differs from the usual strategy of learning global solutions directly. The goal is to achieve a local, physics-informed approximation, where the total numerical flux at a specific interface I + 1/2 is determined only by the adjacent cell averages UI and UI+1 . Thus, our flux representation relies on a compact two-point stencil; this is possible by using the entropy-consistent framework. Due to spatial homogeneity, it is sufficient to use data from a representative local interface (at 1/2 in our case), allowing for a smaller dataset and efficient parallel training. We then utilize the same neural network to compute subgrid fluxes FI+1/2 for I = 0, . . . , Nc − 1 in simulations of the reduced dynamics (13). 8
The model architecture is implemented in JAX/Flax [7] and consists of three neural networks, as shown schematically in Figure 1. First, the Entropy Neural Network parameterizes the convex entropy function η using an Input-Convex Neural Network. It utilizes Softplus activations, nonnegative weight constraints enforced via absolute value operations, and a quadratic skip-connection to guarantee a strictly positive-definite Hessian. Second, the Flux Potential Neural Network utilizes smoothly differentiable Tanh activations to map the dual entropy variables to the flux potential. Finally, the Eddy Viscosity Neural Network employs Swish activation [36] and processes an augmented feature vector containing the local average velocity and an amplified interfacial jump, i.e., UI + UI+1 ξ= , 10|UI − UI+1 | . 2 Its final layer applies a scaled Sigmoid function to restrict the dissipation coefficient Cθ (ξ) strictly within [Cmin , Cmax ] = [0.35, 2.0]. Analysis of the neural network numerical scheme in Section 5 provides analytical guidelines for selecting constants Cmin and Cmax .
Figure 1: Schematic representation of the three branches for the structure-preserving ”Gray-Box” computational pipeline. The flux FI+1/2 (middle) is evaluated using formula (10). Structural corrections (top) are restricted by learning a convex network ηθ (U ) and the flux potential ϕθ (v). State-dependent dissipation (bottom) is implemented by learning the viscosity coefficient Cθ (ξ) and multiplying it by the LLF-type term (equation (15)). Blue boxes are modeled using different neural networks. Green boxes are given by explicit formulas. To train this hybrid architecture, the model is optimized directly on the high-fidelity DNS dataset. Rather than relying on standard mean squared error, which frequently struggles to balance shock capturing with smooth wave preservation, we introduce a Gradient-Weighted Loss function X T 2 L= wI+1/2 ||FI+ with wI+1/2 = 1.0 + γ |UI − UI+1 |, 1/2 − FI+1/2 ||2 batch
where the summation is taken over a batch of samples. The loss evaluates the Mean Squared Error (MSE) between the assembled neural flux prediction and the ”True” coarse flux computed from fine-mesh simulations. To aggressively target high-frequency interfacial jumps without overdissipating smooth regions, the per-sample squared error is scaled by a weight factor wI+1/2 , where 9
γ serves as a tunable hyperparameter controlling the sensitivity to local gradients; the value γ = 0.2 works well in this case. While offline validation metrics naturally favor γ → 0 due to the reduced penalty magnitude, our a posteriori analysis demonstrates that under-regularizing the network removes critical physical constraints, causing under-damped high-frequency errors during dynamic integration. Conversely, excessive weighting (γ ≥ 0.4) introduces severe optimization stiffness. Thus, γ = 0.2 provides the minimum effective regularization necessary for long-term shock stability. Network parameters are optimized via latent-space learning using the Adam optimizer. We implement a cosine decay learning-rate schedule with an initial rate of 10−3 , allowing rapid initial convergence followed by fine-tuning of the network weights. The dynamic cosine decay schedule decreases the learning rate smoothly from an initial value to a prescribed minimum value according to a cosine curve, providing a gradual transition from aggressive parameter updates to smaller corrective updates during training. The model is trained for 200 epochs using mini-batches of size 128, iteratively minimizing the gradient-weighted loss on the training set. Validation is carried out by performing a long simulation of the reduced equation on the coarse grid and comparing spectra and individual solutions. In the next sections, we use a generic notation z (·) and σ(·) to denote hidden layers and activation functions; these design features are different in each section. A comprehensive summary of the physical parameters, architectural capacity, and training configurations for the proposed model is provided in Table 1. Category
Parameter
Value
Physical Setup
Fine Grid Resolution (Nf ) Coarse Grid Resolution (Nc ) Domain Length (L) Time Step (∆t)
512 64 2π 0.001
Dataset Specifications
Training Samples (80%) Validation Samples (20%)
121,832 30,459
Entropy NN (η)
Activations 2 Hidden Layers Activations 1 Hidden Layer Activations 1 Hidden Layer Dynamic Viscosity Bounds
Softplus 16, 16 Tanh 16 Swish, Sigmoid 16 [0.35, 2.0]
Optimizer Learning Rate Schedule Initial Learning Rate Gradient Regularization Weight (γ) Batch Size Epochs
Adam Cosine Decay 1 × 10−3 0.2 128 200
Flux Potential NN (ϕ) Eddy Viscosity NN (ν)
Training & Optimization
Table 1: Summary of physical and neural network hyperparameters for the proposed architecture.
10
A sensitivity analysis on hidden neurons (N ∈ {16, 32, 64}) for the eddy viscosity NN revealed a flat plateau in validation error. This indicates that the underlying physical manifold of the subgrid flux is fully represented by narrower architectures. To maximize computational acceleration relative to the DNS baseline, we selected the parsimonious 16-neuron architecture, strictly minimizing matrix-multiplication overhead during the PDE time-stepping loop. 4.2.1
Input Convex Neural Network (ICNN) for Entropy Formulation
To ensure that the learned discrete conservation law remains hyperbolic, the parameterized entropy function ηθ (u) must be strictly convex with respect to the state variable u. We achieve this by structuring the entropy network as an Input Convex Neural Network (ICNN) [2]. ICNNs enforce overall convexity by satisfying two structural constraints: 1. The nonlinear activation functions must be both convex and non-decreasing. 2. All weight matrices connecting the hidden layers must be strictly non-negative. Since we consider a scalar equation, the input is u ∈ R1 . We satisfy the first constraint by employing the smooth softplus activation function σ(x) = ln(1 + exp(x)). In the first hidden layer, the input undergoes a standard affine transformation followed by the activation: z (1) = σ Wu(1) u + b(1) , (1)
where Wu and b(1) are the standard weights and biases. For the second hidden layer, the network utilizes a skip-connection, directly re-injecting the input u alongside the features from the previous layer. To satisfy the second convexity constraint, (1) the hidden-to-hidden weight matrix Wz is restricted to the non-negative orthant. In our implementation, this is achieved by taking the element-wise absolute value of the trainable parameters, (1) (1) defining the effective weights as W̃z = |Wz |. Thus, the second hidden layer evaluates as: z (2) = σ z (1) W̃z(1) + Wu(2) u + b(2) . The dimension of each hidden layer is 16. To maximize the expressive power of the ICNN, the final output layer is augmented with both linear and quadratic transformations of the original input u. The quadratic weight coefficient, W̃S , is also constrained to be non-negative via an absolute value projection (W̃S = |WS |). The final convex entropy function is computed as: η(u) = z (2) W̃z(2) + b(3) + Wu(3) u + W̃S u2 . Since W̃S ≥ 0 and u2 is a convex function, this quadratic augmentation strictly preserves the overall convexity of the network. This absolute-value bounding strategy provides a robust guarantee of hyperbolicity without the need for complex regularization penalties during training. Note that the standard entropy function for the Burgers equation is quadratic, η(u) = 21 u2 . Since the goal of our framework is to learn an entropy function appropriate for the coarse discretization, the learned entropy function need not coincide exactly with the standard entropy of the Burgers equation. However, including the quadratic term in the output layer provides the network with the exact 11
Burgers entropy as a baseline, allowing it to learn only the additional corrections required by the coarse model. Therefore, this augmentation, most likely, improves the efficiency of training the entropy network. For other hyperbolic systems, it may be advantageous to augment the output layer with a known convex entropy function of the underlying equations. 4.2.2
Flux Potential Network
While the entropy network ηθ (u) guarantees the stability of the system, the dispersive wave-speed corrections are governed by the conservative subgrid flux, defined as the gradient of a scalar potential: ϕ′θ (v). To parameterize this potential function ϕθ , we use a standard Fully Connected Neural Network (FCNN) mapping the entropy variables v to a scalar output. For this specific implementation, the architecture is kept deliberately lightweight. The input v is passed through a single hidden layer of dimension 16 paired with a hyperbolic tangent tanh activation function. The hidden state is computed as: z (1) = tanh Wv(1) v + b(1) . This is followed by a linear output layer to predict the scalar potential: ϕθ (v) = Wv(2) z (1) + b(2) . Taking the exact derivative of this compact network with respect to its inputs provides the conservative flux correction in the reduced equation. 4.2.3
Eddy Viscosity Network
To suppress spurious oscillations without over-damping the physical energy cascade, we introduce an independent Eddy Viscosity network, Cθ . Rather than replacing the baseline numerical dissipation (such as the Local Lax-Friedrichs scheme), this network outputs a dynamic scaling coefficient that modulates the baseline viscosity based on localized feature inputs ξ (e.g., cell averages, jumps, or curvature). The network processes these input features through a single hidden layer of 16 neurons utilizing the swish activation function [36], which ensures smooth, non-vanishing gradient flow during training: (1) (1) z (1) = swish Wf ξ + bf . A critical design element of this network is the strict physical bounding of its output. To prevent artificial energy pileup at the smallest grid scales due to insufficient dissipation, the final layer utilizes a sigmoid activation function σ(x) which is explicitly scaled and shifted: (2) (2) Cθ = 1.65 σ Wf z (1) + bf + 0.35. Since the sigmoid function is strictly bounded by (0, 1), this formulation guarantees that the network output operates exclusively within the range Cθ ∈ [0.35, 2.00]. This architecture ensures that at least 35% of the baseline numerical viscosity is always maintained to suppress spurious high-frequency oscillations.
12
5
Analytical Properties of the Reduced Model
In this section, we describe several analytical properties of the hybrid scheme with flux (14). In particular, we derive conditions under which the numerical scheme is total variation diminishing and discuss entropy stability. The total variation diminishing (TVD) property is a key stability criterion for numerical schemes for hyperbolic conservation laws. A TVD scheme guarantees that the total variation of the numerical solution does not increase in time, thereby suppressing the growth of spurious oscillations near discontinuities. Definition 1 (Discrete total variation). Let {uj }, j = 0, . . . , J be a discrete sequence. Its total variation is defined by J−1 X T V (u) = |uj+1 − uj |. j=0
The following definition introduces the TVD property of a discrete numerical solution un = {un0 , . . . , unJ }, where n denotes the time-step index. Definition 2 (Total variation diminishing scheme). Let unj denote the numerical approximation at time-step n, i.e., unj ≈ u(n∆t, j∆x). A numerical scheme is called total variation diminishing, or TVD, if the total variation of the numerical solution does not increase with time, i.e., T V (un+1 ) ≤ T V (un ) for every n. First, we need to make an assumption about the regularity of the Neural Subgrid Flux. Assumption 1 (Lipschitz regularity of the Neural Subgrid Flux). Consider the Neural Subgrid NN = N (U , U 2 Flux FI+ I I+1 ), where N : R → R is a continuous function. Assume that for each 1/2 admissible interval B = [umin , umax ], N (a, b) is Lipschitz in each argument separately; that is, there is a constant Lnn (possibly Lnn = Lnn (B)) such that N (a, b) − N (a, c) ≤ Lnn |b − c|
N (a, b) − N (c, b) ≤ Lnn |a − c|,
(17)
for all a, b, c ∈ B. Next, we derive TVD conditions for the effective NN scheme and discuss consequences for designing the neural network. The main idea is that the neural corrective flux may introduce additional variation into the numerical solution, while the eddy-viscosity term acts to damp this variation. The TVD conditions derived below quantify the amount of numerical dissipation required to compensate for the corrective flux and guarantee that the total variation does not increase from one time step to the next. Proposition 1 (Sufficient TVD conditions for the neural network-augmented scheme). Let the spatial domain be discretized with uniform mesh spacing ∆x and time step ∆t, and define r = ∆t/∆x. Consider the conservative update UIn+1 = UIn − r FI+1/2 − FI−1/2 , (18) 13
where the composite numerical flux at the interface I + 1/2 is NN NN FI+1/2 = FI+1/2 + FI+ 1/2 (UI , UI+1 ) − DI+1/2 .
(19)
Here
1 f (UI ) + f (UI+1 ) , 2 and the dynamic eddy-viscosity flux is FI+1/2 =
NN FI+ 1/2 = N (UI , UI+1 ),
αI+1/2 (UI+1 − UI ), 2
NN DI+ 1/2 = Cθ,I+1/2
(20)
where Cθ,I+1/2 = Cθ (ξI+1/2 ) is evaluated from the known numerical solution at time level n and treated as an interface coefficient during the update. Assume that 0 < Cmin ≤ Cθ,I+1/2 ≤ Cmax < ∞. Let ∆UI+1/2 = UI+1 − UI and let ΛI+1/2 = max{|f ′ (UI )|, |f ′ (UI+1 )|} be a local wave-speed bound. To avoid degeneracy of the diffusion (20) when ΛI+1/2 = 0, assume that the viscosity uses a positive speed floor, i.e., αI+1/2 = max(ΛI+1/2 , ε),
ε > 0.
Assume also that the neural corrective flux satisfies Assumption 1 on the admissible range, with Lipschitz constant Lnn . The Lipschitz constant Lnn measures the sensitivity of the neural corrective flux to perturbations in its inputs. Consequently, larger values of Lnn require stronger numerical dissipation in order to preserve the TVD property. If, for every interface I + 1/2, Cθ,I+1/2 αI+1/2 ≥ ΛI+1/2 + 2Lnn ,
(21)
and the time step satisfies r Cθ,I+1/2 αI+1/2 + 2Lnn ≤ 1
for every interface I + 1/2,
(22)
then the scheme (18)–(19) satisfies Harten’s TVD conditions. The first condition in (21) requires that the adaptive eddy viscosity be large enough to dominate both the physical wave propagation and the variation introduced by the neural corrective flux. The second condition (22) is a CFL-type restriction that ensures the coefficients in Harten’s incremental representation satisfy the bounds required for the TVD property. The proof proceeds by rewriting the scheme in Harten’s incremental form. The neural corrective flux is first expressed in terms of neighboring solution differences using its Lipschitz continuity, after which the resulting coefficients are compared with Harten’s sufficient conditions for TVD schemes. Proof. The proof consists of expressing each contribution to the numerical flux in terms of neighboring solution differences so that the update can be written in Harten’s incremental form. Once this representation is obtained, the TVD property follows by verifying Harten’s coefficient conditions. 14
Define νI+1/2 = Cθ,I+1/2 αI+1/2 . Then
νI+1/2 ∆UI+1/2 . 2 We now express each component of the numerical flux in terms of the neighboring solution differences ∆UI±1/2 . This allows the update to be written in Harten’s incremental form. First, for the macroscopic central flux, NN DI+ 1/2 =
1 1 FI+1/2 − FI−1/2 = aI+1/2 ∆UI+1/2 + aI−1/2 ∆UI−1/2 , 2 2
(23)
where |aI+1/2 | ≤ ΛI+1/2 from the definition of ΛI+1/2 . Next, the neural corrective flux is treated similarly. Since it depends on two neighboring states, we use the Lipschitz regularity assumption to express its difference in terms of the same neighboring increments. We write NN NN FI+ 1/2 − FI−1/2 = N (UI , UI+1 ) − N (UI−1 , UI ). Adding and subtracting N (UI , UI ) gives NN NN FI+ 1/2 − FI−1/2 = N (UI , UI+1 ) − N (UI , UI ) + N (UI , UI ) − N (UI−1 , UI ) . Thus there exist coefficients qI+1/2 and pI−1/2 such that NN NN FI+ 1/2 − FI−1/2 = qI+1/2 ∆UI+1/2 + pI−1/2 ∆UI−1/2 ,
with qI+1/2 =
N (UI , UI+1 ) − N (UI , UI ) , ∆UI+1/2
pI−1/2 =
N (UI , UI ) − N (UI−1 , UI ) , ∆UI−1/2
(24)
(25)
so that |qI+1/2 | ≤ Lnn ,
|pI−1/2 | ≤ Lnn
(26)
using Assumption 1. The coefficients qI+1/2 and pI−1/2 are difference quotients of N with respect to its second and first arguments, respectively. When the denominator is zero, the associated numerator is also zero and the coefficient may be set equal to zero. The eddy-viscosity contribution satisfies NN NN DI+ 1/2 − DI−1/2 =
νI+1/2 νI−1/2 ∆UI+1/2 − ∆UI−1/2 . 2 2
(27)
Substituting (23), (24), and (27) into (18), we obtain − UIn+1 = UIn + Φ+ I+1/2 ∆UI+1/2 − ΦI−1/2 ∆UI−1/2 .
(28)
Equation (28) has the standard incremental form considered by Harten. The coefficients Φ± determine whether the update is TVD. The coefficients are defined as r (29) Φ+ I+1/2 = 2 νI+1/2 − aI+1/2 − 2qI+1/2 , r Φ− (30) I+1/2 = 2 νI+1/2 + aI+1/2 + 2pI+1/2 . 15
By Harten’s theorem, the incremental scheme (28) is TVD provided that, for every interface, Φ+ I+1/2 ≥ 0,
Φ− I+1/2 ≥ 0,
− Φ+ I+1/2 + ΦI+1/2 ≤ 1.
(31)
It therefore remains to verify that the coefficients Φ± satisfy the three inequalities in (31). We estimate each coefficient using the wave-speed bound and the Lipschitz continuity of the neural corrective flux. Using |aI+1/2 | ≤ ΛI+1/2 and |qI+1/2 | ≤ Lnn , condition (21) gives νI+1/2 − aI+1/2 − 2qI+1/2 ≥ νI+1/2 − ΛI+1/2 − 2Lnn ≥ 0. Therefore Φ+ I+1/2 ≥ 0. Similarly, using |pI+1/2 | ≤ Lnn , νI+1/2 + aI+1/2 + 2pI+1/2 ≥ νI+1/2 − ΛI+1/2 − 2Lnn ≥ 0, and hence Φ− I+1/2 ≥ 0. The final TVD condition concerns the sum of the incremental coefficients. Adding (29) and (30), the macroscopic advective terms cancel: − Φ+ + Φ = r ν 1/2 + pI+1/2 − qI+1/2 . I+ 1 1 I+ /2 I+ /2 Using (26), we obtain − Φ+ I+1/2 + ΦI+1/2 ≤ r νI+1/2 + 2Lnn . Therefore condition (22) implies − Φ+ I+1/2 + ΦI+1/2 ≤ 1.
All three Harten conditions (31) hold, and the scheme is TVD. When ΛI+1/2 > ε, condition (21) can be rewritten as Cθ,I+1/2 ≥ 1 +
2Lnn . ΛI+1/2
Inequalities (21) and (22) have direct consequences for the architecture of the network. Recall that we use a rescaled sigmoid function in the final layer of the dissipation network to ensure that Cmin ≤ Cθ,I+1/2 ≤ Cmax . Thus, the inequality above indicates that it is important to select Cmax > 1 to make the TVD property feasible. At the same time, conditions (21) and (22) are sufficient, not necessary, and hence they should be interpreted as conservative analytical guidelines. The role of the eddy-viscosity term is primarily to stabilize the scheme near shocks and steep gradients; in smooth regions, excessive artificial diffusion may reduce accuracy. In our numerical experiments, we found that imposing Cmin = 1 makes the scheme too diffusive. Therefore, we choose Cmin < 1 and allow the neural network to increase Cθ,I+1/2 adaptively in regions where additional viscosity is needed. The lower bound Cmin is also controlled by the next proposition. Condition (22) can also be checked a posteriori during the numerical simulation to verify that the time step is sufficiently small. We note that strict TVD enforcement is not a necessary condition for a high-resolution scheme to perform well in practice. High-order methods such as WENO [41] and Runge–Kutta Discontinuous Galerkin [12] are not strictly TVD [43, 47], which is consistent with the theoretical barrier established by Godunov’s theorem [19]. The conditions derived in Proposition 1 are sufficient but not 16
necessary: the scheme may satisfy the TVD criterion locally depending on the network output Cθ , but is not required to do so globally in order to remain stable and accurate. The TVD property controls the growth of oscillations in the numerical solution, whereas entropy stability guarantees consistency with the physical entropy inequality satisfied by weak solutions. We next derive a sufficient condition under which the neural corrective flux and the adaptive eddyviscosity term together produce a discretely entropy-stable scheme. Proposition 2 (Discrete entropy stability of the neural network-augmented scheme). Consider the semi-discrete conservative scheme for the Burgers equation 1 2 ut + u =0 2 x with entropy pair 1 1 q(u) = u3 , η(u) = u2 , 2 3 ′ entropy variable v = η (u) = u, and entropy potential 1 ψ(u) = vf (u) − q(u) = u3 . 6 The semi-discrete scheme is
dUI 1 =− FI+1/2 − FI−1/2 , dt ∆x
(32)
NN NN FI+1/2 = FI+1/2 + FI+ 1/2 − DI+1/2 .
(33)
2 UI2 + UI+1 FI+1/2 = 4
(34)
where Here
is the central macroscopic flux, and NN DI+ 1/2 = Cθ,I+1/2
αI+1/2 (UI+1 − UI ), 2
(35)
where ΛI+1/2 = max{|UI |, |UI+1 |},
αI+1/2 = max{ΛI+1/2 , ε},
ε > 0.
Also, assume that the neural corrective flux is constructed using the convex Entropy Neural Network ηθ and the Flux Potential Neural Network ϕθ as follows 1 ′ NN N (UI+1 , UI ) := FI+ ϕθ (ZI ) + ϕ′θ (ZI+1 ) , ZI = ηθ′ (UI ). (36) 1/2 = 2 Note that here we define the ”internal” entropy variables ZI that, in general, are not the same as the entropy variables for the Burgers equation. If, for every interface, |UI+1 − UI | Cθ,I+1/2 αI+1/2 ≥ + 2Lnn , (37) 6 then the semi-discrete scheme satisfies a discrete entropy inequality QI+1/2 − QI−1/2 d η(UI ) + ≤ 0, dt ∆x for a suitable numerical entropy flux QI+1/2 . 17
(38)
The proof follows Tadmor’s entropy stability framework. We first identify the entropy-conservative reference flux and then estimate the entropy production associated with each component of the numerical flux. The Lipschitz regularity assumption is used to bound the entropy production contributed by the neural corrective flux. The adaptive eddy viscosity is finally shown to dominate the remaining positive entropy production. Proof. Following Tadmor’s framework, we compare the composite numerical flux with the entropyconservative flux for Burgers equation. The entropy production associated with each component of the numerical flux is then estimated separately. For the Burgers equation, the physical flux, entropy function, and entropy flux are given by 1 f (u) = u2 , 2
1 η(u) = u2 , 2
1 q(u) = u3 , 3
respectively. The entropy variable is v = u, and the entropy potential is 1 ψ(u) = vf (u) − q(u) = u3 . 6 EC By Tadmor’s condition [44], a two-point numerical flux FI+ 1/2 is entropy-conservative if and only if it satisfies EC VI+1 − VI FI+ (39) 1/2 = ψI+1 − ψI .
Therefore, using entropy variables and the entropy potential for the Burgers equation, VI = UI , and 2 UI2 + UI UI+1 + UI+1 ψ(UI+1 ) − ψ(UI ) EC = . FI+ = 1/2 UI+1 − UI 6 A scheme satisfies (38) if and only if the total entropy production at each interface, EC := EI+1/2 (VI+1 − VI ) FI+1/2 − FI+1/2 ,
(40)
satisfies EI+1/2 ≤ 0. We evaluate the contribution of each component of (33) to EI+1/2 next. Macroscopic Flux. We first consider the contribution of the central macroscopic flux. Since this flux is not entropy conservative, it generates positive entropy production that must later be compensated by numerical dissipation. The difference between the central macroscopic flux and the entropy-conservative flux is EC FI+1/2 − FI+ 1/2 =
2 2 U 2 + UI UI+1 + UI+1 UI2 + UI+1 (UI+1 − UI )2 − I = . 4 6 12
Thus the entropy production associated with the central macroscopic flux is (U 3 I+1 − UI ) EC (UI+1 − UI ) FI+1/2 − FI+ = . 1/2 12
(41)
Neural Corrective Flux. Next, we estimate the entropy contribution of the neural corrective flux. Unlike the macroscopic flux, this contribution depends on the learned neural approximation and is controlled using the Lipschitz regularity assumption. From (36), the Neural Corrective Flux satisfies the neural entropy production equation NN RI+1/2 := (VI+1 − VI ) FI+ 1/2 = (UI+1 − UI )
18
ϕ′θ (ZI ) + ϕ′θ (ZI+1 ) . 2
(42)
Next, using the definition of the neural flux, we can write ϕ′θ (Z(U )) = N (U, U ), and using Assumption 1 we can bound the difference |ϕ′θ (Z(UI+1 )) − ϕ′θ (Z(UI ))| = |N (UI+1 , UI+1 ) − N (UI , UI )| ≤ 2Lnn |UI+1 − UI |. Therefore, |RI+1/2 | ≤ Lnn (UI+1 − UI )2 . (43) Dissipative Flux and Entropy Production. Combining the previous estimates, the total entropy production of the composite numerical flux is bounded by (UI+1 − UI )3 NN − (UI+1 − UI )DI+ 1/2 . 12 Treating the worst-case scenario when the neural entropy production is RI+1/2 > 0, the Tadmor entropy stability condition is therefore RI+1/2 +
(UI+1 − UI )3 + |RI+1/2 |. (44) 12 1 NN 2 Using (35), we obtain (UI+1 − UI )DI+ 1/2 = 2 Cθ,I+1/2 αI+1/2 (UI+1 − UI ) . The eddy-viscosity term therefore compensates both the entropy production of the central flux and the worst-case positive contribution of the neural corrective flux. Thus a sufficient pointwise condition for (44) is NN (UI+1 − UI )DI+ 1/2 ≥
αI+1/2 |UI+1 − UI |3 (UI+1 − UI )2 ≥ + Lnn (UI+1 − UI )2 , 2 12 which is equivalent, for UI+1 ̸= UI , to Cθ,I+1/2
|UI+1 − UI | + 2Lnn . (45) 6 If UI+1 = UI , the inequality is trivially satisfied using the definition of RI+1/2 in (42). Hence the adaptive dissipation dominates all possible positive entropy production, implying the discrete entropy inequality (38). Cθ,I+1/2 αI+1/2 ≥
Condition (37) requires the adaptive eddy viscosity to provide enough dissipation to offset both the entropy generated by the central discretization and the additional entropy production that may arise from the neural corrective flux. Consider the lower bound with RI+1/2 ≤ 0. Then, for ΛI+1/2 ≥ ε, the condition (37) is equivalent to |UI+1 − UI | Cθ,I+1/2 ≥ 6 ΛI+1/2 with ΛI+1/2 = max(|UI+1 |, |UI |). The triangle inequality gives |UI+1 − UI | ≤ |UI | + |UI+1 | ≤ 2ΛI+1/2 . Thus, the standard sufficient condition for discrete entropy stability of the central flux approximation becomes 2ΛI+1/2 1 = . Cθ,I+1/2 ≥ 6ΛI+1/2 3 Thus, condition (37) (same as eq. (45)) is essentially equivalent to Cθ,I+1/2 ≥
1 2Lnn + . 3 ΛI+1/2
Therefore, we select Cmin = 0.35 as the lower bound for the eddy viscosity network. Numerically, as shown in the next section, this bound is sufficient to ensure that no oscillations develop near shocks. 19
6
Numerical Results
In this section, the performance of the neural network-based reduced model is evaluated by comparing its stationary statistical properties with those of the full model. We perform long stationary simulations for the full and reduced models and perform a detailed comparison of several statistical properties. The time step for the simulations of the reduced models is consistent with the time step used in the numerical integration of the full model. For our numerical experiments with the Burgers’ equation, we used the fine-mesh discretization Nf = 512 in simulations of the full model. This high-resolution grid ensures that all relevant physical features are captured. The numerical diffusion in this fine-mesh setup is negligible, meaning further refinement would yield minimal improvements. We generate coarse-grid representations with q = 8. Simulations on the fine-mesh grid are referred to as DNS (Direct Numerical Simulation), serving as a reference for comparison. Reduced model simulations are denoted as LLF-64 for the Local Lax-Friedrichs method with 64 points, and as NN-64 when the neural network subgrid model is used. Since the neural network framework is structure-preserving, there is no risk that the NN-reduced model may generate non-physical behaviors. Energy Spectra. To quantify the model’s statistical accuracy, we evaluate the distribution of kinetic energy across different spatial scales. The energy spectrum is computed in Fourier space by applying the discrete Fourier transform (DFT) to the state variable u(x, t) at each recorded time step, yielding the Fourier coefficients ûk (t) for each wavenumber k. The time-averaged spectral energy is calculated as: Z NT 1 T 1 X Ek = |ûk (t)|2 dt ≈ |ûk (ti )|2 T 0 NT i=1
We generate a long, stationary time series to ensure statistical convergence. The system is integrated for T = 1000 using a time step of ∆t = 0.001. To compute the time-averaged spectra, we sample data with ∆tsample = 0.1, yielding NT = 10, 000 discrete snapshots. The total macroscopic P energy of the system is the sum of the spectral energy over all wavenumbers: Etotal = k Ek . Figure 2 compares the time-averaged energy spectra of the high-resolution DNS, the lowresolution LLF-64 simulation, the TVD-64 scheme by van Leer [49], and the simulation of the reduced NN-64 model. The LLF-64 scheme is clearly too diffusive, artificially draining energy from wavenumbers starting with |k| ≥ 6 and causing a premature spectral decay. The TVD-64 scheme is less diffusive than the LLF-64 simulation; nevertheless, it still underestimates the spectra and the total energy at the resolution Nc = 64. This clearly demonstrates the need for subgrid modeling if the reduced resolution N = 64 is the target. The NN-64 reduced model accurately resolves the inertial range, tracking the DNS energy cascade precisely up to the grid cut-off without exhibiting spurious energy accumulation. We present the relative error for the total energy in each system in Table 2. The NN-64 subgrid modeling framework yields a low relative error compared to the LLF-64 model. This confirms that the localized Eddy Viscosity network correctly balances energy dissipation, ensuring physical conservation principles are maintained at the macroscopic level.
20
Energy Spectra, Ek vs k
100
DNS NN-64 LLF-64 TVD-64
10 1 10 2
Ek
10 3 10 4 10 5 10 6 10 7 10 8
100
101
102
k
Figure 2: Energy spectra EK in simulations of the fully resolved model with Nf = 512, reduced NN64 model, TVD-64 scheme, and coarse LLF-64 simulation. The NN-64 parameterization effectively corrects the excessive numerical dissipation inherent to the coarse LLF-64 and TVD-64 schemes, perfectly recovering the energy cascade. Model DNS NN-64 Model LLF-64 TVD-64
Total Energy (Etotal ) 0.8830 0.8746 0.8466 0.8696
Relative Error with DNS 0.95% 4.12 % 1.5 %
Table 2: Comparison of total energy and the relative error for the different numerical models. The TVD-64 scheme represents the Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL), originally formulated by van Leer [49]. The MUSCL scheme extends the first-order Godunov approach to achieve second-order spatial accuracy via a piecewise linear reconstruction of the state variables at the cell interfaces. In our implementation, we employ the classical van Leer limiter [48], which dynamically enforces the Total Variation Diminishing (TVD) condition by reverting to first-order accuracy in regions of steep gradients or local extrema. The reconstructed interface states are subsequently evaluated using the Local Lax-Friedrichs (LLF) solver to compute the numerical fluxes, ensuring a robust and entropy-stable baseline for performance comparison.
6.1
Comparison with Other Subgrid-Scale Models
To explicitly demonstrate the limitations of traditional closures, we compare the proposed neural architecture against classical static and dynamic Smagorinsky models [17, 27, 42]. Both models are
21
applied as explicit corrective fluxes on top of the Local Lax-Friedrichs (LLF) [39] baseline scheme on the coarse N = 64 grid. In addition, Wasserstein Generative Adversarial Network (WGAN) [1] and Stochastic Mode Reduction (SMR) [15] have been used to model subgrid processes for the same coarse variables (6) in the finite-volume discretization of the forced 1D Burgers equation. Overall, the WGAN and SMR parametrizations perform well in reproducing the spectra and temporal correlations. However, training the WGAN deep learning model was particularly challenging. The SMR approach also performed well, but it slightly overestimated the energy of higher wavenumbers in the spectra. Moreover, we found that when coarse variables are defined as local spatial averages, it may be challenging to compute SMR interaction coefficients for systems of equations and multi-dimensional problems since the SMR approach is a semi-analytical technique that requires knowledge of the statistics of fluctuations yi = ui − UI(i) . Finally, an analytical investigation of both the WGAN and SMR parametrizations is extremely challenging. Static Smagorinsky. The classical static Smagorinsky model [42] computes artificial viscosity directly proportional to the absolute local velocity gradient, given by the formulation νsgs = (Cs ∆x)2 ∂u ∂x . We first evaluate the static Smagorinsky model across a range of fixed coefficients. Figure 3 presents the energy spectra for Cs ∈ {0.01, 0.05, 0.15, 0.20}. For the inviscid Burgers’ equation, this spatial gradient is naturally maximized at the shock front. Consequently, the static model applies its heaviest numerical damping exactly where the wave is steepest. For artificially low values (Cs = 0.01), the subgrid contribution is negligible, leaving the overly dissipative LLF baseline unchanged. However, as Cs approaches standard theoretical values (Cs = 0.15 and 0.20), the excess viscosity severely degrades the wave amplitude.
22
Energy Spectra, Ek vs k
100
DNS NN-64 LLF-64 Smagorinsky-64
10 1
10 3
10 4
10 4
10 5
10 5
10 6
10 6
10 7
10 7 100
101
10 8
102
k
Energy Spectra, Ek vs k
100 10 2
10 4
Ek
10 4 10 5
10 5
10 6
10 6
10 7
10 7 101
k
10 8
102
102
k
Energy Spectra, Ek vs k DNS NN-64 LLF-64 Smagorinsky-64
10 2 10 3
100
101
10 1
10 3
10 8
100
100
DNS NN-64 LLF-64 Smagorinsky-64
10 1
Ek
10 2
10 3
10 8
DNS NN-64 LLF-64 Smagorinsky-64
10 1
Ek
Ek
10 2
Energy Spectra, Ek vs k
100
100
101
k
102
Figure 3: Energy spectra EK in simulations of the fully resolved model with Nf = 512 (blue), reduced NN-64 model (red; nearly overlaps with blue), coarse LLF-64 simulation (green), and Static Smagorinsky model (magenta) with Cs = 0.01 (top left), Cs = 0.05 (top right), Cs = 0.15 (bottom left), and Cs = 0.2 (bottom right). The NN parameterization effectively corrects the excessive numerical dissipation inherent to the coarse LLF-64 scheme, perfectly recovering the energy cascade. Dynamic Smagorinsky Model. We next evaluate the Dynamic Smagorinsky Model (DSM) [17, 27]. Unlike the static variant, DSM dynamically computes the coefficient Cs (x, t) utilizing a b (where typically dual-filtering procedure. By defining a grid-filter scale ∆ and a test-filter scale ∆ b = 2∆), the resolved turbulent stresses at the test-filter level, T , and the grid-filter level, τ , are ∆ c −u bu, b where L denotes the resolved Leonard related exactly by Germano’s identity: L = T − τb = uu stress. Applying the Smagorinsky closure at the grid-filter and test-filter levels gives L ≈ Cs2 M,
where
b 2 ∂x u b ∂x u b + 2∆2 |∂\ M = −2∆ x u| ∂x u.
Following the least-squares argument of Lilly [27], the dynamic coefficient is obtained by minimizing 2 the Germano-identity residual, L − Cs2 M . In the present one-dimensional setting this yields the 23
local scalar estimate Cs2 (x, t) =
L(x, t)M (x, t) , M 2 (x, t)
whenever M ̸= 0. Equivalently, in strictly scalar form, this reduces to Cs2 = L/M . In practice, one may introduce spatial, temporal, or ensemble averaging in the numerator and denominator to regularize this estimate. Figure 4 illustrates the failure of this local DSM formulation for the 1D Burgers test case considered here. The dynamic procedure assumes scale similarity of the modeled SGS stress between the grid-filter and test-filter levels. This assumption is reasonable in inertial-range turbulence but becomes questionable near under-resolved Burgers shocks, where the filtered gradient is dominated by the shock thickness rather than by a self-similar cascade. Across such sharp interfaces, the c2 − u b̄,2 measures the finite resolved variance introduced by the test filter, Leonard stress L = ū while the modeled tensor M is constructed from nonlinear functions of the filtered gradient. The resulting local ratio Cs2 = LM/M 2 can therefore become poorly conditioned or change sign near shocks. When LM < 0, the dynamic procedure yields a negative eddy viscosity, corresponding to local backscatter or anti-diffusion. In a shock-dominated Burgers solution, such anti-diffusive contributions are numerically dangerous because they can oppose the entropy-dissipative regularization needed to obtain the physically admissible solution. Applying ad-hoc clipping removes these negative-viscosity regions, but it does not control large positive values of Cs2 . Consequently, localized peaks in the dynamic coefficient near the shock can produce excessive artificial viscosity, leading to over-damping and distortion of the resolved wave.
24
Energy Spectra, Ek vs k
100
DNS NN-64 LLF-64 Dynamic Smagorinsky-64
10 1 10 2
Ek
10 3 10 4 10 5 10 6 10 7 10 8
100
101
k
102
Figure 4: Energy spectra EK in simulations of the fully resolved model with Nf = 512 (blue), reduced NN-64 model (red; nearly overlaps with blue), coarse LLF-64 simulation (green), and Dynamic Smagorinsky model (magenta). The failure of the Dynamic Smagorinsky Model (DSM) is likely caused by applying the test filter across sharp discontinuities, where the scale-similarity assumption breaks down; the resulting Leonard stress is dominated by the finite shock jump, leading to poorly conditioned dynamic-coefficient estimates and localized excess viscosity that distorts the resolved wave structure. Spatial Correlation Function. To validate how well spatial structures are captured by the reduced model, we analyze the two-point spatial correlation function. This statistical metric quantifies the dependence and spatial coherence between different locations across the domain. For the scalar state variable U (x, t), the time-averaged spatial correlation evaluated at a discrete grid point XI is defined as: N
C(XI ) =
T 1 X U (X0 , tk )U (XI , tk ), NT
(46)
k=1
where NT represents the total number of sampled time snapshots. Because the stochastically forced system is spatially homogeneous and operates within a periodic domain, its statistical properties are translationally invariant. Therefore, computing the correlation with the first spatial grid point X0 is sufficient to characterize the entire domain. Figure 5 (left) presents the normalized spatial correlation function for the resolved flow field. The curve exhibits the characteristic symmetric profile expected in a periodic domain, reflecting the spatial decorrelation length and the recovery of coherency across the periodic boundary. Both the reduced NN-64 and coarse LLF-64 models show agreement with the high-resolution DNS reference. This implies that the NN corrections are not affecting the spatial correlation significantly, and 25
the reduced model accurately reproduces the spatial correlation scales. This confirms that our structure-preserving architecture successfully captures the macroscopic spatial scales and physical coherence of the flow field, despite operating on a truncated grid resolution.
Spatial Correlation Function
1.0 0.8
DNS NN-64 LLF-64
0.8
0.6
0.6
CF( )
0.4
C(x)
Temporal Correlation Function
1.0
0.2 0.0
0.4 0.2
0.2
DNS NN-64 LLF-64
0.4 0.6 0
1
0.0 2
3
x
4
5
6
0
2
4
6
8
10
Figure 5: Left: Spatial Correlation Function for coarse variables UI given by (46), Right: Temporal Correlation Function of coarse variables U0 given by (47). Solid Blue line: high-resolution DNS simulations, Dashed Green line: simulations of the LLF-64 model, Dashed Red line: simulations of the NN-64 reduced model. Temporal Correlation Function. To evaluate the dynamical consistency of the reduced model, we compute the temporal correlation function. This statistical metric quantifies the system’s memory. It indicates how rapidly dynamical variables decorrelate from their prior states. Since the Burgers’ equation is spatially homogeneous, we consider only the spatial location U (X0 , t). The time correlation function is defined as: Z 1 T CF (τ ) = U (X0 , t)U (X0 , t + τ ) dt. (47) T 0 Computationally, the time-correlation function is computed using a moving-window approach. The normalized temporal correlation functions for the coarse variables in NN-64, LLF-64, and DNS are presented in the right part of Figure 5. The NN-64 reduced model demonstrates very good agreement with the high-resolution DNS simulation. For the temporal correlation function, a slight discrepancy is observed between the LLF-64 and DNS curves, indicating that the NN correction terms improve the temporal correlation statistics. The NN parameterization accurately captures the decorrelation timescale and the decay profile of the correlation function. Numerical results with spatial and temporal correlation functions confirm that the reduced-order framework preserves not only the spatial energy cascade but also the accurate statistical evolution and temporal dynamics of the full-scale system. Solution Snapshots. Neural network approximations are often treated as data-driven ”blackbox” models without incorporating physical constraints. As a consequence, neural network approximations, when left uncontrolled, can produce unbounded oscillations or unrealistic solutions. One possible approach is to employ flux limiters [34, 46] to suppress non-physical behavior and spurious oscillations. In this paper, we demonstrate that our approach does not require flux-limiting stabilization techniques. This is a direct consequence of using structure-preserving neural networks, 26
which enforce stability at the architectural level. The structure-preserving neural network design used in this paper guarantees that the subgrid parameterization generates robust, physically consistent solutions, replacing the need for traditional flux-limiting procedures with rigorous architectural constraints. Figure 6 depicts typical solution profiles for four selected times. For comparison, we also include the LLF-64 coarse-model results, which bear little resemblance to the DNS solution and fail to reproduce the relevant dynamics. We can see that the NN model demonstrates excellent agreement with the DNS, successfully capturing shock formations and nonlinear wave propagation. The embedded Eddy Viscosity network is essential in naturally suppressing spurious oscillations near steep gradients, eliminating the need for external flux limiters. t=8.0
t=50.0 DNS NN-64 LLF-64
1
1
0
0
U(x, t)
U(x, t)
2
2
1
1
2
DNS NN-64 LLF-64
2
3 0
1
2
3
x
4
5
6
0
1
2
t=100.0
2
3
x
4
5
6
t=200.0 DNS NN-64 LLF-64
2
1
U(x, t)
U(x, t)
1 0
0
1 1 2 0
1
2
3
x
4
5
DNS NN-64 LLF-64 6
2 0
1
2
3
x
4
5
6
Figure 6: Comparison between velocity snapshots U (x, t) for the Burgers’ equation in simulations of the NN-64 reduced model (red dashed), LLF-64 (green dash-dot), and DNS-512 (blue solid) at times t = 8, 50, 100, 200.
6.2
Extrapolation Results
One important practical concern regarding neural networks is their ability to generalize beyond the training regime. To evaluate this, we apply the NN reduced model without retraining to different forcing regimes and initial conditions, analyzing its robustness and adaptability. Larger Forcing Regimes. To assess the robustness and physical fidelity of our framework, we perform simulations of the reduced model with out-of-distribution forcing. In particular, the forcing amplitude is increased from the baseline training value of A = 1.0 to A = 1.2 and A = 1.4, representing increases of 20% and 40% in the energy injected into the system, respectively. Recall 27
that forcing is given by (2). The forcing is injected into the larger scales and cascades down to the unresolved subgrid scales through nonlinear interactions. Accurate subgrid parameterization must dynamically adapt to this increased energy throughput to maintain stability and spectral accuracy. The energy spectra for simulations with larger forcing are presented in Figure 7. The NN-64 reduced model accurately tracks the DNS energy cascade across the entire resolved wavenumber range. The model achieves this consistency across all forcing regimes without retraining. The localized Eddy Viscosity network, bounded by physical constraints, successfully adjusts the subgrid dissipation dynamically. It applies the precise amount of dissipation required to prevent spurious oscillations while preserving the high-wavenumber tail. We also verified that solution snapshots are reproduced very well in regimes with increased forcing. Figure 8 depicts snapshots in simulations with A = 1.4. The same initial condition and forcing realization are used in simulations with A = 1 depicted in Figure 6. As expected, the overall magnitude of the solution increases with larger forcing, but all networks generalize well in this regime. In particular, the Eddy Viscosity network suppresses spurious oscillations near shocks, while the Structure-Preserving network models the nonlinear wave interactions. Our numerical results confirm that all three networks successfully learned the physics of the nonlinear energy cascade, allowing them to generalize robustly beyond the original training regime.
Energy Spectra, Ek vs k
100 10 1
10 2
10 3
10 3
10 4
10 4
10 5
10 5
10 6
10 6
10 7
10 7
10 8
DNS NN-64 LLF-64
10 1
Ek
Ek
10 2
Energy Spectra, Ek vs k
100
DNS NN-64 LLF-64
100
101
k
10 8
102
100
101
k
102
Figure 7: Energy spectra of DNS, NN-64, and LLF-64 models with increased forcing A = 1.2 (left) and A = 1.4 (right). The NN-64 reduced model is simulated without retraining the NN parametrization.
28
t=8.0
t=50.0 DNS NN-64 LLF-64
3 2
2 1
0
U(x, t)
U(x, t)
1
1 2
0 1 2
3
DNS NN-64 LLF-64
3
4 0
1
2
3
x
4
5
6
0
1
2
t=100.0
3
x
4
5
6
4
5
6
t=200.0 DNS NN-64 LLF-64
3
2
2
1
U(x, t)
U(x, t)
1 0 1
0 1
2 3 0
1
2
3
x
4
5
DNS NN-64 LLF-64 6
2 3
0
1
2
3
x
Figure 8: Comparison between velocity snapshots U (x, t) for the Burgers’ equation in simulations of the NN-64 reduced model (red dashed), LLF-64 (green dash-dot), and DNS-512 (blue solid) with increased forcing A = 1.4 at times t = 8, 50, 100, 200. Compare with simulations with A = 1 in Figure 6. Note that the NN parameterization was trained on the baseline regime A = 1 and is evaluated here without retraining. Step Initial Conditions. To further evaluate how well our framework generalizes to out-ofsample regimes, we consider discontinuous initial conditions. While the Neural Network parameterization was trained exclusively on smooth initial wave profiles, we evaluate its performance here on a discontinuous step function (a Riemann-type problem) without any retraining. The specific initial condition considered here is ( 2π 4π 0.5, 3 ≤x< 3 u(x, 0) = −0.5, otherwise Handling sharp discontinuities can be particularly challenging for machine learning parametrizations. Black-box models trained exclusively on smooth data typically suffer from severe, nonphysical Gibbs oscillations or instability when encountering shocks (see e.g., [13, 16, 24, 34, 46]). In contrast, our structure-preserving NN parametrization remains remarkably stable and does not require any adjuSMRents for discontinuous data. Figure 9 compares the DNS, NN-64, and LLF-64 solution snapshots for the discontinuous initial condition at t = 0.1, 10, 40, 50. At early simulation times (e.g., t = 0.1), the reduced-order model accurately captures the sharp, near-vertical gradients of the initial step without overshoots. As the simulation evolves (t = 10, 40, 50), the discontinuities propagate, interact, and decay. Through29
out this evolution, the NN-64 reduced model consistently reproduces the correct shock speed and amplitude of the high-resolution DNS reference. This successful extrapolation highlights the robustness of our structure-preserving NN parametrization approach and the embedded Eddy Viscosity network in particular. Because the architecture enforces physical bounds and localized dissipation dynamically, this network successfully recognizes and stabilizes steep gradients even in a physical regime absent from the training data. Overall, the NN-64 parametrization provides a substantial improvement not only for averaged quantities but also for individual solution trajectories. t=0.1 0.4
1
U(x, t)
U(x, t)
0.2 0.0 0.2
t=10.0
2
DNS NN-64 LLF-64
0 1
0.4
DNS NN-64 LLF-64
2
0.6 0
1
2
3
x
4
5
6
0
1
2
t=40.0
x
4
5
6
5
DNS NN-64 LLF-64 6
t=50.0
1.5
1.0
1.0
0.5
U(x, t)
0.5
U(x, t)
3
0.0 0.5
0.0 0.5 1.0
1.0
1.5
DNS NN-64 LLF-64
1.5 2.0 0
2.0 1
2
3
x
4
5
6
0
1
2
3
x
4
Figure 9: Comparison between velocity snapshots U (x, t) for the Burgers’ equation in simulations of the NN-64 reduced model (red dashed), LLF-64 (green dash-dot), and DNS-512 (blue solid) with step initial conditions at times t = 0.1, 10, 40, 50. Note that the NN parameterization was trained on solutions with smooth initial conditions and is evaluated here without retraining. The training time for the neural network model is approximately 12 minutes. This fast training time reflects the design (relatively small depth and width) of all three neural networks. Table 3 reports the inference cost for the fully resolved DNS model with Nf = 512, the coarse LLF-64 and NN-64 models using the same time step as the DNS, and the NN-64 model using a time step four times larger, 4∆t. The inference time of the coarse NN-64 model using the smaller time step is greater than that of the DNS because evaluating the neural networks requires additional operations. However, the coarser spatial resolution allows the use of a larger time step, which substantially reduces the inference cost of the NN-64 model. We verified that simulations performed with the larger time step reproduce the same statistical quantities and trajectory predictions. We expect the relative computational savings to become more significant for multidimensional 30
problems, where the cost of fully resolved simulations increases rapidly with spatial resolution. Table 3: Inference wall-clock time for simulation with T = 1000 for the DNS (Nf = 512) and simulations on a uniform coarse grid (Nc = 64). Scheme
Inference Time (s)
DNS LLF-64 NN-64 NN-64 with 4∆t
7
243.46 165.59 329.90 80.87
Conclusions
In this paper, we develop a structure-preserving neural network (NN) approach for parametrizing subgrid physical processes in coarse simulations of hyperbolic conservation laws. This structurepreserving approach relies on the underlying analytical properties of hyperbolic systems, but does not require explicit knowledge of the entropy function. Instead, the entropy function is learned using an input-convex neural network. Subgrid fluxes are then estimated using another network to learn the flux potential. A third neural network approximates the Eddy Viscosity and has design features similar to traditional eddy viscosity models. The main contribution of this work is the formulation and numerical validation of a structure-preserving neural-network framework for subgrid closure. The novelty lies in adapting this methodology to represent unresolved fluxes in a prescribed coarse discretization and in combining the learned entropy function and flux potential with a neural-network eddy-viscosity closure. We demonstrate that our NN reduced model performs very well and reproduces the statistical properties of solutions with high accuracy. In particular, there is no over-damping at higher wavenumbers in the energy spectra of coarse variables UI (see Figure 2), which is usually difficult to achieve. We also demonstrate that the NN parametrization is structurally stable and that the solutions of the NN reduced model do not develop spurious oscillations. Consequently, the proposed NN parametrization approach does not require additional stabilization techniques. Numerical results in section 6.2 show that the NN reduced model generalizes well outside the training regime. In simulations with larger forcing, both the total energy and the magnitude of solutions increase, and the NN parametrization operates outside of the training regime. The NN reduced model reproduces the energy spectra and individual solutions in simulations with increased forcing very well. Numerical simulations with discontinuous initial conditions also demonstrate that the NN reduced model generalizes well to regimes with discontinuous initial data. Results with larger forcing are particularly relevant for potential applications in atmospheric fluid dynamics, where it is important to consider increased atmospheric forcing on the ocean due to climate change scenarios. The numerical comparisons with the static and dynamic Smagorinsky models further demonstrate the advantages of the proposed formulation. For the Burgers-equation test cases considered here, which include solutions with shocks, the static model introduces excessive dissipation, whereas the dynamic model can produce excessively large stresses because its scale-invariance assumption is not well suited to these cases. The comparison with the second-order TVD method also provides 31
a conventional shock-capturing baseline. These results show that the proposed closure improves the representation of unresolved dynamics relative to the classical models considered in this study. The analytical results in Section 5 demonstrate that the resulting NN parametrization is entropy stable and provide useful guidelines for selecting the range for the diffusion constant in the Eddy Viscosity Network. This highlights the advantage of a structure-preserving approach to data-driven flux estimation in the reduced model. In contrast to black-box closure strategies, the proposed formulation preserves key stability properties while retaining the flexibility of neural-network-based parametrizations. This combination of analytical tractability, numerical robustness, and improved predictive accuracy makes the approach particularly attractive for reduced-order modeling of hyperbolic conservation laws. Overall, the NN parametrization yields a considerable improvement in the accuracy of both the energy spectra and individual solution trajectories. The model also extrapolates well to regimes with stronger forcing, suggesting that the learned closure captures robust coarse-scale dynamics rather than merely interpolating within the training distribution. A further advantage is that accurate performance is obtained using relatively small neural networks, making the approach computationally lightweight and straightforward to train compared, for instance, with WGANbased closure models. Finally, because the parametrization is formulated at the level of the fluxes, it remains amenable to analytical investigation, providing a useful bridge between data-driven modeling and structure-preserving numerical analysis. The structure-preserving approach allows the use of a small number of hidden layers and neurons in each layer. For instance, the Eddy Viscosity network has only one hidden layer. In addition, we expect that our approach can be extended to other hyperbolic problems, but applications to more complex equations, and especially to systems of equations, introduce additional network-design considerations. For instance, it is important to properly design the Eddy Viscosity network to avoid over-damping at higher wavenumbers of the coarse variables, which is a common problem in many closure models. Therefore, it is necessary to select an appropriate vector of feature variables ξ in (15), and possibly increase the number of layers or split the Eddy Viscosity network into several networks with scalar outputs to avoid inaccurate estimation of the vector C(ξ). This requires a separate investigation and will be addressed in future work. We do not foresee conceptual issues in extending our approach to systems of equations and multidimensional problems. Therefore, extending the structure-preserving framework to systems of conservation laws and multidimensional problems represents a natural and highly promising direction for future work. We expect that the overall architecture of the Entropy Neural Network and the Flux Potential Neural Network would remain similar. However, it might be necessary to increase the number of layers and/or the number of neurons to incorporate multidimensional spatial stencils. Developing an efficient and accurate Eddy Viscosity network might require a separate investigation, as discussed in the previous paragraph. For multidimensional systems, the entropy function would depend on the full vector of state variables, and the neural-network architecture would need to account for fluxes in multiple spatial directions. The Eddy Viscosity Network would also need to incorporate multidimensional spatial derivatives and interactions between different directions. The resulting increase in the number of state variables, spatial stencils, and model parameters may also lead to more demanding and potentially less well-conditioned optimization problems. Determining suitable architectures and training strategies that preserve entropy stability and provide adequate dissipation without introducing excessive numerical diffusion will therefore require a dedicated investigation. Moreover, extending the entropy-stability analysis would likely
32
require additional structural assumptions on the neural networks, potentially imposing further constraints on their architecture.
Acknowledgments The authors thank Dr. Lu Zhang for helpful discussions.
References [1] J. Alcala and I. Timofeyev, Subgrid-scale parametrization of unresolved scales in forced Burgers equation using generative adversarial networks (GAN), Theoretical and Computational Fluid Dynamics, 35 (2021), pp. 875–894. [2] B. Amos, L. Xu, and J. Z. Kolter, Input convex neural networks, in International conference on machine learning, PMLR, 2017, pp. 146–155. [3] P. Benner, S. Gugercin, and K. Willcox, Survey of projection-based model reduction methods for parametric dynamical systems, SIAM Review, 57 (2015), pp. 483–531. [4] M. R. Bin Mizan, M. Olshanskii, and I. Timofeyev, A parametric tensor rom for the shallow water dam break problem, Computers & Fluids, 309 (2026), p. 107005. [5] M. R. Bin Mizan, I. Timofeyev, and M. Olshanskii, Reduced-order modeling of parameterized visco-plastic shallow flows, arXiv preprint arXiv:2605.06526, (2026). [6] T. Bolton and L. Zanna, Applications of deep learning to ocean data inference and subgrid parameterization, Journal of Advances in Modeling Earth Systems, 11 (2019), pp. 376–399. [7] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, Y. Katariya, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. WandermanMilne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs; http://github.com/jax-ml/jax, 2018. [8] N. D. Brenowitz and C. S. Bretherton, Spatially extended tests of a neural network parametrization trained by coarse-graining, Journal of Advances in Modeling Earth Systems, 11 (2019), pp. 2728–2744. [9] J. M. Burgers, A mathematical model illustrating the theory of turbulence, Advances in applied mechanics, 1 (1948), pp. 171–199. [10] M. Chantry, H. Christensen, P. Dueben, and T. Palmer, Opportunities and challenges for machine learning in weather and climate modelling: hard, medium and soft ai, Philosophical Transactions of the Royal Society A, 379 (2021), p. 20200083. [11] Z. Chen, A. Gelb, and Y. Lee, Learning the dynamics for unknown hyperbolic conservation laws using deep neural networks, SIAM Journal on Scientific Computing, 46 (2024), pp. A825– A850.
33
[12] B. Cockburn and C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convectiondominated problems, Journal of Scientific Computing, 16 (2001), pp. 173–261. [13] E. J. Coutinho, H. Zhao, A. M. Cunha, S. Afonso, and S. Mukherjee, Physicsinformed neural networks with adaptive localized artificial viscosity, Journal of Computational Physics, (2023). [14] C. M. Dafermos, Hyberbolic conservation laws in continuum physics, Springer, 2005. [15] S. I. Dolaptchiev, U. Achatz, and I. Timofeyev, Stochastic closure for local averages in the finite-difference discretization of the forced Burgers equation, Theor. Comput. Fluid Dyn., 27(3-4) (2013), pp. 297–317. [16] O. Fuks and H. A. Tchelepi, Limitations of physics informed machine learning for nonlinear two-phase transport in porous media, Journal of Machine Learning for Modeling and Computing, 1 (2020). [17] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Physics of fluids a: Fluid dynamics, 3 (1991), pp. 1760–1765. [18] E. Godlewski and P.-A. Raviart, Numerical approximation of hyperbolic systems of conservation laws, vol. 118, Springer, 2021. [19] S. K. Godunov, A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations, Matematicheskii Sbornik, 47 (1959), pp. 271–306. Translated from Russian. Cornell Aeronautical Laboratory translation, 1969. [20] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu, Strong stability preserving Runge-Kutta and multistep time discretizations, World Scientific, 2011. [21] J. S. Hesthaven, Numerical methods for conservation laws: From analysis to algorithms, SIAM, 2017. [22] G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted eno schemes, Journal of computational physics, 126 (1996), pp. 202–228. [23] V. M. Krasnopolsky, M. S. Fox-Rabinovitz, and D. V. Chalikov, New approach to calculation of atmospheric model physics: Accurate and fast neural network emulation of longwave radiation in a climate model, Monthly Weather Review, 133 (2005), pp. 1370–1383. [24] A. S. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney, Characterizing possible failure modes in physics-informed neural networks, Advances in Neural Information Processing Systems, 34 (2021), pp. 26548–26560. [25] M. Lesieur, O. Métais, and P. Comte, Large-eddy simulations of turbulence, Cambridge University Press, 2005. [26] R. J. LeVeque, Numerical methods for conservation laws, vol. 132, Springer, 1992. [27] D. K. Lilly, A proposed modification of the germano subgrid-scale closure method, Physics of Fluids A: Fluid Dynamics, 4 (1992), pp. 633–635. 34
[28] M. Lino, S. Fotiadis, A. A. Bharath, and C. D. Cantwell, Current and emerging deep-learning methods for the simulation of fluid dynamics, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479 (2023). [29] L. Liu, L. Zhang, and A. Gelb, Neural entropy-stable conservative flux form neural networks for learning hyperbolic conservation laws, Journal of Computational Physics, 553 (2026), p. 114719. [30] A. J. Majda, I. Timofeyev, and E. Vanden-Eijnden, A mathematics framework for stochastic climate models, Comm. Pure Appl. Math., 54 (2001), pp. 891–974. [31] A. V. Mamonov and M. A. Olshanskii, Interpolatory tensorial reduced order models for parametric dynamical systems, Computer Methods in Applied Mechanics and Engineering, 397 (2022), p. 115122. [32]
, Tensorial parametric model order reduction of nonlinear dynamical systems, SIAM Journal on Scientific Computing, 46 (2024), pp. A1850–A1878.
[33]
, A priori analysis of a tensor ROM for parameter dependent parabolic problems, SIAM Journal on Numerical Analysis, 63 (2025), pp. 239–261.
[34] M. A. H. Mojamder, Z. Xu, M. Wang, and I. Timofeyev, Parametrization of subgrid scales in long-term simulations of the shallow-water equations using machine learning and convex limiting, Fluids, 11 (2026), p. 75. [35] A. Quarteroni, A. Manzoni, and F. Negri, Reduced Basis Methods for Partial Differential Equations: An Introduction, Springer, 2015. [36] P. Ramachandran, B. Zoph, and Q. V. Le, Searching for activation functions, arXiv preprint arXiv:1710.05941, (2017). [37] S. Rasp, M. S. Pritchard, and P. Gentine, Deep learning to represent subgrid processes in climate models, Proceedings of the National Academy of Sciences, 115 (2018), pp. 9684–9689. [38] C. W. Rowley and S. T. Dawson, Model reduction for flow analysis and control, Annual Review of Fluid Mechanics, 49 (2017), pp. 387–417. [39] V. V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, Journal of Computational Mathematics and Mathematical Physics, 1 (1961), pp. 267–279. [40] P. Sagaut, Large eddy simulation for incompressible flows: an introduction, Springer, 2006. [41] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Lecture Notes in Mathematics, 1697 (1998), pp. 325–432. [42] J. Smagorinsky, General circulation experiments with the primitive equations: I. the basic experiment, Monthly weather review, 91 (1963), pp. 99–164. [43] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM Journal on Numerical Analysis, 21 (1984), pp. 995–1011.
35
[44] E. Tadmor, The numerical viscosity of entropy stable schemes for systems of conservation laws. I, Mathematics of Computation, 49 (1987), pp. 91–103. [45]
, Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems, Acta Numerica, 12 (2003), pp. 451–512.
[46] I. Timofeyev, A. Schwarzmann, and D. Kuzmin, Application of machine learning and convex limiting to subgrid flux modeling in the shallow-water equations, Mathematics and Computers in Simulation, 238 (2025), pp. 163–178. [47] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 3rd ed., 2009. [48] B. Van Leer, Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme, Journal of computational physics, 14 (1974), pp. 361–370. [49]
, Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov’s method, Journal of computational Physics, 32 (1979), pp. 101–136.
[50] M. Zacharuk, S. I. Dolaptchiev, U. Achatz, and I. Timofeyev, Stochastic subgridscale parametrization for one-dimensional shallow-water dynamics using stochastic mode reduction, Quarterly Journal of the Royal Meteorological Society, 144 (2018), pp. 1975–1990.
36