{'fraction_non_alphanumeric': 0.030380181365677834, 'fraction_numerical': 0.009611808793660782, 'mean_word_length': 4.907063197026022, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Subgrid-scale parametrization of unresolved scales in forced Burgers equation using Generative Adversarial Networks (GAN)

Stochastic subgrid-scale parametrizations aim to incorporate effects of unresolved processes in an effective model by sampling from a distribution usually described in terms of resolved modes. This is an active research area in climate, weather and ocean science where processes evolved in a wide range of spatial and temporal scales. In this study, we evaluate the performance of conditional generative adversarial network (GAN) in parametrizing subgrid-scale effects in a finite-difference discretization of stochastically forced Burgers equation. We define resolved modes as local spatial averages and deviations from these averages are the unresolved degrees of freedom. We train a Wasserstein GAN (WGAN) conditioned on the resolved variables to learn the distribution of subgrid flux tendencies for resolved modes and, thus, represent the effect of unresolved scales. The resulting WGAN is then used in an effective model to reproduce the statistical features of resolved modes. We demonstrate that various stationary statistical quantities such as spectrum, moments, autocorrelation, etc. are well approximated by this effective model.


Introduction
Many physical processes in mathematics, engineering and other applied sciences are represented by differential equations. In practical numerical applications partial differential equations are discretized in space and time. Often finite-difference or finite-volume numerical schemes are used for spatial discretization which results in a system of ordinary differential equations or difference equations. However, in doing so, model error are introduced and some of the physical processes may not be fully resolved on the defined spatial grid. Either, some processes evolve on a much smaller spatial (subgrid) scale compared to the chosen mesh size or there are inherent properties of the system (e.g. nonlinearity, chaotic behavior) that cannot be captured by the chosen discretization. The problem of incorporating the approximate effects of unresolved processes or sub-grid scale processes in an effective model is called parametrization (or sub-grid scale parametrization, in our context).
Developing accurate and fast parametrizations has been an active research area in ocean, climate and atmospheric science for many decades. For deterministic parametrization schemes a function relating unresolved modes with every state of the resolved modes is constructed either by using known physics or existing high-resolution data. This leads to incorporating the mean impact of subgrid processes on resolved modes. To be realistic, some measure of variability of the impacts of these processes must be taken into account as proposed by Hasselman [21] in a visionary paper in 1976. Introduction of randomness leads to stochastic parametrizations.
In recent decades, there have been a considerable number of studies focusing on stochastic parametrizations. Many methods have been developed which aim to parametrize effects of sub-grid scales by incorporating stochastic terms with the goal of reproducing the statistical behavior of resolved scales. To mention a few, there are techniques which use regression models (e.g. [39,25]), add stochastic perturbations to the mean parametrization terms [34,14], utilize data-driven methods such as NARMAX fitting [9,27], or use conditional Markov chains [10,24,26]. There are also more rigorous approaches (see recent review paper [33]) including stochastic mode reduction strategy exploiting scale separation between the resolved and sub-grid scale processes [28,29,30,16,15], reduction techniques based on the Zwanzig-Mori formalism [6,7,8,22], turbulence closure techniques [17,18,36]. In particular, stochastic mode-reduction has been applied recently to reproduce the statistical behavior of local spatial averages in the truncated Burgers-Hopf model [12,11] and the Shallow-Water equations [42]. A more relaxed scheme that relies on weak coupling and uses ideas from response theory has been developed in [40]. Another approach uses maximum entropy principle for estimating probability density of unresolved scales [35,37]. Stochastic parametrizations have been proven to be extremely useful in weather and climate science for analyzing prediction uncertainty, representing sub-grid processes, capturing response to forcing and inducing regime transitions (see e.g. [13]).
Machine learning has become increasingly popular in recent years offering possible data-driven solutions for problems in many areas of physical and natural sciences. Machine learning also received increased attention in applied and computational fluid dynamics as well. One advantage of most machine learning methods is that they do not necessarily require assumptions about the model equations and/or structure of the data. Recently, neural networks (NNs) became one of the most popular area of machine learning and there are several applications of NNs in various models relevant to atmospheric research. In particular, NNs were used to represent clouds [32] and convection [31] in a global circulation model. NN was also used to develop a parametrization for residual heating and moistening in a global circulation model [3]. Various deep learning methods such as reservoir computing and LSTM neural networks were used to predict trajectories of chaotic dynamics (e.g. [38,5]). Methods used in these papers are classified as supervised learning and results in deterministic parametrizations. Here we use generative models to develop NN-based stochastic parametrizations.
Application of generative models is one of the relatively new areas in machine learning. These are primarily techniques from unsupervised learning in which a network is trained to produce samples from a given distribution. One popular generative model is generative adversarial networks (GAN) introduced in [19]. This model consists of two networks -generator and discriminator. Given a training set, these two networks compete with each other -the generator attempts to produce new samples consistent with the training set, and the discriminator evaluates them with the goal of distinguishing whether these samples came from the original distribution. Thus, GAN aims to learn the underlying distribution for the training data in order to produce new samples that a critic (i.e. discriminator) qualifies to be coming from that distribution. It is easy to see that using GAN would result in a stochastic parametrization. A modified version of this approach was applied to Lorenz 96 model [23]. Recently, an alternative training metric for generative adversarial networks has been proposed. It particular, it was demonstrated that Wasserstein GAN (WGAN) might have certain advantages over the traditional GAN during training [2,20].
In this paper a conditional WGAN is trained to learn the distribution of subgrid processes and then used to generate possible values of unresolved modes based on the current state of the resolved variables. We apply this approach for parametrizing local residuals in a flux discretization of stochastically forced Burgers' equation. Specifically, a WGAN conditioned on resolved modes is trained to produce samples from the distribution of subgrid fluxes to be used in an effective model for local averages.
The rest of the paper is organized as follows. In section 2.1 we discuss general methodology for defining coarse variables and representing the effective equations using flux discretization. In section 2.2 we discuss application of this methodology to the forced Burgers-Hopf model. In section 2.3 we discuss configuration and training of the neural network used to represent the subgrid fluxes. In section 3 we present numerical results comparing stationary statistics of the full and reduced models and discuss performance of the WGAN parametrization with respect to changes in forcing and resolution. Final remarks are presented in section 4.

Subgrid flux parametrization
Assume that a system is described by a high dimensional state vector u which can be decomposed as u = (U, y) where U are the resolved variables and y are unresolved (subgrid) degrees of freedom. In the context of this paper, vector u of time-dependent variables is obtained by a suitable finite-difference or finite-volume discretization of a partial differential equation. Typically, vector u is high-dimensional, while dim(U) dim(u).
The complete dynamics of u can be recast as a coupled system for resolved and unresolved variables where F(U) represents the interactions among the resolved modes, G(U, y) represents interactions between the resolved and unresolved variables, etc. In the context of subgrid scale parametrization developed in this paper, the first equation (1) describes the time evolution of large-scale variables on a coarser grid. The term G(U, y) represents subgrid scale tendencies arising from the interaction among the resolved and unresolved degrees of freedom. In the most general sense, the aim of parametrization is to approximate G(U, y) in terms of the resolved modes only so that the resulting closed system approximates, in some sense, the behavior of resolved variables U in the full dynamics. This paper is focused on stochastic parametrizations where approximations to G are drawn from the distribution of subgrid-scale tendencies. In addition, it is necessary to condition the distribution of G on the current values of the large-scale variables, U.
The method developed in this paper applies to system based on flux difference discretization of partial differential equations that can be written in conservative form. More precisely, consider a state field u(x, t) which can be discretized on a uniform fine mesh M f = {0, 1, 2, ..., N − 1} resulting in a system of ODEs where ∆x = x i+1/2 − x i−1/2 . The value u i is the numerical approximate for either the value of u(x, t) at the center of interval [x i−1/2 , x i+1/2 ] or for the average 1 ∆x x i+1/2 x i−1/2 u(x, t)dx depending on whether the numerical method used is a finite difference or a finite volume scheme, respectively. More applicable to finite volume methods, F i+1/2 and F i−1/2 are discrete approximations for the flux at the right and left ends of the interval i, respectively. Here, ρ represents forcing. Equation above represents the full dynamics. Typically, direct numerical simulations (DNS) of the equation (3) (same as (1), (2)) on a small enough mesh are computationally expensive and it is desirable to obtain a closed-form system for variables U alone.
We define the resolved modes as the mean of u i over n neighboring high resolution grid points, that is, for each cell I in a coarse mesh M c = {0, 1, 2, N/n − 1}. The unresolved degrees of freedom are then the residuals or deviations y i = u i − U I(i) , where index i for y i and u i refers to the coarse cell I where fine cell i is located. The idea of using averaging operator is similar to large-eddy simulation (LES) modelling in turbulence. Evidently, some degree of scale separation can be expected between U I and y i by this construction, but it might depend on the number of averaged points, n. Our goal is to develop a closed effective equation for variables U I on a coarser mesh of size n∆x. To this end, we consider coupled dynamical system for U = {U I , I ∈ M c } and y = {y i , i ∈ M f }. There are various ways to decompose the dynamics in (3) and define fluxes F(U) and G(U, y). However, a natural way is to use the same flux discretization for F(U) as for the full dynamics, but now applied to the coarse mesh. This results in the system where ρ U I and ρ yi are the projections of the forcing on the respective modes. We will refer to the coupled system (5), (6) as the full model. Note that F n(I+1)−1/2 is the flux at the right boundary of the coarse cell I and it is a function of both the resolved and unresolved modes. We can express this flux as where G I+1/2 is referred to as the subgrid flux at the right boundary of cell I. Term F I+1/2 (U) describes self-interactions of local averages on the coarse mesh using the same discrete flux used in (3). Physical processes which cannot be adequately resolved with the mesh size n∆x are represented by G I+1/2 . Also, note that although we write the entire vector U and y as predictors, only few neighboring 1 U I and y i are actually used in the right-hand size of the equation for U I .
We would like to point out that decomposition (5), (6), and (7) is a particular case of more general form of the resolved -subgrid variables decomposition presented in (1), (2). Decomposition in (5), (6), and (7) imposes a particular structure on the right-hand side of the full model, and thus, on the righthand side of the effective model for the resolved variables. In particular, the model in (5), (6), (7) involves subgrid fluxes in an explicit form and, thus, the goal is to develop parametrization for the subgrid fluxes G I+1/2 (U, y). Formulation in (5), (6), (7) has several potential advantages. In particular, (i) subgrid fluxes defined in (7) may contain features that can be lost when net flux tendencies in (1) are considered instead, (ii) approximating the distribution of subgrid flux might present less complexity compared with the distribution of the net flux, (iii) decomposition (5) with (7) enforces mass conservation (at least on average) that is also achieved in application of flux difference scheme for U in the full system, (iv) same flux difference scheme is used for the effective equation as in the full model with parametrized subgrid fluxes considered as corrections; this allows for better understanding of relative strengths of flux terms and potential introduction of empirical corrections.
Replacing the sudgrid fluxes G I+1/2 with a suitable parametrizationG I+1/2 , the effective (or reduced) equation for large-scale variables becomes 1 Number depends on the stencil used to represent the numerical flux Here, we assume thatG I+1/2 for each I is a random variable, possibly conditioned on some neighboring resolved variables. ParametrizationG I+1/2 represent mechanism for effectively sampling the distribution of G I+1/2 and, thus, allowing for an efficient numerical simulation of the effective equation in (8). In general, there are no restrictions on the form ofG I+1/2 . In this paper we utilize a neural network to representG I+1/2 . In particular, we use generative adversarial network (GAN) to generate samples from the distribution of G I+1/2 . Since the parametrization is inherently stochastic, we focus our attention of reproducing the statistical properties of the resolved variables. Emphasis on statistical properties is, for instance, motivated by necessity to develop effective model for ensemble forecasts in climate and weather predictions.

Forced Burgers' equation
Viscous Burgers equation is one of the well-studied nonlinear partial differential equations and it is often used as a prototype model for turbulence [4]. We consider a stochastically forced Burgers equation over a periodic domain [0, L] given by where ν is the viscosity coefficient and ρ is a large-scale stochastic forcing. To produce a system admitting the form given in (3) we adopt the flux discretization used in [41] and the discrete approximation for the flux becomes Large-scale forcing is applied to low wavenumbers k ∈ K. The forcing term is chosen such that ρ yi = 0 for all i, and where ∆t is the integration time step of the full model, A is the amplitude, α k and β k are independent random coefficients sampled from N (0, 1).
For the Burgers' equation in (9) with the flux discretization (10) the self-interactions of the resolved variables and the subgrid flux in (7) become where and the subgrid flux depends on resolved variables U I and U I+1 , as well as on the unresolved modes y r I and y l I+1 . Here y r I and y l I+1 represent residuals located the right boundary of cell I and left boundary of cell I +1, respectively. Thus, for simpler notations, instead of using fine mesh index, we use l and r, signifying lef t and right, to refer to location of the residuals with respect to its coarse cell index subscript. One essential thing to note about applying the averaging operator used in (4) is that the resulting flux difference scheme admits a viscosity coefficient n times larger than that in the full model. This is the direct consequence of applying the discrete averaging operator (4) to the full dynamics (3) with flux definition in (10).
The subgrid flux term can further be divided into subgrid contributions to the nonlinear part, u 2 /2, and to the viscosity term, ν∂ x u. That is, we can represent the subgrid flux as There are several reasons for such decomposition. First, terms on the right-hand side of (12) have different complexity -the first term involves quadratic interaction between resolved and unresolved modes, and the second term only involves linear relationship between unresolved modes. Second, we can easily investigate scalability of the stochastic parametrization with respect to varying n.
Given the discrete flux and forcing above, we use Runge-Kutta 3 time-stepping scheme for integrating the full model. We use the following parameters in our simulations We perform direct numerical simulations (DNS) of the full model in (9) and generate long stationary timeseries of u i (t). We then compute the resolved variables using (4) and compute numerically their statistical properties, such as spectra, distribution, correlation function, and kurtosis. We then train a conditional WGAN to reproduce subgrid fluxes and perform simulations of the effective model in (8). We then compare statistical behavior of resolved variables in the DNS and simulations of the effective model.

WGAN parametrization
In this section, we present the structure of the network and describe the training procedure. we also describe data used for training and verification.

Network
We use the conditional variant of WGAN as our network with the aim of generating samples dependent on the state of local averages. Moreover, our goal is to achieve a local parametrization, that is, produce approximate subgrid flux at each I + 1/2 conditioned only on immediate neighboring local averages U I and U I+1 . Since the resolved and unresolved modes exhibit spatial homogeneity, it is sufficient to construct and train a generator using data for a particular cell I. In our case we consider I = 0 and we use time series of U 0 , U 1 , G 1/2 , and G 1/2 for training and testing.
The generator accepts a vector Z consisting of two random number independently drawn from the uniform distribution on [−1, 1], and the conditions U 0 and U 1 at a specific time t. It outputs two values -G 1/2 . It then decides how likely are the generated values drawn from the distribution of subgrid flux tendencies by giving "scores" to all the inputs. (Note that in vanilla GAN, discriminator gives probabilities of being "real" to inputs. Here we use "scores" because our training uses Wasserstein GAN. More details is presented in the next subsection.) Both the generator and discriminator are fully-connected feed-forward networks with 3 hidden layers, 16 1/2 , and G 1/2 ) are saved until 6 × 10 5 samples for each are obtained. The first 10 5 is used for training and the rest is for validation. In training, we use Wassertein GAN with gradient penalty (WGAN-GP) since its loss function has more interpretability than that of the game-theoretic loss of vanilla GAN. The gradient penalty is included in the loss function to impose Lipschitz condition required for convergence of Wassertein distance between distributions. We refer to [2,20] for thorough discussion on this matter. The code for training and training data set is available at github.
Stochastic gradient descent method is implemented for optimization with minibatch size of 400. Initial learning rate in Adam optimizer is set to 2 × 10 −5 . During training, Wassertein distance between the distribution of "true" subgrid flux tendencies and the generated subgrid flux tendencies using the validation set is computed every epoch. Optimized weights and biases at epoch 100 are used to generate samples of subgrid flux tendencies in the reduced model (11).

