Modelling and numerical approximation of a 2.5D set of equations for mesoscale atmospheric processes

# Modelling and numerical approximation of a 2.5D set of equations for mesoscale atmospheric processes

Dante Kalise Dipartimento di Matematica, Università di Roma “La Sapienza”, P. Aldo Moro 2, 00185 Roma, Italy. Ivar Lie Research and Development Department, StormGeo AS, Universitetsgata 8, 0164 Oslo, Norway
###### Abstract

The set of 3D inviscid primitive equations for the atmosphere is dimensionally reduced by a Discontinuous Galerkin discretization in one horizontal direction. The resulting model is a 2D system of balance laws where with a source term depending on the layering procedure and the choice of coupling fluxes, which is established in terms of upwind considerations. The “2.5D” system is discretized via a WENO-TVD scheme based in a flux limiter centered approach. We study four tests cases related to atmospheric phenomena to analyze the physical validity of the model.

###### keywords:
primitive equations, layering, discontinuous Galerkin, upwind flux, WENO-TVD schemes, test cases for dynamical cores
journal:

## 1 Introduction

The so-called primitive equations used to model atmospheric and ocean dynamics are not well-posed for any reasonable set of boundary conditions. This important, and relatively old result, was obtained by Oliger and Sundström in Oliger and Sundström (1978). The problem of formulating a well-posed set of “primitive” equations is of obvious importance both from the mathematical and numerical side, and has recently attracted substantial research activity. One example of such research was undertaken in Chen et al. (2008), where the authors have formulated a type of dimensionally-reduced equations, called 2.5D equations, where well-posedness was analyzed. However, this paper uses linearized primitive equations and the layering procedure was performed via a Continuous Galerkin approach using orthogonal piecewise linear functions. This imposes some limitations regarding the linear character of the system and the number of elements to be considered in the expansion. What we want to communicate here is an alternative derivation of a 2.5D set of full nonlinear equations and its approximation by a high-order finite volume scheme. We consider the analysis of well-posedness for these equations to be a far too difficult task. Instead we perform numerical experiments with well-known test cases normally used to test dynamical cores of atmospheric models Giraldo and Restelli (2008); Klemp and Skamarock (1994). Good results from such experiments are a strong indication that our formulation of a 2.5D set of equations is well-posed. The dimensional reduction procedure that we propose is based on a Discontinuous Galerkin approach (DG) Cockburn and Shu (2001); this class of schemes have proved to be successful in the resolution 2D atmospheric models Giraldo and Restelli (2008). We use the same basic principles to divide the domain into 2D layers where the equations are locally approximated and coupled by suitable fluxes. This approach leads to a system of 2D balance laws, which is solved via a WENO-TVD scheme Titarev and Toro (2005), using a flux limiter centered approach Billett and Toro (2000) for space discretization and a Runge-Kutta TVD scheme for time discretization. Since the primitive equations are not well-posed, one may wonder what approach is taken in the many atmospheric and oceanic models that exist. The answer is that a pragmatic, and hence not strictly rigorous, approach is taken. In practice this means that one imposes boundary conditions that are not boundary conditions in the ordinary sense. An example is the use of relaxation zones near the boundary, where the solution in the inner of the domain is relaxed towards a certain boundary value. While this is stable numerically (in most cases), it can hardly be said to be satisfactory in the physical sense. And not in the mathematical sense either.

The remaining part of the paper is organized as follows: In section 2 we present the full 3D governing equations in conservative form for a dry atmosphere, and its ”semi-discretization” in a spatial dimension so the equations are formulated in 2.5D as explained above. The finite volume-based numerical scheme is presented in detail in section 3, and in section 4 we present a number of numerical experiments using both well-known tests for atmospheric model dynamical cores, as well as new tests. The results are discussed for each test, and we summarize the results along with directions for further work in section 5.

## 2 A dimensionally-reduced system for atmospheric dynamics

This section is concerned with the derivation of a dimensionally-reduced set of equations describing the atmosphere. We first present the original three-dimensional system of primitive equations which is then layered via a DG ansatz in order to obtain a 2.5D model.

### 2.1 A three-dimensional system of balance laws

Our starting point is the set of equations modelling the evolution in time of a dry atmosphere. We impose conservation of mass, momentum and energy, and consider effects of rotation and gravity, together with neglecting friction; this leads to a set of 3D inviscid primitive equations cast in conservative form

 ∂tQ+∂xF+∂yG+∂zH=S, (1)

where

 Q=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρρuρvρwρθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,F=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρuρu2+Pρuvρuwρuθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,G=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρvρuvρv2+Pρvwρvθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (2)
 H=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρwρwvρwvρw2+Pρwθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,S=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0fρv−fρu−ρg0⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (3)

