Stochastic subgrid-scale parameterization for one-dimensional shallow water dynamics using stochastic mode reduction
Matthias Zacharuk, Stamen I. Dolaptchiev, Ulrich Achatz, Ilya Timofeyev
SStochastic subgrid-scale parameterization for one-dimensionalshallow water dynamics using stochastic mode reduction ∗ Matthias Zacharuk † Stamen I. Dolaptchiev ‡ Ulrich Achatz ‡ Ilya Timofeyev § August 17, 2018
Abstract
We address the question of parameterizing the subgrid scales in simulations of geophysicalflows by applying stochastic mode reduction to the one-dimensional stochastically forced shallowwater equations. The problem is formulated in physical space by defining resolved variables aslocal spatial averages over finite-volume cells and unresolved variables as corresponding residuals.Based on the assumption of a time-scale separation between the slow spatial averages and thefast residuals, the stochastic mode reduction procedure is used to obtain a low-resolution modelfor the spatial averages alone with local stochastic subgrid-scale parameterization coupling eachresolved variable only to a few neighboring cells. The closure improves the results of the low-resolution model and outperforms two purely empirical stochastic parameterizations. It is shownthat the largest benefit is in the representation of the energy spectrum. By adjusting only asingle coefficient (the strength of the noise) we observe that there is a potential for improvingthe performance of the parameterization, if additional tuning of the coefficients is performed.In addition, the scale-awareness of the parameterizations is studied.
Atmospheric processes encompass a large spectrum of spatial and temporal scales. These rangefrom several millimeters and seconds for boundary layer turbulence up to 10 meters and severalweeks (and even longer) for planetary wave dynamics. Due to limited computer resources numericalatmospheric models cannot describe all these processes on all scales simultaneously. However, thedifferent scales are interacting in a complex manner and this leads to the challenging problemof parameterizing the effect of the unresolved subgrid-scale (SGS) processes onto the resolvedones. Examples include the parameterization of synoptic and mesoscale eddies in planetary scaleatmospheric models (e.g. [47, 60]), momentum and temperature fluxes in the atmospheric boundarylayer [54] or SGS Reynold stresses in large eddy simulations (e.g. [49]).In this context, stochastic elements have become increasingly popular. Stochastic parameteri-zations can reduce a systematic model error, represent uncertainty in predictions, or trigger regimetransitions (e.g. [42, 6]). Typically some ad-hoc SGS model is assumed and the correspondingcoefficients are optimized (tuned) so as to obtain the best possible agreement, in some sense, with ∗ Submitted to Q.J.R. Meteorol. Soc. † Corresponding author; Institut f¨ur Atmosph¨are und Umwelt, Fachbereich Geowissenschaften/Geographie,Goethe-Universit¨at ([email protected]). ‡ Institut f¨ur Atmosph¨are und Umwelt, Fachbereich Geowissenschaften/Geographie, Goethe-Universit¨at. § Department of Mathematics, University of Houston, Houston, TX 77204-3008, ( [email protected] ). a r X i v : . [ phy s i c s . f l u - dyn ] A ug bservations or high-resolution simulations. Examples in comprehensive climate and weather mod-els are stochastically perturbed parameterization tendencies [7, 43] or stochastic kinetic energybackscatter [52, 5]. Empirical Ornstein-Uhlenbeck (OU) processes have been used in some studiesof low-frequency and large-scale atmospheric variability (e.g. [61, 41, 46]), which can be extendedto include quadratic nonlinearities as well as a time correlated stochastic forcing [32, 30].With regard to SGS parameterizations in climate models issues can arise from the fact thatthey are typically tuned to optimally represent the statistics of the present-day climate. If climatechanges due to some external forcing, it is not guaranteed that the tuned parameters are stilloptimal. The fluctuation-dissipation theorem might be able to provide corrections [3, 48] in somecases, but such an approach relies on the perturbations being sufficiently weak. Moreover, thereis a need for scale-aware parameterizations in atmosphere modeling, as model resolution increasescontinuously and mesh refinement techniques become widely used. In addition, the consistencybetween particular SGS parameterizations and the numerical discretization becomes important.These considerations motivate the development of other approaches where the SGS parameter-ization is derived from first principles, if possible without any empirical parameter optimization.The direct interaction approximation (DIA) introduced by [31] allowed to successfully apply sta-tistical dynamical closure theory in relevant geophysical flows [17, 20, 18]. In the presence of timescale separation, the asymptotic method of averaging has been applied [23, 25, 4, 40]. This methodrequires an estimation of the invariant measure of the fast scales conditioned on the slow scales,which might limit its applicability when going to high-dimensional systems. Another promising ap-proach, without any empirical component, is based on the maximum entropy principle [56, 58, 57].Recently, [63, 64] have introduced a new method originating from response theory. This methodrelies on a weak coupling between resolved and unresolved scales and it has been applied to simpleand more complex settings [62, 11, 59].The DIA parameterization has been successfully applied to barotropic [17] and primitive equa-tions model [18], it has been extended to include the effects of mean flow and topography [20].The DIA closure is derived in spectral space by considering the evolution of second-order cumulantand response function. Next, the nonlinear damping rate and nonlinear noise are introduced. Thisresults in a globally coupled SGS model in spectral space. However, techniques have been proposedto simplify the equations and obtain locally coupled models reproducing the spectra from directnumerical simulations [17, 18].Another nearly self-consistent possibility that exploits a separation of time scales betweenresolved and unresolved scales is the stochastic mode reduction (SMR) procedure proposed by[37, 35, 38]. The SMR is a homogenization technique for multiscale systems ([27, 26, 33, 44] anda recent overview [45] and references therein) and it is supplemented by an empirical step, wherethe fast SGS self-interactions in the evolution equation for the unresolved modes are replaced byan OU process. Following this step, an analytical derivation of a stochastic parameterization forthe fast modes is possible, rigorously valid in the limit of infinite time scale separation.So far the SMR procedure has already been applied to balanced models, such as barotropic [16]and quasi-geostrophic [15] dynamics. The separation between resolved and unresolved scales hasbeen performed by using empirical orthogonal functions (EOFs). However, EOFs are sometimes notable to guarantee a sufficient separation of the underlying time scales (e.g. Figure 3 in [15]). TheSMR carried out in spectral space [16, 15] is quite similar to the DIA closure approach. In particular,the main goal of both techniques is to represent the subgrid processes by a nonlinear damping anda state-dependent noise and both techniques have been utilized successfully in geophysical flows.In applications of the SMR to spectral space the resulting reduced model is globally coupledwith linear, quadratic and even cubic terms. This hampers the applicability of the technique whenhigh-dimensional systems with large number of resolved modes are considered. However, the latter2roblem can be avoided by applying the SMR in physical space to a finite-volume discretization ofthe equations. Such discretization does not per se include global coupling as in spectral discretiza-tions, since grid cells interact directly only with a small number of neighbors. Finite-volume schemesare traditionally applied in ocean models (e.g. [22]), regional atmospheric modeling (e.g. [53]) andrecently even for global atmospheric models as well [50, 39, 51]. In the examples above complexboundaries, such as continental boundaries in ocean modeling or non-periodic lateral boundaries inregional area atmospheric modeling, necessitate the use of discretizations and SGS parameteriza-tions formulated in physical space. This motivated [12, 13] to consider a local approach where theresolved variables are defined by local spatial averages and the SGS flow by deviations from these av-erages, a configuration typically encountered in large-eddy turbulence parameterization (e.g. [49]).The local definition leads to a local SGS parameterization, coupling only near neighbors, as shownfor the Burgers equation [12, 13]. The efficient local stochastic SGS parameterization allows toconsider large numbers of resolved scales. In addition, the clear gap of spatial scales between theresolved and unresolved variables enables a more pronounced time-scale separation.Obviously Burgers equation represents a highly idealized prototype model for testing variousstatistical and closure methods and it is necessary to verify the applicability of the SMR for localspatial averages for more realistic fluid-dynamical models. One step in this direction is performed inthis work by applying the approach to a stochastically forced one-dimensional shallow water layer(1DSW). It incorporates at least two issues of relevance in the general context. First, in contrast tothe Burgers equation the 1DSW allows for gravity waves. Secondly, if formulated in flux form, theshallow-water flow dynamical equations entail non-polynomial nonlinearities. This problem is ofbroader relevance, since such highly nonlinear terms appear in the general compressible fluid flowequations as well, in the pressure-gradient acceleration.The work presented here can be summarized as follows. Based on a high-resolution finite-volumediscretization of the shallow-water equations we use in Sec. 2 local spatial averages to define coarseand slow (resolved) variables and, via corresponding residuals, fine and fast (unresolved) variables.The assumed time-scale separation is verified numerically. The SMR theory for obtaining an SGSparameterization of the unresolved modes is then introduced and applied to the specific problem.In Sec. 3 we discuss the practical implementation, and also introduce, for comparison, two purelyempirical approaches. Results from model simulations with the various SGS parameterizations arethen compared in Sec. 4. Here we also investigate the scale awareness of the approaches, i.e. theirability to be applied at different resolutions without any re-tuning. Conclusions are finally drawnin Sec. 5. We consider a stochastically forced one dimensional shallow water layer with periodic boundaries,using as variables the height of the fluid h and the momentum hu , where u is the velocity. Thegoverning equations with plane topography (e.g. [55]) read in flux form ∂ t (cid:18) hhu (cid:19) = − ∂ x (cid:32) hu − ν∂ x h ( hu ) h + g h − ν∂ x hu (cid:33) + (cid:37)(cid:37)(cid:37), (1)with a large-scale stochastic forcing (cid:37)(cid:37)(cid:37) (see Sec. 3) and a mass weighted diffusion with the constantparameter ν . 3or a high-resolution spatial discretization the domain of length L is divided into N fine intervals∆ x = L/N , labelled by a small index i ∈ { , , ..., N − } . With this the equations in (1) can bediscretized by a symmetric finite-volume scheme ddt (cid:18) h i hu i (cid:19) = − x (cid:16) F i + − F i − (cid:17) + (cid:37)(cid:37)(cid:37) i , (2)with the discrete forcing (cid:37)(cid:37)(cid:37) i and the flux at the boundary given by F i + = 12 ( hu ) i +1 + ( hu ) i − ν h i +1 − h i ∆ x ( hu ) i +1 h i +1 + ( hu ) i h i + g h i +1 + g h i − ν ( hu ) i +1 − ( hu ) i ∆ x . (3)The discrete flux form (2) conserves total mass N (cid:80) N − k =0 h k and total momentum N (cid:80) N − k =0 hu k in the absence of forcing. Given our choice of forcing and dissipation, the number of fine cells ischosen large enough so as to resolve all processes occurring. Hence in the following simulationsusing (2) , with N large enough, will be called direct numerical simulation (DNS). As a representation of the typical situation of atmospheric models with insufficient resolution, weintroduce a second discretization, with N c = N/n coarse cells, each consisting of n fine cells of theinitial discretization, and labelled by the capital index I ∈ { , , ..., N c − } . Associated with thecoarse grid coarse variables H and HU , also called resolved variables, are defined by local spatialaverages inside a coarse box (cid:18) H I HU I (cid:19) = 1 n n ( I +1) − (cid:88) k = nI (cid:18) h k hu k (cid:19) . (4)Further, fine variables h (cid:48) and hu (cid:48) , referred to as unresolved or SGS variables, are defined usingthe deviations of the initial variables from the corresponding coarse variables (cid:18) h (cid:48) i hu (cid:48) i (cid:19) = (cid:18) h i hu i (cid:19) − (cid:18) H I [ i ] HU I [ i ] (cid:19) . (5)Here I [ i ] denotes here the index of the coarse cell with the i -th fine cell placed inside. The coarseand fine variables can be used to express (2) as ddt (cid:18) H I HU I (cid:19) = − F ( I +1) n − − F nI − n ∆ x + (cid:37)(cid:37)(cid:37) I , (6) ddt (cid:18) h (cid:48) i hu (cid:48) i (cid:19) = − F i + − F i − ∆ x + F ( I [ i ]+1) n − − F I [ i ] n − n ∆ x , (7)where we assume that the forcing (cid:37)(cid:37)(cid:37) I acts only onto the coarse variables. By collecting all resolvedvariables H I , HU I in one vector x ∈ R M c and all SGS variables h (cid:48) i , hu (cid:48) i in another vector z ∈ R M with M c = 2 N c and M = 2 N , (6) and (7) can be rewritten as˙ x i = (cid:37) xi + a xi ( x ) + b xzi ( x , z ) , (8)˙ z i = b zi ( x , z ) + c zi ( z ) . (9)4ere (cid:37) xi results from the forcing and terms have been regrouped so as to identify the coarse variableself-interactions a xi ( x ), the coupling terms b xzi ( x , z ), b zi ( x , z ) and the fine variable self-interactions c zi ( z ). Complete neglect of the SGS variables yields the bare-truncation model˙ x i = (cid:37) xi + a xi ( x ) . (10)This low resolution model is defined on the coarse grid with N c grid cells and it lacks an SGSparameterization. A difficulty in the application of SMR is caused by the terms involving 1 /h in the flux (3) , asthey represent a nonlinearity of an arbitrary order. So far SMR has only been applied to systemswith quadratic nonlinearities. It can, however, handle nonlinearities of arbitrary polynomial form.Hence a solution could be expanding everywhere but in a xI ( x ) terms with 1 /h in a finite Taylor seriesaround the mean fluid height H . It turns out sufficient, however, to simply replace 1 /h ≈ / H in order to reproduce the statistics of the DNS. Thus, the bare truncation part of the model iscomputed exactly but in all other terms involving SGS modes this approximation is used, leadingto an approximation of (8) and (9), where SGS-variable nonlinearities take a quadratic form˙ x i = (cid:37) xi + a xi ( x ) + (cid:0) L xzij z j + B xxzijk x j z k + B xzzijk z j z k (cid:1) , (11)˙ z i = (cid:0) L zxij x j + B zxxijk x j x k + B zxzijk x j z k (cid:1) + (cid:0) L zzij z j + B zzzijk z j z k (cid:1) . (12)Here and in the following we make use of Einstein’s summation convention and the summationindex is running up to either M c or M depending on if x or z is involved. The linear and quadraticinteraction coefficients L xzij , L zxij , L zzij , B xxzijk , B xzzijk , B zxxijk , B zxzijk , B zzzijk are given in Appendix A. The SGS variables z i are not independent, since the corresponding local spatial average over a coarsecell vanishes by definition. Thus, one degree of freedom for each coarse cell has to be eliminated.This is achieved by Fourier-transforming h (cid:48) and hu (cid:48) locally inside each coarse cell. The Fourieramplitude of the zero-wavenumber is equal to the vanishing local average inside the coarse cell.Hence by discarding this wavenumber component one degree of freedom can be eliminated. Thisdefines the new independent SGS variables θ i θ i = ˆ T ij z j (13) z i = ˆ R ij θ j , (14)where the matrices ˆ R ∈ R M × M f and ˆ T ∈ R M f × M are constructed from the inverse and forwardFourier transformation and M f = M − M c is the number of independent SGS variables. Withthis preparation one can move to the next step of SMR, i.e. replacing the SGS self-interactions L zzij z j + B zzzijk z j z k by an empirical OU process. The SGS equation (12) becomes dθ i = ˆ T ik (cid:16) L zxkj x j + B zxxkjl x j x l + B zxzkjl x j ˆ R lm θ m (cid:17) dt + Γ ij θ j dt + σ i dW i , (15)where Γ ij denote the coefficients of the negative-definite OU drift, σ ij = σ i δ ij those of a diagonaldiffusion tensor, and dW i Wiener increments. Note that no sum over i is taken in the Wiener term5n (15). We assume that Γ couples only SGS modes corresponding to the same coarse cell. Underthis assumption Γ has a block-diagonal form and the resulting SGS closure is local, coupling onlyneighbors and next-neighbors of a coarse cell. Since SMR assumes further that the OU process isthe dominant term in the SGS equation (see below), the OU drift Γ can be estimated from thelagged covariance of θθθ [21] θθθ ( t ) θθθ T ( t + τ ) = C ( τ ) = C (0) e Γ T τ , (16)where ( · ) denotes a time average. By integrating over time one can solve for Γ (cid:0) Γ T (cid:1) − = − C (0) − (cid:90) ∞ C ( τ ) dτ . (17)The computation of the time integral in (17) is performed using the numerically efficient Cooper-Haynes algorithm [34]. Note that despite the block-diagonal form, Γ still allows for a couplingbetween both SGS variables h (cid:48) i and hu (cid:48) i inside each coarse cell. Because of the spatial homogeneityof the considered shallow-water model the coefficients of Γ are the same for each matrix block andcan be obtained by averaging over the estimates from the different coarse cells. Using the steadyLyapunov equation Γ C (0) + C (0)Γ T = − σσ T , the diagonal diffusion coefficients are found from σ i = (cid:112) − ik C (0) ki . (18)Moreover, it has turned out to be useful to observe that in typical applications Γ is diagonalizableand has distinct eigenvalues (e.g. [9]). We hence introduce new variables y i y i = U − ij θ j , (19)where the real invertible matrix U is from the real Jordan canonical form decomposition of ΓΓ = U Λ U − . (20)By applying the transformation matrices R = ˆ RU and T = U − ˆ T , the model equations (11)and (15) can be written in terms of the new variables as˙ x i = (cid:37) xi + a xi ( x ) + b xi ( x , y ) , (21)˙ y i = b yi ( x , y ) + Λ ij y j + Σ i ˙ W i . (22)Here we use the notation b xi ( x , y ) = L xzij R jk y k + B xxzijk R kl x j y l + B xzzijk R jl R km y l y m = L xyij y j + B xxyijk x j y k + B xyyijk y j y k , (23) b yi ( x , y ) = T ij (cid:0) L zxjk x k + B zxxjkl x k x l + B zxzjkl R lm x k y m (cid:1) = L yxij x j + B yxxijk x j x k + B yxyijk x j y k , (24)and we have also introduced effective drift coefficients Σ i = (cid:113) U − ij U − ij σ j with pairwise identicalnoise parameters for pairs of complex eigenvalues, see Appendix C in [12] for the details.6 .3.3 Homogenization The remaining step is the derivation of an effective equation for the coarse variable x alone, usingthe homogenization technique [37, 45], with terms taking SGS effects into account. The mainassumption of the SMR is the presence of distinct time scales in the considered variables. So farthe model is spatially separated into coarse and fine variables. This does not necessarily imply aseparation of the underlying time scales. However, as it will be shown later, for the consideredregime the separation in space also induces a separation in time, with the resolved variable x acting on a slower time scale than the SGS variable y . Hence a time-scale separation factor (cid:15) (cid:28) b x → b x /(cid:15) , b y → b y /(cid:15) , and Λ ij y j +Σ i ˙ W i → Λ ij y j /(cid:15) +Σ i ˙ W i /(cid:15) ,where then b x , b y , Λ ij y j , and Σ i ˙ W i are all O (1), obtaining˙ x i = (cid:37) xi + a xi ( x ) + 1 (cid:15) b xi ( x , y ) , (25)˙ y i = 1 (cid:15) b yi ( x , y ) + 1 (cid:15) Λ ij y j + 1 (cid:15) Σ i ˙ W i . (26)The above scaling implies that the bare truncation part (cid:37) x + a x ( x ) acts on the slowest, thecoupling terms b x and b y on a faster and the SGS self-interactions on the fastest time scale. In thefollowing the corresponding backward Fokker-Planck equation (FPE) ∂ t p = L p + 1 (cid:15) L p + 1 (cid:15) L p , (27)for the probability density function (PDF) p is considered, in the limit of an infinite time scaleseparation (cid:15) →
0. The operators on the right hand side are defined as L = − ( (cid:37) xi + a xi ( x )) ∂ x i , (28) L = − b xi ( x , y ) ∂ x i − b yi ( x , y ) ∂ y i , (29) L = − Λ ij y j ∂ y i − Σ i ∂ y i . (30)Next, the PDF in (27) is expanded in terms of (cid:15) : p = p + (cid:15)p + (cid:15) p + ... , which leads to thefollowing set of equations O ( (cid:15) − ) : 0 = L p , (31) O ( (cid:15) − ) : 0 = L p + L p , (32) O ( (cid:15) ) : p = L p + L p + L p . (33)The leading order O ( (cid:15) − )-equation shows that p is in the null space of L and therefore itdoes not depend on the fast variables: p = p ( x ). The O ( (cid:15) − )-equation can be solved for p : p = − L − L p if the solvability condition P L p = 0 (34)is satisfied, where the projection operator P , projecting onto the null space of L is utilized. Withthis result the last O ( (cid:15) )-equation can be written as an effective FPE for p only7 t p = L p − P L L − L p . (35)This is the backward FPE of a low-resolution model, for the coarse variables alone, whichconsists of the bare truncation part L p and an SGS parameterization P L L − L p . The null-space projection P and inverse of the OU operator L are detailed in Appendix B. Using these,one finds that the stochastic differential equation corresponding to the effective FPE (35) can bewritten as dx i = [ (cid:37) xi + a xi ( x ) + β i ( x )] dt + dξ i ( x ) . (36)Here β i represents the deterministic part and dξ i the stochastic part of the SGS parameterization,containing both additive and multiplicative noise terms. One finds that the deterministic closure β i is β i = (cid:90) ∞ dτ (cid:28) b xj ( x , y ) ∂b xi ( x , y ( τ )) ∂x j (cid:29) + (cid:104) yy T (cid:105) − jm (cid:90) ∞ dτ (cid:68) y m b yj ( x , y ) b xi ( x , y ( τ )) (cid:69) − (cid:90) ∞ dτ (cid:42) ∂b yj ( x , y ) ∂y j b xi ( x , y ( τ )) (cid:43) , (37)where the expectations (cid:104)·(cid:105) are taken over the OU statistics, and y represents an OU trajectory withinitial condition y = y (0). The stochastic closure dξ i takes the form dξ i = √ B ij dW j , (38)where the matrix elements B ij are obtained from the decomposition B ik B jk = (cid:90) ∞ dτ (cid:10) b xi ( x , y (0)) b xj ( x , y ( τ )) (cid:11) . (39)With these results one can also show that back-transforming b x /(cid:15) → b x , b y /(cid:15) → b y , andΛ ij y j /(cid:15) + Σ i ˙ W i /(cid:15) → Λ ij y j + Σ i ˙ W i leaves β and dξ unchanged so that (36) is the desired low-resolution model.Finally getting back to the specific case, (21) – (24) , Appendix C shows that the solvabilitycondition (34) is satisfied. Moreover, inserting (23) and (24) for b x and b y yields β i = (cid:16) L xymj + B xxymlj x l (cid:17) B xxyimk ( C S ) jk + (cid:0) L yxmo x o + B yxxmop x o x p (cid:1) (cid:16) L xyik + B xxyijk x j (cid:17) (cid:104) yy T (cid:105) − mn ( C S ) nk + B xyyijk B yxympo x p (cid:104) yy T (cid:105) − mn ( C T ) onjk , (40)where the tensors C S and C T are given by( C S ) jk = (cid:90) ∞ dτ (cid:104) y j (0) y k ( τ ) (cid:105) , (41)( C T ) onjk = (cid:90) ∞ dτ ( (cid:104) y o (0) y j ( τ ) (cid:105)(cid:104) y k ( τ ) y n (0) (cid:105) + (cid:104) y o (0) y k ( τ ) (cid:105)(cid:104) y j ( τ ) y n (0) (cid:105) ) . (42)With this the decomposition (39) becomes B ik B jk = B xyyimn ( C T ) mnkl B xyyjkl + (cid:0) L xyin + B xxyikn x k (cid:1) ( C S ) nm (cid:16) L xyjm + B xxyjlm x l (cid:17) . (43)8he prescription would be to perform this decomposition every time step. However, this wouldbe very expensive. Therefore we neglect cross-correlations between the dξ i in different cells andapproximate them by dξ i ≈ (cid:112) B ij B ij dW i , = (cid:113) B xyyijn ( C T ) jnkl B xyyikl dW i + (cid:114) (cid:16) L xyin + B xxyijn x j (cid:17) ( C S ) nk (cid:0) L xyik + B xxyilk x l (cid:1) dW i , (44)which we call effective stochastic forcing. In each coarse cell the approximated stochastic term hasthus the same variance as its exact counterpart. The stochastic term (44) consists of an additivepart, which acts on both variables, and a multiplicative part, that acts only on HU . An importantfeature of the SGS parameterization, with deterministic part (40) and stochastic component (44),is that it couples a volume cell only to its neighbors and next neighbors. This allows applicationof the approach to large systems. For the validation of our approach we consider a stochastically forced periodic shallow-water layerof horizontal extent L = 10 km and mean height H = 10 km. The diffusion constant is ν = 10 km day -1 . A large-scale stochastic forcing [8] is applied to the momentum equation (cid:37)(cid:37)(cid:37) I = (cid:32) (cid:80) k =1 µα k √ k ∆ t cos (cid:16) π (cid:16) kIn ∆ xL x + ψ k (cid:17)(cid:17) (cid:33) . (45)Normally distributed random numbers α k and ψ k are used, the amplitude parameter µ is 10 km / day and the forcing acts onto the leading Fourier modes 1 ≤ k ≤
3. Various model set-upshave been chosen as follows, a summary is given in Table 2. In all cases the integrations have beendone over 10 days with 10 outputs per day. Reference is provided by direct numerical simulations, integrating (2) with N = 512 volume cells.A 4th-order Runge-Kutta-scheme is used, with a time step ∆ t = 10 − days. In two intermediate steps in the application of the SMR, first the nonlinearities affected by SGSdynamics have been kept quadratic by replacing 1 /h → / H , and then the SGS nonlinear self-interactions have been replaced by an empirical OU process, leading to the system (21) and (22).Direct integration of these equations, henceforth termed OU-DNS, thus appears as a useful check ofthe validity of the SMR approach. However, directly using the OU parameters estimated from (17)and (18) turned out not to be stable enough. Therefore following [2] an additional scale-selectivedamping has been supplemented to the OU drift in each coarse cell. This has been done in thespectral representation of the latter, see (15), by replacing9 → Γ + γ . . . . . . γ , with γ = − α . . . . . . ( n − ∈ R ( n − × ( n − . The diagonal matrix γ represents damping of the Fourier modes inside each coarse cell with anamplitude proportional to the squared wave number, with here α = 90 day − . As also in all otherstochastic integrations outlined below, the time integration has been done by a split-step methodwith a 4th-order Runge-Kutta step for the deterministic part and an Euler-Mayurama step for thestochastic part. The time step is ∆ t = 10 − day, as in the deterministic DNS. With the high-resolution simulations as reference we can validate the SMR approach for providingan SGS parameterization for low-resolution models that only use the coarse cells. We comparethe performance of this approach also to that of more simple purely empirical parameterizations.Considered approaches are as follows where, unless otherwise stated, the number of coarse cellsemployed was always N c = 64 with an averaging interval of n = 8. Two slightly different approaches have been chosen to obtain low-resolution models. The first isdefined by the original discretized equations (2), but with a lower spatial resolution of N c = N/n .This is henceforth referred to as low-resolution model (LRM) . The second variant is the bare-truncation model (BRT) defined in (10), with a resolution of N c as well. The difference betweenBRT and LRM is the diffusion in the models. In the BRT it is proportional to ν/ ( n ∆ x ), and inthe LRM to ν/ ( n ∆ x ) , implying that the BRT has an effective diffusion by a factor n stronger,when compared to the LRM. Three types of low-resolution simulations with stochastic SGS parameterization have been tested.
SMR parameterization.
The low-resolution model (36) to be validated is the BRT supple-mented by the SMR SGS parameterization consisting of the deterministic and stochastic compo-nents (40) and (44), it is referred to as BRT-SMR. For stability reasons the BRT-SMR diffusivityhad to be increased in corresponding simulations to ν = 2 · km day − . The time step employedwas ∆ t = 2 · − day. Empirical OU parameterizations for BRT and LRM.
As a quality measure for the SMRapproach we also consider low-resolution simulations with an empirical OU SGS parameterization,denoted by BRT-OU or LRM-OU, depending on the low-resolution dynamical core used togetherwith the empirical OU SGS parameterization. As in the SMR parameterization only coupling toneighbors and next neighbors is taken into account. The BRT-OU, for example, can be written as dx i = (cid:16) (cid:37) xi + a xi ( x ) + ˜Γ ij ˆ x Ij (cid:17) dt + ˜ σ i dW i . (46)10here the vector ˆx I has 10 components and encompasses the values of H and HU in the fivecoarse cell from I − I + 2, where I is the cell index corresponding to the variable x i . TheOU parameters ˜Γ and ˜ σ have been estimated using a standard maximum likelihood approach [24],yielding ˜Γ ij = b xi ˆ x k (cid:16) ˆxˆx T − (cid:17) kj , (47)˜ σ i = ∆ t (cid:104) b xi − ˜Γ ij ˆ x Ij (cid:105) , (48)with the superscript of ˆx I suppressed in (47). For the LRM-OU the corresponding parametershave been determined in the same manner. In contrast to the estimation of the OU processesin the SMR by (17) and (18), here the integrated lagged covariance function could not be used,because the SGS effects b xi ( x , y ) as such do not satisfy a prognostic equation dominated by an OUprocess. Replacing (17) and (18) by maximum-likelihood estimates would have been an option aswell. Corresponding tests have shown a slightly deteriorated performance, however. In the following we show step by step the essential results from our various simulation experiments.Autocorrelations and spectra turned out to be qualitatively similar for momentum and surface-height. We therefore focus below on the latter.
Fig. 1 (left) displays the time dependence of the autocorrelation of the resolved variable H andof the SGS variable h (cid:48) in the DNS. A slowly decaying oscillation is visible in the autocorrelationof H . The period of this oscillation is nearly equal to the time τ = L/ √ g H ≈ .
37 day, requiredfor gravity waves to pass once through the domain. This shows that the model has some intrinsicdynamics and is not dominated by forcing and diffusion. The autocorrelation of the SGS variable h (cid:48) decays much faster to zero than that of H . The large difference in the correlation time betweenSGS and resolved variables indicates that the assumption of time-scale separation between x and y is met to a good agreement.The spatial distribution of the variance of h (cid:48) is displayed in Fig. 1 (right). The variance islowest in the middle of a coarse cell, and gradually increases towards the cell boundaries. Thisspatial shape is explained in Appendix E as being due to a spatially decreasing autocorrelation of h . The potential-energy spectrum from the DNS is displayed in the left panel of Fig. 2. With theconsidered forcing and diffusion parameters one obtains an inertial range with spectral index 2 upto around wavenumber kL/ π = 64. There is a small kink in the spectrum after wavenumber 3,due to the forcing acting only onto the first three modes.The deviations from a Gaussian in the fourth order moments of H and HU are less than 4%and 2%, respectively. In addition we find nearly vanishing odd moments (not shown). We concludethat the statistics of the resolved variables are close to Gaussian. As described above, the replacement 1 /h → / H in the SGS nonlinearities leads to the system (11)and (12) with strictly quadratic nonlinearities, as required for the application of the SMR method.11imulations with this model reproduce the DNS data nearly perfectly (not shown).Replacing the SGS self-interactions by an empirical OU process leads to the system (21) and(22). The corresponding OU-DNS reproduces the correlations of the DNS with minor differences(not shown). The energy spectrum from the OU-DNS, projected onto the coarse grid is displayedin Fig. 2 (right). It follows the DNS spectrum for the first 7 wavenumbers and then drops belowit. This indicates a too strong damping at high wavenumbers, which seems to be due to theintroduction of the deterministic part in the OU-process. The spatial variance of h (cid:48) from the OU-DNS model is presented in Fig. 1 (right). It follows with small deviations the structure from theDNS model. Before considering results from low-resolution simulations with the various SGS parameterizations,we first address low-resolution simulations without any parameterizations, to provide a useful ref-erence. The time dependence of the autocorrelation function of H from these simulations is shownin Fig. 1 (left). The amplitude of the oscillation of the auto-correlation from the LRM simulationsis significantly weaker than from the DNS and has a relative error of 6 . . HU even by 111 . N = 512 grid points, it is possible to perform low-resolution DNS with at least256 spatial points. However, further reducing the number of points in the low-resolution DNSsignificantly corrupts the spectra of the first 32 wavenumbers. Thus, this demonstrates the needfor SGS parameterizations if one wants to reduce the spatial resolution beyond N = 256. Energy spectrum.
The potential-energy spectra obtained from low-resolution simulations withthe different parameterizations are shown in Fig. 3 (right). The overall qualitative behavior ofthe shallow-water layer can be reproduced with the SMR parameterization but there is too muchenergy in scales up to around wavenumber 12 and too little energy in higher wavenumbers. On theother hand, the bare truncation model with empirical stochastic corrections (BRT-OU) and the low-resolution model (LRM-OU) do not reproduce all the details of the spectra sufficiently accurately.In particular, the LRM-OU spectrum contains significantly too little energy in all wave numbers.The BRT-OU simulation can reproduce the true spectrum well in the first 15 wavenumbers butfails completely at higher wavenumbers.
Moments.
The statistical moments from the various simulations are summarized in Table 1.In general all closures show high relative errors, the empirical OU parameterizations underestimateand the SMR parameterization overestimates the moments. Errors in the HU -moments are smallestfor the low-resolution simulation using the SMR parameterization. In the H -moments BRT-SMR12as lower errors than LRM-OU but larger than BRT-OU. The BRT-SMR model overestimates thefluctuations in the coarse H -variable, which is consistent with the result for potential energy spectradiscussed above. This implies that the SGS stochastic forcing representing the energy backscatteris too strong. Improvement of Energy Balance.
To further improve the performance of our SGS parametriza-tion using the stochastic mode reduction, we consider BRT-SMR model with reduced SGS stochasticforcing. This is motivated by the fact that there are several assumptions in the SMR approach(e.g. time-scale separation, representing the fast variables to the leading order by the OU process,polynomial form approximation of the interaction terms in the equation for momentum). Thus, weconsider BRT-SMR models where the stochastic part in the SMR parameterization is reduced by40% or completely neglected. Stochastic terms in the SMR parametrization represent the energybackscatter of small scales onto resolved large scales. It has been recognized that proper model-ing of this phenomena is particularly important in the context of geophysical turbulence (see e.g.[42, 43, 5]). Therefore, we study how well the SMR parametrization reproduces this process andwhether various approximations introduced in the context of applying the SMR to the shallowwater equation impose additional sensitivity of the BRT-SMR model to the stochasticity of theclosure.The reduction of the stochastic part by 40% (defining BRT-SMR-0.6) significantly reduces theerror in the variance of H to 2 .
9% and of HU to − . H and HU can be reduced by a different percentage, thus further optimizing theperformance of the BRT-SMR model. With the choice of the 40% reduction of SGS noise, theperformance in the energy spectrum can be improved for the first 8 wave numbers, but for higherwave numbers the energy content drops, see Fig. 4. However, we wold like to emphasize that thefirst 8 wavenumbers contain approximately 97% of the potential energy. From Fig. 4 it is alsovisible that already the deterministic part of the closure can significantly improve the spectrum ascompared to BRT. Correlation.
The time autocorrelations from low-resolution simulations (BRT or LRM) witheither the empirical OU or the SMR parameterization are depicted in Fig. 3 (left). One can seethat application of the SMR parameterizations leads reproducing the autocorrelation from the DNSwith small differences in amplitude. The relative error of the correlation is 3 . .
5% for the LRM-OU simulation and6 .
6% for the BRT-OU simulation. The SGS noise reduction in the BRT-SMR-0.6 model does notsignificantly affect the correlation function (not depicted for this model). The correlation functionfor the BRT-SMR-0.6 overlaps with the correlation function computed using the BRT-SMR. Thiscan be intuitively understood since the correlation function in many stochastic models is determinedprimarily by the strength of the deterministic terms.
The advantage of the SMR parameterization is that it can be adapted easily to changes in the modelsetup, and in many situations it does not have to be recalculated. This has been investigated byconsidering larger averaging intervals of n = 16 ,
32 , resulting in different spatial resolutions N c = N/n = 16 ,
32. To adjust the SMR closure to the changed resolution, we use (74) and (76)from the Appendix D. Note that no re-determination of the model is necessary. In contrast to this,13o modification rule exists for the empirical closures in the BRT-OU and the LRM-OU model. Wekeep those parameterizations unchanged for the considered cases. Whereas the LRM-OU remainsstable, the BRT-OU is unstable in both cases.The potential energy spectra from integrations of the resulting stable models are displayed inFig. 5. For comparison the corresponding DNS projection is shown as well. In both cases integrationof the low-resolution models with SGS parameterization yield less energy in the resolved flow thanthe DNS. However, application of the SMR SGS parameterization leads to better agreement withthe DNS, especially for n = 16. Both low-resolution simulations can capture the time correlationwell (not shown). The applicability of subgrid-scale (SGS) parameterizations to a wide range of parameters of adynamical system such as the atmosphere, and their ability to be easily used at different modelconfigurations, requires that they are based on first principles as much as possible. Stochasticmode reduction (SMR) as suggested by [37, 35, 38], i.e. homogenization applied to a systemwith its nonlinear fast-variable self-interactions replaced by an empirical Ornstein-Uhlenbeck (OU)process, is a promising option in this direction. Geophysical applications of the SMR so far wereperformed always in spectral space [16, 15]. However, in many applications, such as ocean modelingor regional climate modeling, SGS parameterizations in physical space are required. In order toconstruct such parameterization we use the local approach suggested by [12, 13], and tested withinthe framework of the Burgers equation. A central aspect of this approach is the discrimination,within a finite-volume formulation of the high-resolution dynamics, between slowly varying spatialaverages, that are resolved explicitly, and more rapidly varying deviations from those, that are tobe parameterized.As a next step towards the application of this technique to real atmospheric flows, our workvalidates the applicability of the SMR in the context of one-dimensional shallow-water (1DSW)flow. This introduces two general features. (1) Gravity waves are included as well as (2) high-ordernon-polynomial nonlinearities that generally affect compressible flows. After the validation of therequired time-scale separation between local averages and small-scale flow, the latter issue hasbeen handled by replacing, in all dynamical terms affecting or affected by the fast small-scale flow,the inverse of the water-column height by the inverse of its global equilibrium value. This limitsthe corresponding nonlinearities to quadratic. Further replacing all small-scale self-interactions byan empirical OU process yields a representation of the dynamics that allows model simulationsin rather good agreement with simulations of the unmodified 1DSW equations. We could henceproceed and apply the homogenization technique to obtain an explicit low-resolution model for thelocal averages, with an SMR SGS parameterization of the small-scale flow coupling only a smallnumber of neighboring cells.This model has been validated against data from high-resolution simulations of 1DSW flow.It is shown that the SMR SGS parameterization improves the energy spectrum at the smallerresolved scales in comparison with both simulations without SGS parameterizations and simulationsusing an empirical OU SGS parameterization. In the error of some statistical moments no clearbenefit of the SMR SGS parameterization is present. However, we demonstrated that the errorcan be considerably lowered by diminishing the stochasticity in the SMR closure. In particular,the variance error of the SMR SGS model can be reduced to 2 .
9% in H . We also found that thiscomes along with less energy at high wavenumbers. We conjecture that the performance of theBRT-SMR model can be improved further by empirically adjusting the coefficients of the SMR14GS parametrization. Finally, we also show that the closure can easily be adapted to changesin the model parameters. This enables a scale awareness, which allows to utilize the SMR SGSparameterization for different spatial resolutions and leads to improvements compared to empiricalSGS schemes.In a related study within the framework of the Lorenz 96 model, [59] have recently demonstratedparameter-awareness of a parameterization derived using response theory [63] by changing thetime-scale separation between resolved and unresolved scales. How far this extends to our setting,where scale-adaptivity is considered with regard to the number of resolved modes, remains to beinvestigated. Moreover, as shown by [62] for comparatively simple models, the SGS scheme of[63] does outperform the SMR parameterization at smaller time lags, but on longer time scalesit converges to the SMR result in the limit of infinite time scale separation. Still, it would beinteresting to extend the work of [10] by comparing both approaches in more complex applications.Motivated by the DIA closure of [17, 19] applied a stochastic modeling approach accountingfor memory effects of the turbulent eddies and constructed SGS parameterization from a high-resolution simulation. The resulting SGS model is local in spectral space and includes lineareddy drain viscosity and stochastic backscatter viscosity. The same approach was successfullyused by [28, 29] to construct scale-aware SGS parameterizations, which reproduce the spectraexactly. Interestingly, the SMR provides additional nonlinear deterministic correction terms andmultiplicative noise terms. Such terms might become important in situations where effects due totopography, intermittency or large-scale flow are relevant. In addition, the scaling laws found by[28, 29] for the eddy viscosities suggest that similar scaling laws might be valid for the parametersof the OU process in Fourier space, used in the SMR approach. This might improve further theresults on scale-adaptivity presented here.Potentially an issue is that we had to increase diffusivity in order to stabilize the low-resolutionmodel with SMR SGS parameterization. This is a well known issue with purely empirical SGSparameterizations (e,g, [1]). In the present semi-analytical approach it might be overcome by usingthe energy conserving discretization of [14]. As pointed out by [36] there is a connection betweenenergy conservation by the discretized nonlinearities and the cubic damping term in the SMR SGSparameterization. Indeed, in the studies of [16, 15, 12] the discrete treatment of the nonlinear termsconserves energy and the resulting SMR SGS schemes are stable.Even from the present results, however, we conclude that it appears worthwhile further movingtowards the application of the SMR SGS parameterizations to low-resolution simulations in generalcompressible flows. Next step would be to increase the complexity by considering two-dimensionalshallow-water flow and by including rotational effects. Such system contains dispersive inertialgravity waves as well as geostrophic balanced flow. One interesting question in this regard is if theeffect of high-frequency, small-scale gravity waves on the large-scale gravity waves and geostrophicflow can be parameterized using the present local SMR approach. Acknowledgments
We want to thank the reviewers for their comments and suggestions which helped to improve thedraft version of the manuscript. The code for the SGS parameterizations is available upon request.MZ, SD and IT thank the German Research Foundation (DFG) for partial support through grantDO 1819/1-1. IT thanks for partial support by the grant ONR N00014-17-1-2845. UA thanks DFGfor partial support through grant AC 71/7-1. 15
Interaction coefficients
To define the interaction coefficients in (11), (12), the following notation is used( L xzlm z m , B xxzlmn x m z n , B xzzlmn z m z n ) = (cid:40)(cid:0) L Hz , B Hxz , B
Hzz (cid:1) if x l denotes H I (cid:0) L HUz , B
HUxz , B
HUzz (cid:1) if x l denotes HU I , (49)( L zzlm z m , L zxlm x m ) = (cid:16) L h (cid:48) z , L h (cid:48) x (cid:17) if z l denotes h (cid:48) i (cid:16) L hu (cid:48) z , L hu (cid:48) x (cid:17) if z l denotes hu (cid:48) i , (50)( B zxxlmn x m x n , B zxzlmn x m z n , B zzzlmn z m z n ) = (cid:16) B h (cid:48) xx , B h (cid:48) xz , B h (cid:48) zz (cid:17) if z l denotes h (cid:48) i (cid:16) B hu (cid:48) xx , B hu (cid:48) xz , B hu (cid:48) zz (cid:17) if z l denotes hu (cid:48) i , (51)with I ∈ { , , ..., N c − } and i ∈ { , , ..., N − } . The linear interaction coefficients read L Hz = − n ∆ x (cid:16) hu (cid:48) n ( I +1) + hu (cid:48) n ( I +1) − − hu (cid:48) nI − hu (cid:48) nI − (cid:17) (52)+ νn ∆ x (cid:16) h (cid:48) n ( I +1) − h (cid:48) n ( I +1) − − h (cid:48) nI + h (cid:48) nI − (cid:17) L HUz = νn ∆ x (cid:16) hu (cid:48) n ( I +1) − hu (cid:48) n ( I +1) − − hu (cid:48) nI + hu (cid:48) nI − (cid:17) (53) L h (cid:48) x = − x (cid:0) HU I [ i +1] − HU I [ i − (cid:1) + ν ∆ x (cid:0) H I [ i − − H I [ i ] + H I [ i +1] (cid:1) (54)+ 12 n ∆ x (cid:0) HU I [ i ]+1 − HU I [ i ] − (cid:1) − νn ∆ x (cid:0) H I [ i ]+1 − H I [ i ] + H I [ i ] − (cid:1) (55) L hu (cid:48) x = ν ∆ x (cid:0) HU I [ i − − HU I [ i ] + HU I [ i +1] (cid:1) (56) − νn ∆ x (cid:0) HU I [ i ]+1 − HU I [ i ] + HU I [ i ] − (cid:1) (57) L h (cid:48) z = − x (cid:0) hu (cid:48) i +1 − hu (cid:48) i − (cid:1) + ν ∆ x (cid:0) h (cid:48) i − − h (cid:48) i + h (cid:48) i +1 (cid:1) − L Hz (58) L hu (cid:48) z = ν ∆ x (cid:0) hu (cid:48) i − − hu (cid:48) i + hu (cid:48) i +1 (cid:1) − L HUz , (59)where in (58), (59) the terms L Hz and L HUz are given in (52), (53), but with the index I replaced by I [ i ]. The nonlinear interaction coefficients read16 Hxz = B Hzz = B h (cid:48) xx = B h (cid:48) xz = B h (cid:48) zz = 0 (60) B HUxz = − n ∆ x (cid:34) H (cid:16) HU I +1 hu (cid:48) n ( I +1) + HU I hu (cid:48) n ( I +1) − (61) − HU I hu (cid:48) nI − HU I − hu (cid:48) nI − (cid:17) + g (cid:16) H I +1 h (cid:48) n ( I +1) + H I h (cid:48) n ( I +1) − − H I h (cid:48) nI − H I − h (cid:48) nI − (cid:17) (cid:35) B HUzz = − n ∆ x (cid:34) H (cid:16) ( hu (cid:48) n ( I +1) ) + ( hu (cid:48) n ( I +1) − ) − ( hu (cid:48) nI ) − ( hu (cid:48) nI − ) (cid:17) (62)+ g (cid:16) ( h (cid:48) n ( I +1) ) + ( h (cid:48) n ( I +1) − ) − ( h (cid:48) nI ) − ( h (cid:48) nI − ) (cid:17) (cid:35) B hu (cid:48) xx = − x (cid:18) H (cid:16) HU I [ i +1] − HU I [ i − (cid:17) + g (cid:16) H I [ i +1] − H I [ i − (cid:17)(cid:19) (63)+ 12 n ∆ x (cid:20) H (cid:16) HU I [ i ]+1 − HU I [ i ] − (cid:17) + g (cid:16) H I [ i ]+1 − H I [ i ] − (cid:17)(cid:21) (64) B hu (cid:48) xz = − x (cid:20) H (cid:0) HU I [ i +1] hu (cid:48) i +1 − HU I [ i − hu (cid:48) i − (cid:1) + g (cid:0) H I [ i +1] h (cid:48) i +1 − H I [ i − h (cid:48) i − (cid:1)(cid:21) − B HUxz (65) B hu (cid:48) zz = − x (cid:18) H (cid:0) ( hu (cid:48) i +1 ) − ( hu (cid:48) i − ) (cid:1) + g (cid:0) ( h (cid:48) i +1 ) − ( hu (cid:48) i − ) (cid:1)(cid:19) − B HUzz , (66)where in (65), (66) the terms B HUxz and B HUzz are given in (61), (62), but with the index I replaced by I [ i ]. B Null-space projection and inverse of the OU backward Fokker-Planck operator
To determine the projection operator P and the inverse operator L − one can consider an auxiliaryprocess described by the backward FPE ∂ t χ = L χ with the conditional PDF χ (˜ y , | y , τ ). Theinvariant measure of this process p s (˜ y ) = lim τ →−∞ χ (˜ y , | y , τ ) defines the projection operator( P g ) ( x ) = (cid:90) d y g ( x , y ) p s ( y ) = (cid:104) g ( x , y ) (cid:105) , (67)with the expectation (cid:104)·(cid:105) of g ( x , y ) with respect to the invariant measure p s . The inverse operator L − applied onto a function f ( x , y ) is found to be given by (cid:0) L − f (cid:1) ( x , y ) = (cid:90) ∞ dτ (cid:90) d y f ( x , ˜ y ) χ (˜ y , τ | y , . (68)Both operators applied consecutively yield 17 P gL − f (cid:1) ( x ) = (cid:90) d y g ( x , y ) p s ( y ) (cid:90) ∞ dτ (cid:90) d ˜ y f ( x , ˜ y ) χ (˜ y , τ | y , (cid:90) ∞ dτ (cid:90) d y g ( x , y ) p s ( y ) (cid:90) d ˜ y f ( x , ˜ y ) χ (˜ y , τ | y , (cid:90) ∞ dτ (cid:104) g ( x , y ) f ( x , ˜ y ( τ )) (cid:105) , (69)where in the lagged covariance in the last line ˜ y ( τ ) is understood to be an OU trajectory withinitial condition ˜ y (0) = y . C Solvability condition
The solvability condition (34) can be rewritten to
P b xi ∂ x i p = 0 since p = p ( x ). It is fulfilled ifthe even stronger condition P b xi = 0 holds. This is the case here. Using (23), defining κ jk = R jl R km (cid:104) y l y m (cid:105) (70)and inserting the transformation z = R y , one obtains with (cid:104) y i (cid:105) = 0 P b xi = B xyyijk (cid:104) y j y k (cid:105) = B xzzijk κ jk = B xzzijj κ jj . (71)We have used in the last step that B xzzijk = B xzzijk δ jk , as can be seen from Appendix A. Since B xzzijj κ jj is a difference between fluxes at the right and at the left of a cell, the total sum over i vanishes (cid:88) i B xzzijj κ jj =0 . (72)The homogeneity of (cid:104) y l y m (cid:105) implies that each element in the sum is identical and thus must vanish.This proofs the solvability condition. D The SMR SGS parameterization for changed resolution
The SMR SGS parameterization can be adapted to different coarse grid resolutions without anyrecalculation. This can be achieved by collecting interaction coefficients proportional to differentpowers of n . For example the linear coefficients L yx can be written as L yx = ˜ L yx + n ˆ L yx , whereˆ L yx results from all terms in L yx multiplied by n and ˜ L yx from terms independent of n . Similarlyall other coefficients can be split in this way, in particular we have L xy = n ˆ L xy , B xxy = n ˆ B xxy ,etc.. This separation leads to the deterministic closure β i = (cid:18) n ˆ L xymj + 1 n ˆ B xxymlj x l (cid:19) n ˆ B xxyimk ( C S ) jk + (cid:18)(cid:20) ˜ L yxmo + 1 n ˆ L yxmo ) (cid:21) x o + (cid:20) ˜ B yxxmop + 1 n ˆ B yxxmop (cid:21) x o x p (cid:19) n (cid:16) ˆ L xyik + ˆ B xxyijk x j (cid:17) (cid:104) yy T (cid:105) − mn ( C S ) nk + 1 n ˆ B xyyijk (cid:20) ˜ B yxympo + 1 n ˆ B yxympo (cid:21) x p (cid:104) yy T (cid:105) − mn ( C T ) onjk (73)= 1 n ˜ β i + 1 n ˆ β i (74)18here in the last steps the results are summarized with respect to the power of n . The effectivestochastic closure analogously yields dξ i ≈ n (cid:113) B xyyijn ( C T ) jnkl ˆ B xyyikl dW i + 1 n (cid:114) (cid:16) ˆ L xyin + ˆ B xxyijn x j (cid:17) ( C S ) nk (cid:16) ˆ L xyik + ˆ B xxyilk x l (cid:17) dW i (75)= 1 n d ˆ ξ i . (76) E Spatial shape of the SGS variance
In continuous space the surface-height mean in the first coarse cell, with length L c , and the corre-sponding SGS deviations can be written as H = 1 L c (cid:90) L c h ( x ) dx , h (cid:48) ( x ) = h ( x ) − H . (77)This leads to spatial dependence in the variance of h (cid:48) h (cid:48) ( x ) − h (cid:48) = h ( x ) − h ( x ) H + H − h (cid:48) , (78)Due to spatial homogeneity of h ( x ) and h (cid:48) = 0, only the second term can be assumed to bespatially dependent. Now by assuming an exponentially decaying spatial correlation h ( x ) h ( x (cid:48) ) ∼ exp ( − α | x − x (cid:48) | ), with a decay rate α >
0, the middle term becomes h ( x ) H = 1 L c (cid:90) L c h ( x ) h ( x (cid:48) ) dx (cid:48) ∼ L c (cid:26)(cid:90) x exp (cid:2) α ( x (cid:48) − x ) (cid:3) dx (cid:48) + (cid:90) L c x exp (cid:2) α ( x − x (cid:48) ) (cid:3) dx (cid:48) (cid:27) = 1 αL c (cid:110) − e − αx − e α ( x − L c ) (cid:111) = 2 αL c (cid:26) − e − αLc cosh (cid:20) α (cid:18) x − L c (cid:19)(cid:21)(cid:27) , (79)describing the characteristic U-shape of the fine variable variance in the first coarse cell x ∈ [0 , L c ]. The same considerations hold for all other coarse cells. References [1] U Achatz and G Schmitz. On the Closure Problem in the Reduction of Complex AtmosphericModels by PIPs and EOFs: A Comparison for the Case of a Two-Layer Model with ZonallySymmetric Forcing.
Journal of the Atmospheric Sciences , 54(20):2452–2474, 1997.[2] Ulrich Achatz and Grant Branstator. A Two-Layer Model with Empirical Linear Correctionsand Reduced Order for Studies of Internal Climate Variability.
Journal of the AtmosphericSciences , 56(17):3140–3160, 1999. 19able 1: First column in both tables shows the spatially averaged second ( k = 2) and fourth ( k = 4)order centered statistical moments values of DNS. The other columns contain the relative errors ofdifferent low-resolution simulations. Upper table for H and lower table for HU . H DNS LRM BRT LRM-OU BRT-OU BRT-SMR BRT-SMR-0.6 BRT-SMR-0.0 k = 2 2.846 (km) k = 4 23.56 (km) HU DNS LRM BRT LRM-OU BRT-OU BRT-SMR BRT-SMR-0.6 BRT-SMR-0.0 k = 2 2082 (10 (km) / d) k = 4 1.277E+07 (10 (km) / d) N cells˙ x i = (cid:37) xi + a xi ( x ) + b xzi ( x , z ); ˙ z i = b zi ( x , z ) + c zi ( z )BRT Bare-Truncation Model with spatial resolution of N c = N/n ˙ x i = (cid:37) xi + a xi ( x )LRM Low-Resolution ModelDNS with spatial resolution of N c = N/n .OU-DNS DNS with h -approximation and SGS self-interactions replaced by an OU process˙ x i = (cid:37) xi + a xi ( x ) + b xi ( x , y ), ˙ y i = b yi ( x , y ) + Λ ij y j + Σ i ˙ W i .BRT-OU BRT with empirical Ornstein-Uhlenbeck parameterization dx i = (cid:16) (cid:37) xi + a xi ( x ) + ˜Γ ij ˆ x Ij (cid:17) dt + ˜ σ i dW i .LRM-OU LRM with empirical Ornstein-Uhlenbeck parameterizationBRT-SMR BRT with Stochastic Mode Reduction parameterization dx i = [ (cid:37) xi + a xi ( x ) + β i ( x )] dt + dξ i ( x )BRT-SMR-0.6 BRT-SMR with stochastic forcing reduced to 60% dx i = [ (cid:37) xi + a xi ( x ) + β i ( x )] dt + 0 . dξ i ( x ).BRT-SMR-0.0 BRT-SMR without stochastic forcing dx i = [ (cid:37) xi + a xi ( x ) + β i ( x )] dt [3] Ulrich Achatz, Ulrike L¨obl, Stamen I. Dolaptchiev, and Andrey Gritsun. FluctuationDissipa-tion Supplemented by Nonlinearity: A Climate-Dependent Subgrid-Scale Parameterization inLow-Order Climate Models. Journal of the Atmospheric Sciences , 70(6):1833–1846, jun 2013.[4] Ludwig Arnold, Peter Imkeller, and Yonghui Wu.
Reduction of deterministic coupledatmosphere-ocean models to stochastic ocean models: A numerical case study of the Lorenz-Maas system , volume 18. 2003.[5] J. Berner, G. J. Shutts, M. Leutbecher, and T. N. Palmer. A Spectral Stochastic KineticEnergy Backscatter Scheme and Its Impact on Flow-Dependent Predictability in the ECMWFEnsemble Prediction System.
Journal of the Atmospheric Sciences , 66(3):603–626, 2009.20 τ τ DNSLRMBRT 0 100 200 300 400 500 60000.020.040.060.080.1 Variance h’ [(km) ]Distance [km] DNSOU−DNS
Figure 1: Left: spatially averaged time autocorrelation of the resolved variable H and SGS variable h (cid:48) . Results from the DNS and two low-resolution simulations without SGS parameterizations (LRMand BRT). On the time axis the characteristic gravity wave time τ is marked, see Sec. 4 4.1. Inthe upper right corner a shorter time interval is presented, which resolves the decay of h (cid:48) . Right:the variance of the SGS variable h (cid:48) for 32 fine grid points and an averaging interval n = 8. Resultsfrom the DNS and OU-DNS.[6] Judith Berner, Ulrich Achatz, Lauriane Batt´e, Lisa Bengtsson, Alvaro De La C´amara,Hannah M. Christensen, Matteo Colangeli, Danielle R.B. Coleman, Daaaan Crommelin,Stamen I. Dolaptchiev, Christian L.E. Franzke, Petra Friederichs, Peter Imkeller, HeikkiJ¨arvinen, Stephan Juricke, Vassili Kitsios, Fran¸cois Lott, Valerio Lucarini, Salil Mahajaa-jaajan, Timothy N. Palmer, C´ecile Penland, Mirjajana Sakradzijaja, Jin Song Von Storch,Antje Weisheimer, Michael Weniger, Paul D. Williams, and Jun Ichi Yano. Stochastic pa-rameterization toward a new view of weather and climate models. Bulletin of the AmericanMeteorological Society , 98(3):565–587, 2017.[7] R Buizza, M Miller, and TN Palmer. Stochastic representation of model uncertainties in theECMWF Ensemble Prediction System.
Quarterly Journal of the Royal Meteorological Society ,125(560):2887–2908, oct 1999.[8] Alexei Chekhlov and Victor Yakhot. Kolmogorov turbulence in a random-force-driven Burgersequation: Anomalous scaling and probability density functions.
Physical Review E , 52(5):5681–5684, 1995.[9] Timothy Delsole. Stochastic models of quasigeostrophic turbulence.
Surveys in Geophysics ,25(2):107–149, 2004.[10] Jonathan Demaeyer and S Vannitsem. Comparison of stochastic parameterizations in theframework of a coupled ocean-atmosphere model. 01 2018.[11] Jonathan Demaeyer and St´ephane Vannitsem. Stochastic parametrization of subgrid-scaleprocesses in coupled ocean-atmosphere systems: benefits and limitations of response theory.
Quarterly Journal of the Royal Meteorological Society , 143(703):881–896, jan 2017.21 −10 −8 −4 WavenumberPotential Energy Spectrum [(km) /d ] k −2 DNS 1 2 4 8 16 3210 −9 −8 −4 WavenumberPotential Energy Spectrum [(km) /day ] DNSOU−DNSLRMBRT
Figure 2: Left: the potential energy spectrum computed from h in the DNS. Right: the poten-tial energy spectrum computed from H in the DNS, OU-DNS and two low-resolution simulationswithout SGS parameterizations (LRM and BRT). τ τ DNSLRM−OUBRT−OUBRT−SMR 1 2 4 8 16 3210 WavenumberPotential Energy Spectrum [(km) /day ] DNSLRM−OUBRT−OUBRT−SMR
Figure 3: The spatially averaged time autocorrelation (left) and the potential energy spectrum(right) from DNS and three low-resolution simulations with SGS parameterizations (LRM-OU,BRT-OU, BRT-SMR).[12] S. I. Dolaptchiev, U. Achatz, and I. Timofeyev. Stochastic closure for local averages in thefinite-difference discretization of the forced Burgers equation.
Theoretical and ComputationalFluid Dynamics , 27(3-4):297–317, 2013.[13] Stamen I. Dolaptchiev, Ilya Timofeyev, and Ulrich Achatz. Subgrid-scale closure for theinviscid Burgers-Hopf equation.
Communications in Mathematical Sciences , 11(3):757–777,2013.[14] Ulrik S. Fjordholm, Siddhartha Mishra, and Eitan Tadmor. Well-balanced and energy stableschemes for the shallow water equations with discontinuous topography.
Journal of Computa-tional Physics , 230(14):5587–5609, 2011. 22 WavenumberPotential Energy Spectrum [(km) /day ] DNSBRTBRT−SMRBRT−SMR−0.6BRT−SMR−0.0
Figure 4: The potential energy spectrum in DNS, BRT, BRT-SMR, BRT-SMR with a dampedstochastic forcing dξ → . dξ (BRT-SMR-0.6) and BRT-SMR with neglected stochastic forcing dξ → WavenumberPotential Energy Spectrum n=16 [(km) /day ] DNSLRM−OUBRT−SMR 1 2 4 810 WavenumberPotential Energy Spectrum n=32 [(km) /day ] DNSLRM−OUBRT−SMR
Figure 5: The Spectrum for an averaging interval of n = 16 (left) and n = 32 (right) in DNS, LRM-OU and BRT-SMR. The simulations with LRM-OU and BRT-SMR have a resolution of N c = 32(left) and N c = 16 (right).[15] Christian Franzke and Andrew J Majda. Low-Order Stochastic Mode Reduction for a Proto-type Atmospheric GCM. Journal of the Atmospheric Sciences , 63(2):457–479, 2006.[16] Christian Franzke, Andrew J Majda, and Eric Vanden-Eijnden. Low-Order Stochastic ModeReduction for a Realistic Barotropic Model Climate.
Journal of the Atmospheric Sciences ,62(6):1722–1745, jun 2005.[17] J. S. Frederiksen and A. G. Davies. Eddy Viscosity and Stochastic Backscatter Parameteri-zations on the Sphere for Atmospheric Circulation Models.
Journal of Atmospheric Sciences ,54:2475–2492, October 1997. 2318] J. S. Frederiksen, M. R. Dix, and A. G. Davies. The effects of closure-based eddy diffusion onthe climate and spectra of a gcm.
Tellus A , 55(1):31–44, 2003.[19] J. S. Frederiksen and S. M. Kepert. Dynamical Subgrid-Scale Parameterizations from DirectNumerical Simulations.
Journal of Atmospheric Sciences , 63:3006–3019, November 2006.[20] JS Frederiksen. Subgrid-scale parameterizations of eddy-topographic force, eddy viscosity,and stochastic backscatter for flow over topography.
JOURNAL OF THE ATMOSPHERICSCIENCES , 56(11):1481–1494, JUN 1 1999.[21] Crispin Gardiner.
Stochastic Methods . Springer-Verlag Berlin Heidelberg, 4 edition, 2009.[22] Dale B. Haidvogel and Aike Beckmann.
Numerical Ocean Circulation Modeling.
ImperialCollege Press, London, 1999.[23] K. Hasselmann. Stochastic climate models Part I. Theory.
Tellus , 28(6):473–485, jan 1976.[24] J. Honerkamp.
Stochastic Dynamical Systems: Concepts, Numerical Methods, Data Analysis .Wiley-VCH, 1994.[25] Peter Imkeller and Jin-Song von Storch, editors.
Stochastic Climate Models . Birkh¨auser Basel,Basel, 2001.[26] R. Z. Khasminsky. A limit theorem for the solutions of differential equations with randomright-hand sides.
Theory Prob. Applications , 11:390–406, 1966.[27] R. Z. Khasminsky. On stochastic processes defined by differential equations with a smallparameter.
Theory Prob. Applications , 11:211–228, 1966.[28] V. Kitsios, J. S. Frederiksen, and M. J. Zidikheri. Subgrid Model with Scaling Laws forAtmospheric Simulations.
Journal of Atmospheric Sciences , 69:1427–1445, April 2012.[29] V. Kitsios, J. S. Frederiksen, and M. J. Zidikheri. Scaling laws for parameterisations of subgrideddy-eddy interactions in simulations of oceanic circulations.
Ocean Modelling , 68:88–105,August 2013.[30] Dmitri Kondrashov, S. Kravtsov, A. W. Robertson, and M. Ghil. A hierarchy of data-basedENSO models.
Journal of Climate , 18(21):4425–4444, 2005.[31] Robert H. Kraichnan. The structure of isotropic turbulence at very high reynolds numbers.
Journal of Fluid Mechanics , 5(4):497–543, 1959.[32] S. Kravtsov, D. Kondrashov, and M. Ghil. Multilevel Regression Modeling of Nonlinear Pro-cesses: Derivation and Applications to Climatic Variability.
Journal of Climate , 18(21):4404–4424, nov 2005.[33] T. G. Kurtz. A limit theorem for perturbed operator semigroups with applications to randomevolution.
J. Funct. Anal. , 12:55–67, 1973.[34] Nicholas J. Lutsko, Isaac M. Held, and Pablo Zurita-Gotor. Applying the FluctuationDis-sipation Theorem to a Two-Layer Model of Quasigeostrophic Turbulence.
Journal of theAtmospheric Sciences , 72(8):3161–3177, 2015.2435] A. Majda, I. Timofeyev, and E. Vanden-Eijnden. A priori tests of a stochastic mode reductionstrategy.
Physica D: Nonlinear Phenomena , 170(3-4):206–252, 2002.[36] Andrew J Majda, Christian Franzke, and Daan Crommelin. Normal forms for reduced stochas-tic climate models.
Proceedings of the National Academy of Sciences of the United States ofAmerica , 106(10):3649–53, 2009.[37] Andrew J. Majda, Ilya Timofeyev, and Eric Vanden Eijnden. A mathematical framework forstochastic climate models.
Communications on Pure and Applied Mathematics , 54(8):891–974,2001.[38] Andrew J. Majda, Ilya Timofeyev, and Eric Vanden-Eijnden. Systematic Strategies forStochastic Mode Reduction in Climate.
Journal of the Atmospheric Sciences , 60(14):1705–1722, 2003.[39] Detlev Majewski, D¨orte Liermann, Peter Prohl, Bodo Ritter, Michael Buchhold, ThomasHanisch, Gerhard Paul, Werner Wergen, and John Baumgardner. The Operational GlobalIcosahedralHexagonal Gridpoint Model GME: Description and High-Resolution Tests.
MonthlyWeather Review , 130(2):319–338, 2002.[40] Adam H. Monahan and Joel Culina. Stochastic averaging of idealized climate models.
Journalof Climate , 24(12):3068–3088, 2011.[41] Matthew Newman, Prashant D. Sardeshmukh, Christopher R. Winkler, and Jeffrey S.Whitaker. A Study of Subseasonal Predictability.
Monthly Weather Review , 131(8):1715–1732, 2003.[42] T. N. Palmer. A nonlinear dynamical perspective on model error: A proposal for non-localstochastic-dynamic parametrization in weather and climate prediction models.
Quarterly Jour-nal of the Royal Meteorological Society , 127(572):279–304, 2001.[43] T N Palmer, R Buizza, F Doblas-Reyes, T Jung, M Leutbecher, G J Shutts, M Steinheimer,and A Weisheimer. Stochastic Parametrization and Model Uncertainty.
ECMWF TechnicalMemoranda , 598:1–42, 2009.[44] G. Papanicolaou. Some probabilistic problems and methods in singular perturbations.
RockyMountain J. Math , 6:653–673, 1976.[45] Grigorios A. Pavliotis and Andrew M. Stuart.
Multiscale Methods , volume 53 of
Texts Appliedin Mathematics . Springer New York, New York, NY, 2008.[46] Kathy Pegion and Prashant D. Sardeshmukh. Prospects for Improving Subseasonal Predic-tions.
Monthly Weather Review , 139(11):3648–3666, 2011.[47] V. Petoukhov, A. Ganopolski, V. Brovkin, M. Claussen, A. Eliseev, C. Kubatzki, and S. Rahm-storf. CLIMBER-2: A Climate System Model of Intermediate Complexity. Part I: ModelDescription and Performance for Present Climate.
Climate Dynamics , 16:1–17, 2000.[48] Martin Pieroth, Stamen I. Dolaptchiev, Matthias Zacharuk, Andrey Gritsun, and UlrichAchatz. Climate-Dependence in Empirical Parameters of Subgrid-Scale Parameterizationsusing the Fluctuation-Dissipation Theorem.
Journal of the Atmospheric Sciences , 2018. sub-mitted. 2549] Stephen B. Pope.
Turbulent Flows . Cambridge University Press, Cambridge, 2000.[50] P. R´ıpodas, A. Gassmann, J. F¨orstner, D. Majewski, M. Giorgetta, P. Korn, L. Korn-blueh, H. Wan, G. Z¨angl, L. Bonaventura, and T. Heinze. Icosahedral Shallow Water Model(ICOSWM): results of shallow water test cases and sensitivity to model parameters.
Geosci-entific Model Development Discussions , 2(1):581–638, 2009.[51] M. Satoh, T. Matsuno, H. Tomita, H. Miura, T. Nasuno, and S. Iga. Nonhydrostatic icosahe-dral atmospheric model (NICAM) for global cloud resolving simulations.
Journal of Compu-tational Physics , 227(7):3486–3514, 2008.[52] Glenn Shutts. A kinetic energy backscatter algorithm for use in ensemble prediction systems.
Quarterly Journal of the Royal Meteorological Society , 131(612):3079–3102, 2005.[53] W.C. Skamarock and Coauthors. A Description of the Advanced Research WRF Version 3.Technical Report NCAR/TN-475+STR, NCAR, 2008.[54] Roland B. Stull, editor.
An Introduction to Boundary Layer Meteorology . Springer Netherlands,Dordrecht, 1988.[55] Geoffrey K Vallis.
Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-ScaleCircu- lation . 2006.[56] W. T M Verkley. A maximum entropy approach to the problem of parametrization.
QuarterlyJournal of the Royal Meteorological Society , 137(660):1872–1886, 2011.[57] W. T M Verkley, P. C. Kalverla, and C. A. Severijns. A maximum entropy approach to theparametrization of subgrid processes in two-dimensional flow.
Quarterly Journal of the RoyalMeteorological Society , 142(699):2273–2283, 2016.[58] W. T M Verkley and C. A. Severijns. The maximum entropy principle applied to a dynamicalsystem proposed by Lorenz.
European Physical Journal B , 87(1), 2014.[59] Gabriele Vissio and Valerio Lucarini. A proof of concept for scale-adaptive parameterizations:the case of the Lorenz ’96 model.
Quarterly Journal of the Royal Meteorological Society , 2017.[60] Andrew J. Weaver, Michael Eby, Edward C. Wiebe, Cecilia M. Bitz, Phil B. Duffy, Tracy L.Ewen, Augustus F. Fanning, Marika M. Holland, Amy MacFadyen, H. Damon Matthews, Ka-trin J. Meissner, Oleg Saenko, Andreas Schmittner, Huaxiao Wang, and Masakazu Yoshimori.The UVic earth system climate model: Model description, climatology, and applications topast, present and future climates.
Atmosphere-Ocean , 39(4):361–428, 2001.[61] Christopher R. Winkler, Matthew Newman, and Prashant D. Sardeshmukh. A linear model ofwintertime low-frequency variability. Part I: Formulation and forecast skill.
Journal of Climate ,14(24):4474–4494, 2001.[62] Jeroen Wouters, Stamen Iankov Dolaptchiev, Valerio Lucarini, and Ulrich Achatz. Parame-terization of stochastic multiscale triads.
Nonlinear Processes in Geophysics , 23(6):435–445,2016.[63] Jeroen Wouters and Valerio Lucarini. Disentangling multi-level systems: averaging, correla-tions and memory.
Journal of Statistical Mechanics: Theory and Experiment , 2012(03):P03003,2012. 2664] Jeroen Wouters and Valerio Lucarini. Multi-level Dynamical Systems: Connecting the RuelleResponse Theory and the Mori-Zwanzig Approach.