Results
Aside from bare truncation model (BTR) which integrates the reduced equation (8) in time without any subgrid flux corrections, we also compare the performance of GAN parametrization with a polynomial regression with additive noise (POLY) similar to one used in [39]. This reduced model parametrizes subgrid scales by fitting a third-order polynomial with Gaussian additive noise estimated from regression residuals. The deterministic part of POLY parametrization is a multivariate polynomial for each G (i) 1/2 (i = 1, 2) regressed against U 0 and U 1 . Coefficient estimation was done in MATLAB. The same polynomial parametrization (with independent noise realizations) is applied for each coarse cell.
Time integration: For BTR, Runge-Kutta of order 3 (RK3) is used. For the other two reduced models (POLY and GAN), a splitting scheme is implemented. That is, the deterministic part of the subgrid flux tendency is integrated using RK3 while the stochastic term (additive noise or GAN-generated samples) is then added to the result in the manner similar to how Euler-Maruyama handles stochastic term in simulating SDEs. Time step used in simulations of reduced models is equal to the time step in the numerical integration of the full model.

Stationary statistics
Statistical properties of trajectories are often used to access the performance of reduced models. To this end, we performed five long simulations with different realizations of forcing. All statistics are taken as ensemble and spatial averages (since the system is spatially homogeneous).
Energy spectrum and probability density for the resolved variables in the full model and different reduced models are presented in Figure 2. Stationary variance and forth moment for U I are also presented in Table 1. Energy spectra clearly demonstrates that a sub-grid parametrization is necessary, since the bare-truncation (BTR) underestimates the spectrum at high wavenumbers. The role of the parametrization is then to introduce corrections which mimic the sub-grid fluxes and reproduce the statistical properties more accurately. Both polynomial (POLY) and neural-network (GAN) reduced models reproduce energy spectrum quite accurately with GAN reduced model providing a slightly better approximation at the highest wavenumbers. Figure 2 shows that the pdf of the reduced variables is approximately normal. It is evident that the shape is best resembled by the GAN parametrization. In addition, the GAN parametrization also reproduces one-point statistics (see Table 1) slightly more accurately compared to the polynomial reduced model. Mean values and third moments of reduced variables are omitted as all of them are very close to zero, similar to the DNS results. Because of higher viscosity and absence of any corrections, bare truncation is expected to underestimate the variance and hence results in a PDF which has higher values near the mean. Polynomial model provides some improvement for one-point statistics compared to the bare truncation, but still the pdf is slightly narrower than the results of the DNS. Typically, it is more challenging to reproduce two-point statistical properties, compared to the one-point stationary statistics. Therefore, we also compute two-time stationary statistics and compare performance of the three reduced models with the results of DNS. Two-time stationary correlation function and kurtosis  Table 1: Variance and fourth moments of resolved modes in full and reduced models. are depicted in Figure 3. Kurtosis is a measure of non-Gaussianity and it is given by The value of K(s) is exactly one for Gaussian models. All three reduced models reproduce stationary auto-correlation very well. The GAN parametrization predicts the non-Gaussian features better compared to the polynomial and bare-truncation reduced model. Relative errors for K(s) for the POLY model are approximately 8-10% compared to approximately 3% relative errors for K(s) in the GAN model. This is particularly evident for short lags. Clearly, GAN reduced models performs very well in terms of capturing long-term statistical properties of local averages. Polynomial model also performs reasonably well. However, recall that polynomial model includes additive Gaussian noise. Stationary statistics indicates that local averages are almost Gaussian -we can see that the pdf is very close to a Gaussian distribution and kurtosis does not deviate from significantly from one. We can expect that polynomial reduced model will perform worse for stronger non-Gaussian models. On the other hand, due to the general nature of the noise in the GAN reduced model, we expect that this model should perform adequately for non-Gaussian processes as well. This will be investigated in subsequent papers. Reduced models significantly accelerate computations. In particular, reduced GAN model is approximately three times faster compared with the DNS.