The dynamics are expressed in terms of the density of the fluid , the 3D velocity field , the pressure and the potential temperature , which relates to the thermodynamic temperature via

 (4)

The system is closed by the equation of state for an ideal gas,

 P=C0(ρθ)γ,C0=RγdPRd/cv0. (5)

The model parameters are: the Coriolis parameter , the gravitational acceleration , the atmospheric pressure at sea level , the gas constant for dry air , the specific heat of dry air at constant pressure and volume , the specific heat of dry air at constant volume and its ratio . Moreover, by defining the Exner pressure

 π≡(PP0)R/cp, (6)

the expression for the total energy of the system is given by,

 e=cvθπ+12(u2+v2+w2)+gz. (7)

The system (1) governs the time evolution of the conserved variables in by the interaction between the physical fluxes and the source term . A more detailed representation of the atmospheric dynamics could include an equation for conservation of moisture that will affect the thermodynamics and hence the dynamics (expressed by the fluxes). Physical parameterizations in a weather model, like surface parameterizations and condensation processes, together with coordinate transformations, like the extensively used terrain-following coordinates will generate terms appearing in the source . In the next subsection we present a dimensionally-reduced version of the original 3D set of equations where we characterize the behavior of the variables in one spatial dimension, in a way such that the action of its associated flux is transferred to the source term. Although this can be seen as a step in the numerical discretization of the original set of 3D equations, our study will be focused on the physical significance of the dimensionally-reduced (or 2.5D) set of equations as a system of balance laws on its own.

### 2.2 Towards a 2.5D model

In Chen et al. (2008), the authors obtained a 2.5D model from the linearized primitive equations via a Continuous Galerkin expansion of the variables in the -direction. It has to be emphasized that the derivation made considerable use of the linearity of the system, and that the basis functions were a set of orthogonal piecewise linear functions; removing any of those aspects leads to a model whose complexity is not clearly related those of the underlying physics. On the other hand, there is an obvious limitation related to the trade-off between orthogonality and locality of the basis functions; the aforementioned article considered two elements with a basis of three orthogonal piecewise functions, which is just a slight modification of the classical hat functions. An increment in the number of elements, together with the orthogonality condition, will produce a set of basis function with increasing support, significantly different from the analogous set of hat functions. The approach that we present pays a cost by removing the continuity requirement on the solution (in the classical Galerkin framework), and switches to a DG formulation that allows, in a straightforward and tractable manner, to deal with nonlinearity and locality at the same time. We start by defining our set of equations over . Initial and boundary conditions will be specified later. We consider a uniform partition in the y-direction into elements, i.e.

 Y=[0,Ly]=Ny⋃i=1Yi,Yi=[yi−1/2,yi+1/2]yi+1/2−yi−1/2=Δy,∀i. (8)

Next, in every element we consider a local classical Galerkin ansatz by multiplying eq. (1) by a test function and integrating with respect to , leading to

 ∫yi+1/2yi−1/2VT{∂tQ+∂xF+∂yG+∂zH−S}dy=0. (9)

Mapping every element into the canonical element by means of

 y=1−ξ2yi−1/2+1+ξ2yi+1/2, (10)

and, as usual, performing integration by parts, leads to

 Δy2∫1−1VT{∂tQ+∂xF+∂zH−S}dξ+VTG|1−1=∫1−1∂ξVTGdξ. (11)

We consider a local set of basis functions and expand restricted to :

 Qi(x,ξ,z,t)=Ni∑j=0qij(x,z,t)ϕj(ξ). (12)

The choice of the set of basis functions is an open question and will depend on the application and the complexity of the computations. In this particular case, in order to be able to carry on in a tractable manner the calculations for a nonlinear model, it is important that we chose a set of orthogonal functions, and therefore selecting a set of orthogonal polynomials (as the Legendre polynomials for instance) will allow us to preserve a reasonable model but also to include high-order approximations. Throughout this article however, we will restrict ourselves to the most basic case in this framework, i.e., to consider Legendre polynomials up to degree 0, which is nothing but to consider a piecewise constant approximation of the variables inside every element. In Kalise et al. (2011), the piecewise linear extension of this procedure was performed in a similar framework, and poses no additional difficulties except for the increase in the computational cost associated with the amount of local expansion coefficients. When we consider 2 elements in the -direction, and and the first Legendre polynomial , the 2.5D system of equations reads,

 ∂tQ1+∂xF(Q1)+∂zH(Q1) =Δy{G(y=0.5Ly)−G(y=0)}+S(Q1), (13) ∂tQ2+∂xF(Q2)+∂zH(Q2) =Δy{G(y=Ly)−G(y=0.5Ly)}+S(Q2). (14)

At this point we must close the system by giving computable expressions for the evaluation of at the boundaries and the intermediate value. We first address the problem of the intermediate value, and then we expand this procedure to the boundary fluxes. Following the same ideas of the finite volume framework, we seek to express the intermediate state as a function of adjacent cells, i.e.

 G(y=0.5Ly)≈G(Q1,Q2). (15)

