[April 13, 2013] A Conservative Integration of the Pseudo-Incompressible Equations with Implicit Turbulence Parameterization [Monthly Weather Review] Tweet (Monthly Weather Review Via Acquire Media NewsEdge) ABSTRACT Durran's pseudo-incompressible equations are integrated in a mass and momentum conserving way with a new implicit turbulence model. This system is soundproof, which has two major advantages over fully compressible systems: the Courant-Friedrichs-Lewy (CFL) condition for stable time advancement is no longer dictated by the speed of sound and all waves in the model are clearly gravity waves (GW). Thus, the pseudo-incompressible equations are an ideal laboratory model for studying GW generation, propagation, and breaking. Gravity wave breaking creates turbulence that needs to be parameterized. For the first time the adaptive local deconvolution method (ALDM) for implicit large-eddy simulation (ILES) is applied to non- Boussinesq stratified flows.ALDMprovides a turbulence model that is fully merged with the discretization of the flux function. In the context of non-Boussinesq stratified flows this poses some new numerical challenges-the solution of which is presented in this text. In numerical test cases the authors show the agreement of the results with the literature (Robert's hot-cold bubble test case), they present the sensitivity to the model's resolution and discretization, and they demonstrate qualitatively the behavior of the implicit turbulence model for a 2D breaking gravity wave packet. (ProQuest: ... denotes formulae omitted.) 1. Introduction A main effect of acoustic waves in the atmosphere is a rapid adjustment to a balanced state with respect to pressure-density perturbations. As long as one is not interested in this adjustment process itself but rather in atmospheric dynamics on comparatively longer time scales, it is often attractive to filter acoustic waves from the dynamic equations. Thus, simplified (i.e., soundproof) equations enable a focus on the dynamics of more significant mesoscale processes, while the absence of fast waves also allows for larger time steps in numerical integrations. A classic example is the Boussinesq system, where the prognostic continuity equation is replaced by a diagnostic divergence constraint on the wind field. This yields a good approximation of atmospheric dynamics on spatial scales smaller than the atmospheric-scale height. Deep convection and wave propagation, however, experience effects of the vertical decrease of ambient density that cannot be reproduced in the Boussinesq approximation. This density effect causes, for example, gravity waves-typically generated in the troposphere-to grow in amplitude until they get unstable, break, and finally dissipate and interact with the large-scale flow in the middle atmosphere (Lindzen 1981; Fritts and Alexander 2003). To include such effects, two more general approaches have been proposed. One is the well-known anelastic equations first introduced by Batchelor (1953) and Ogura and Phillips (1962) and further developed by Lipps and Hemler (1982) and Lipps (1990). They use a more general divergence constraint, which takes the decreasing background density into account as aweighting factor for the wind. A potential drawback of the anelastic system is that it requires the leading-order background atmosphere to be close to isentropic. Some indications exist that corresponding errorsmight often be negligible even for realistic tropospheric stratifications (Smolarkiewicz and Dörnbrack 2008; Klein et al. 2010; Smolarkiewicz and Szmelter 2011). However, an approach that does not constrain the background potential temperature as much as in the anelastic equations is offered by the pseudoincompressible system (Durran 1989; Durran and Arakawa 2007; Durran 2008). By introducing the socalled pseudo-incompressible density, the compressible coupling between pressure and density fluctuations and thus the transfer from elastic potential to kinetic energy is cut off. Smolarkiewicz and Dörnbrack (2008) underlined that the Durran system retains the full momentum equation and, consequently, admits unabbreviated baroclinic production of vorticity, as also observed by Klein (2009). Achatz et al. (2010) demonstrated a multiscale asymptotic consistency between the pseudo-incompressible system and the compressible Euler equations to leading order, which cannot be observed for the anelastic system. The pseudo-incompressible equations are thus suitable for studying full-scale gravity wave dynamics from their generation until their breaking. In this problem the distances covered bywaves of interest, both in the horizontal and in the vertical, are typically far too large to allow an explicit direct simulation of the generated turbulence. For its parameterization an implicit large-eddy-simulation (ILES) method as developed by Adams et al. (2004) and Hickel et al. (2006, 2007) seems attractive. This socalled adaptive local deconvolution method (ALDM) represents a merger of numerical method and subgridscale (SGS) turbulence parameterization. Instead of keeping numerical truncation errors small, which can be quite costly or unfeasible, and later dissipate the solution with an explicit turbulence model, numerical truncation is deliberately exploited and tuned to act as turbulence model. This holistic approach results in a particularly reliable and efficient method, while turbulence spectra are matchedwell.ALDMhas been calibrated to be consistent with eddy damped quasi-normal Markovian (EDQNM) spectral turbulence theory (Hickel et al. 2006). Indeed, comparisons with direct numerical simulation (DNS) data of turbulent flows showed better spectral behavior than standard turbulence models. ALDM was originally developed for incompressible flows (Hickel et al. 2006, 2007) and subsequently extended to fully compressible turbulence (S. Hickel 2011, unpublishedmanuscript; Hickel and Larsson 2008). Although its already widespread application in engineering applications (Klar et al. 2011; Hickel et al. 2011) had been limited to the description of the effects of unstratified turbulence, recent results by Remmler and Hickel (2012a,b) demonstrate that the method can also be applied to stratified turbulence governed by the Boussinesq equations. Here we report the first application of ALDM to full-scale atmospheric flows beyond the Boussinesq approximation. For this purpose a novel fully conservative discretization method for Durran's pseudoincompressible system has been developed, as previously achieved, to the best of our knowledge, only by Smolarkiewicz and Dörnbrack (2008) and Smolarkiewicz and Szmelter (2011). As the implicit turbulence model of ALDMis fully merged with the spatial discretization, an extension to another system of partial differential equations can be challenging. This paper shows how these challenges have been resolved within the framework of our newly developed atmospheric flow solver Pseudo Incompressible Flow Solver with Implicit Turbulence (pincFloit) model. We emphasize that many numerical issues discussed here, such as the convergence criterion for the pseudo-incompressible equations and a divergence correction for reducing the number of Poisson-solver iterations, are relevant to any kind of finite-volume method. The paper is organized as follows. In section 1 we recall the pseudo-incompressible equations and present them in a conservative and scaled form suitable for the discretization by a conservative finite-volume model. The general approach to the discretization is given in section 2. Some special issues concerning pitfalls in the application of conservative finite-volume methods to stratified fluids with the pseudo-incompressible divergence constraint are discussed in section 3. Numerical test cases in section 4 close the text, while section 5 gives a summary and conclusions. a. Governing equations The pseudo-incompressible equations were first derived by Durran (1989). If we restrict ourselves to adiabatic nonrotating dynamics, then the Durran system in conservative flux form reads ...(1) ...(2) ...(3) with the Mach number Ma, the Froude number1 Fr, and the Reynolds number Re; the velocity u = (u, y, w) with its zonal, meridional, and vertical components, respectively; and the potential temperature u and the pseudoincompressible (effective) density r. Further quantities are the Exner pressure ...(4) with k = R/cp, where cp is the specific heat at constant pressure and R is the specific gas constant of air; a constant reference pressure p0; and the pressure p. Following the notation of Klein (2009), we introduce P defined by ...(5) with the ratio of specific heat at constant pressure and constant volume g = cp/cV set to 1.4 in all our calculations. The pseudo-incompressible (effective) density is defined through the equation of state ...(6) This effective density is only coupled to the potential temperature and no longer to the (total) Exner pressure, which was replaced by the background-state Exner pressure p in the equation of state. Through this approximation the system of equations becomes soundproof. Since the background state-given by ...-satisfies the equation of state, the pseudo-incompressible (effective) density2 satisfies ...(7) For better numerical accuracy and to avoid the wellbalancing problem it is sensible to work with the fluctuation quantities for Exner pressure ...(z) and potential temperature u0 ...(z). The viscous stress tensor P is given by ...(8) with the Kronecker symbol dij and the dynamic shear viscosity coefficient h. The effect of volume viscosity is not considered. Note that the viscosity term is responsible for the resolved grid-scale dissipation and is part of the model even if an implicit LES is used since the latter treats the subgrid-scale dissipation. b. Transport of potential temperature Substituting the density of the thermodynamic relation (5) into the effective continuity equation (2) and using the divergence constraint (3) we obtain a thermodynamic equation for the transport of potential temperature ...(9) which is equivalent to D(1/?)/Dt=Du/Dt=0. This shows that without heat sources the potential temperature is simply advected with the flow. If we combine effective continuity [(2)] and transport of potential temperature [(9)] we obtain the relation ...(10) for the mass-weighted potential temperature. Using the divergence constraint ... (??u) we see that ?? is constant in time. Note that for diabatic flows a heat source enters the RHS of the divergence constraint and the transport equation of the potential temperature. 2. Discretization with ILES We present the implicit LES proposed by S. Hickel (2011, unpublished manuscript), Hickel et al. (2006, 2007), and Hickel and Larsson (2008) named ALDM. Originally introduced for incompressible and fully compressible flows, we adopt the approach to the pseudo-incompressible equations. We present general ideas in the main part and put some of themore technical aspects into the appendices. a. General setup All variables are stored in a C-grid fashion. Figure 1 shows a finite-volume cell of the density with the velocities defined at the cell interfaces. Each momentum component has its own cell shifted with respect to the mass cell by half a cell in the corresponding direction. To calculate the fluxes across the cell surfaces, point values of velocity and scalars are reconstructed there. To denote the reconstructed value at the right wall of cell ijk we use the symbol ... and the superscripts L, F, B, U, and D for reconstruction to the left, forward, backward, upward, and downward cell interface, respectively; please see also Fig. 2. The equations for these reconstructions terms are given in the appendices. b. Adaptive deconvolution (reconstruction) To obtain the values at the cell interface, ideas from numerical gas dynamics are borrowed and extended. In the weighted essentially nonoscillatory (WENO)method (Shu 1997) data within a cell is reconstructed using the cell averages from neighboring cells as depicted in Fig. 3 on the example of a fifth-order WENO scheme with third-order reconstruction polynomials. Note that the parabolas do not interpolate the data points since they are finite-volume cell averages and not finite-difference point values. At the cell interface xi11/2 three reconstructed values exist from the three stencils. In WENO the reconstructed values are averaged using solution adaptive weights. These weights are chosen in such a way that the order of accuracy in smooth flows approaches the maximum for the given stencil width while unphysical oscillations are avoided by keeping the total variation in nonsmooth regions bounded. ALDM differs from WENO in two ways. First, with ALDM the weights in the convex combination of the reconstruction polynomials are used to optimize the nonlinear spectral numerical dissipation that acts as an implicit SGS model rather than maximizing the formal order of accuracy. Second, in contrast to WENO, where polynomials of a certain degree n are used to obtain a convergence order 2n 2 1, in ALDM first-, second-, and third-order reconstructions are blended. The blending weights, again, are degrees of freedomthat are used for SGSmodel tuning. Formore details please see the appendices or Hickel et al. (2006, 2007). c. Numerical flux function To calculate the update of a quantity in cell (i, j, k) we need to discretize the divergence operator, which is done in a conservative finite-volume way: ...(11) where f, g, and h indicate the fluxes in the three spatial directions as depicted in Fig. 4 for the flux of zonal momentum. The numerical flux function fi+1/2,j,k is the second ingredient of ALDM for ILES. It consists of a central term of high order and an artificial viscosity term. Mass and momentum have slightly different flux functions. 1) MASS FLUX FUNCTION For scalar transport it was shown in (Hickel et al. 2007) that the following flux function leads to an implicit SGS model that is consistent with turbulence theoretical constraints: ...(12) where ...(13) The first part of the flux function leads to a high-order central difference. Here it is important to work with the filtered (cell centered, volume averaged) velocity u, which satisfies the pseudo-incompressible divergence constraint. The second term is a dissipative term, with a diffusion coefficient given in (13), which itself depends on the roughness of the velocity field and a tuning parameter sC. In a smooth velocity field a density distribution is (almost) not diffused numerically with this SGS ansatz. Throughout our calculations we use for the tuning parameter ...(14) as proposed in Hickel et al. (2007) for air. 2) MOMENTUM FLUX FUNCTION The location of the fluxes on a C-grid is depicted in Fig. 4. To showthe principle of the numerical flux function for the momentum we give the meridional transport of zonal momentum ...(15) where ...(16) At first we note that all quantities used in (15) are reconstructed to the forward interface center of the momentum cell of rui,j,k. The interpolation of the density ... is-at first glance-not unique but it can be constrained so as to make momentum and effective continuity equation consistent-an aspect discussed in section 3a. The first part of the flux function (15) is a central interpolation leading to a highorder central difference for the divergence operator. The second term is a numerical diffusion term that dissipates momentum according to the magnitude of the second derivative of the reconstructed (deconvolved in LES terms) velocity u. The SGS dissipation coefficient ... given in (16) is proportional to the roughness of the volume-averaged (filtered in LES terms) velocity u. A list of all flux components can be found in appendix A. Also note that su depends on the location if the grid is not uniform-that is, with varying mesh size (Hickel et al. 2006). In pincFloit a uniform Cartesian mesh is used so that this parameter is a constant ...(17) 3. Numerical issues The implicit LES was devised for incompressible flows by Hickel et al. (2006) and for fully compressible flows by Hickel and Larsson (2008) and S. Hickel (2011, unpublished manuscript). The application of ALDM to a model with atmospheric background stratification and a background density rapidly decaying in the vertical lead to a number of problems-the solution of which we discuss in this section. a. Flux function: Consistency between continuity and momentum equation So far we have not devised a rule how to interpolate the densities p needed in the momentum fluxes. We require that a pure density distribution with homogeneous background wind should solely advect the density while maintaining the constant background wind for all times. This goal can be achieved by assuring that effective continuity and momentum equation do-in this special setting-exactly the same, namely pure transport of mass. Note that inconsistent interpolation of p leads to the appearance of artificial velocity fluctuations. For simplicity, we assume at first a constant wind in zonal direction. We list all assumptions for the state at time t = tn: ... The following discussion applies most directly to multistage time-stepping schemes, such as Runge-Kutta, which can be seen as a series of appropriately weighted forward Euler steps, here called predictor step. With only zonal transport of momentum the equation for velocity predictor ... reads: ... For ease of reading we omitted the spatial indices j and k. The density at the zonal-velocity position is taken to be pi+1/2 = (pi + pi+1)/2. Note that for constant u the correction term containing s is absent. The position of density, velocity, and fluxes can be seen in Fig. 5. We now enforce ... (i.e., no change of velocity in time), apply the assumption ui = u = const, and obtain ...(18) for the momentum equation. With the same assumptions, the continuity equation for cell i reads ... and for cell i + 1 ... The arithmetic mean of both continuity equations yields ...(19) Note that this arithmetic averaging agrees with the interpolation ...(20) used to obtain the density for the momentum ri11/2,j,kui,j,k. We now require that momentum transport [(18)] and mass transport [(9)] are equal: ... The two differences are depicted in Fig. 6. It can easily be checked that this equality is obtained if the density ^r is interpolated like ...(21) where we have added the two missing indices j and k for completeness. This interpolation rule, which is the only possible solution of the form ..., is illustrated in Fig. 7. This reasoning can be generalized to 3D. A list of all interpolation rules can be found in appendix B. b. Temporal discretization 1) CHOOSING THE TIME STEP PincFloit can be run with a fixed time step or with a variable time step chosen according to several stability criteria to obtain an efficient and stable numerical model. The stability criterion by Courant, Friedrichs, and Levy (CFL) limits the time step in proportion to the gridpoint spacing but inversely proportional to the advecting wind speed. We thus choose the time step in the following way for our model: ...(22) The CFL number n is set to 0.9 in most of the calculations in conjunction with third-order Runge-Kutta schemes. The CFL condition for a thermal bubble-initially at rest in a windless atmosphere-would allow an infinitely large first time step. Therefore the acceleration a must be considered as well (S. Remmler, TU Munich, 2010, personal communication): ...(23) Solving for the positive root of ?t leads to ...(24) To avoid division by zero, we follow the algorithm, here in 1D: first obtain umax=max|u| and amax5maxjaj, then calculate the time step according to the standard CFL condition ?tCFL = min(v?x/|u|, ?tmax) and if the speedup ?u = amax?tCFL is essential (i.e., if ?u > eu max with e.g., e = 10-2), then compute ...(25) We want to resolve the oscillations related to gravity waves (GW) in time. The highest possible GW frequency is the Brunt-Väisäläfrequency N. A sensible time step limit is therefore given by ?tGW=1/N. For our standard GWP test case with an isothermal background of T = 300 K and N = 0.018 s-1 we obtain a time step limitation due to gravity wave oscillations of ?tGW = 55 s. The time step to be used by the scheme is the minimum of the time steps dictated by the various time step restrictions and a Dtmax to avoid too-large numerical errors. 2) LOW-STORAGE RUNGE-KUTTA AND DIVERGENCE CONSTRAINT We have implemented (among other time schemes) the low-storage Runge-Kutta method of third order by Williamson (1980): ...(26) ...(27) ...(28) We derive the pressure correction procedure applied within each Runge-Kutta stage and show how this additional projection step changes the convergence properties of the overall scheme. Without heat source, the divergence constraint on the velocity ...(29) can be ensured by correcting the velocity field after the predictor step. Care has to be taken because this pressure correction has to be applied within each Runge- Kutta stage m. The equation for the tendency and the predictor of the momentum are given by ...(30) ...(31) ...(32) where ui denotes a velocity component and rm11 is known from the update of the effective continuity equation; FM comprises all fluxes, forces, and sources except for the pressure gradient; and am, bm are the weights for the low-storage Runge-Kutta method by Williamson (1980). The old pressure p0m does not guarantee the divergence constraint on u*. Therefore, a new pressure p0m11 is sought that ensures that the velocity ... at the next Runge-Kutta stage satisfies (29): ...(33) Taking the difference between the last two equations (32) and (33) and applying the discrete divergence operator d/dxi we obtain for the corrector procedure ...(34) ...(35) ...(36) where in the second line we sum over double indices giving the discrete divergence operator on the LHS and a Laplacian on the RHS. Equation (35) is a discrete Poisson equation for the pressure correction ?p' with the right-hand side given as the residual divergence of the predicted velocity field u* weighted with P. Note that the predicted tendency ?u* must also be corrected because it is used in the following Runge-Kutta stage; that is, we have to supplement the equation ...(37) with ... . Alternatively to the correction of velocity, the momentum tendency could be corrected ...(38) It is known that a fractional step method deteriorates the temporal convergence order of the overall numerical scheme. To quantify this error we tested three third-order Runge-Kutta schemes with the 1D gravity wave packet test case: the low-storage Runge-Kutta by Williamson (1980) (LS-Will-RK3), a low-storage total-variation-diminishing (TVD) Runge-Kutta (LS-TVD-RK3), and a classical, nonlow-storage TVD Runge-Kutta scheme (CL-TVDRK3) (Gottlieb and Shu 1998). The results are presented in Fig. 8. If the projection step is switched offall third-order schemes show the expected third-order convergence as Dt / 0. If the projection is switched on the convergence drops down to first order as a simple forward Euler. Note, however, that the overall error of the third-order schemes is one order of magnitude smaller so that it pays offto work with the higher-order schemes. c. ALDM with background stratification It turned out that the turbulence model should not be applied to the full density but only to the perturbations. Otherwise the model would ''see'' roughness in the background, which it would try to smooth out (i.e., it would try to establish a linear background profile). This can also be seen on the equation level: since the divergence in the effective continuity [(2)] must be discretized with the ALDM flux function it cannot be consistent with the divergence constraint (3). We consider an unperturbed isentropic background, for simplicity, so that the divergence constraint becomes ...(39) and the effective continuity equation with p=p simplifies to ...(40) The straightforward discretization of (2) does not satisfy (40). The discrete divergence ...(41) with ...(42) does not vanish because the background density has a nontrivial profile leading to reconstruction values that no longer satisfy the divergence constraint ...(43) A simple way out is to split the total density into the background p and a fluctuation part p' and apply the reconstruction only to p'. The ILES damping term in the flux will only depend on ..., while the density in the convective term will be composed of background and reconstructed fluctuation: ...(44) The flux in the above situation takes the form ...(45) which leads to ...(46) Accordingly, the flux function (12) is replaced by ...(47) In Fig. 9 the effect of introducing the fluctuation density is shown with isolines of potential temperature for the hot bubble test case at t = 20 min (see section 4 for the setup). On the leftthe discretization with the total density leads to oscillations in the solution, while the introduction of density fluctuation (right) produces a smooth solution. The small-scale structures on the leftare on the grid scale and are numerical artifacts. Especially the oscillations at the top of the domain are unphysical since the solution should remain smooth away from the bubble-induced vortices. The smooth solution on the right compares well with other simulations of the test case in the literature (e.g., Robert 1993; Klein 2009; Mendez-Nunez and Carroll 1994). A similar positive effect can be observed for the 1D gravity wave packet (see section 4 for the setup) shown in Fig. 10. If the total density is reconstructed, severe oscillations appear in the solution (left). ALDM applied only to the density fluctuations leads to a smooth solution (right), which is the correct solution since the gravity wave should not break and produce turbulence in this regime. d. Solving the Poisson problem Throughout our calculations we use the biconjugate gradient stabilized method (BiCGSTAB) (Meister 1999) algorithm to solve the Poisson problem. A small parameter e is introduced in the abort criterion, which will be discussed in section 2. 1) SCALING OF THE POISSON EQUATION If we solve the Poisson equation to satisfy ... , we obtain a height-dependent error in the velocity since ... depends exponentially on height. In Fig. 11 this is shown for a (initially) uniform, isothermal atmosphere at rest with a domain ranging up to zmax = 150 km with a Poisson-solver tolerance of e = 10-7 and a fixed time step of 1 s. One way out of this would be to tighten the tolerance-a very inefficient way.We consider a different approach: scaling of the Poisson equation. To see the necessity we analyze the problem by looking at the influence of the divergence error on the effective continuity equation and, consequently, via the momentum equation on the velocity field. The effective continuity equation ... with the density split as ... reads ...(48) The second term is related to the divergence constraint in the following way. Assuming a Poisson-solver tolerance of e and using ... , then ...(49) and, consequently, ...(50) The divergence error in the density therefore has the following height dependency, assuming the second term of the RHS of (50) to be negligible in this order estimate: ...(51) An error in the density leads to a buoyant acceleration in the vertical momentum equation ...(52) ...(53) ...(54) ...(55) where we used relation (51) for the last step. Now it is obvious that the velocity perturbation w0 created because of the divergence error is height dependent and will grow exponentially like ... . To obtain height independency for w0 we see from (54) that the divergence error in the density should scale as ...(56) Then the divergence error in the effective continuity equation [i.e., first term of the RHS of (50)] should satisfy ...(57) which is equivalent to the following scaling of the divergence constraint ...(58) If we apply this pressure scaling ... while keeping the Poisson-solver tolerance at e = 10-7 we obtain the results presented in Fig. 12. With the same computational effort, velocity and potential temperature perturbations were reduced by four orders of magnitude. Note for the anelastic divergence constraint a scaling with p would be in order. 2) PHYSICAL CONVERGENCE CRITERION Smolarkiewicz et al. (1997) discussed stopping criteria for the iterative solution of the Poisson problem in the context of the Boussinesq system. Their main argument is a physical one: an iteration should be stopped if the relative density change ?p/p due to the divergence error within one time step ?t is smaller than some e. A convergence to machine precision would be a waste of time since other discretization errors are then orders of magnitudes larger. ... Here we discuss it in the context of the pseudoincompressible system. From (48) and (50) we recall that the divergence constraint enters the effective continuity equation in the form ...(59) The density change due to the divergence error in one time step is ...(60) and the relative change with respect to the background density is ...(61) To ensure that the relative change of density is small throughout the computational domain, a sensible abort criterion is ...(62) Note that the scaling of the divergence is consistent with the scaling we introduced in the section above. This is because we also assume in both cases that density errors are proportional to the background density. 3) CORRECTING THE DIVERGENCE ERROR Smolarkiewicz and Margolin (1994) showed that there can be an enormous difference between the Eulerian and the semi-Lagrangian approach: the flux-formformulation needs a much smaller Poisson-solver tolerance e in the iterative elliptic solver to avoid unphysical solutions than the advective formulation.3 This leads to a higher number of iterations and increase in computing time. In our model we also observed this inefficiency of the Poisson solver and propose the following remedy, which, on the one hand, violates the conservation property of the effective continuity equation. On the other hand this error is of the order of the Poisson-solver tolerance and is therefore acceptable-depending on the application. To understand the modification let us write once more the effective continuity equation ... in a different form: ...(63) so that the term containing the divergence error becomes visible. If the tolerance e is very small the term can be neglected in all other cases it has an (unphysical) effect. To overcome this discrepancy, we simply subtract this term from the equation, which is the same as adding it as a source term in the effective continuity equation: ...(64) In the following we study the effect of the proposed ''divergence-error correction'' for the 1D gravity wave packet, which is a very sensitive test case. In Fig. 13 the 1D gravity wave packet described in section 4 is shown. If the effective continuity equation is not changed (leftcolumn) we can see that for a Poisson-solver tolerance of e = 10-5 (top) the scheme is not even stable and for e = 10-7 (bottom) there are oscillations in the higheraltitude regions. In contrast, the introduction of the source term in the effective continuity equation stabilizes the scheme and produces correct results for tolerances as low as e = 10-5. In Table 1 we see how the average number of iterations per call to BiCGSTAB4 is reduced by introducing the divergence-error correction. The reason for this behavior is clear: if the error is not produced in the effective continuity equation, no unphysical dynamic is induced in the momentum equation-and an initially unperturbed, divergence-constraint satisfying statemaintains this divergence constraint. Note that for e = 10-8 the divergence-error correction no longer has influence. Also note that in some instances BiCGSTAB is not called at all, leading to average numbers smaller than 1. 4. Test case results a. Hot bubble test case For the hot bubble test case described by Mendez- Nunez and Carroll (1994) we assume a neutrally stratified atmosphere with u0 = 300 K. The initial perturbation by the hot thermal is given by ...(65) with the radial distance given by ...(66) The radius is set to r0 = 2.5 km and the initial height is zc = 2.75 km. The bubble is placed horizontally in the middle of the domain. This test case was chosen to demonstrate the importance of transporting only the density perturbation instead of the full density; see Fig. 9 and the discussion in section 3c. The test shows interfacial instabilities as already reported in Grabowski and Clark (1991) for a thermal with moist physics and a high-resolution (direct simulation) approach. We now consider a variant of the hot bubble test case with a more realistic atmosphere. The troposphere is assumed isentropic with utr5300 K, the tropopause is set at ztr = 12 km, and the stratosphere is assumed isothermal with a constant temperature given by the temperature at the tropopause ...(67) with the pressure at the tropopause ...(68) The potential temperature profile in the stratosphere is given by ...(69) The result of the simulation with the implicit turbulence model at a resolution of 100 × 100 at 20-min model time is shown in Fig. 14. Depicted are the isolines of potential temperature perturbation u0 5u2u ranging from 23.2 to 2.2 Kin steps of 0.5 K. Positive perturbations are marked with a ''1'' and negative perturbations with a ''2.'' As for a convective cell in moist atmosphere, the tropopause acts as a natural border for the vertical movement of the hot bubble and its spreading looks similar to experiments with a solid wall at z5ztr.Nevertheless, the updraftproduced by the hot bubble induces perturbations at the tropopause, which travel as gravity waves into the stratosphere. Wave crests and troughs are visible in (14) as elliptic isolines with positive and negative values of potential temperature. b. Bubble test case by Robert Robert (1993) discussed the test case of two colliding bubbles: a large hot bubble with a small cold bubble placed horizontally offcentered above. The background is an isentropic atmosphere with ?=300K. Both bubbles have the following structure: ...(70) with ...(71) The warm bubble has the data ?u0 = 0.5 K, r0 = 150 m, and s = 50 m, and is positioned at x0 = 500 m, z0 5 300 m. The smaller cold bubble has the data ?Du0 5 20.15 K, r0 = 0 m, and s = 50 m, and is positioned at x0 = 560 m, z0 = 640 m. The domain spans [2500 m, 500 m] × [0, 1000 m]. No-slip boundary conditions are applied at all four walls of the domain. The CFL number is set to 0.9, leading to a time step of 3.0 s for 100 3 100 and 1.5 s for the 200 × 200 resolution. An interesting aspect is the existence of two scales introduced by the large and small bubble and the fact that the cold bubble is horizontally ''offcentered,'' which creates a highly unsymmetric solution. In Fig. 15 we compare the isolines of density perturbation r0 5r2r after 10 min for two resolutions: 100 × 100 on the leftcolumn and 200 × 200 on the right. The solution for the Euler equations with an upwind scheme by Robert is presented in the top line of the figure; while the solution with pincFloit- that is, for the pseudo-incompressible equations with ALDM as implicit turbulence model-is shown in the bottom line. The resemblance of the solutions is remarkable. A major difference is the vortex at the lower-right corner: pincFloit produces more structure in this flow region. As it turns out this structure is also obtained with an upwind scheme [pincFloit with a monotone upstream-centered scheme for conservation laws (MUSCL) as flux solver] for a higher resolution (400 3 400). This implies that ALDM has less numerical dissipation than standard (second order) upwind schemes at a comparable resolution. A variant of Robert's test case with a uniform horizontal background wind u0 = 1.67 m s21 is shown in Fig. 16 with a resolution of 200 3 200. The velocity is set so that the bubbles fully traverse the periodic domain within 10 min. The plotted levels for the isolines of density perturbation r0 5r2r are the same for the experiment without background wind (leftcolumn) and with uniform background wind (right column). The results with and without wind are not identical but very close to each other. The small deviation can be attributed to numerical truncation errors, which are different if the fluid system is shifted along the grid and exponentially grow in time owing to the nonlinear transport processes. c. 1D gravity wave packet On an isothermal reference state with T00 = 300 K a wave packet is superimposed with the initial buoyancy amplitude ...(72) with the vertical wavenumber m0 = 2p/lz,0 and the normalized amplitude a = 0.1, the isothermal Brunt- Väisäläfrequency for T00 = 300 K of N = 0.018 s21, the half-width s = 5 km, and the center of the wave packet at zc = 30 km. Using the polarization relations we set the following initial fields: ...(73) ...(74) ...(75) ...(76) For the wave lengths we set lx5lz51 km, defining the initial wavenumber k0 = (k0, m0) and the intrinsic frequency ... . The domain of the full model has the size (lx, lz) = (1 km, 60 km) at a resolution of nx × nz = 64 × 3840. In both test cases the domain is periodic in the horizontal and has a no-slip solid wall boundary condition at the bottom and top. A sponge layer at the top can be switched on to avoid spurious reflections of the gravity wave packet. In Rieper et al. (2012, manuscript submitted to J. Fluid Mech.) the extended, weakly nonlinear Wentzel-Kramers-Brillouin (WKB) theory by Achatz et al. (2010) is verified-that is, a WKB model and a full model5 are compared with respect to the predictions of the WKB theory. Here we recall the major results in short: there is very good agreement for the propagation of wave 1 (initial wavenumber k0). The fullmodel shows a wave 2 (wavenumber 2k0) oscillating about a smoothly increasing amplitude predicted by theWKB theory. For smaller l (i.e., for smaller e in the multiscale asymptotic ansatz), these oscillations vanish. Higher harmonics should not be inducedwithin the orders of magnitude considered in the WKB theory (Achatz et al. 2010) and the full model only produces higher harmonics of a distinctly smaller order of magnitude. In this text we analyze the behavior of full and WKB model with respect to numerical issues: influence of resolution, flux function, Poisson-solver tolerance, and the divergence-error correction introduced in section 3. Throughout these calculations we used a fixed (small) time step ?t = 1.0 s in order to fully resolve the buoyancy oscillations. The corresponding CFL numbers depend on the chosen resolution and are stated at due place. In Fig. 17 it is shown how the various harmonics evolve in time in the full model and according to the WKBtheory depending on the spatial resolution. Shown are the amplitude maxima of buoyancy for the wavenumbers ak0 with a = 2, 3, 4, and 5. The resolution nl is given with respect to the number of points per wavelength of the initial wave with k0. The time step of 1 s leads to a CFL number of 1/ 60 for the high-resolution case and 1/ 120 for the low-resolution case. For the coarse resolution with nl = 16-corresponding to a domain resolution of nx 3 nz = 16 3 960-the wave 2 amplitude falls behind the WKB-predicted value. On the other hand, the higher harmonics (a . 2), which should be absent according to WKB, are of comparable magnitude as wave 2. Obviously, a poor resolution leads to a feed in of energy from wave 2 to the higher harmonics. If the resolution is doubled (right of Fig. 17) the higher harmonics are clearly of a smaller order of magnitude. If the resolution is doubled again (not shown here) the higher harmonics practically vanish and the WKB predictions are well met. Note that the wave 1 amplitude of the full model (also not shown) always agrees very well with the WKB prediction. At the resolution nl = 32 the influence of the flux function (not shown here) is only marginal. We compare three different flux functions: central difference scheme (CDS), MUSCL, and ALDM produce a second harmonic of comparable magnitude. This amplitude is smoothly varying for MUSCL and ALDM but shows a rough evolution for CDS. The differences appear in the higher harmonics. CDS practically does not produce any higher harmonics. MUSCL has very small oscillations in the odd harmonics of O(1%) of the second harmonic, while ALDM has O(10%)-oscillations in the third harmonic, and O(1%)-oscillations in the harmonics higher than 3 for a Poisson-solver tolerance of e = 1029. The tolerance of the Poisson solver also influences the spectrum of the solution as can be seen in Fig. 18. At a tolerance of e5 1025 the full model produces spurious amplitudes in the higher harmonics: a.1. Only for small tolerances e , 1026 do these harmonics disappear. Note that wave 1 (not shown here) is in good agreement to the WKB predictions for tolerances as large as e = 1025. These calculations were done at a high resolution of nl 5 64 with a fixed time step of 1 s and a corresponding CFL number of 1/ 30. The tolerance of the Poisson solver can be relaxed if the divergence-error correction as described in section 3 is applied. If it is switched offwe obtain the result shown in Fig. 19. The higher harmonics do not agree withWKB predictions in contrast to the results shown on the right side of Fig. 18 with the divergence-error correction switched on. Obviously, even wave 1 (left) is completely wrong; see Fig. 19 left. These calculations were also done with nl = 64, a fixed time step of 1 s, and a corresponding CFL number of 1/ 30. d. Breaking of a 2D gravity wave packet With this test we evaluate how the implicit turbulence model handles the breaking of a gravity wave packet and how the induced turbulence is modeled qualitatively. Since this text focuses on numerical aspects we do not analyze the turbulence characteristics, which are part of future work. The wave packet is initialized as described in section 4 albeit with an amplitude having a Gaussian profile in two spatial directions with the initial buoyancy amplitude given by ...(77) with the wave lengths lx = lz = 4 km and a normalized amplitude factor of a = 0.9. The domain is 80 km 3 80 km and resolved with nx 3 nz = 640 3 640. The time step was fixed to 1 s leading to a CFL number of 1/ 6. The background is isothermal with T00 = 300 K. In Fig. 20 the isolines of the total potential temperature u5u1u0 are shown at the initial state (top left), after 10 min (top right), and in according order after 20, 30, 40, and 50 min. Already after 10 min the steepening of the gradients becomes visible. The first overturning of the gravity waves occurs after about 40 min. In this text we do not give an analysis of the turbulence modeling but we can see-at least-that the solution remains smooth and physical, grid-scale oscillations are avoided by the scheme. 5. Summary and conclusions We have presented a conservative way of discretizing Durran's pseudo-incompressible equations, which was-to the best of our knowledge-until now only done by Smolarkiewicz and Dörnbrack (2008) and by Smolarkiewicz and Szmelter (2011) on unstructured meshes. New is the implementation of a turbulence model right into the numerical flux function. The adaptive local deconvolution method (ALDM) by S. Hickel (2011, unpublished manuscript), Hickel et al. (2006, 2007), and Hickel and Larsson (2008) is for the first time used in the context of non-Boussinesq flows with a stratified background. We analyze numerical difficulties and propose ways around them: * To avoid spurious oscillations and to avoid the activation of the turbulence model on an unperturbed atmosphere, ALDM must only be applied to the deviation r0 5r2r from the background. * The physically motivated abort criterion proposed by Smolarkiewicz et al. (1997) for solving the elliptic Helmholtz equation in the context of the Euler equation is adapted to the pseudo-incompressible system. * To avoid oscillations in a uniform atmosphere at high altitude the Poisson equation has to be scaled with ru. This causes a height-independent divergence error in the velocity field, allowing us to reduce the Poissonsolver tolerance to acceptable values with respect to the overall computational effort. * A correction term in the effective continuity equation allows us to reduce the number of iterations of the Poisson solver even more. It reduces the effect of the residual divergence error present in the velocity field after the projection step. Note that these methods are applicable to any kind of finite-volume method and are not restricted to the ALDM discretization. Also note that the implicit turbulence scheme can be switched off. The resulting high-order central difference scheme could then be combined with an explicit turbulence model for testing and comparing purposes. In both cases-implicit and explicit-the modeling introduces unresolved subgrid- scale dissipation. The resolved dissipation on the grid scale is done by the discretization of the viscous stress tensor. Note that this does not lead to ''double counting'' since turbulence model and physical viscosity act on different scales. With the hot-cold bubble test case byRobert (1993)we validate the code for convective problems. The similarity between the solutions shows three things: * Euler and Durran's equation produce similar results if applied to convective problems. * Upwind schemes and ALDM produce similar results if applied to well-resolved laminar flows and ALDM has less numerical dissipation on coarse grids. * The solution only weakly depends on a uniform mean background wind. With the hot bubble test case by Mendez-Nunez and Carroll (1994) we demonstrate the relevance of applying ALDM only to the density perturbations and the capability of the model to cope with convection in the presence of a stratified background. The 1D gravity wave packet (GWP) test case is used to demonstrate the influence of the numerics on the spectral development of a GWP. We compare the results with those predicted by the extended WKB theory (Achatz et al. 2010): * If the resolution is too low, higher harmonics are excited-which is not of physical but only of numerical origin. * The tolerance of the Poisson solver has a direct impact on the spectrum, especially on the higher harmonics. * The introduction of a divergence-error correction improves the quality of the spectrum for lower tolerance of the Poisson solver. In a 2D GWP breaking test case we visualize how the ALDM turbulence model keeps a smooth solution without creating unphysical oscillations on the grid scale. Since this text focuses on numerical issues we do not go further into analyzing the turbulence characteristics, which is planned for future work. As turbulence in real GW breaking events is a threedimensional phenomenon, it remains open to show how a 3D GWP breaks and how the ALDM scheme behaves in such a more physical 3D context. Since calculations in 3D are time consuming, this part of the research has to be postponed until a parallel version of pincFloit is available. So far the model is restricted to the dry atmosphere. We plan to include moist processes in order to be able to study the interaction between gravity waves and clouds. With the implicit turbulence scheme the moist quantities would be transported as a passive tracer-an application for which ALDM has already been tested (Hickel et al. 2007). The latent heat term appears on the right-hand side of the divergence constraint. Acknowledgments. Wethank the Leibniz-Gemeinschaft(WGL) for support within their PAKT program. S.H. and U.A. thank Deutsche Forschungsgemeinschaftfor partial support through the MetStröm Priority Research Program (SPP 1276) and through Grants HI 1273/1-1 Ac71/4-1. 1 Please note that in the present paper the Froude number is defined as follows: ... with reference velocity uref and reference length lref. 2 In the following we only use the term ''density'' but always refer to the pseudo-incompressible (effective) density. 3 It is also shown that this problem can be circumvented by transporting only the perturbation of the potential temperature instead of the full quantity. 4 Note that each iteration within BiCGSTAB has two calls of the linear operator. 5 Since the flow remains laminar in the analyzed test case in that paper no turbulencemodel was needed and-without loss of generality-a MUSCL scheme was used. REFERENCES Achatz, U., R. Klein, and F. Senf, 2010: Gravity waves, scale asymptotics, and the pseudo-incompressible equations. J. Fluid Mech., 663, 120-147. Adams, N., S. Hickel, and S. Franz, 2004: Implicit subgrid-scale modeling by adaptive deconvolution. J. Comput. Phys., 200, 412-431, doi:10.1016/j.jcp.2004.04.010. Batchelor, G. K., 1953: The conditions for dynamical similarity of motions of a frictionless perfect-gas atmosphere. Quart. J. Roy. Meteor. Soc., 79, 224-235. Durran, D. R., 1989: Improving the anelastic approximation. J. Atmos. Sci., 46, 1453-1461. _____, 2008: A physically motivated approach for filtering acoustic waves fromthe equations governing compressible stratified flow. J. Fluid Mech., 601, 365-379, doi:10.1017/S0022112008000608. _____, and A. Arakawa, 2007: Generalizing the Boussinesq approximation to stratified compressible flow. C. R. Mec., 355, 655-664, doi:10.1016/j.crme.2007.08.010. Fritts, D. C., and M. J. Alexander, 2003: Gravity wave dynamics and effects in the middle atmosphere. Rev. Geophys., 41, 1003, doi:10.1029/2001RG000106. Gottlieb, S., and C.-W. Shu, 1998: Total variation diminishing Runge-Kutta schemes. Math. Comput., 67, 73-85, doi:10.1090/ S0025-5718-98-00913-2. Grabowski, W. W., and T. L. Clark, 1991: Cloud-environment interface instability: Rising thermal calculations in two spatial dimensions. J. Atmos. Sci., 48, 527-546. Hickel, S., 2008: Implicit turbulence modeling for large-eddy simulation. Ph.D. dissertation, Technische Universität München, 204 pp. _____, and J. Larsson, 2008: An adaptive local deconvolution model for compressible turbulence. Proc. 2008 Summer Program, Stanford, CA, Stanford University Center for Turbulence Research, 85-96. _____, N. A. Adams, and J. Domaradzki, 2006: An adaptive local deconvolution method for implicit LES. J. Comput. Phys., 213, 413-436, doi:10.1016/j.jcp.2005.08.017. _____,_____, and N. Mansour, 2007: Implicit subgrid-scale modeling for large-eddy simulation of passive-scalar mixing. Phys. Fluids, 19, 095102, doi:10.1063/1.2770522. _____, M. Mihatsch, and S. J. Schmidt, 2011: Implicit large eddy simulation of cavitation in micro channel flows. Proc. Third Int. Cavitation Forum, Warwick, United Kingdom, WIMRC. [Available online at https://files.warwick.ac.uk/shengcaili/ browse/Proceedings12011.] Klar, J.-U., C. Breitsamter, S. Hickel, and N. A. Adams, 2011: Integrated experimental-numerical analysis of high agility aircraftwake vortex evolution. J. Aircr., 48, 2050-2058. Klein, R., 2009: Asymptotics, structure, and integration of soundproof atmospheric flow equations. Theor. Comput. Fluid Dyn., 23, 161-195. _____, U. Achatz, D. Bresch, O. M. Knio, and P. K. Smolarkiewicz, 2010: Regime of validity of sound-proof atmospheric flow models. J. Atmos. Sci., 67, 3226-3237. Lindzen, R. S., 1981: Turbulence and stress owing to gravity wave and tidal breakdown. J. Geophys. Res., 86 (C10), 9707-9714. Lipps, F., 1990: On the anelastic approximation for deep convection. J. Atmos. Sci., 47, 1794-1798. _____, and R. Hemler, 1982: A scale analysis of deep moist convection and some related numerical calculations. J. Atmos. Sci., 39, 2192-2210. Meister, A., 1999: Numerik Linearer Gleichungssysteme: Eine Einführung in Moderne Verfahren (Numerical Methods for Linear Systems of Equations: An Introduction to Modern Methods). Vieweg, 222 pp. Mendez-Nunez, L. R., and J. J. Carroll, 1994: Application of the MacCormack scheme to atmospheric nonhydrostatic models. Mon. Wea. Rev., 122, 984-1000. Ogura, Y., and N. A. Phillips, 1962: A scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci., 19, 173-179. Remmler, S., and S. Hickel, 2012a: Direct and large-eddy simulation of stratified turbulence. Int. J. Heat Fluid Flow, 35, 13-24. _____, and _____, 2012b: Spectral structure of stratified turbulence: Direct numerical simulations and predictions by large eddy simulation. Theor. Comput. Fluid Dyn., doi:10.1007/S00162- 012-0259-9, in press. Robert, A., 1993: Bubble convection experiments with a semiimplicit formulation of the Euler equations. J. Atmos. Sci., 50, 1865-1873. Shu, C.-W., 1997: Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. ICASE Rep. 97-65, 78 pp. Smolarkiewicz, P. K., and L. G. Margolin, 1994: Variational solver for elliptic problems in atmospheric flows. Appl. Math. Comput. Sci., 4, 527-551. _____, and A. Dörnbrack, 2008: Conservative integrals of adiabatic Durran's equations. Int. J. Numer. Methods Fluids, 56, 1513-1519. _____, and J. Szmelter, 2011: A nonhydrostatic unstructured-mesh soundproof model for simulation of internal gravity waves. Acta Geophys., 59, 1109-1134, doi:10.2478/s11600-011-0043-z. _____, V. Grubisi c, and L. G. Margolin, 1997: On forward-in-time differencing for fluids: Stopping criteria for iterative solutions on anelastic pressure equations. Mon. Wea. Rev., 125, 647-654. Williamson, J., 1980: Low-storage Runge-Kutta schemes. J. Comput. Phys., 35, 48-56, doi:10.1016/0021-9991(80)90033-9. FELIX RIEPER Institut fur Atmospheuroare und Umwelt, Goethe-Universiteuroat Frankfurt, Frankfurt, Germany STEFAN HICKEL Lehrstuhl fur Aerodynamik und Streuroomungsmechanik, TU Meurounchen, Munich, Germany ULRICH ACHATZ Institut fur Atmospheuroare und Umwelt, Goethe-Universiteuroat Frankfurt, Frankfurt, Germany (Manuscript received 25 January 2012, in final form 17 July 2012) Corresponding author address: Felix Rieper, Deutscher Wetterdienst, Frankfurter Str. 135, D-63067 Offenbach, Germany. E-mail: [email protected] (ProQuest: Appendix omitted.) (c) 2013 American Meteorological Society [ Back To TMCnet.com's Homepage ]