Sensitivity to forcing
It is important to develop parametrizations which can be easily adapted to different parameter and forcing regimes without additional re-training. Parametrization's ability to perform adequately with respect to changes in parameters/forcing/model resolution is an important practical factor. Practical significance of parametrization is significantly enhanced if it can perform well in a relatively wide parameter regime without re-fitting or re-training. Therefore, we test the performance of the GAN parametrization with respect to increase in the magnitude of the forcing. In particular, we increase the magnitude of the forcing, A, by 10% and 20% and compare the performance of the GAN reduced model with DNS. We would like to emphasize that the GAN reduced model is not re-trained, i.e. we use GAN reduced model trained using data generated with the original forcing amplitude.
Polynomial model is excluded from the comparison here because the polynomial reduced model produced cases when the system's energy explodes. With increased forcing amplitude, there are occurrences of local averages higher in magnitude than those in the data generated at the original forcing. This leads to the incorrect estimation of sub-grid fluxes by the polynomial regression and the polynomial model often diverges. The most likely reason for a poor performance of the polynomial reduced model is inadequate estimation of sub-grid fluxes in regions with high gradients. When an under-or overestimation of subgrid flux corrections occurs in polynomial model, near-discontinuties in regions with high gradients may not be sufficiently "tamed" by viscosity which leads to more instances of non-physical shocks or oscillations. Although this event can be hard to disprove in the case of GAN, it was observed that GAN did not encounter this problem.
Numerical results for comparison of the GAN reduced model with DNS in the regime with increased forcing are presented in Figure 4 and Table 2. The GAN reduced model performs very well in both cases. We would like to point out that increase in forcing magnitude by 10% and 20% results in the increase in variance of the resolved variables by approximately 12% and 26%, respectively. One-point statistical quantities are reproduced perfectly by the GAN reduced model. Two-point statistics is also reproduced very-well, including the kurtosis K(s) and these results are not presented here only for the brevity of presentation.   Results in this section demonstrate that the GAN reduced model is able to reproduce the behavior of resolved variables in DNS with different forcing amplitude. We would like to point out that changes in forcing amplitude resulted in substantial increase in the variability of resolved variables. Therefore, the GAN reduced model proves capable of extrapolating subgrid flux tendencies outside of the parameter regime where it was trained.