As it is extensively described in the literature, the two main choices for such procedure are centered (averaged) values, or states obtained via upwind considerations. The typical centered fluxes can be interpreted as particular averaging operators in space and time, which does not fit in a proper manner in our scheme, as time integration will be considered only for the resulting 2.5D model; therefore we opt for upwind fluxes to determine the coupling at the interface. We recall that the flux in the Euler equations has the homogeneity property

 G(Q)=C(Q)Q, (16)

with

 C(Q)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣00100−uvvu00−v202v0c2sγθ−wv0wv0−θv0θ0v⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,c2s=∂ρP=C0(ρθ)γ−1γθ, (17)

where the speed of sound in the fluid. This matrix admits a representation of the form

 C(Q)=RΛR−1, (18)

with

 R=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣11001u010uv−cs√γv00v+cs√γw001wθ000θ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,Λ=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣v−cs√γ00000v00000v00000v00000v+cs√γ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (19)

We decompose into positive and negative parts, by defining

 v+=max(v,0),v−=min(v,0), (20)

and we obtain

 Λ+=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣000000v+00000v+00000v+00000v++cs√γ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,Λ−=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣v−−cs√γ00000v−00000v−00000v−000000⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (21)

Note that in the last expression we have assumed a low Mach number for the fluid (), as expected in atmospheric dynamics. Thus, the intermediate state is defined as

 G(y=0.5Ly)≈G(Q1,Q2)=G+(Q1)+G−(Q2), (22)

where

 G+(Q1) = RΛ+R−1Q1=12(v+cs√γ)⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρρuρ(v+cs√γ)ρwρθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (23) G−(Q2) = RΛ−R−1Q2=12(v−cs√γ)⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ρρuρ(v−cs√γ)ρwρθ⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (24)

The same argument can be applied at the boundaries; the tests performed in section 4 considered solid wall boundary conditions, which were implemented by defining ghost cells values and where,

 ρ0 = ρ1,u0=u1,v0=−v1,w0=w1,θ0=θ1, (25) ρ3 = ρ2,u3=u2,v3=−v2,w3=w2,θ3=θ2. (26)

This final definition allows us to write a direct expression for the 2.5D model:

 ∂tQ1+∂xF(Q1)+∂zH(Q1) = S1(Q1,Q2), (27) ∂tQ2+∂xF(Q2)+∂zH(Q2) = S2(Q1,Q2), (28)
 S1(Q1,Q2) = Δy{G+(Q1)+G−(Q2)−G+(Q0)−G−(Q1)}+S(Q1), (29) S2(Q1,Q2) = Δy{G+(Q2)+G−(Q3)−G+(Q1)−G−(Q2)}+S(Q2). (30)

As we previously pointed out, the resulting dimensionally-reduced model maintains the same structure for the physical fluxes in the and -directions, but has transformed the flux into a “reactive” source term that generates a coupling between the layers. The next section deals with the numerical approximation of this coupled system of balance laws.

## 3 The numerical scheme

In this section we present a finite volume scheme for system of balance laws derived in the previous section. For the sake of simplicity, we develop the scheme for a single set of equations of the form,

 ∂tQ+∂xF(Q)+∂zH(Q)=S(Q), (31)

being eqns. (27)-(28) a particular case with the augmented state . We first indicate that our strategy will be based in a splitting scheme, as it suggested in Chen and Toro (2004) given the flux choice that we will make. Thus, we will first establish a numerical scheme for the system of conservation laws,

 ∂tQ+∂xF(Q)+∂zH(Q)=0, (32)

to be combined with a procedure for the resolution of the source term dynamics

 ∂tQ=S(Q). (33)

In order to approximate eq. (32), we begin by meshing the spatial domain into uniform control volumes of size . Inside every control volume we average with respect to and leading to the semi-discrete scheme

 dQi,j(t)dt=−1Δx(Fi+1/2,j−Fi−1/2,j)−1Δz(Hi,j+1/2−Hi,j−1/2)≡Li,j(Q), (34)

where,

 Qi,j=1Δx1Δz∫xi+1/2xi−1/2∫zj+1/2zj−1/2Q(x,z,t)dzdx, (35)
 Fi+1/2,j =1Δz∫zj+1/2zj−1/2F(Q(xi+1/2,z,t))dz, (36) Hi,j+1/2 =1Δx∫xi+1/2xi−1/2H(Q(x,zj+1/2,t))dx. (37)

We approximate the expressions in eqns. (36)-(37) by suitable Gaussian quadrature formulas,

 Fi+1/2,j≈12∑zGpwzGpF(Q(xi+1/2,zGp,t))dz, (38) Hi,j+1/2≈12∑xGpwxGpH(Q(xGp,zj+1/2,t))dz, (39)

where and are prescribed Gauss points with corresponding weights and respectively. The computation of eqns. (38)-(39) is performed via a high-resolution approach that makes use of a WENO reconstruction procedure; after this step is completed, a polynomial of prescribed order is obtained at every cell, and therefore, at every cell interface, accurate flux calculations can be performed by taking extrapolated boundary values. We briefly describe the WENO reconstruction procedure that is used in this article; we opted for the technique described in Altmann et al. (2007); Balsara et al. (2009) in its third order (quadratic reconstruction) version. This technique makes extensive use of the structure of the reconstruction procedure in one dimension, adding some additional mixed terms (“cross terms”) that are efficiently computed by reduced stencils. It is an optimal and easy way to implement the algorithm for achieving high-order reconstructions in 2 and 3 dimensions; it also defines an unique polynomial in every cell, which is particularly useful when space dependent source terms such as viscosity are considered. At a given time (the subscript indicating time is omitted throughout this derivation), given the set of averaged values for the whole domain, at every cell, the reconstruction procedure seeks a quadratic expansion upon a linear combination of Legendre polynomials rescaled in local coordinates expressed in the form

 Q(x,z)=Q0+QxP1(x)+QxxP2(x)+QzP1(z)+QzzP2(z)+QxzP1(x)P1(z), (40)
 P1(x)=xP2(x)=x2−112. (41)

Except for the last term in (40), every coefficient can be computed by performing a dimension-by-dimension reconstruction, which we now illustrate. We assign the subscript ”0” to the cell where we are computing the coefficients, other values indicating location and direction with respect to (note that the notation is coherent with the fact that the first coefficient in the expansion , holds , i.e., the centered value). Next, for this particular problem we define three stencils

 S1={Q−2,Q−1,Q0},S2={Q−1,Q0,Q1},S3={Q0,Q1,Q2}, (42)

and in every stencil we compute a polynomial of the form

 Q(i)(x)=Q(i)0+Q(i)xP1(x)+Q(i)xxP2(x)i=1,2,3. (43)

The coefficients are given by

 S1 : Q(1)x=−2Q−1+Q−2/2+3Q0/2,Q(1)xx=(Q−2−2Q−1+Q0)/2, (44) S2 : Q(2)x=(Q1−Q−1)/2,Q(2)xx=(Q−1−2Q0+Q1)/2, (45) S3 : Q(3)x=−3Q0/2+2Q1−Q2/2,Q(3)xx=(Q0−2Q−1+Q2)/2. (46)

For every polynomial we calculate a smoothness indicator defined as

 IS(i)=(Q(i)x)2+133(Q(i)xx)2, (47)

leading to the following WENO weights:

 ω(i)=α(i)∑3i=1α(i),α(i)=λ(i)(ϵ+IS(i))r, (48)

where is a parameter introduced in order to avoid division by zero; usually . The scheme is rather insensitive to the parameter , which we set . The parameter is usually computed in an optimal way to increase the accuracy of the reconstruction at certain points; we opt for a centered approach instead, thus , while . The 1D reconstructed polynomial is given by

 Q(x)=ω(1)Q(1)(x)+ω(2)Q(2)(x)+ω(3)Q(3)(x). (49)

Next, a 1D reconstruction in the direction is performed in a totally analogous way. Finally, we address the computation of the mixed term , which is calculated in a 2D fashion. Keeping the same convention regarding location subscripts as in 1D, Altmann et al. (2007) considers 4 formulas for the cross term upon taking all the moments around the cell. The expressions for the cross term are:

 Q(1)xz = Q1,1−Q0,0−Qx−Qz−Qxx−Qzz, (50) Q(2)xz = −Q1,−1+Q0,0+Qx−Qz+Qxx+Qzz, (51) Q(3)xz = −Q−1,1+Q0,0−Qx+Qz+Qxx+Qzz, (52) Q(4)xz = Q−1,−1−Q0,0+Qx+Qz−Qxx−Qzz, (53)

and the corresponding smoothness indicators are given by

 IS(i)=4(Q(i)xx)2+4(Q(i)zz)2+(Q(i)xz)2. (54)