Scalability
Another important practical property of parametrizations is scalability. Scalability reflects the ability to utilize the same parametrization without re-training in a reduced model with different spatial (and possibly temporal) resolution. Typically, if the spatial resolution is changed (e.g. coarser or finer mesh than originally designed is used), we would like to use the same parametrization scheme (possibly properly scaled) to perform simulations of reduced models on a new mesh without re-training or re-fitting. Resolution decrease is motivated by necessity to run fast crude simulations in, for example, ensemble predictions, and increase in resolution is motivated by technological advances and future increase in computational capabilities. Thus, in the future it should be possible to run reduced models with higher spatial resolution, and it is important to verify that stochastic parametrizations utilized in reduced models can adequately respond without retraining.
Here we demonstrate that simple linear scaling in GAN parametrization is sufficient to adapt GAN reduced model with respect to changes in mesh size and perform simulations with a different spatial resolution. Of course, we cannot expect that the same parametrization would perform adequately if the spatial resolution is changed drastically. Therefore, we test several cases -when the mesh becomes twice finer and two and four times coarser. In particular, we consider the following setup. Recall, that we denote the averaging window as n (see (4)). If the averaging window is changed to n , then we use a linear scaling and the generated subgrid fluxes (output of GAN) is multiplied by a factor n/n . This linear scaling is motivated by the application of a systematic stochastic model-reduction strategy to finite-difference discretization of the Burgers-Hopf equation [12,11] and the shallow-water model [42].
The reference setup in section 3.1 corresponds to n = 16 and here we consider cases n = 4, n = 8, and n = 32. Second and fourth moments in simulations with n = 16 and n = 4, 8, 32 are presented in Table 3. In addition, in this Table we also present results for simulations with non-scaled (NS) GAN parametrizations. We can see that the scaling introduced earlier improves the performance of the reduced model for all cases n = 4, 8, 32. Numerical results for n = 4 and n = 8 are very close to the results for n = 16 and, thus, scaling of the GAN parametrization does not produce a large effect. On the other hand, numerical results for n = 32 clearly indicate that scaling of the GAN parametrization improves the performance of the reduced model. This is particularly important for simulations with the coarser resolution n = 32. We would like to point out that variance of reduced variables does not increase significantly for cases n = 4 and n = 8 (only 2.7% and 1.5% relative increase, respectively), but for n = 32 variance is reduced by approximately 8% compared with n = 16. This is likely the reason why the reduced model for n = 32 is more sensitive to the scaling of GAN parametrization compared to simulations with n = 4 and n = 8. Two point stationary statistics are also reproduced very well by reduced models with scaled GAN parametrization. Overall, there is a good agreement between the direct numerical simulations and simulations of the reduced model with scaled GAN parametrization. In addition to results presented in Table 3, utilizing scaled GAN parametrization also improves the shape of the spectrum for high wavenumbers for n = 4 and n = 8 (see [1]). Simulations with unscaled GAN parametrization with n = 32 completely fail to reproduce stationary spectrum and probability distribution of resolved variables [1]. These results are not presented here only for the brevity of presentation.