Note that in the first part of the reconstruction, when the weights were computed, a larger suboptimal weight was assigned to the central stencil, which is a way to ensure stability and robustness of the algorithm by sacrificing additional order in the approximation (for more details, see Dumbser and Käser (2007)). However, for this term, the numerators assigned to the corresponding ’s remains the same for every expression. The computation of this term concludes the reconstruction procedure, and now we have at our disposal one polynomial per cell that can be used to calculate values at the boundaries or inside the cell. The next step in our numerical scheme consists of the calculation of the numerical fluxes (38)-(39), which will use extrapolated boundary values of the reconstructed polynomials. Rather than the use of the classical WENO scheme (as in Chan et al. (1994) or Jiang and Shu (1996)), which performs this calculation via a first order flux, we opt for the WENO-TVD approach described in Titarev and Toro (2005). We make use of the 2D extension of the flux-limiter-centered scheme (FLIC) approach presented in Billett and Toro (2000); Toro (2009), which is a second-order, centered and non-oscillatory flux. In our case, it consists of a flux-limited version of a generalized Lax Wendroff flux, using as a low-order flux the GFORCE (generalized first order centered) flux Dumbser et al. (2009), which can be interpreted as a convex combination of Lax-Friedrichs and Lax-Wendroff-type of fluxes:

 FFLICi+1/2,j=FGFORCEi+1/2,j+ψi+1/2,j(FLWi+1/2,j−FGFORCEi+1/2,j), (55)

where

 FGFORCEi+1/2,j =FGFORCEi+1/2,j(QLi+1/2,j,QRi+1/2,j)=ωFLWi+1/2,j+(1−ω)FLFi+1/2,j, (56) FLFi+1/2,j =12(F(QLi+1/2,j)+F(QRi+1/2,j)−12ΔxΔt(QRi+1/2,j−QLi+1/2,j)), (57) FLWi+1/2,j =F(Q∗i+1/2,j), (58) Q∗i+1/2,j =12(QLi+1/2,j+QRi+1/2,j)−ΔtΔx(F(QRi+1/2,j)−F(QLi+1/2,j)). (59)

The parameter varies between 0 and 1, and is chosen in a compatible manner with the CFL number in order to ensure monotonicity; throughout the numerical experiments presented in this article, we will restrict ourselves to the usual FORCE flux, i.e., . We have omitted the formulas for the remaining cell boundaries, but they can be derived in a straightforward manner. Also note that even though the formulas are written along the boundary ’’, the use of the Gaussian quadrature formula will replace the axes’’ by Gauss points and therefore this subscript must be understood in that sense. It is important to notice that so far we are deriving expressions for the semi-discrete approximation of the system of conservation laws, however, the fluxes include the parameter which arises from the averaging operators that originate these fluxes. Thus, in the spatial discretization of the system, the time stepping enters just as a parameter. At the end of the derivation of the scheme, when we present the time discretization of eq. (34), will be considered as “marching parameter” in the sense that its inclusion in the formulas will generate an updated state in time. The function is a flux limiter; a slight variation of the usual limiters has to be considered in this context since we use a centered flux instead of an upwind approach (the reader can refer to (Toro, 2009, Ch 13.) for more details); in our case we mainly use the VanLeer limiter, which on its centered version reads:

 ψ(r)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0if r≤0,2r1+rif 0≤r≤1,ϕg+2r(1−ϕg)1+rif r≥1,ϕg=1−|c|1+|c|, (60)

where corresponds to the Courant number which depends on the problem. To implement the limiter in the ”conservation paradigm”, the flow parameter is defined via the total energy of the system,

 e=cvθπ+12(u2+v2+w2)+gz. (61)

The left and right flow parameters are then defined as,

 rLi+1/2,j=eRi−1/2,j−eLi−1/2,jeRi+1/2,j−eLi+1/2,j,rRi+1/2,j=eRi+3/2,j−eLi+3/2,jeRi+1/2,j−eLi+1/2,j, (62)

and finally,

 ψi+1/2,j=min(ψ(rLi+1/2,j),ψ(rRi+1/2,j)). (63)

The above described procedure starts with a set of averaged values and ends with a numerical approximation of the space operators involved in eq. (32). The resulting scheme is still continuous in time, and we conclude this section by discretizing this operator in a manner that is consistent with the choices that we have made in the generation of the space discretization operator. At a given starting time , we begin by considering the semi-discrete scheme

 dQi,j(t)dt=Li,j(Q), (64)

bringing the system to a final state with a time stepping . In order to preserve high-order and non-oscillatory properties in time, we consider the well-known family of explicit TVD Runge-Kutta schemes Gottlieb and Shu (1998), in particular its third order version

 Qn+13i,j = Qni,j+ΔtLi,j(Qni,j), (65) Qn+23i,j = 34Qni,j+14Qn+13i,j+14ΔtLi,j(Qn+13i,j), (66) Qn+1i,j = 13Qni,j+23Qn+23i,j+23ΔtLi,j(Qn+23i,j). (67)

We conclude this section with the inclusion of the source term. The source term appearing in eq. (33) will not depend on space nor space derivatives, and therefore it can be averaged in space and solved in the same manner as the above presented time discretization, by replacing by . If we denote by the the fully discrete operator that brings the system of conservation laws (32) units ahead in time, and by the fully discrete operator that updates the source term (33) in units, we preserve, at least, second order accuracy in time by implementing a Strang splitting Strang (1968) in the form

 Qn+1i,j=S(Δt/2)L(Δt)S(Δt/2)Qni,j. (69)