Conclusions
A stochastic subgrid flux parametrization is developed here using conditional GAN. We demonstrated that Generative Adversarial Networks can be successfully trained using the Wasserstein metric to reproduce subgrid corrections for fluxes of resolved variables. The WGAN version of the Adversarial Networks significantly improves conversion during the training process. In addition, to estimate the subgrid flux between two neighboring cells, the GAN is conditioned on the current value of the solution in these two cells. This conditioning is essential for reproducing the subgrid fluxes. Resolved variables are defined here as local spatial averages. This makes the reduced model local in space, i.e. only two neighboring cells are coupled by a GAN parametrization which estimates the flux between these two cells. Thus, if the model is non-homogenious in space, the GAN parametrization can potentially be trained separately for different spatial regions.
Our aim here is to reproduce statistical properties of long-term turbulent simulations. In particular, we compare spectra, moments, and two-point temporal statistical quantities. We would like to point out that we simulate the reduced model for a long time and the reduced model does not blow up; all runs of the reduced model produce stable trajectories. The reduced model with GAN parametrization reproduces statistical properties of the full model quite well. We illustrate our approach using the forced Burgers-Hopf model. In addition, we investigate numerically the performance of the GAN parametrization with respect to changes in forcing and different spatial resolutions. To this end, we utilize the original parametrization to simulate reduced models with 10% and 20% increase in forcing. We also introduce linear scaling for the GAN parametrization to compute reduced models with different spatial resolutions. We demonstrate that this linear scaling is essential for the reduced model to accurately reproduce statistical properties of solutions. We would like to emphasize that the GAN parametrization is not re-trained on the new data (e.g. simulations with increased forcing or changed resolution). Therefore, we demonstrate that the reduced models with the GAN paramterization can be used in a wide range of regimes without re-training the neural network.
The approach developed here is particularly suitable for finite-volume (or finite-difference) discretizations of partial differential equations since resolved variables are defined as local spatial averages and subgrid fluxes are represented using corresponding averages of the right-hand side of the equation. We demonstrated applicability of this approach for a prototype nonlinear model which exhibits complex turbulent behavior. Thus, this approach stands ready to be applied in more complex models of fluid dynamics. This will be investigated in subsequent papers.