In this way we have derived a fully discrete, high-order and total variation diminishing numerical scheme for treating a system of balance laws arising from dimensional reduction of the original set of 3D primitive equations. In the next section, numerical experiments are performed with the aim of understanding the consequences of the dimensional reduction and the coupling effects, and its ability to represent atmospheric phenomena in a plausible way. The test are motivated by standard study cases for non-hydrostatic model development (see Klemp and Skamarock (1994); Robert (1993); Jebens et al. (2011)), and therefore a qualitative comparison with published model outputs is possible.

## 4 Numerical experiments

The purpose of the numerical experiments is to use well-known test cases that are used extensively to test dynamical cores of atmospheric models so that our results can easily be compared with results in the literature. There are quite a few such test cases published, and we have selected three. In addition we have constructed a new test case that is used to show the capabilities of the dimensionally-reduced model.

### 4.1 Convective bubble in a neutral atmosphere

This first test case, that has been previously studied in Robert (1993); Pudykiewicz and Smolarkiewicz (1992); Giraldo and Restelli (2008) (among others), is performed by using the scales and settings as in Jebens et al. (2011). The domain is , with and chosen according to,

 Δt=CFLΔxmaxi,j(|vel−cs|,|vel+cs)|, (70)

where the number is set to , and is the module of the velocity field (this setting of the time step will be kept for the remaining tests). Both layers are initialized at rest, and we consider the following potential temperature profiles:

 θ1=¯θ1+θ′1,θ2=¯θ2, (71)

with,

 (72)

The reference states and are in hydrostatic balance at their respective layers, and therefore the density is initialized via the hydrostatic balance law,

 cpθdπdz=−g, (73)

together with the equation of state for an ideal gas,

 ρ=P0RdθπcvRd, (74)

both evaluated at the corresponding reference state . We implement solid wall boundary conditions along the domain. The underlying physics of this test dictate that, as the temperature perturbation is warmer than the background state, a buoyancy force will push the perturbation upwards. As it starts rising, it also starts experiencing horizontal expansion because of the same buoyancy effect, eventually generating a mushroom-shaped cloud. This has been observed in many tests based on 2D sets of Euler equations. We observe a similar behavior in the layered model, as shown in figure 1. Figure 2 shows the evolution of the potential temperature residual between both layers, which preserves horizontal symmetry but decreases in magnitude; recall that the second layer is initialized without a perturbation and therefore the system tends to an equilibrium state between both layers. Velocity fields in figure 3 are consistent with results obtained in 2D tests performed in the references mentioned in the description of the test; vector plots in figure 4 illustrate the deformation process to which the temperature perturbation is subjected. Finally, 5 exhibits conservation of energy of the total system, and also confirms the equilibrium reached by the layers.

### 4.2 Interaction between hot and cold bubbles

This test is a variation of the test presented in Robert (1993), which studied the interaction between positive and negative potential temperature perturbations, modified to keep the same scale of our first test case; all the domain parameters are set as in test 4.1. Horizontal velocity is initialized in both layers with a value of in order to make the test more stringent. A periodic boundary condition is set in the lateral -direction while solid boundary walls are kept at the bottom and the top of the domain. Vertical velocity in both layers is set to zero. Thermodynamic variables are computed in exactly the same way as in the first test, but we add a cold perturbation to the second layer,

 θ′2={−15cos(πL2)L≤1,0i.o.c.,L=12000√x2+(z−8000)2. (75)

As we have previously described in the first test case, in a neutral atmosphere, a warm temperature perturbation is expected to raise; by the same principles, a cold perturbation is expected to fall. Figure 6 shows the initial temperature for both layers. Figure 7 illustrates the experiment: placed in the same vertical axis, the warm bubble raises while the cold bubble falls, and at certain instant both bubbles collide starting an eddy interaction which is governed by the same buoyancy effects. In this experiment we also included horizontal velocity and we can observe that the perturbations are horizontally translated with a proper speed (both layers were initialized at the same velocity) and preserving symmetry all along the experiment. Figure 8 shows a decrease in the magnitude of the potential temperature residual and figures 10 and 9 exhibit the associated velocity field, where we can observe the formation of eddies. Figure 11 illustrates the evolution of potential temperature extreme values in every layer, and it can be seen that for both maximum and minimum, the layers have a tendency to reach equilibrium values. Total energy is preserved as shown in figure 12 in a similar manner as in the first test.

### 4.3 Adjustment of shear between two layers

This third test case aims to study the level of interaction and adjustment between the two layers in the model. We again consider the same domain parameters, except for a final simulation time of . However, variables are initialized in a different way. First of all, horizontal velocity is set to a value of in the first layer, while in the second layer it has a value of . The horizontal velocity is initialized in the first layer with a logarithmic profile

 u1=50(log(zLZ+1))1/2 (76)

and it is set to zero in the second layer. Vertical velocity in both layers is set to zero. Thermodynamic reference states are different: no potential temperature perturbation is considered in any case, the first layer having the same neutral atmosphere as in the previous test cases, while the second is considered to be a stable atmosphere with

 ¯θ2=θ0exp(N2gz), (77)

where and is the Brunt-Väisälä frequency with a value of . In both cases is computed as in the previous tests. Boundary conditions are set again to be periodic in the lateral -direction and solid wall boundary conditions in .

In figure 13 we present the velocity profile for the velocity component in the x-direction. We see that the profiles are adjusted to some logarithmic profile from above in layer 1, and from below in layer 2. This is as expected. The same adjustment is seen for the potential temperature in figure 14, the adjusted atmosphere becomes stable with a stability which is a ”mean value” between the two layers.

In figure 15 we present the time evolution of the velocity component in the y-direction, in a cross section at . So here we expect to see a decay of an initial shear, produced by the initial shear in the x-direction. And this is what figure 15 shows, a decay from a positive value in layer 1, and from a negative value in layer 2.

Note the interesting ”undershoot” in layer 1, and the ”overshoot” in layer 2: The profiles maintains their shape, before the final decay towards zero. Recall that this test is without explicit friction or parameterizations of the shear layer, so it is basically thermodynamics () that controls the process. That the adjustment works well is also illustrated in figure 16, which shows the conservation of total energy and the adjustment of the energy in the two layers.

This test case can be extended to have a rising bubble in one of the layers as well as other perturbations to simulate front behavior; one can also introduce parameterizations of the shear layer, for example by considering a simple turbulence model (a Smagorinsky model is one choice).

### 4.4 Inertia-gravity waves through a periodic channel

The last case that we study is a variation of the inertia-gravity waves test case proposed by Klemp and Skamarock (1994). The purpose of this test is twofold: first, we qualitatively compare the behavior of our model for the generation of wave motion in a nonhydrostatic scale, but we also study the nonlinear effects of the coupling. We set , with . In a first model run, the horizontal velocity is set to a value of for the first layer, while any other velocity variable is set to zero. Both layers are initialized with the same stable atmosphere profile as in the previous example, and a temperature perturbation is added to the first layer, in the form,

 θ′1=θpsin(πzLZ)1+(x−xpa)2, (78)

with , , and . In a second model run, we leave the first layer at rest and add a velocity to the second layer. We remove the temperature perturbation from the first layer and assign it to the second, only changing the parameter , thus expecting a similar behavior as in the first model run but in the opposite direction. A third and final run considers a combination of the previous model runs, i.e., one perturbation in every layer, the first layer with a positive horizontal displacement and the second layer with a negative horizontal displacement. Figure 18 shows the results of the first and second model runs, initial condition and result at , first run in the two upper panels and second run in the two lower panels. These results are in very good agreement with the results in Klemp and Skamarock (1994) and other test results in the literature. In figure 19 we present the results of the third model run. We see immediately that these results are not superpositions of the results in figure 18, so we see the effects of a nonlinear interaction between the two inertia-gravity waves. In figure 20 we present the velocity fields for the nonlinear interaction case. We see time-dependent fan-like structures typical of nonlinear wave interaction and almost perfect symmetry of the fields. Figure 17 shows the evolution of extremal values for the horizontal velocity in every layer. It can be seen that both layers tend to an equilibrium, while an interaction between maximal and minimal values occurs. It is not easy to explain all the details in the results, for example the development from to , since the interaction is purely nonlinear. However, it is possible to see effects of focussing, which is a purely nonlinear wave phenomenon. In the solutions the symmetry is well preserved indicating that the is performing in a good manner; conservation of energy is shown in figure 21 .

## 5 Summary and concluding remarks

We have derived a set of dimensionally-reduced fully nonlinear primitive equations to model atmospheric dynamics. Our goal was to show that the 2.5D model is able to reproduce well-known atmospheric phenomena properly. A DG approach was used for the dimensional reduction, whereas a WENO-TVD method was used for discretization. We think that finite volume methods are the most promising discretization method for atmospheric models since it ensures discrete conservation. If one uses, as we have done, equations in conservative form with conserved variables, the discretized equations should model the conserved quantities in the atmosphere with high accuracy. The numerical experiments with extensively used test cases for atmospheric model cores have shown that the dicretization is accurate, stable, energy-conserving and do not produce spurious oscillations. No artificial diffusion or any other form of smoothing is used. We think that this is a strong indication that our set of equations are well-posed, and that a higher order discretization is indeed valuable. Our approach can be extended in several directions: we could use higher order elements in the dimensional reduction and/or even higher order finite volume methods. The purpose of this would be to investigate the accuracy versus the computational costs. This can easily be done for well-known test cases where we have results to compare with, but we would also like to use more complicated tests cases from real weather situations. The latter test cases could be used to investigate how one can model certain weather phenomena where one, from a physical viewpoint, have a certain insight in the behavior. There is also an option of using other methods than the DG method for the dimensional reduction, but our opinion is that this method is good natural choice, because it poses no difficulty for handling nonlinear equations. Also on the computational side, the actual computational costs has not been explicitly addressed in this paper. However, it is worth to mention that every block of the code can be completely vectorized, while the dimensionally-reduced equations have a natural form of parallelism by parallelizing over the layers. This is an important topic since the existing dynamical cores of atmospheric models which we would like to compare with in term of computational cost, are parallelized.

## References

• Oliger and Sundström (1978) Oliger, J., Sundström, A.. Theoretical and practical aspects of some initial boundary value problems in fluid dynamics. SIAM Journal on Applied Mathematics 1978;35(3):419–446.
• Chen et al. (2008) Chen, Q., Temam, R., Tribbia, J.. Simulations of the 2.5d inviscid primitive equations in a limited domain. Journal of Computational Physics 2008;227(23):9865–9884.
• Giraldo and Restelli (2008) Giraldo, F., Restelli, M.. A study of spectral element and discontinuous galerkin methods for the navier-stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases. Journal of Computational Physics 2008;227(8):3849 – 3877.
• Klemp and Skamarock (1994) Klemp, J., Skamarock, W.. Efficiency and accuracy of the klemp-wilhelmson time-splitting technique. Monthly Weather Review 1994;122(11):2623–2630.
• Cockburn and Shu (2001) Cockburn, B., Shu, C.W.. Runge-kutta discontinuous galerkin methods for convection-dominated problems. Journal of Scientific Computing 2001;16:173–261.
• Titarev and Toro (2005) Titarev, V., Toro, E.. Weno schemes based on upwind and centred tvd fluxes. Computers & Fluids 2005;34(6):705–720.
• Billett and Toro (2000) Billett, S., Toro, E.. Centred tvd schemes for hyperbolic conservation laws. IMA Journal of Numerical Analysis 2000;20(1):47–79.
• Kalise et al. (2011) Kalise, D., Lie, I., Toro, E.. High-order finite volume schemes for layered atmospheric models. Paper I of this thesis; 2011.
• Chen and Toro (2004) Chen, G.Q., Toro, E.. Centered difference schemes for nonlinear hyperbolic equations. Journal of Hyperbolic Differential Equations 2004;1(3):531–566.
• Altmann et al. (2007) Altmann, C., Balsara, D.S., Dumbser, M., Munz, C.D.. A sub-cell based indicator for troubled zones in rkdg schemes and a novel class of hybrid rkdg+hweno schemes. Journal of Computational Physics 2007;226(1):586–620.
• Balsara et al. (2009) Balsara, D.S., Dumbser, M., Munz, C.D., Rumpf, T.. Efficient, high accuracy ader-weno schemes for hydrodynamics and divergence-free magnetohydrodynamics. Journal of Computational Physics 2009;228(7):2480–2516.
• Dumbser and Käser (2007) Dumbser, M., Käser, M.. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems. Journal of Computational Physics 2007;221(2):693–723.
• Chan et al. (1994) Chan, T., Liu, X.D., Osher, S.. Weighted essentially non-oscillatory schemes. Journal of Computational Physics 1994;115(1):200 – 212.
• Jiang and Shu (1996) Jiang, G.S., Shu, C.W.. Efficient implementation of weighted eno schemes. Journal of Computational Physics 1996;126(1):202–228.
• Toro (2009) Toro, E.F.. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer-Verlag; third ed.; 2009.
• Dumbser et al. (2009) Dumbser, M., Hidalgo, A., Toro, E.. FORCE schemes on unstructured meshes I: Conservative hyperbolic systems. Journal of Computational Physics 2009;228(9):3368–3389.
• Gottlieb and Shu (1998) Gottlieb, S., Shu, C.W.. Total variation diminishing runge-kutta schemes. Mathematics of Computation 1998;67(221):73–85.
• Strang (1968) Strang, G.. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 1968;5(3):pp. 506–517.
• Robert (1993) Robert, A.. Bubble convection experiments with a semi-implicit formulation of the euler equations. Journal of the Atmospheric Sciences 1993;50(13):1865–1873.
• Jebens et al. (2011) Jebens, S., Knoth, O., Weiner, R.. Partially implicit peer methods for the compressible euler equations. Journal of Computational Physics 2011;230(12):4955–4974.
• Pudykiewicz and Smolarkiewicz (1992) Pudykiewicz, J., Smolarkiewicz, P.. A class of semi-lagrangian approximations for fluids. Journal of the Atmospheric Sciences 1992;49(22):2082–2096.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters