@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Earth, Ocean and Atmospheric Sciences, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Wong, May Wai San"@en ; dcterms:issued "2014-04-22T15:05:55Z"@en, "2014"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description "Traditional semi-Lagrangian dynamical solvers are widely used in current global numerical weather prediction (NWP) and climate models, but are known to lack inherent mass-conserving properties. Some newer approaches utilize a cell-integrated (conservative) semi-Lagrangian (CISL) semi-implicit solver, which is inherently mass-conserving. However, existing CISL semi-implicit solvers lack consistent formulation among the discrete continuity equation and other discrete conservation equations for scalar tracers such as water vapour and air pollutants. Such inconsistency can lead to spurious generation or removal of scalar mass. In this dissertation, a new cell-integrated semi-Lagrangian (CISL) semi-implicit nonhydrostatic solver is presented with consistent discrete mass conservation equations for air and all tracers, and which preserves the shape of tracer-mass distribution. The discretization does not depend on a mean reference state, but maintains the same framework as typical semi-implicit CISL solvers, where a linear Helmholtz equation is constructed and a single application of the cell-integrated transport scheme is needed for scalar transport. Tests of this new solver are made for a series of increasingly complex flow scenarios. The initial testbed utilizes the hydrostatic, incompressible, shallow-water equations, for which the new solver is shown to be numerically stable. It maintains accuracy comparable to other existing solvers even for a highly nonlinear unstable jet. The second suite of tests are for nonhydrostatic two-dimensional (x-z) fully compressible flows in the atmosphere as governed by the moist Euler equations, which compare well to several idealized benchmark test cases from the literature. The third flow scenario is complex orography, where the nonhydrostatic equations are transformed to use a terrain-following height coordinate. Results from test cases of dry and moist flows over idealized mountain shapes are presented. In summary, the prototype development work presented in this dissertation shows that the proposed CISL nonhydrostatic solver with conservative and consistent transport may be a desirable candidate for a dynamical core in comprehensive global NWP and climate models."@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/46519?expand=metadata"@en ; skos:note "A fully-compressible nonhydrostatic cell-integrated semi-Lagrangianatmospheric solver with conservative and consistent transportbyMay Wai San WongB. Sc. Atmospheric Science, The University Of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Atmospheric Science)The University Of British Columbia(Vancouver)April 2014© May Wai San Wong, 2014AbstractTraditional semi-Lagrangian dynamical solvers are widely used in current global numerical weatherprediction (NWP) and climate models, but are known to lack inherent mass-conserving proper-ties. Some newer approaches utilize a cell-integrated (conservative) semi-Lagrangian (CISL) semi-implicit solver, which is inherently mass-conserving. However, existing CISL semi-implicit solverslack consistent formulation among the discrete continuity equation and other discrete conservationequations for scalar tracers such as water vapour and air pollutants. Such inconsistency can lead tospurious generation or removal of scalar mass.In this dissertation, a new cell-integrated semi-Lagrangian (CISL) semi-implicit nonhydrostaticsolver is presented with consistent discrete mass conservation equations for air and all tracers, andwhich preserves the shape of tracer-mass distribution. The discretization does not depend on a meanreference state, but maintains the same framework as typical semi-implicit CISL solvers, where alinear Helmholtz equation is constructed and a single application of the cell-integrated transportscheme is needed for scalar transport.Tests of this new solver are made for a series of increasingly complex flow scenarios. The initialtestbed utilizes the hydrostatic, incompressible, shallow-water equations, for which the new solveris shown to be numerically stable. It maintains accuracy comparable to other existing solvers evenfor a highly nonlinear unstable jet. The second suite of tests are for nonhydrostatic two-dimensional(x-z) fully compressible flows in the atmosphere as governed by the moist Euler equations, whichcompare well to several idealized benchmark test cases from the literature. The third flow scenario iscomplex orography, where the nonhydrostatic equations are transformed to use a terrain-followingheight coordinate. Results from test cases of dry and moist flows over idealized mountain shapesare presented. In summary, the prototype development work presented in this dissertation showsthat the proposed CISL nonhydrostatic solver with conservative and consistent transport may be adesirable candidate for a dynamical core in comprehensive global NWP and climate models.iiPrefaceThe main body of this dissertation is composed of work from two accepted journal papers (Chapters2 and 3) and one submitted (Chapter 4).Chapter 1I wrote this chapter with editing provided by Prof. Stull.Chapter 2Wong, M., W.C. Skamarock, P.H. Lauritzen, and R.B. Stull, 2013: A cell-integrated semi-Lagrangian semi-implicit shallow-water model (CSLAM-SW) with conservative and consistenttransport. Mon. Wea. Rev. 141:2545-2560• I designed the proposed formulation for consistent transport in a cell-integrated semi-Lagrangiansemi-implicit shallow-water model with Dr. W.C. Skamarock.• I analyzed the effectiveness of the consistent transport scheme in CSLAM-SW.• I provided most of the machinery and tools needed in the CSLAM-SW solver (includingdiscretizing the system, the elliptic solver, trajectory computations, and marriage of solvercomponents to the existing CSLAM transport scheme).• Dr. P.H. Lauritzen provided the CSLAM transport scheme.• Dr. J.B. Klemp provided the Eulerian solver that was used for comparison purposes (whichI modified from using a flux-form leapfrog explicit time-stepping scheme to an advectiveleapfrog semi-implicit scheme to provide a fair comparison with CSLAM-SW).• I wrote the manuscript with editing provided by the co-authors and four anonymous reviewers.I prepared most of the figures for this publication; Dr. P.H. Lauritzen provided Fig. 2.1 andDr. W.C. Skamarock provided Fig. 2.2.Chapter 3Wong M., W.C. Skamarock, P.H. Lauritzen, J.B. Klemp, and R.B. Stull, 2014: A compressiblenonhydrostatic cell-integrated semi-Lagrangian semi-implicit solver (CSLAM-NH) with consistentand conservative transport. Mon. Wea. Rev. doi:10.1175/MWR-D-13-00210.1, in press.iii• I extended the design of CSLAM-SW for a nonhydrostatic system to formulate CSLAM-NH.• I analyzed and compared the results from CSLAM-NH with an existing Eulerian nonhydro-static solver (see below).• As in Chapter 2, I provided most of the machinery needed in the nonhydrostatic solver (in-cluding the coupling of an existing microphysics scheme to the solver and other componentsthat I adapted from the code I developed in CSLAM-SW).• Dr. J.B. Klemp provided the Eulerian nonhydrostatic solver, used for comparison purposes,and the Kessler microphysics scheme.• I wrote the manuscript with editing provided by the co-authors and three reviewers. I pre-pared all the figures (except for Fig. 3.1, which is the same as Fig. 2.1 that was provided byDr. P.H. Lauritzen) for this publication, with labelling touch-ups by Dr. W.C. Skamarock.Chapter 4A version of Chapter 4 has been submitted for publication:Wong, M., W.C. Skamarock, P.H. Lauritzen, J.B. Klemp, and R.B. Stull. Testing of a cell-integrated semi-Lagrangian semi-implicit nonhydrostatic atmospheric solver (CSLAM-NH) withidealized orography. (Submitted, February 2014.)• I reformulated CSLAM-NH to use terrain-following height coordinates and an implicit Rayleighdamping layer.• I provided the modifications to CSLAM-NH for it to work with idealized terrain.• I designed CSLAM-NH to use an iterative scheme to handle the buoyancy terms implicitly,and in the testing process, I found that recomputation of the trajectories are needed to allowfor large time steps over terrain.• I analyzed and compared the results from CSLAM-NH with those from the literature, andin one case, I modified the Eulerian solver (provided by Dr. J.B. Klemp for Chapter 3) forcomparison with CSLAM-NH.• I wrote the manuscript and prepared all the figures for this manuscript, with editing providedby the co-authors.Chapter 5I drew all the conclusions and wrote this chapter, with editing provided by Prof. Stull.ivTable of contentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 The role of a dynamical core in an atmospheric model . . . . . . . . . . . . . . . . 11.2 Previous related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Eulerian vs. Lagrangian advection schemes . . . . . . . . . . . . . . . . . 21.2.2 Traditional semi-Lagrangian semi-implicit method . . . . . . . . . . . . . 41.2.3 Mass conservation in semi-Lagrangian schemes . . . . . . . . . . . . . . . 51.3 Dissertation contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.4 Dissertation outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 A two-dimensional shallow-water equations solver with conservative and consistenttransport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Mass conservation and consistency in SLSI solvers . . . . . . . . . . . . . . . . . 142.2.1 CSLAM — a CISL transport scheme . . . . . . . . . . . . . . . . . . . . 142.2.2 A discrete semi-implicit continuity equation in velocity-divergence formusing CSLAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.3 Numerical inconsistency in semi-implicit continuity equations in a velocity-divergence form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17v2.3 A consistent and mass-conserving semi-implicit SW solver . . . . . . . . . . . . . 182.4 Test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4.1 A radially-propagating gravity wave . . . . . . . . . . . . . . . . . . . . . 232.4.2 Bickley jet – Ro = 0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.4.3 Gaussian jet – Ro = 5.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.4.4 Gaussian jet – Ro = 5.0 with shape-preservation . . . . . . . . . . . . . . . 312.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353 Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver . . . . . . . 383.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.3 A consistent and mass-conserving nonhydrostatic solver . . . . . . . . . . . . . . 413.3.1 CSLAM — a conservative transport scheme . . . . . . . . . . . . . . . . . 413.3.2 Trajectory algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3.3 Discretization of the momentum equations . . . . . . . . . . . . . . . . . 433.3.4 Discretization of the thermodynamic equation . . . . . . . . . . . . . . . . 453.3.5 Helmholtz equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.6 Discretization of the continuity equation . . . . . . . . . . . . . . . . . . . 493.3.7 Discretization of moisture conservation equations . . . . . . . . . . . . . . 493.3.8 Diabatic processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.3.9 Diagnostic equation of state . . . . . . . . . . . . . . . . . . . . . . . . . 503.3.10 Consistency and shape-preservation . . . . . . . . . . . . . . . . . . . . . 513.4 Desirable properties of CSLAM-NH . . . . . . . . . . . . . . . . . . . . . . . . . 513.5 Idealized test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.5.1 Density current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.5.2 Gravity wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.5.3 2D (x-z) squall line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664 Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed verticalcoordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.2.2 CSLAM — a cell-integrated semi-Lagrangian transport scheme . . . . . . 72vi4.2.3 Discretized momentum equations . . . . . . . . . . . . . . . . . . . . . . 744.2.4 Conservative and consistent flux-form equations . . . . . . . . . . . . . . 764.2.5 Helmholtz equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.2.6 Iterative centered-implicit time-stepping scheme . . . . . . . . . . . . . . 784.2.7 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.2.8 Implicit Rayleigh damping . . . . . . . . . . . . . . . . . . . . . . . . . . 804.3 Idealized test cases: results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.3.1 Linear mountain waves over bell-shaped mountain . . . . . . . . . . . . . 814.3.2 Downslope windstorm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.3.3 Moist flow over a mountain in a nearly neutral environment . . . . . . . . 874.3.4 Computational performance (This section is not in the manuscript.) . . . . 934.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.1 Summary of contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.2 Impact on the broader community . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.3 Further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.3.1 Dimensionality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.3.2 Enhancement of CSLAM-NH stability . . . . . . . . . . . . . . . . . . . . 1065.3.3 Cost-effectiveness in a parallel environment . . . . . . . . . . . . . . . . . 1075.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116A.1 Numerical schemes for comparison . . . . . . . . . . . . . . . . . . . . . . . . . . 116A.1.1 A two-time-level traditional semi-Lagrangian semi-implicit model . . . . . 116A.1.2 An Eulerian leap-frog semi-implicit advective model . . . . . . . . . . . . 117viiList of tablesTable 3.1 Statistics for the density current simulations at time 15 min using CSLAM-NHat various grid resolutions and time step sizes. Comparison values from the non-hydrostatic x-z solver using SLICE in Melvin et al. (2010) are also presented.REFC25 are values taken from Straka et al. (1992). θ ′sampled are solutions sam-pled at 200 m for comparison with values in Straka et al. (1992). . . . . . . . . 56viiiList of figuresFigure 1.1 Numerical domain of dependence the upstream finite difference scheme [adaptedfrom Durran (2010), Fig. 2.1] . . . . . . . . . . . . . . . . . . . . . . . . . . 2Figure 1.2 Schematic diagram of a Lagrangian backward trajectory used in semi-Lagrangianadvection schemes [adapted from Durran (2010), Fig. 6.1]. The subscript A de-notes the arrival grid point, and subscript D denotes departure grid point. Blackdots indicate the Eulerian grid as in Fig. 1.1. . . . . . . . . . . . . . . . . . . 3Figure 1.3 Reconstruction of the subgrid-cell distribution using the piecewise quasi-biparabolicreconstruction method of Rancic (1992). The columns represent known Eule-rian cell-averaged values and the smooth surface represents the quasi-biparabolicreconstruction function (elevated in the figure for clarity). . . . . . . . . . . . 7Figure 2.1 (a) Exact departure cell area (δA∗, dark grey region) and the correspondingarrival grid cell (∆A, light grey region). (b) Departure cells in CSLAM (δA) arerepresented as polygons defined by the departure locations of the arrival gridcell vertices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.2 Definition of an Eulerian arrival grid cell, and its associated velocities at the cellfaces (ul,ur,vt ,vb) and cell corners (uc,vc)i for i=1, 2, 3, 4. . . . . . . . . . . . 17Figure 2.3 Comparison of the height field L2 error norms for the radially-propagating gravity-wave solutions. Errors are plotted at time T = 1× 105 s for the (a) linear(∆h = 10 m and h0 = 990 m) and (b) non-linear (∆h = 500 m and h0 = 1000 m)test cases computed on a 500 m mesh. Note the different scales in the plots. . . 24Figure 2.4 Specific concentration error (q− q0) in LKM for a divergent flow initializedwith a constant q0 = 1 in the (a) linear (∆h = 10 m and h0 = 990 m) and (b)nonlinear (∆h = 500 m and h0 = 1000 m) height perturbation cases. Note thedifferent scales in the plots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.5 Variation of specific concentration error (q− q0) (maximum absolute error,mean absolute error, and root-mean-squared error) with time step size in LKM(solid line) and CSLAM-SW (dashed line) for the (a) linear height perturbationand (b) nonlinear height perturbation cases. . . . . . . . . . . . . . . . . . . . 27ixFigure 2.6 Initial mean height h0 (top) and velocity v0 (bottom) profiles for the Bickley jet(∆h = 1 m, ∆v = 1 m s−1) and Gaussian jet (∆h = 50 m, ∆v = 56 m s−1). . . . 28Figure 2.7 Solutions of the Bickley jet at time T = 5×106 s (after 2500 time steps) for Ro= 0.1, Fr = 0.1 and Cradv = 0.2. Plotted are positive (solid line) and negative(dashed line) vorticity between −1×10−5 s−1 and 1×10−5 s−1 with a contourinterval of 5×10−7 s−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.8 Solutions of the Gaussian jet for Ro = 5.0 and Cradv = 0.56 at time T = 1.8×105s (after 1800 time steps). Plotted are positive (solid line) and negative (dashedline) vorticity between −5×10−4 s−1 and 5×10−4 s−1 with a contour intervalof 5×10−5 s−1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 2.9 CSLAM-SW solutions of the Gaussian jet for Ro = 5.0 at three different times(left to right on each row) of the simulation: at time T = 5×104 s, 1.0×105 s,and 1.4× 105 s. (a - c) Solutions using a Cradv of 0.56 (same simulation as inFig. 2.8) (d - f) Solutions using a larger Cradv of 2.5. Plotted are positive (solidline) and negative (dashed line) vorticity between −5×10−4 s−1 and 5×10−4s−1 with a contour interval of 5×10−5 s−1. . . . . . . . . . . . . . . . . . . . 33Figure 2.10 Specific concentration error (q− q0) in LKM for the Gaussian jet at time T =1.8×105 s, initialized with a constant q0 = 1 field. . . . . . . . . . . . . . . . 34Figure 2.11 Specific constituent concentration q at time T = 1.8× 105 s. Initial minimumand maximum q are 0.1 and 1.0 respectively. Regions with unphysical over-shooting (red) and undershooting (purple) are highlighted. . . . . . . . . . . . 37Figure 3.1 (a) Exact departure cell area (δA∗, dark grey region) and the correspondingarrival grid cell (∆A, light grey region). (b) Departure cells in CSLAM (δA) arerepresented as polygons defined by the departure locations of the arrival gridcell vertices. (Wong et al., 2013) . . . . . . . . . . . . . . . . . . . . . . . . . 41Figure 3.2 Geometric representation of the Lagrangian flux divergence, defined as the flux-area difference between the Eulerian arrival grid cell (solid square) and the de-parture cell (dashed polygon) in one time step. Velocities associated with theEulerian grid cell at the cell faces (ul,ur,wt ,wb) and cell vertices (uc,wc)i fori = 1,2,3,4 are also shown. White arrows indicate the computed trajectories ofeach departure grid cell vertex. . . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 3.3 Potential temperature perturbation (K) after 15 min. Contour intervals are every1 K, starting at 0.5 K. Mean wind U = 0 m s−1. . . . . . . . . . . . . . . . . . 55xFigure 3.4 Potential temperature perturbation (K) after 15 min. Contoured as in Fig. 3.3.Solution is translated using a mean wind U = 20 m s−1. . . . . . . . . . . . . 57Figure 3.5 Potential temperature perturbation (K) from CSLAM-NH after 15 min for gridspacing ∆x = ∆z = 100 m using time step sizes ∆t = 3 and 4 s. Mean windU = 0 m s−1. The Eulerian split-explicit scheme (not plotted) was numericallyunstable for these time steps, as it required ∆t ≤ 2 s for numerical stability ofthis gravity current. Contoured as in Fig. 3.3. . . . . . . . . . . . . . . . . . . 58Figure 3.6 Potential temperature perturbation (K) after 50 min. Contour intervals are every5×10−4 K (zero contour line not plotted). Solid lines indicate positive contoursand dashed lines indicate negative contours. Solution is translated using a meanwind U = 20 m s−1. Horizontal axis has also been translated with the meanwind so the line of symmetry remains at x = 0. . . . . . . . . . . . . . . . . . 59Figure 3.7 Potential temperature perturbation (K) solutions of the gravity wave case usingincreasingly large CSLAM-NH time steps (∆x = ∆z = 1 km) where (a)-(c) U =0 m s−1 and (d)-(f) U = 20 m s −1. Contoured as in Fig. 3.6. . . . . . . . . . . 60Figure 3.8 Vertical cross-sections of vertical velocity (color shading in m s−1) and solidcontour of the convective cloud structure (qc = 0.1 g kg−1) at times 0.6, 0.8,1.3, 1.7 h of the simulation for the 500 m grid-spacing runs with a time step of5.0 s from (left) CSLAM-NH, (middle) 5th-order split-explicit Eulerian model,and (right) 2nd-order split-explicit Eulerian model. . . . . . . . . . . . . . . . 63Figure 3.9 Moisture statistics including surface precipitation rate (kg s−1), accumulatedsurface precipitation (kg), and condensation rate (kg s−1) from the microphysicsusing CSLAM-NH, Eulerian 5th-order horizontal advection, and Eulerian 2nd-order horizontal advection at ∆x = ∆z = 500 m. . . . . . . . . . . . . . . . . . 64Figure 3.10 Updraft intensities using CSLAM-NH (red) and Eulerian 5th-order horizontaladvection (black) at ∆x = ∆z = 500 m and ∆t = 5 s (solid), and ∆x = ∆z = 250m and ∆t = 2.5 s (dashed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 3.11 Timing and intensity of the maximum vertical updraft using ∆x =∆z = 500 m atdifferent CSLAM-NH time step sizes (solid lines), as compared to the Eulerian5th-order horizontal advection vertical velocity (dashed lines). (Only first houris plotted.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66xiFigure 3.12 Mixing ratio errors (g kg −1) due to numerical inconsistency associated with(3.47). The passive tracer is initialized with a uniform mixing ratio field of1.0 g kg−1. The consistent formulation in CSLAM-NH (which does not use(3.47)) ensures mixing ratio constancy of the same passive tracer up to machineroundoff of order 10−14 (not shown). . . . . . . . . . . . . . . . . . . . . . . . 67Figure 4.1 Discrete departure cells in CSLAM-NH are approximated using straight edges(shaded in grey). The departure cell vertices (black circles) are computed usingbackward-in-time trajectories (arrows) from the vertices (white circles) of theEulerian arrival grid cell (white box). The CSLAM transport scheme is stableas long as the discrete departure grid cells are (a) non-self-intersecting, andbecomes problematic if (b) the departure cell self-intersects since the scheme isno longer mass-conserving. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.2 Analytical steady-state solution for the linear hydrostatic mountain waves ob-tained following Smith (1980). . . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 4.3 CSLAM-NH simulation of a linear hydrostatic wave (U/Na ≈ 0.1) for a lowwide mountain showing vertical velocity w (in m s−1) after T = 15000 s using(top) ∆t = 10 s (Cr = 0.2) and (bottom) ∆t = 75 s (Cr = 1.5). Contour intervalis 5×10−4 m s−1. Mean wind is from left to right. . . . . . . . . . . . . . . . 82Figure 4.4 Analytical steady-state solution for the linear nonhydrostatic mountain wavesobtained following Smith (1980). . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.5 CSLAM-NH simulation of a linear nonhydrostatic wave (U/Na = 1) for a nar-row mountain showing vertical velocity w (in m s−1) after T = 9000 s using(top) ∆t = 5 s (Cr = 0.125) and (bottom) ∆t = 60 s (Cr=1.5). Contour intervalis 6×10−4 m s−1. Mean wind is from left to right. . . . . . . . . . . . . . . . 83Figure 4.6 CSLAM-NH simulation for the 1972 Boulder windstorm case (a) potential tem-perature θ (K) (with contour interval of 8 K), (b) horizontal velocity U (m s−1)(with contour interval of 8 m s−1), and (c) Richardson number sgn(Ri)|Ri|0.5(−5 ≤ Ri ≤ 1 are plotted with contour interval of 0.5, following Doyle et al.(2000); grey contour shows negative values) at T= 3 h using a time step ∆t = 10 s. 85Figure 4.7 Time series of the simulated maximum CSLAM-NH downslope wind speeds (ms−1) for the 1972 Boulder windstorm case using different horizontal smoothingcoefficients, KD (m4 s−1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87xiiFigure 4.8 A self-intersecting departure cell (highlighted in red with vertices marked byblack circles) in CSLAM-NH when a large time step size of 25 s is used for thestrongly sheared flow in the 1972 Boulder downslope windstorm case. Blackcircles indicate departure grid cell vertices, and white circles the Eulerian arrivalgrid cell vertices. Arrows symbolize the computed backward-in-time trajecto-ries. Trajectories and the arrival grid cell associated with the self-intersectingdeparture cell are highlighted in red. . . . . . . . . . . . . . . . . . . . . . . . 88Figure 4.9 CSLAM-NH cloud-water mixing ratio (g kg−1) at time 5 h from an initiallysaturated nearly neutral flow (with an initial qc = 0.05 g kg−1) over a 700 mhill. White region above ground indicates subsaturated air (qc ≡ 0). . . . . . . 91Figure 4.10 (a) CSLAM-NH cloud-water mixing ratio qc (g kg−1) at time 5 h 10 mins froman initially saturated nearly neutral flow (with an initial qc = 0.05 g kg−1) overa 2 km mountain. White region above ground indicates subsaturated air (qc ≡0 g kg−1) (b) CSLAM-NH rain-water mixing ratio qr (g kg−1) at the samesimulation time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 4.11 Computational times for the moist neutral flow over a 2 km mountain usingdifferent CSLAM-NH time-step sizes (∆t). Times are normalized to that using∆t = 20 s for comparison purposes. The solid line indicates computational timesusing the iterative centered-implicit scheme. Dashed line indicates using thenon-iterative CSLAM-NH with ∆t = 5 s. . . . . . . . . . . . . . . . . . . . . 94Figure 4.12 Comparison of the computational times using CSLAM-NH and an Euleriansplit-explicit solver. (a) CPU times (in seconds) for the two solvers for a rangeof tracer counts. (b) Projected computational times using best linear fits. (c)Ratio of CPU times of CSLAM-NH to the Eulerian solver. . . . . . . . . . . . 96xiiiAcknowledgmentsThis dissertation would not be possible without the help of many people. First and foremost, I wouldlike to thank my research supervisor, Prof. Roland Stull, for letting me join his team in the earlydays as an undergraduate research student, for introducing me to scientific research and numericalweather prediction, and for providing invaluable advice and endless support through all these years.I am also indebted to my supervisory committee. I wish to thank Dr. William Skamarock for hisgenerous guidance and great mentorship throughout this dissertation. His contagious enthusiasmhas helped me bring this project to fruition, and, in a way, has shown me how to become a betterscientist. Dr. Peter Lauritzen provided expert advice on developing semi-Lagrangian schemes, andI would like to thank him for the great work experience. I wish to thank Dr. Henryk Modzelewskifor all that he has taught me since the very beginning, and Prof. Neil Balmforth for his time andguidance.I would also like to extend my gratitude to Dr. Joseph Klemp for his constructive feedback atvarious stages of the model development and manuscript preparation.It was truly an honour learning from and working with you all.Funding for this research was provided by Prof. Stull’s Discovery Grant from the CanadianNatural Sciences and Engineering Research Council and Prof. Stull’s BC Hydro contract. Muchof the research was done during summer visits to the National Center for Atmospheric Research inBoulder, Colorado, through the Advanced Study Graduate Visitor Program (2011, 2012) as well asa month-long visit co-sponsored by the Mesoscale and Microscale Meteorology Division and theClimate Global Division.I would like to thank all my teammates (past and present) on the Weather Research ForecastTeam and all the graduate students, post-doctorates, and other scientists I have met along the way(both in EOAS and at NCAR) for all the great memories and stimulating discussions. I thank myparents and my siblings for always supporting me, despite my being so far away from home for somany years. To all my dear friends, for always being there for me and for making me laugh on daysI most needed it. To all my running partners, for keeping me balanced. And last but not least, specialthanks to my very supportive fiance´, my best friend (and a wonderful cook), for always believing inme.xivDedicationDedicated to my mah mah (grandmother), Low LamxvChapter 1Introduction1.1 The role of a dynamical core in an atmospheric modelFluid motion in the atmosphere is described by a system of equations that includes the equations ofmotion, the continuity equation, the thermodynamic equation, water vapour conservation equation,and the ideal gas law. These equations depict how features in the atmosphere evolve and affectweather conditions at the earth’s surface. Analytical solutions are typically not available due tothe nonlinearity of the system, so numerical solutions are sought by discretizing the equations intime and space in numerical weather prediction (NWP) computer models. The output from thesecomputer models is a useful tool for meteorologists in the generation of weather forecasts.An NWP model can be broadly divided into two main components: the model dynamics andthe model physics. The “dynamics” include physical processes resolvable at the grid scale of anNWP model and are approximated using a fluid flow solver (also termed the ‘dynamical core’ ofthe model). These processes include wind evolution, and some of its driving forces such as horizon-tal and vertical advection of momentum, pressure-gradient force, Coriolis force, and gravitationalforce. The dynamical core is also responsible for the advection of all scalars (moisture, tracers, etc.)by the Reynolds-averaged winds, and these advective processes are carried out by the transportscheme of the fluid flow solver.The “physics” component refers to the subgrid-scale physics that are approximated using pa-rameterization schemes. These schemes parameterize physical processes unresolved at the grid cell,such as cloud microphysics, cumulus convection, turbulence, radiative fluxes, and modification offluxes due to land-surface characteristics. The validity of these physical parameterization schemes,as well as how they are coupled to the dynamics, are ongoing research topics that will not be ad-dressed here. The work presented in this dissertation will focus on the dynamics component of anNWP model, in particular on the development of a prototype solver serving as a viable alternativedynamical core for an NWP model.1Chapter 1: Introduction1.2 Previous related work1.2.1 Eulerian vs. Lagrangian advection schemesThe advection (or transport) scheme of a fluid flow solver is responsible for computing the timetendencies due to advection. A simplified one-dimensional passive advection equation for mixingratio q in a flow with constant velocity U is given below:∂q∂ t +U∂q∂x = 0. (1.1)As an example of an Eulerian advection scheme, the upstream method is used to discretize (1.1):qt+∆ti −qti∆t+Uqti−qti−1∆x= 0, (1.2)where ∆t and ∆x are the time step and grid spacing, respectively. The superscripts represent thetime level at which the terms are evaluated, and the subscripts denote the position of the grid point.Assuming that values for q from the previous time step are known, we can solve for qt+∆ti .The stencil of the scheme (1.2) is shown as grey squares, and the numerical domain of depen-dence is shaded in light grey in Fig. 1.1. If the advection velocity U translates the parcel beyondthis stencil, the scheme violates a numerical stability condition called the Courant-Friedrichs-Lewy(CFL) condition. Qualitatively, the CFL condition “requires that the numerical domain of depen-dence of a finite-difference scheme include the domain of dependence of the associated partialdifferential equation” (Durran, 2010). In an Eulerian scheme such as (1.2), the discrete advectionscheme poses a limit on the maximum time step by the CFL stability condition.qi∆xqi-1qi-2xttt+∆t qiqi-1qi-2∆tFigure 1.1: Numerical domain of dependence the upstream finite difference scheme [adaptedfrom Durran (2010), Fig. 2.1]2Chapter 1: Introduction∆xxttt+∆t qA∆tqDU∆tFigure 1.2: Schematic diagram of a Lagrangian backward trajectory used in semi-Lagrangianadvection schemes [adapted from Durran (2010), Fig. 6.1]. The subscript A denotes thearrival grid point, and subscript D denotes departure grid point. Black dots indicate theEulerian grid as in Fig. 1.1.The simplified passive advection equation in (1.1) can be rewritten as the Lagrangian derivative(also known as the material derivative or the total derivative),dqdt= 0, (1.3)which can be solved numerically usingqt+∆tA −qtD∆t= 0, (1.4)where subscript A denotes the arrival position and subscript D denotes the departure location, cal-culated using the wind speed U . In other words, the mixing ratio q is conserved along the trajectory(arrow connecting from qD to qA in Fig. 1.2). The Lagrangian transport method allows the elimina-tion of the advection terms by numerically solving for the total derivative. Moreover, the departurelocation defines the domain of dependence of the continuous problem. Since this point is calcu-lated explicitly, it will always be included in the numerical domain of dependence and therefore, theLagrangian method is not subject to the CFL condition.Although fully Lagrangian schemes allow for larger time steps than Eulerian schemes, theyhave the disadvantage that the parcels will eventually become unevenly distributed in deformedflows. Features in parts of the domain where parcels are sparse will become under-resolved. Tomaintain an even resolution, traditional semi-Lagrangian schemes were devised. Semi-Lagrangianschemes use a combination of two reference frames: Eulerian and Lagrangian. For schemes that3Chapter 1: Introductionuse backward trajectories, a new set of parcels is traced back a time step. The arrival grid pointsare chosen to be on a pre-defined regular Eulerian grid (Fig. 1.2). However, the departure locationtypically does not coincide with the Eulerian grid (such as the departure point in Fig. 1.2). Sincethe solution from the previous time step is only known on the Eulerian grid, some interpolation isneeded to obtain the value at the departure point. Staniforth and Coˆte´ (1991) provides a review ontraditional semi-Lagrangian schemes.1.2.2 Traditional semi-Lagrangian semi-implicit methodSemi-Lagrangian schemes are the most widely used space-time numerical integration scheme incurrent global NWP and climate models. The popularity of these schemes can be mainly attributedto their lenient numerical stability condition for advection. Larger stable advection time steps helpmake these numerical models more computationally efficient. The gain in efficiency can then beused to offset any increase in computational time needed in making longer forecasts, using moresophisticated physics schemes, and/or making higher resolution forecasts. In general, terrain andland surface features are better resolved using higher model resolution.The early work of semi-Lagrangian transport schemes was pioneered by Fjørtoft (1952, 1955),Wiin-Nielsen (1959), Krishnamurti (1962) and Sawyer (1963) for a more accurate representationof advection. Robert (1969) proposed combining an implicit treatment of the terms responsible forfast propagating waves in the atmosphere with an explicit treatment of the advection terms. Byimplicitly solving the gravity and acoustic wave terms, larger stable time step sizes are allowedas compared to using fully explicit schemes. Later, Robert (1981, 1985) suggested combining thetraditional semi-Lagrangian method for the advection terms with the implicit integration of the fastpropagating terms, which combines the extended stability of both schemes. This method has sincebeen widely adopted at numerous weather research centres.At the time of writing this dissertation, the following weather research centres are using semi-Lagrangian transport schemes in their numerical models: the Canadian Meteorological Centre (Coˆte´et al., 1998), the European Centre of Medium-Range Weather Forecasts (ECMWF) (Ritchie et al.,1986; Ritchie et al. 1995), the United Kingdom Meteorological Office (Davies et al., 2005), the Dan-ish Meteorological Institute (Nair and Machenhauer, 2002), the China Meteorological Association(Yang et al., 2008), Me´te´o-France (Bubnova´ et al., 1995), Hydrometcentre of Russia (RosHydromet,2011), and the Japan Meteorological Agency (JMA, 2013).As mentioned, traditional semi-Lagrangian schemes often require some grid-point interpolationto retrieve the value at a departure point. A widely known caveat of traditional semi-Lagrangianschemes is that it does not ensure mass conservation due to grid-point interpolation. The distri-bution of atmospheric pressure is the main driver for the flow features we see in our weather. As4Chapter 1: Introductionatmospheric pressure is the weight of air mass above, errors in the transport of the latter will likelyhave an impact on the prediction of meteorological features. Hence, mass-conserving advectionschemes are desirable, which is the motivation for this research.The following section describes some of the past modelling efforts in ensuring mass conserva-tion in models.1.2.3 Mass conservation in semi-Lagrangian schemesMass fixersNot all transport schemes inherently conserve mass, but this property can be enforced using a priorior a posteriori constraints (Rasch and Williamson, 1990a). A priori mass-conserving schemes areoften formulated with the desirable property inherent in the design of the numerical method, whereasa posteriori schemes typically attain mass conservation through an adjustment procedure, and areknown as ‘fixers’, following Rasch and Williamson (1990a).Different a posteriori adjustment procedures have been devised to overcome the lack of massconservation in semi-Lagrangian schemes. Priestley (1993) developed a mass restoration schemethat minimizes the local mass change in one time step by constraining the interpolation coefficientsin the solution of a monotonic semi-Lagrangian transport scheme. Moorthi et al. (1995) appliedan adjustment formula such that the entire model surface pressure field (equivalent to mass) ismultiplied by a constant factor that preserves the global mean surface pressure at its initial value.Williamson and Rasch (1994) enforced a conservation fixing step after the advection calculation,and restores global mass conservation selectively over regions more strongly affected by the advec-tion scheme. These schemes, albeit carried out in a rather ad-hoc manner, can all maintain globalmass conservation.Inherently conservative transportRecently, stricter approaches of ensuring mass conservation in semi-Lagrangian schemes have beendevised. These schemes are referred to as conservative semi-Lagrangian transport schemes, alsoknown as cell-integrated semi-Lagrangian (CISL) schemes (Rancic, 1992; Nair and Machenhauer,2002; Zerroukat et al., 2002; Lauritzen et al., 2010). In essence, these schemes are finite-volumemethods that utilize the fact that the integral mass over a control volume is conserved along itstrajectory. The problem then reduces to a remapping problem from an Eulerian arrival grid to aLagrangian departure grid (for a backward-in-time scheme). The three main steps to a CISL schemeare: (i) departure cell approximation, (ii) subgrid-cell reconstruction of the field on the Eulerian grid,5Chapter 1: Introductionand (iii) remapping of the field onto the Lagrangian grid. These are explained as follows.(i) There are multiple ways to define the discrete Lagrangian departure cell over which theintegration of the mass field is carried out. A rigorous approach is to trace the vertices of theEulerian arrival grid cell back to their departure locations, which define the corners of the departurecell (e.g. Rancic, 1992; Lauritzen et al., 2010). The departure cell edges can then be approximatedas straight lines, forming irregular quadrilaterals. This set of irregular quadrilateral cells form agrid that typically does not coincide with the Eulerian arrival grid (similar idea to that shown fortraditional schemes in Fig. 1.2).Another approximation is to compute the mass over an approximated area of the discrete depar-ture cell. The area of the discrete departure cell is approximated with edges parallel to the coordinateaxes to simplify the remapping procedure (e.g. Nair and Machenhauer, 2002). The horizontal posi-tions of the vertical cell edges coincide with the midpoints between the computed departure pointsof the arrival grid cell vertices. The vertical positions of the horizontal cell edges are defined bythose of the departure cell vertices. This approximation leads to a computational departure cell withmore sides (generally eight, but possibly fewer in some special circumstances), but the departurecell is now made up of rectangular areas that align with the Eulerian grid.(ii) Since the Lagrangian departure grid cells typically do not coincide with the Eulerian gridcells, a reconstruction of the subgrid-cell distribution of the transported variable is needed. Asubgrid-cell reconstruction function provides a best fit continuous function based on known Eule-rian cell-averaged values. Fig. 1.3 shows an example of a quasi-biparabolic fit (the smooth surface,elevated in the figure for clarity) of the Eulerian cell-averaged values (columns). The function isthen integrated over the Lagrangian departure cell (as opposed to using interpolation as in traditionalsemi-Lagrangian schemes).Various subgrid-cell reconstruction methods have been reported in the literature: piecewise con-stant method, piecewise linear method, piecewise parabolic method, and even higher-order methodssuch as the piecewise cubic method and piecewise quartic method. The lower the order of accuracyof the method, the more inherently damping the scheme. However, more reconstruction coefficientsare required in higher-order schemes, which can become expensive to compute. To reduce the costassociated with the computation of the coefficients, Nair and Machenhauer (2002) introduced anextension of the quasi-biparabolic method of Rancic (1992) (shown in Fig. 1.3) and used only fiveof the nine coefficients in a fully two-dimensional piecewise biparabolic fit. The method has beenshown to still give accurate representation despite the simplification, and also has the benefit of astraightforward implementation of shape-preserving filters (next section).(iii) Given the departure cell position and the reconstruction function, the final step is to carry outthe remapping of the transported variable from the Eulerian arrival grid to the Lagrangian departure6Chapter 1: IntroductionFigure 1.3: Reconstruction of the subgrid-cell distribution using the piecewise quasi-biparabolic reconstruction method of Rancic (1992). The columns represent knownEulerian cell-averaged values and the smooth surface represents the quasi-biparabolicreconstruction function (elevated in the figure for clarity).grid. One approach is to do the remapping of the mass field from one grid to another using acascade approach, i.e. the integral of the reconstruction function over the Lagrangian departure cellis computed incrementally one dimension at a time, as in Rancic (1992), Nair and Machenhauer(2002), and Zerroukat et al. (2002). Others have shown that a more fully 2D remapping procedurecan also be done, as in Lauritzen et al. (2010), who used the approach in Dukowicz and Kodis (1986)and simplified the 2D remapping procedure by converting the areal integral into a line integral usingthe Gauss-Green theorem. As long as the integral is exact and there is no gap or overlap betweenthe discrete departure cells, local (and therefore, global) mass conservation is ensured (Lauritzenet al., 2011).Shape-preserving transportRelated to mass conservation is the desirable property of shape-preservation in tracers. Schemesthat ensure shape-preservation prevent new minima and/or maxima from forming in the advectedfield. In the broader literature, they are sometimes referred as monotonic limiters, essentially non-oscillatory limiters, or positive-definite limiters for those that eliminate negative values. In atmo-spheric science, shape-preserving filters are often applied in moisture and tracer transport as theminimum and maximum specific concentrations are physically bounded, e.g. the lower end must be7Chapter 1: Introductionnon-negative and the upper end is bounded by the initial maximum value for pure transport cases.Spurious minimum and maximum can be generated through an interpolation or remapping proce-dure, and the need of shape-preserving limiters is dependent on the advection scheme operators.Shape-preserving schemes can be applied either a priori or a posteriori. A priori enforcementof the shape-preserving property can be applied for example within the interpolation procedure byconstraining the derivative estimates in an interpolation scheme (Rasch and Williamson, 1990b;Bermejo and Staniforth, 1991). A posteriori fixers reclaim the initial shape of the field by apply-ing adjustments afterwards. For example, in the shape-preserving fixer described in Williamsonet al. (1987), local restoration of moisture at negative specific humidity points can be made by in-creasing the values to zero there. But to ensure global moisture conservation, specific humiditiesat the neighbouring grid points are then reduced accordingly. If there is insufficient moisture in thelocal neighbourhood, then the global moisture field in the domain is reduced by a constant factorto compensate for the moisture restoration at the negative points (and to preserve total moisturecontent).In CISL schemes, shape-preserving filters can be applied within the subgrid-cell reconstruction.Since these schemes advect flux-form variables such as moisture mass, care should be taken to onlyapply the shape-preserving limiter on the specific humidity. Often, new extrema in the moisturemass are physically acceptable, e.g. in a divergent flow where there is convergence/divergenceof mass. Nair and Lauritzen (2010) show one way of implementing an a priori shape-preservinglimiter on only the specific concentration in the subgrid-cell reconstruction of the tracer mass.Consistent transport between air density and tracersAs mentioned earlier, the maximum time step of an explicit time-stepping scheme with semi-Lagrangian transport is restricted by fast-propagating gravity and acoustic waves, which led to thedevelopment of traditional semi-Lagrangian semi-implicit schemes. Following the devising of CISLtransport schemes, Machenhauer and Olk (1997) demonstrated for a one-dimensional shallow-watersystem that their CISL advection scheme, like traditional semi-Lagrangian schemes, can also be im-plemented with a semi-implicit time-integration scheme. Lauritzen et al. (2006) later derived andextended their discrete CISL continuity equation for a two-dimensional shallow-water system. Sim-ilar effort has also been carried out by Zerroukat et al. (2009) using the inherently mass-conservingsemi-Lagrangian transport scheme SLICE of Zerroukat et al. (2002).In existing semi-implicit CISL continuity equations, a divergence term is implicitly coupled tothe pressure-gradient terms in the momentum equations. This divergence term will be referred toas the semi-implicit correction to the explicit transport (of air density). Unlike air density however,atmospheric tracers (including moisture) are not physically coupled to the wind field, and their8Chapter 1: Introductionmasses are passively (and explicitly) transported by the advection scheme in the solver. As a result,many models use a different numerical method for tracer transport (e.g. explicit transport) than forsolving the continuity equation for air mass (Lauritzen et al., 2011).If the discrete tracer conservation equation is numerically consistent with the continuity equa-tion, the former should be identical to the continuity equation for a constant tracer mixing ratio ofone. As discussed in Lauritzen et al. (2011), air mass and tracer mass may be conserved individuallyusing a mass-conserving transport scheme (enforced either a priori or a posteriori), but numericalerrors due to inconsistent transport schemes may appear when the mixing ratio is diagnosed fromthe prognostic air and tracer mass variables. These numerical errors will lead to spurious gener-ation or loss of tracer mass, despite using an inherently conserving transport scheme. The needfor numerically consistent transport of tracers and air density motivated the first step of developinga semi-implicit CISL shallow-water equations solver in Lauritzen et al. (2006) in the first place.However, they had found that the use of a mean reference state in the semi-implicit correction termrendered it difficult to ensure numerical consistency in the transport.1.3 Dissertation contributionsIn this dissertation, a new cell-integrated semi-Lagrangian (CISL) semi-implicit nonhydrostaticsolver with conservative and consistent transport of air and tracers is presented. The achievementsencompassed by this dissertation are fourfold:1. A fully-compressible nonhydrostatic prototype solver is developed, in anticipation of the ef-fective resolution of current global NWP and climate models extending into the nonhydro-static regime in the foreseeable future.2. To ensure conservative transport, a CISL transport scheme is chosen, and the use of thisCISL transport scheme in solving the fully-compressible nonhydrostatic atmospheric systemis demonstrated for the first time.3. To ensure consistent transport for all scalars, a new discrete semi-implicit CISL continuityequation is devised.4. To the best of my knowledge, at the time of writing this dissertation, coupling of a CISL non-hydrostatic atmospheric solver to diabatic effects from a simple cloud microphysics schemeis published for the first time.The design of the proposed nonhydrostatic solver considers the following. The goal is to furtherimprove the numerical consistency in existing semi-Lagrangian transport methods used for global9Chapter 1: IntroductionNWP and climate modelling purposes. An alternative method in ensuring consistency betweenair density and scalar transport is devised. Other aspects of the model utilize existing numericalmethods, such as the traditional semi-Lagrangian scheme for momentum transport and a conjugate-residual solver as the elliptic solver for the semi-implicit time-integration scheme. As mentioned inthe previous section, traditional semi-Lagrangian semi-implicit schemes are widely used in currentoperational global NWP models, such as the Canadian Global Environmental Multiscale (GEM)model and the ECMWF Integrated Forecasting System (IFS).To demonstrate that the CISL semi-implicit nonhydrostatic solver is a feasible option for imple-mentation in a more comprehensive NWP model, a two-dimensional x-z prototype solver is devel-oped and applied to a suite of idealized test cases. The solver is tested using simplified boundaryconditions: periodic in x and rigid top and bottom boundaries. The Arakawa C grid (Arakawa andLamb, 1977) is used for a more accurate representation of gradient and divergence terms, whichhave generally been regarded as the dominant terms in fine-scale weather.To ensure conservative properties for dry air mass, tracer mass, and potential temperature, thecontinuity equation, tracer mass conservation equation, and the thermodynamic equation are allcast in the flux-form. The approach is based on the Advanced Research Weather Research andForecasting (WRF-ARW) model (Skamarock and Klemp, 2008), and is also used in the recentlydeveloped variable-resolution Model Prediction Across Scales-Atmosphere model (MPAS-A) (Ska-marock et al., 2012). All explicit scalar advection is carried out by the conservative semi-Lagrangiantransport scheme. A new semi-implicit discretization scheme for the continuity equation is devisedand used to maintain numerical consistency among the CISL scalar conservation equations.1.4 Dissertation outlineThe model development work is documented as three main stages (in the following three chapters):1. A two-dimensional shallow-water equations (SWE) solver with conservative and consistenttransport (Chapter 2).Overview: Shallow-water systems are often used as initial testbeds because of their resem-blance to the fully-compressible Euler equations. The SWE represent the horizontal dynamicsin the three-dimensional hydrostatic primitive equations, i.e. they permit Rossby and inertia-gravity waves, while excluding acoustic waves. The new consistent semi-implicit CISL con-tinuity equation is tested first on the two-dimensional shallow-water equations. The goal isto test that the new SWE solver is stable and accurate for this simplified system. The relativeimportance of numerical consistency in tracer mass conservation is illustrated using compar-isons with the discrete semi-implicit CISL continuity equation presented in Lauritzen et al.10Chapter 1: Introduction(2006). In addition, this stage validates the coding of the solver components, such as theLagrangian trajectory computations, Lagrange polynomial interpolation schemes, the ellipticsolver, and the marriage of existing conservative semi-Lagrangian transport subroutines to theSWE solver.2. Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver with conservativeand consistent transport (Chapter 3).Overview: The fully-compressible nonhydrostatic (Euler) equations allow a wide range ofatmospheric motions. In addition to gravity waves, compressibility of the atmosphere allowsfast acoustic waves to propagate. To test the new solver in these motions, the SWE solveris extended for the two-dimensional vertical-slice nonhydrostatic equations. The nonhydro-static equations require solving the vertical momentum equation and accounting for buoyancyeffects, which may impose a constraint on the numerical stability of the scheme. The semi-implicit nonhydrostatic solver is constructed using an implicit treatment of the fast acousticwaves, the proposed semi-implicit CISL continuity equation, and consistent formulations forall other scalar transport. The goal is to assess this solver in a set of idealized benchmark testscommonly used for nonhydrostatic atmospheric flows. Subgrid-scale latent heat processes,such as evaporation and condensation, are also included in the nonhydrostatic solver using asimple cloud microphysics scheme. The importance of numerical consistency is shown usinga 2D idealized convective storm simulation.3. Extension to a two-dimensional (x-z) nonhydrostatic solver with a transformed vertical coor-dinate (Chapter 4).Overview: Idealized mountain flows are simulated using the proposed solver. This stage ofthe development validates that the new CISL nonhydrostatic solver works for more realisticbottom boundary conditions, and demonstrates that only a little modification to the mechanicsof the nonhydrostatic solver is required for the change in the vertical coordinate system. Thesolver is applied for a suite of linear and nonlinear mountain flows, including an orographiccloud formation simulation.A summary of the dissertation, the broader impact of the effort, and an outline of future researchdirections for this work to be extended into a more comprehensive weather and climate model aregiven in the Conclusions (Chapter 5).11Chapter 2A two-dimensional shallow-waterequations solver with conservative andconsistent transport2.1 IntroductionSemi-Lagrangian semi-implicit (SLSI) schemes have been widely used in climate and numericalweather prediction (NWP) models since the pioneering work of Robert (1981) and Robert et al.(1985). The more lenient numerical stability condition in these schemes allows larger time stepsand thus increased computational efficiency. Traditional semi-Lagrangian schemes are not inher-ently mass-conserving due to their use of grid-point interpolation, and the lack of conservation canlead to accumulation of significant solution errors (Rasch and Williamson, 1990a; Machenhauer andOlk, 1997). To address this issue, conservative semi-Lagrangian schemes, also called cell-integratedsemi-Lagrangian (CISL) transport schemes (Rancic, 1992; Laprise and Plante, 1995; Machenhauerand Olk, 1997; Zerroukat et al., 2002; Nair and Machenhauer, 2002; Lauritzen et al., 2010), havebeen developed. Although CISL transport schemes allow for locally (and thus globally) conserva-tive transport of total fluid mass and constituent (i.e. tracer) mass, an issue related to conservationremains when they are applied in fluid flow solvers: the lack of consistency between the numeri-cal representation of the total mass continuity and constituent mass conservation equations (Jo¨ckelet al., 2001; Zhang et al., 2008). The lack of numerical consistency between the two can lead to theunphysical generation or removal of model constituent mass, which can introduce significant errorsin applications such as chemical tracer transport (Machenhauer et al., 2009).Our testbed for developing and testing CISL-based fluid flow solvers are the shallow-water (SW)equations on an f -plane:∂u∂ t +u∂u∂x + v∂u∂y − f v−g′ ∂h∂x = 0, (2.1)12Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport∂v∂ t +u∂v∂x + v∂v∂y + f u−g′ ∂h∂y = 0, (2.2)∂h∂ t +∇ · (hv) = 0, (2.3)∂ (hq)∂ t +∇ · (hqv) = 0 (2.4)where v = (u,v) is the horizontal velocity vector, f is the Coriolis parameter, g′ is the reducedgravity (defined as g′ = g∆ρ/ρ , where g is gravitational acceleration, ∆ρ is the density differencebetween the two layers, and ρ is the mean density), h is the total fluid depth (a surrogate for totalfluid mass), and hq is the depth portion (mass fraction) of an arbitrary constituent, where q is itsspecific concentration. Numerical consistency is satisifed if, for q0 = 1, the discretization schemeof the constituent equation (2.4) collapses to that for the continuity equation (2.3), also known asfree-stream preservation.The difficulty in maintaining consistency, as will be discussed in more detail, can partly be at-tributed to the conventional linearization around a constant mean reference state in the semi-implicitform of a CISL continuity equation. To eliminate the reference state, Thuburn (2008) developed afully-implicit CISL-based scheme for the shallow-water equations that requires solving a nonlin-ear Helmholtz equation at every time step. The solution of the Helmholtz equation is potentiallyproblematic and expensive (Thuburn et al., 2010). To reduce the dependence of their semi-implicitscheme on a reference state, Thuburn et al. (2010) used an alternative iterative approach to solve thenonlinear system, but it requires multiple calls to a Helmholtz solver per time step, again makingthe scheme potentially expensive.In addition to consistency and mass conservation, another desirable property is that the newscheme should be shape-preserving. A shape-preserving scheme ensures that no new unphysicalextrema are generated in a field due to the numerical scheme (e.g. Machenhauer et al., 2009). Forexample, specific concentrations of a passive constituent should not go outside the range of its initialminimum and maximum values. Non-shape-preserving schemes may generate unphysical specificconcentrations, such as negative concentration values due to undershooting.In this paper, using a shallow-water system, we present a new SLSI formulation that uses aCISL scheme for mass conservation and ensures numerical consistency between the total massand constituent-mass fields. The new scheme is based on the CISL transport scheme called theConservative Semi-LAgrangian Multi-tracer (CSLAM) transport scheme developed by Lauritzenet al. (2010). Like other typical conservative SLSI solvers, the algorithm requires a single linearHelmholtz equation solution and a single application of CSLAM. To ensure shape-preservation, thescheme is further extended to use existing shape-preserving filters.13Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport(b)(a)Figure 2.1: (a) Exact departure cell area (δA∗, dark grey region) and the corresponding arrivalgrid cell (∆A, light grey region). (b) Departure cells in CSLAM (δA) are represented aspolygons defined by the departure locations of the arrival grid cell vertices.The paper is organized as follows. In section 2.2, the conservative semi-Lagrangian schemeCSLAM is described and a discussion of the issue of consistency between total-mass and constituent-mass conservation in its semi-implicit formulation is provided. A new consistent semi-implicit dis-cretization of the CSLAM continuity equation, including the implementation of the shape-preservingschemes, is proposed in section 2.3. Results from four test cases are presented in section 2.4, high-lighting the stability and accuracy of the new scheme for linear and highly-nonlinear flows, as wellas showing the shape-preserving ability of the scheme. And finally, in section 2.5, a summary of theresults and a potential extension of the new scheme are given.2.2 Mass conservation and consistency in SLSI solvers2.2.1 CSLAM — a CISL transport schemeThe CSLAM transport scheme is a backward-in-time CISL scheme, where the departure grid cellarea δA∗ is found by tracing the regular arrival grid cell area ∆A back in time one time-step ∆t (Fig.2.1a). The CSLAM discretization scheme for (2.3) is given byhn+1exp ∆A = hn∗δA∗, (2.5)where the superscript denotes the time level, hn+1exp is the explicit cell-averaged height solution com-puted by integrating the height field hn over δA∗, which gives departure cell-averaged height value14Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transporthn∗. The departure cell area δA∗ in CSLAM is found through iterative trajectory computations fromthe four vertices of an arrival grid cell (unfilled circles in Fig. 2.1b) to their departure points (filledcircles in Fig. 2.1b). The departure cell area is then approximated using straight lines as cell edges(dark grey region δA in Fig. 2.1b). To integrate the height field over δA, CSLAM implements aremapping algorithm that consists of a piecewise biparabolic subgrid-cell reconstruction of the hnfield, and then the integration of the reconstruction function over the departure cell area. The areaintegration in CSLAM is transformed into a series of line integrals using the Gauss-Green theorem,and involves solving for a set of weights that depends only on the departure cell boundary. The useof line integrals greatly enhances the transport scheme’s computational efficiency for multi-tracertransport as the weights can be reused for all tracer species in the model. For full details on thetransport scheme, see Lauritzen et al. (2010).2.2.2 A discrete semi-implicit continuity equation in velocity-divergence form usingCSLAMLauritzen et al. (2006) (which we will refer to as LKM) developed an SLSI SW equations solverusing the explicit CISL transport scheme of Nair and Machenhauer (2002). For the momentumequations (2.1) and (2.2), they used a traditional SLSI discretization [(A.1) and (A.2) in the Ap-pendix but without time-off-centering]. Their momentum equations are then implicitly coupled toa velocity divergence correction term in the continuity equation. In this paper we follow the con-struction of the SW equations solver described in LKM, but we use CSLAM as the explicit CISLtransport scheme. The discrete semi-implicit CISL continuity equation given in LKM [eq.(31) inLKM] ishn+1 = hn+1exp −∆t2H0[∇eul ·vn+1−∇lag · v˜n+1]+∆t2H0[∇eul ·vn−∇lag ·vn]δA∗∆A, (2.6)where hn+1exp is as described above, ∆t is the model time step, H0 is the constant mean referenceheight, vn+1 is the velocity field implicitly coupled to the momentum equations, v˜n+1 = 2vn−vn−1is the velocity field extrapolated to time-level n+1, and vn is the velocity field at time-level n. Theirsemi-implicit correction term [first term in brackets in (2.6)] is the correction to the explicit solutionhn+1exp from the CSLAM scheme, and the second term in brackets in (2.6) is a predictor-correctorterm (where the overbar denotes the departure cell-averaged value). The implicit linear terms areobtained, as in the traditional approach [e.g., Kwizak and Robert (1971); Machenhauer and Olk(1997)], by linearizing the height field around a constant mean reference state, and hence (2.6)15Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportresults in a velocity-divergence form. The notations ∇eul and ∇lag denote discretized divergenceoperators based on the Eulerian and Lagrangian forms respectively. Using notations in Fig. 2.2, theEulerian divergence operator is given by∇eul ·v =1∆x(ur−ul)+1∆y(vt − vb). (2.7)The Lagrangian divergence operator [eq.(25) in LKM] is given by∇lag ·v =1∆A∆A−δA∆t, (2.8)and is computed as the change in cell area in one time step.The form of the semi-implicit correction term in (2.6) is due to the split-divergence approxima-tion [eq. (26) in LKM]∇ ·vn+1/2 ≈12[∇ · v˜n+1 +∇ ·vn](2.9)applied to the linearized divergence term of the semi-implicit continuity equation. The split-divergenceapproximation is used to evaluate the linear divergence term at the mid-point trajectory (at time-leveln+1/2). As explained in Lauritzen et al. (2006), this approximation stems from their trajectory al-gorithm, where the trajectory is approximated as two segments: (i) from the departure point tothe trajectory mid-point (computed iteratively), and (ii) from the mid-point to the arrival grid point(computed using extrapolated winds; see Fig. 1 in LKM). Since the Lagrangian divergence is cal-culated based on the change of cell area over time, and departure cell areas are computed usingthe split-trajectory algorithm, the split-approximation can also be applied to the divergence term(Lauritzen et al., 2006).Ideally, to be consistent, the implicit and the extrapolated divergences would both be solvedin a Lagrangian fashion; however, this would lead to a nonlinear elliptic equation instead of astandard Helmholtz equation (Lauritzen, 2005). To retain a linear elliptic equation, Lauritzen et al.(2006) implemented a predictor-corrector approach to correct for the Eulerian discretization of theimplicit divergence term, and found that this step was necessary to maintain stability in their model.In our implementation of the LKM solver using CSLAM, we follow the approach of Lauritzenet al. (2006), where the predictor-corrector term [second term in brackets in (2.6)] is evaluated byintegrating the departure cell-averaged value over δA∗.16Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportul urvbvt(uc ,vc )1 (uc ,vc )2(uc ,vc )3(uc ,vc )4Figure 2.2: Definition of an Eulerian arrival grid cell, and its associated velocities at the cellfaces (ul,ur,vt ,vb) and cell corners (uc,vc)i for i=1, 2, 3, 4.2.2.3 Numerical inconsistency in semi-implicit continuity equations in avelocity-divergence formNumerical consistency between total mass and constituent mass is difficult to maintain in semi-implicit CISL schemes such as LKM. The prognostic constituent mass variable hq is typically solvedexplicitly using(hq)n+1 = (hq)n+1exp , (2.10)where (hq)n+1exp is the CISL explicit solution, h is the shallow-water height (analogous to total airmass in a full model), and q is the specific concentration of an arbitrary constituent. (From hereon, the constituent mass will be written without the parentheses for simplicity.) The cell-integratedtransport equation in its flux-form helps conserve constituent mass, analogous to the amount ofwater vapour and other passive tracers in an atmospheric model — an important constraint especiallyfor long simulations. Since the departure cell areas are the same for both total fluid mass and theconstituent mass, the weights of the line integrals in CSLAM will need to be computed only onceper time step, and represents one of the advantages of this scheme.If the discrete constituent equation is consistent with the discrete continuity equation, the formershould reduce to the latter when q = 1, and an initially spatially uniform specific concentration fieldshould remain so. For a divergent flow, however, the semi-implicit correction term in (2.6) maybecome large enough such that (2.10), in its explicit form, is no longer consistent (Lauritzen et al.,2008).Alternatively, one can formulate the discrete constituent equation by including similar semi-17Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportimplicit correction and predictor-corrector terms as in (2.6), i.e.hqn+1 = hqn+1exp −∆t2(HQ)0[∇eul ·vn+1−∇lag · v˜n+1]+∆t2(HQ)0[∇eul ·vn−∇lag ·vn]δA∗∆A, (2.11)where (HQ)0 is a constant mean reference constituent mass, the velocities vn+1 are solutions fromthe Helmholtz solver, and v˜n+1 and vn are the same velocities as in (2.6).However, (2.11) is only strictly numerically consistent when q = Q0, where Q0 is the constantmean reference specific concentration, and that (HQ)0 is the product of H0 and Q0. Moreover, thedependence on a constant mean reference constituent mass (HQ)0 may create a source of numericalerrors for regions with little constituent mass. For example, in regions where hqn+1exp = 0, if the flowis highly-divergent such that the terms in square brackets in (2.11) are non-zero, spurious constituentmass will be erroneously generated due to a non-zero constant mean constituent mass. Similarly, inareas where hqn+1exp is a non-zero constant, spurious deviation from constancy can be generated bythe correction terms.The issue with an inconsistent constant mean reference state for the total fluid mass and con-stituent mass fields can be resolved with the formulation we present in the next section.2.3 A consistent and mass-conserving semi-implicit SW solverOur new scheme ensures numerical consistency between the continuity and constituent equationsby formulating the discrete equations, specifically the semi-implicit correction and the predictor-corrector terms, in flux form instead of a velocity-divergence form. The goal is to avoid the use ofa constant reference state, such as (2.6). We test this approach for the SW equations, and refer tothe model using the flux-form scheme as CSLAM-SW. We formulate the semi-implicit flux-formcontinuity equation ashn+1 = hn+1exp −∆t2[∇eul · (hn+1exp vn+1)−∇lag · (hn+1exp v˜n+1)]+∆t2[∇eul · (hnvn)−∇lag · (hnvn)]δA∗∆A, (2.12)and use the explicit CSLAM solution hn+1exp as the reference state in the semi-implicit correctionterm. The shallow-water model CSLAM-SW, like the LKM model, couples the semi-implicit heightcontinuity equation with the traditional semi-Lagrangian momentum equations, as described in theAppendix, and solves the resulting elliptic system with a conjugate-gradient Helmholtz solver.18Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportTo ensure consistency, we simply express the constituent equation ashqn+1 = hqn+1exp −∆t2[∇eul · (hqn+1exp vn+1)−∇lag · (hqn+1exp v˜n+1)]+∆t2[∇eul · (hqnvn−∇lag · (hqnvn)]δA∗∆A, (2.13)where hqn+1exp is the explicit CSLAM update to the constituent mass, the velocities vn+1 in ∇eul ·(hqn+1exp vn+1) are from the SLSI solution, and hqn and vn are the constituent mass and velocityat time-level n respectively. This scheme also resolves the problem of spurious generation of con-stituent mass for regions with near-zero specific concentration (as described in the previous section).The specific concentration q is diagnosed by decoupling the constituent mass usingqn+1 =hqn+1hn+1. (2.14)We note that to ensure numerical consistency, we must eliminate machine-roundoff and convergenceerrors in the Helmholtz solver. In solving for hqn+1, we substitute the solutions of vn+1 derived fromthe Helmholtz solution hn+1 into (2.13). Prior to diagnosing q using (2.14), we must correct thesolution hn+1 by substituting solutions of vn+1 back into (2.12); otherwise, the values of hn+1 canbecome inconsistent with hqn+1. The consistent hn+1 solution is then used to solve for q using (2.14)and in the next time step. To compute hqn+1exp , we follow Nair and Lauritzen (2010) in separating thesubgrid-cell reconstructions for h and q, and then compute hq(x,y) usinghq(x,y) = hq+q(h−h), (2.15)where h = h(x,y) and q = q(x,y) are the reconstruction functions, and (h,q) are cell averages.The new flux-form conservation equations (2.12) and (2.13) involve the computation of an Eu-lerian flux-divergence and a Lagrangian flux-divergence using extrapolated velocities. Using themesh described in Fig. 2.2, the discrete Eulerian flux-divergence is given as∇eul · (hv) =1∆x[(hxu)r− (hxu)l]+1∆y[(hyv)t − (hyv)b], (2.16)where ∆x and ∆y are the grid spacing in the x- and y-directions, and each of the fluxes are evaluatedas hxrur, hxl ul , hyt vt , and hybvb, respectively (see Appendix for operator definitions).The Lagrangian flux-divergence in (2.13) needs to be consistent with the Lagrangian velocity-divergence (2.8). To derive the new operator, we begin by computing the Lagrangian backward-trajectories of the arrival grid cell vertices given in Fig. 2.2. We define the arrival cell corner points19Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportto be at (~x1,~x2,~x3,~x4), i.e. (x1,y1), (x2,y2), (x3,y3), (x4,y4), and the departure cell corner points as~xd1 =~x1−∆t · (uc, vc)1, (2.17)~xd2 =~x2−∆t · (uc, vc)2, (2.18)~xd3 =~x3−∆t · (uc, vc)3, (2.19)~xd4 =~x4−∆t · (uc, vc)4, (2.20)where (uc,vc)i = (uy,vx)i denote the x- and y-velocity components at the ith vertex, where i =1,2,3,4.The area of the departure cell is computed asδA = 12[~xd21×~xd41 +~xd43×~xd23], (2.21)where~xd21 =~xd2−~xd1;~xd41 =~xd4−~xd1;~xd43 =~xd4−~xd3; and~xd23 =~xd2−~xd3. We can then rewritethe departure cell area asδA = ∆x∆y−∆t[Fr−Fl +Ft −Fb]. (2.22)whereFr = uryy∆y− (uc2vc3−uc3vc2)∆t/2, (2.23)Fl = ulyy∆y− (uc1vc4−uc4vc1)∆t/2, (2.24)Ft = vtxx∆x− (uc3vc4−uc4vc3)∆t/2, (2.25)Fb = vbxx∆x− (uc2vc1−uc1vc2)∆t/2 (2.26)(see Appendix for operator definitions).Using (2.22), the velocity divergence can be written as:D=1∆x∆y[Fr−Fl +Ft −Fb], (2.27)which is identical to the Lagrangian divergence (2.8). The first flux term in each of Fr, Fl , Ft ,and Fb is identical to the Eulerian velocity-divergence using velocities averaged to the cell edgesand the remaining terms give the geometric correction for a Lagrangian representation (see Fig.9 in Lauritzen, 2005). Using this velocity divergence, we now approximate the Lagrangian flux-20Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportdivergence term in equation (2.12) as:∇lag · (hv) =1∆x∆y[Frhxr−Flhxl +Fthyt −Fbhyb]. (2.28)Using (2.16) and (2.28) and replacing h with hq, we can further combine each of the terms inbrackets of the constituent equation (2.13), which becomeshqn+1 = hqn+1exp −∆t2[∇eul · (hn+1exp qn+1exp∗v′n+1)]+∆t2[∇eul · (hnqn∗v′n)]δA∗∆A, (2.29)where∇eul · (hq∗v′) =1∆x[hxrq∗xr(ur−Fr/∆y)−hxl q∗xl (ul−Fl/∆y)]+1∆y[hyt q∗yt (vt −Ft/∆x)−hybq∗yb(vb−Fb/∆x)]. (2.30)The corrective velocity v′ is defined as the difference between the velocity field used in the Eulerianflux divergence (2.16) and that derived from the Lagrangian flux areasFr,Fl ,Ft , andFb, dividedby the cell face length. The corrective velocity v′n+1 in (2.29) is computed using vn+1 from theHelmholtz solver and the Lagrangian flux areas based on extrapolated winds divided by the cellface length. The velocity v′n used in the predictor-corrector term in (2.29) is computed using thevelocity field vn at time-level n and the Lagrangian flux areas based on vn, and again divided bythe cell face length. Shape-preserving schemes, e.g. the first-order upwind scheme, or higher-ordermethods such as flux-corrected transport schemes or flux-limiter schemes, can then be applied to thefluxes in (2.29). The first-order upwind scheme is used here, where the upstream values (denotedby the asterisks) qn+1exp∗and qn∗ at each cell face are determined by the directions of v′n+1 and −v′nrespectively [see e.g., Durran (2010), eq.(5.109)]. The first-order upwind scheme is numericallydiffusive (Durran, 2010), but the damping effect on the correction and predictor-corrector termsshould be minimal as the corrective velocities v′n+1 and v′n are typically very small. To ensureshape-preservation in the explicit CSLAM solution, we implement a simple 2D monotonic filter(Barth and Jespersen, 1989) that searches for new local minima and maxima in the reconstructionfunction of q, and scales the function if these values exceed those in the neighbouring cell.Testing of the CSLAM-SW model [based on (2.12) for h, and (A.1) and (A.2) for the velocitycomponents] revealed an instability related to the averaging of the C-grid velocities to the cell corner21Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportpoints in the continuity equation and its interaction with the rotational modes. Following Randall(1994), we can write a generalized discretized dispersion relation for the linearized shallow-waterequations asω3−ω(c2lvlh + c2kukh + fu fv)− ic2( fukhlv− fvkulh) = 0, (2.31)where the terms fu and fv are the discrete Coriolis operators, ku and lv are the discrete height-gradient operators, kh and lh are the discrete velocity-divergence operators in the continuity equa-tion (the letter subscripts refer to the equations in which they appear), and c2 = gH. In the lin-earized shallow-water dispersion relation for C-grid, the last two terms on the L.H.S. of (2.31),fukhlv and fvkulh, cancel and thus there are no numerical frequencies ω with imaginary parts thatamplify in time. Although the CSLAM-SW model uses the C-grid, we have found that the dis-cretization of the linearized Lagrangian divergence is equivalent to taking an average of the u and vvelocities to the corners of the grid cell followed by an averaging back to the cell-faces, i.e. the dis-cretization is equivalent to using a 1-2-1-averaging of the u velocities in the y-direction, and of the vvelocities in the x-direction, at the Eulerian grid cell faces. This averaging leads to non-cancellationof fukhlv and fvkulh, and growing modes. We have found that using the averaging operators vxyxxand uxyyy(see Appendix for operator definitions) on the Coriolis terms in the x- and y-momentumequations, respectively, recovers the cancellation and eliminates the unstable mode.2.4 Test casesWe present four test problems involving divergent flows: a radially-propagating gravity wave (withtwo different initial perturbations), and two highly-nonlinear barotropically unstable jets [the Bick-ley and the Gaussian jets from Poulin and Flierl (2003)]. The gravity-wave problem (section 2.4.1)is a simple case to assess the stability and accuracy of the new SLSI solver (CSLAM-SW) withrespect to an imposed mean flow speed and the gravity-wave propagation speed. We also use thistest case to highlight the issue of numerical inconsistency in the constituent transport scheme ofLKM. The nonlinearity of the unstable jet in the second problem is particularly useful in testingthe stability limits of the new scheme. The Bickley jet (section 2.4.2) has a moderate gradient inthe initial height profile, while the steeper profile in the Gaussian jet (section 2.4.3) drives a moreunstable jet. These strong gradients provide a severe test for advection schemes. In addition to thosefrom LKM, solutions from a traditional semi-Lagrangian formulation and an Eulerian formulation(see Appendix) are also presented for comparison. We use the highly-divergent Gaussian jet caseto compare the solutions between the shape-preserving CSLAM-SW solver described by (2.29) andthe LKM with a shape-preserving explicit transport scheme (section 2.4.4).22Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport2.4.1 A radially-propagating gravity waveA non-rotating ( f = 0) 2D radially-propagating gravity wave is initiated by a circular height pertur-bation h′ and advected by a mean background flow:u(x,y, t = 0) = u0 = 1.2 m s−1, (2.32)v(x,y, t = 0) = v0 = 0.9 m s−1, (2.33)h(x,y, t = 0) = h0 +h′, (2.34)whereh′ =12∆h[1+ cos(pir10 km)], if r ≤ 10 km,0, otherwise,(2.35)and h0 is the initial background height, ∆h is the magnitude of the initial height perturbation, r =√(x− xc)2 +(y− yc)2, and (xc,yc) is the center of a 200 km × 200 km domain. We performtests for two different initial height perturbations: a linear case with ∆h = 10 m and h0 = 990 m;and a nonlinear case with ∆h = 500 m and h0 = 1000 m. A reduced gravitational acceleration ofg′ ≈ 0.0204 m s−2 is used, giving an initial gravity wave speed c =√g′h0 ≈ 4.5 m s−1 for the twocases. The mean advection speed(√u20 + v20 = 1.5 m s−1)is chosen to emulate the speed ratio ofthe fastest advection of sound waves (≈ 300 m s−1) in the atmosphere to the speed of the jet stream(≈ 100 m s−1). The background flow velocities u0 6= v0 are also chosen to ensure that the flow doesnot align with the mesh.The model domain consists of 400 × 400 grid cells, with a grid spacing of ∆x = ∆y = 500 m,and is periodic in both x- and y-directions. Since there is no analytical solution to the test problem, toevaluate CSLAM-SW, we produce a fine-resolution Eulerian reference solution with a grid spacingof ∆x = ∆y = 100 m and a time step of ∆t = 10 s. The center of the gravity wave disturbance in thereference solution is stationary (i.e. u0 = v0 = 0 m s−1), and we compare the solutions by translatingthe gravity wave disturbance in CSLAM-SW to the center of the domain.In addition to CSLAM-SW, we also run the two initial perturbation cases using LKM, the tra-ditional semi-Lagrangian formulation, and an Eulerian formulation. We use the l2-norm of error asthe error measure, which for a uniform mesh isl2 =√∑i, j[h(i, j)−href(i, j)]2√∑i, j[href(i, j)]2, (2.36)23Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport101 102 10310−610−510−410−3(a) L2 error norms, linear case l 2∆t (s) EulerianCSLAM−SWLKMTRAD−SL101 102 10310−310−210−1(b) L2 error norms, nonlinear case l 2∆t (s) EulerianCSLAM−SWLKMTRAD−SLFigure 2.3: Comparison of the height field L2 error norms for the radially-propagating gravity-wave solutions. Errors are plotted at time T = 1×105 s for the (a) linear (∆h = 10 m andh0 = 990 m) and (b) non-linear (∆h = 500 m and h0 = 1000 m) test cases computed on a500 m mesh. Note the different scales in the plots.where i, j are the grid indices, h(i, j) is the model solution, and href(i, j) is the Eulerian high-resolution reference solution. The l2-norm of error in the height field (at time T = 1× 105 s) fordifferent time step sizes is shown in Fig. 2.3 for all four models. Results from both the linear andnonlinear initial perturbations are plotted. The time truncation error in CSLAM-SW is very com-parable to those in the other two semi-Lagrangian models for both cases. Except for the Eulerianmodel, all model solutions converge as the time step size is reduced to less than ∆t = 50 s. Atthis point, differences between the errors are mainly due to the spatial discretization schemes (more24Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportnoticeably in the nonlinear case). The Eulerian model and the traditional semi-Lagrangian modelhave a commonality that they both use a ‘true’ C-grid divergence operator in the continuity equa-tion; whereas as discussed in section 2.3, the CISL computation of divergence in both CSLAM-SWand LKM consists of an extra averaging operator. For this reason, one may see a smaller spatialdiscretization error in the traditional semi-Lagrangian model and “coarse” Eulerian model whencompared to an Eulerian high-resolution reference solution than those in the CISL models, as is thecase in Fig. 2.3.To evaluate the consistency in CSLAM-SW and LKM, a constituent with an initially constantspecific concentration distribution (q0 = 1) is initialized in each model. The CSLAM explicit trans-port scheme conserves constituent mass in both models; however, as discussed in section 2.2.3,when numerical consistency is violated, constancy of the specific concentration is not guaranteed,and generation or removal of constituent mass is possible. The specific concentration is diagnosedby decoupling the constituent mass variable using (2.14). A time step of ∆t = 100 s is used. Fig.2.4 shows an example of the specific concentration error in LKM at time T = 1× 105 s for boththe linear and nonlinear perturbation cases. The error is largest near the leading edge of the gravitywave, where the flow is most divergent and the semi-implicit correction term is non-zero. Fig. 2.5shows the variation in error with time step size for both the linear and nonlinear perturbations atthe same simulation time as in Fig. 2.4. The error measures used are the maximum absolute error,the mean absolute error, and the root-mean-squared error. Errors in the solutions from LKM andCSLAM-SW are shown in solid and dashed lines respectively. Since the inconsistent semi-implicitcorrection in (2.6) is proportional to ∆t, errors in the scalar field grow with time step size, which canbecome a major issue for semi-Lagrangian models that take advantage of larger stable time steps.For the nonlinear test, the maximum absolute error from LKM is in the order of 10−2 to 10−1, andis significant for constituents like water vapour which has a typical mixing ratio of roughly 0.1% to3% in air. On the other hand, CSLAM-SW using a consistent formulation is free-stream preserving(up to machine roundoff) for both cases and all time-step sizes tested.2.4.2 Bickley jet – Ro = 0.1The stability of CSLAM-SW is further evaluated with two perturbed jets; we begin with the Bick-ley jet from Poulin and Flierl (2003). The Bickley jet is simulated at the Rossby number, Ro =U/ f L = 0.1, where U is the flow velocity scale, f is the Coriolis parameter and L is the lengthscale of the jet width. We choose the Froude number, Fr = ( f L)2/g′H = 0.1. The jet is character-ized by greater heights to the left of the channel and dropping off to smaller heights to the right,geostrophically-balanced by a mean flow velocity down the channel (Fig. 2.6). A height pertur-bation is superimposed at the initial time, causing wave amplification and eventual breaking of the25Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxyDeviation (q−q0) from the initially constant q0 = 1in the nonlinear perturbation case 100 200 300 400100200300400−1 0 1x 10−3xyDeviation (q−q0) from the initially constant q0 = 1in the linear perturbation case 100 200 300 400100200300400−1 0 1x 10−5q−q0q−q0x10−3(b) nonlinear case(a) linear case−5x10Figure 2.4: Specific concentration error (q−q0) in LKM for a divergent flow initialized with aconstant q0 = 1 in the (a) linear (∆h = 10 m and h0 = 990 m) and (b) nonlinear (∆h = 500m and h0 = 1000 m) height perturbation cases. Note the different scales in the plots.26Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport101 102 10310−1610−1410−1210−1010−810−610−410−2100102∆t (s)Error measures(a) Linear perturbation Max Abs ErrorMean Abs ErrorRMSE101 102 10310−1610−1410−1210−1010−810−610−410−2100102∆t (s)Error measures(b) Nonlinear perturbation Max Abs ErrorMean Abs ErrorRMSEFigure 2.5: Variation of specific concentration error (q−q0) (maximum absolute error, meanabsolute error, and root-mean-squared error) with time step size in LKM (solid line) andCSLAM-SW (dashed line) for the (a) linear height perturbation and (b) nonlinear heightperturbation cases.jet into vortices, and formation of a vortex street along the channel. These vortex streets consist ofthin filaments of vorticity with strong horizontal velocity shear, making it a good test because it ischallenging for all numerical schemes. A more detailed description of the evolution of these jetscan be found in Poulin and Flierl (2003).The initial geostrophically balanced mean state (u0, v0, and h0) and height perturbation h′ of theBickley jet is given by:u(x,y, t = 0) = u0 = 0, (2.37)v(x,y, t = 0) = v0 =−g′∆hf asech2( xa), (2.38)h(x,y, t = 0) = h0 +h′, (2.39)27Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transport−1000 km x = 0 km 1000 kmHeight, h 0 (m)Profile of the initial geostrophically−balanced mean state∆h = 1 m, 50 m h = 100 m| |a = 100 km−1000 km x = 0 km 1000 kmVelocity, v 0 (m s−1)∆v = 1 m/s, 56 m/sv = 0 m/sFigure 2.6: Initial mean height h0 (top) and velocity v0 (bottom) profiles for the Bickley jet(∆h = 1 m, ∆v = 1 m s−1) and Gaussian jet (∆h = 50 m, ∆v = 56 m s−1).whereh0 = 100−∆h tanh( xa), (2.40)h′ = 0.1∆h sech2( xa)sin(2piyYn). (2.41)The parameter ∆h is the maximum amplitude of the height perturbation and depends on Ro, g′ is thegravitational acceleration, a is the jet width, Y is the length of the channel, and n is the wavenumbermode of the height perturbation. In our simulations, L = a = 1× 105 m, X (width of channel)=Y = 2×106 m, f = 1×10−4 s−1, and g′ = 10 m s−2. For the specified scale of the jet width and aflow with Fr = 0.1, the mean height of h0 is 100 m. The amplitude of the height perturbation ∆h = 1m is determined by the scale of the initial geostrophically balanced flow speed (U ∼ 1 m s−1) forRo = 0.1. We choose the most unstable mode of wavenumber n = 3 (Poulin and Flierl, 2003) forall of our jet simulations.Each grid domain has 202 × 202 grid cells and a grid spacing of ∆x = ∆y = 9950 m, withsolid boundary conditions at x = −X/2 and x = X/2 and periodic boundary conditions in y wherey ∈ [−Y/2,Y/2 ]. A time step of ∆t = 2000 s was used in all simulations. Based on the initialgravity-wave speed c ≈ 32 m s−1 and initial flow speed |v| = 1 m s−1, the Courant numbers areCrgw = 6.4 and Cradv = 0.2 respectively.28Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106(d) Eulerian(c) LKM(b) TRAD­SL(a) CSLAM−SWFigure 2.7: Solutions of the Bickley jet at time T = 5×106 s (after 2500 time steps) for Ro =0.1, Fr = 0.1 and Cradv = 0.2. Plotted are positive (solid line) and negative (dashed line)vorticity between −1× 10−5 s−1 and 1× 10−5 s−1 with a contour interval of 5× 10−7s−1.29Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportTo maintain numerical stability in the Eulerian model, we implemented a second-order ex-plicit diffusion term with a numerical viscosity parameter βx = βy = ν∆t/∆x2 = 0.02 (where νis analogous to the physical viscosity). This value corresponds to the numerical Reynolds number,Re = UL/ν = 102, a factor of 10 smaller than that used in the forward-in-time Eulerian model ofPoulin and Flierl (2003). Explicit diffusion was not applied to any of the semi-Lagrangian mod-els because the schemes have sufficient inherent damping to maintain numerical stability. For thetraditional semi-Lagrangian model, however, we found that time-off-centering in the semi-implicitscheme was needed to maintain stability.Fig. 2.7 shows the solutions from CSLAM-SW and the three comparison models. Althoughthe exact form of the initial height perturbation was not provided in Poulin and Flierl (2003), wewere able to reproduce results very similar to theirs [cf. Fig. 4c of Poulin and Flierl (2003)].The most noticeable difference among the different model solutions is in the shape and magnitudeof the relative vorticity maxima and minima. CSLAM-SW showed very similar vortex shapes tothose from LKM and TRAD-SL. The vortices in the Eulerian results are similar to those from theEulerian model of Poulin and Flierl (2003). The difference between the Eulerian solution and thesemi-Lagrangian solutions can be attributed to the inherent damping in the reconstruction step ofthe CISL schemes and the grid-point interpolation in the traditional semi-Lagrangian scheme.2.4.3 Gaussian jet – Ro = 5.0The third test case is the Gaussian jet with Ro = 5.0. Similar to the Bickley jet, the Gaussian jet hasFr = 0.1, and has an initially geostrophically-balanced mean-state with greater heights to the left ofthe channel and dropping off to smaller heights to the right (Fig. 2.6). The main difference betweenthe two jets is that the Gaussian jet has a slightly steeper height profile at the center of the channel,and therefore, produces a more pronounced nonlinear flow, especially at larger Ro. The initial meanstate and height perturbation for the Gaussian jet is given as:u(x,y, t = 0) = u0 = 0, (2.42)v(x,y, t = 0) = v0 =−2g′∆h√pi f aexp(−(x/a)2), (2.43)h(x,y, t = 0) = h0 +h′, (2.44)30Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportwhereh0 = 100−∆h erf( xa), (2.45)h′ = 0.1∆h( 2√piexp(−(x/a)2)sin(2piyyLn), (2.46)and the notation is as before. All the parameters remain the same, except ∆h = 50 m for Ro =5.0, and ∆t = 100 s is used. With an initial gravity-wave speed and maximum flow speed of 38 ms−1 and 56 m s−1 respectively, Crgw = 0.4 and Cradv = 0.56. We note that U > c, i.e. the flow issupercritical. Despite the existence of supersonic waves in the solution, CSLAM-SW is stable evenat larger Courant numbers.As pointed out in Poulin and Flierl (2003), jets in this Rossby regime are highly unstable andof particular interest is the formation of an asymmetric vortex street with triangular cyclones andelliptical anticyclones. As the vortex street is advected towards the deeper water, a strong cut-offcyclone develops due to vortex stretching (adjacent to the main anticyclonic feature). All of ourmodels, including CSLAM-SW, were able to reproduce these features [Fig. 2.8, cf. Fig. 10e ofPoulin and Flierl (2003)]. As in the Bickley jet case, we find that CSLAM-SW produced solutionssimilar to the other two semi-Lagrangian models (LKM and TRAD-SL).In addition to comparing solutions of CSLAM-SW at time steps allowable by the Eulerianscheme, we also tested the stability of CSLAM-SW at a much larger Cradv = 2.5. Figs. 2.9a-c show solutions at various times from the previous CSLAM-SW simulation (Cradv = 0.56), andFigs. 2.9d-f show solutions at each of the corresponding time for Cradv = 2.5, using the largest timestep allowable by the Lipschitz condition for this flow. The solution from the Cradv = 2.5 simulationis almost identical to the solution using Cradv = 0.56.The CSLAM-SW is numerically stable for the highly nonlinear flow in the Gaussian jet and atCourant numbers much greater than unity. To check that consistency and shape-preservation in sucha highly divergent flow can be maintained, we repeat the Gaussian jet case using CSLAM-SW andthe shape-preserving extensions described in section 2.3.2.4.4 Gaussian jet – Ro = 5.0 with shape-preservationThe shape-preserving CSLAM-SW solver (2.29) is tested using the divergent flow of the Gaussianjet as described in section 2.4.3. We also test the LKM solver with the Barth and Jespersen (1989)filter implemented in the explicit scalar transport scheme of hqn+1exp . All parameters are as describedin section 2.4.3, and a time step of ∆t = 100 s is used for results in Figs. 2.10 and 2.11.To test for numerical consistency in the two solvers, we repeat the consistency test described31Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106 (d) Eulerian(c) LKM(a) CSLAM−SW (b) TRAD­SLFigure 2.8: Solutions of the Gaussian jet for Ro = 5.0 and Cradv = 0.56 at time T = 1.8×105s (after 1800 time steps). Plotted are positive (solid line) and negative (dashed line)vorticity between −5× 10−4 s−1 and 5× 10−4 s−1 with a contour interval of 5× 10−5s−1.32Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 106xy−1 −0.5 0 0.5 1x 106−1−0.500.51 x 1065 4 (e) T = 1.0 x 10 s(d) T = 5.0 x 10 s (f) T = 1.4 x 10 s5 5 (a) T = 5.0 x 10 s4 (b) T = 1.0 x 10 s (c) T = 1.4 x 10 s5 Figure 2.9: CSLAM-SW solutions of the Gaussian jet for Ro = 5.0 at three different times (leftto right on each row) of the simulation: at time T = 5×104 s, 1.0×105 s, and 1.4×105s. (a - c) Solutions using a Cradv of 0.56 (same simulation as in Fig. 2.8) (d - f) Solutionsusing a larger Cradv of 2.5. Plotted are positive (solid line) and negative (dashed line)vorticity between −5× 10−4 s−1 and 5× 10−4 s−1 with a contour interval of 5× 10−5s−1.in section 2.4.1 by initializing a constant specific concentration field q0 = 1. The shape-preservingCSLAM-SW solution is able to maintain numerical consistency between h and hq up to machineroundoff for this highly divergent flow and the result is independent of time-step size. As for LKM,despite the shape-preserving transport scheme in the solver, numerical inconsistency is still an issuewith a maximum absolute error (defined as the deviation from q0 = 1) of 6.79× 10−3, a meanabsolute error of 4.82×10−4, and a root-mean-squared error of 1.06×10−3 at time T = 1.8×105s (Fig. 2.10), and as in section 2.4.1, the error is a function of the time-step size (not shown).To compare the shape-preservation ability between CSLAM-SW and LKM, we initialize aspecific-concentration distribution that varies only in the x-direction and has a sharp gradient that33Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxyDeviation from the initial constant q0 = 1 field(q−q0) 100 200100200−1 0 1x 10−3q−q0Specific concentration errorx10−3Figure 2.10: Specific concentration error (q− q0) in LKM for the Gaussian jet at time T =1.8×105 s, initialized with a constant q0 = 1 field.coincides with the center of the initial jet:q(x,y, t = 0) =1.0, if −X/2≤ x < 0.0.1, if 0≤ x < X/2.(2.47)Solutions of q diagnosed from hq from the non-shape-preserving CSLAM-SW, LKM with shape-preserving transport, and the shape-preserving CSLAM-SW are presented in Figs. 2.11a-c. Thesimulation time T = 1.8×105 s in the figure corresponds to the vorticity field shown in Fig. 2.8.For the non-shape-preserving CSLAM-SW solver (Fig. 2.11a), q reaches an unphysical peakvalue of 1.233 and an unphysical minimum value of -0.145 (specific concentrations cannot be neg-ative). The LKM solver with shape-preserving transport (Fig. 2.11b) has less severe errors thanthe non-shape-preserving CSLAM-SW, but loses its shape-preserving ability due to numerical in-consistency. The minimum and maximum q values are 0.09997 and 1.0063 respectively at timeT = 1.8×105 s. The overshooting of q (which may generate spurious constituent mass) appears tobe greater in amplitude than the undershooting for this flow. Overshooting occurs mostly within the34Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportstrongest anticyclones (negative vorticity centers on the left side of the channel, highlighted in solidblack lines in Fig. 2.11b). Using the shape-preserving CSLAM-SW solver (Fig. 2.11c), minimumand maximum values of q are kept within its physical limits (0.1 and 1.0 respectively, up to machineroundoff) and shape-preservation is ensured.2.5 ConclusionA conservative and consistent semi-Lagrangian semi-implicit solver is constructed and tested forshallow-water flows (CSLAM-SW). The model uses a new flux-form discretization of the semi-implicit cell-integrated semi-Lagrangian continuity equation that allows a straight-forward imple-mentation of a consistent constituent transport scheme. Like typical conservative semi-Lagrangiansemi-implicit schemes, the algorithm requires at each time step a single Helmholtz equation solutionand a single application of CSLAM.Specifically, our new discretization uses the flux divergence as opposed to a velocity diver-gence that requires linearization about a constant mean reference state. For traditional semi-implicitschemes, the dependence on a constant mean reference state makes it difficult to ensure consistencybetween total fluid mass and constituent mass. When numerical consistency is not maintained, con-stituent mass conservation can be violated even for solvers that use inherently-conservative transportschemes. More unacceptably, constituent fields may no longer preserve their shapes, e.g. losingconstancy or positive-definiteness.We have shown an example of a traditional discrete cell-integrated semi-Lagrangian semi-implicit continuity equation (LKM), in which inconsistency can generate significant numerical er-rors in the specific constituent concentration. The inconsistent semi-implicit correction term inLKM causes errors to grow proportionally with time step size and with the nonlinearity of the flow.The ideal radially-propagating gravity wave tests using the LKM solver showed a maximum ab-solute error in an initially constant specific concentration (q0 = 1) field ranging from an order of10−7 to 10−3 in the linear case, and an order of 10−4 to 10−1 in the nonlinear case. The orders ofmagnitude of these errors are significant relative to the specific concentration of tracers and watervapour in the atmosphere. The consistent formulation in the new CSLAM-SW on the other handeliminates these errors (up to machine roundoff).The new flux-form solver (CSLAM-SW) is tested for a range of flows and Courant numbersfor the shallow-water system, and is stable and compares well with other existing semi-implicitschemes, including a two-time-level traditional semi-Lagrangian scheme and an Eulerian leap-frogscheme. The Gaussian jet test (the more nonlinear jet of the two presented) showed that CSLAM-SW remains numerically stable when large time steps are used.35Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportWe have also identified and eliminated a computational unstable mode in CSLAM-SW andLKM, using the discrete dispersion relation of the linearized shallow-water equations. The numer-ical instability, associated with the Lagrangian divergence operator on a C-grid, can be eliminatedby introducing a new averaging operator on the Coriolis terms in the momentum equations.Shape-preservation in CSLAM-SW is ensured by applying a 2D shape-preserving filter in theCSLAM transport scheme and the first-order upwind scheme to compute the predictor-correctorand flux-form correction terms. As shown in the Gaussian jet case, without any shape-preservingfilter, unphysical negative and unreasonable positive specific concentrations may develop due toundershoots and overshoots. For inconsistent formulations such as that in LKM, the use of a shape-preserving explicit transport scheme cannot guarantee shape-preservation either due to numericalconsistency errors. CSLAM-SW, on the other hand, allows for straightforward implementationof existing shape-preserving schemes and filters and ensures shape-preservation (up to machineroundoff).The initial testing of the semi-implicit formulation in CSLAM-SW shows promising results.We are currently implementing the extension of CSLAM-SW to a 2D (x-z) non-hydrostatic, fully-compressible atmospheric solver. The desirable properties of mass conservation, consistency, andshape-preservation for moisture variables and tracers will likely be important for both short- andlong-term meteorological applications.36Chapter 2: A two-dimensional shallow-water equations solver with conservative and consistent transportxySpecific concentration, q −1 −0.5 0 0.5 1x 106−1−0.500.51 x 1060.1 1xySpecific concentration, q −1 −0.5 0 0.5 1x 106−1−0.500.51 x 1060.1 1xySpecific concentration, q −1 −0.5 0 0.5 1x 106−1−0.500.51 x 1060.1 1Specific concentration, qSpecific concentration, q Specific concentration, q(b) Shape−preserving LKM (c) Shape−preserving CSLAM−SW(a) Non−shape−preserving CSLAM−SWFigure 2.11: Specific constituent concentration q at time T = 1.8×105 s. Initial minimum andmaximum q are 0.1 and 1.0 respectively. Regions with unphysical overshooting (red)and undershooting (purple) are highlighted.37Chapter 3Extension to a two-dimensional (x-z)nonhydrostatic atmospheric solver3.1 IntroductionSemi-Lagrangian semi-implicit (SLSI) schemes have been widely used in climate and numericalweather prediction (NWP) models since the pioneering work of Robert (1981) and Robert et al.(1985). The more lenient numerical stability condition in these schemes allows larger time stepsand thus increased computational efficiency. Traditional semi-Lagrangian schemes are not inher-ently mass-conserving due to their use of grid-point interpolation, and the lack of conservation canlead to accumulation of significant solution errors (Rasch and Williamson, 1990a; Machenhauer andOlk, 1997). To address this issue, conservative semi-Lagrangian schemes, also called cell-integratedsemi-Lagrangian (CISL) transport schemes (Rancic, 1992; Laprise and Plante, 1995; Machenhauerand Olk, 1997; Zerroukat et al., 2002; Nair and Machenhauer, 2002; Lauritzen et al., 2010), havebeen developed. Although CISL transport schemes, when applied in fluid flow solvers, allow forlocally (and thus globally) conservative transport of total fluid mass and constituent (i.e. tracer)mass, a lack of consistency arises between the numerical representation of the total dry air massconservation, to which we will refer as the continuity equation, and constituent mass conservationequations (Jo¨ckel et al., 2001; Zhang et al., 2008; Wong et al., 2013). Numerical consistency in theflux-form equation for a tracer requires the equation for a constant tracer field to correspond numer-ically to the mass continuity equation; this consistency ensures that an initially spatially uniformpassive tracer field will remain so. The lack of numerical consistency between the two can lead tothe unphysical generation or removal of model constituent mass, which can introduce significanterrors in applications such as chemical tracer transport (Machenhauer et al., 2009).Recently, Wong et al. (2013) introduced a new flux-form formulation of the semi-implicit CISLheight conservation equation for the shallow-water equations (SWE) solver. They showed thatthe scheme is accurate and stable even for highly-nonlinear barotropically-unstable jets and largeCourant numbers. They also found that the use of a shape-preserving filter in an inconsistent for-38Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solvermulation of the continuity equations is ineffective, highlighting the importance of numerical con-sistency in these models.In this paper, the flux-form semi-implicit SWE formulation is extended to the fully-compressibletwo-dimensional (x-z) moist nonhydrostatic equations for the atmosphere. We refer to this newconservative and consistent nonhydrostatic solver as CSLAM-NH. A nonhydrostatic model permitsfast-moving internal gravity and acoustic waves. Here, we integrate the terms responsible for theacoustic waves in a semi-implicit manner to allow large time steps while maintaining stability forthese waves. As in Wong et al. (2013), our nonhydrostatic solver is based on the ConservativeSemi-LAgrangian Multi-tracer (CSLAM) transport scheme, a CISL transport scheme developedby Lauritzen et al. (2010) that has been implemented in NCAR’s High-Order Methods ModelingEnvironment [HOMME; Erath et al. (2012)].The semi-implicit CISL nonhydrostatic solver has six main advantages and desirable properties.As we will show, our nonhydrostatic cell-integrated semi-Lagrangian solver is (1) inherently mass-conserving, (2) shape-preserving, and, with the new formulation, (3) has numerically consistenttransport. The discretization (4) does not depend on a mean reference state, but maintains the sameframework as typical semi-implicit CISL solvers, where (5) a single linear Helmholtz equation issolved and (6) a single application of CSLAM is needed per time step.The paper is organized as follows. The governing equations of the two-dimensional fully-compressible nonhydrostatic system are first described in section 3.2. We then present the proposeddiscretization of the governing equations, including a consistent formulation of the moisture conser-vation equations (section 3.3). The desirable properties of the nonhydrostatic solver are discussedin section 3.4. We test the nonhydrostatic solver with three idealized test cases and compare resultswith an Eulerian split-explicit time-stepping scheme (section 3.5). A fourth test case on numericalconsistency is also presented in section 3.5 to demonstrate the shape-preserving ability of the solverwith additional passive tracers. A summary is given in section 3.6.3.2 Governing equationsThe model governing equations are the two-dimensional (x-z) moist Euler equations in Cartesiangeometry:∂u∂ t +u∂u∂x +w∂u∂ z =−piρmγRd∂Θ′m∂x +Fu, (3.1)∂w∂ t +u∂w∂x +w∂w∂ z =−piρmγRd∂Θ′m∂ z +gρm[ρdpi ′pi −ρ′m]+Fw, (3.2)∂Θm∂ t +∇ · (Θmv) = FΘ, (3.3)39Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver∂ρd∂ t +∇ · (ρdv) = 0, (3.4)∂Q j∂ t +∇ · (Q jv) = FQ j , (3.5)p = p0(RdΘmp0)γ, (3.6)where pi = (p/p0)κ is the Exner function, κ = Rd/cp, γ = cp/cv = 1.4, Rd = 287 J kg−1 K−1,cp = 1003 J kg−1 K−1, and g = 9.81 m s−2. Perturbation variables from a time-independent hy-drostatically balanced background state are used to reduce numerical errors in the calculations ofthe pressure gradient terms (Klemp et al., 2007). The hydrostatically balanced background stateis defined as dp(z)/dz = −ρd(z)g. Perturbation variables are defined as Θm = ρd(z)θ(z) +Θ′m,pi = pi+pi ′, ρd = ρd(z)+ρ ′d , and the moist density ρm = ρd(1+qv +qc +qr), where qv, qc, and qrare the mixing ratios for water vapor, cloud, and rainwater, respectively. The F(·) terms representdiffusion, and any diabatic effects and parameterized physics when moisture is present.As in Klemp et al. (2007), fluxes are coupled to the dry density ρd . The flux variables are givenasΘm = ρdθm and Q j = ρdq j, (3.7)where θm is the modified potential temperature θm = θ(1 + a′qv) where a′ ≡ Rv/Rd ' 1.61 andq j = (qv,qc,qr).The momentum equations are cast in their advective form, and all other equations, i.e., fordensity, potential temperature, and moist species, are cast in their conservative flux-form. Pressureis a diagnostic variable given by the equation of state. The governing equations are based on Klempet al. (2007); the pressure gradient terms in (3.1) and (3.2) have been recast in terms of Θ′m using(3.6) to derive the relation∇p = γRdpi∇Θm, (3.8)and enables us to form an implicit equation for Θ′ (section 3.3). The equations are still exact andno approximations have been applied. The only difference from the governing equations in Klempet al. (2007) is that their momentum equations are cast in the conservative flux-form, whereas theadvective form is used here to facilitate the use of the traditional semi-Lagrangian method.40Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver(b)(a)Figure 3.1: (a) Exact departure cell area (δA∗, dark grey region) and the corresponding arrivalgrid cell (∆A, light grey region). (b) Departure cells in CSLAM (δA) are represented aspolygons defined by the departure locations of the arrival grid cell vertices. (Wong et al.,2013)3.3 A consistent and mass-conserving nonhydrostatic solver3.3.1 CSLAM — a conservative transport schemeTo ensure mass conservation, we utilize an inherently conservative semi-Lagrangian transport schemecalled CSLAM (Lauritzen et al., 2010). The CSLAM transport scheme is a backward-in-time CISLscheme1, where the departure grid cell area δA∗ is found by tracing the regular arrival grid cell area∆A back in time one time-step ∆t (Fig. 3.1a). The CSLAM discretization scheme for the lhs of(3.3), (3.4), and (3.5) is given byφ n+1exp ∆A =∫δA∗φ ndA = φ n∗ δA∗ (3.9)where φ = Θm,ρd , or Q j. The superscript denotes the time level, and φ n+1exp is the explicit cell-averaged transport term computed by integrating the field φ n over the departure cell area δA∗,which gives the cell-averaged departure value φ n∗ .The departure cell area δA∗ in CSLAM is found through iterative trajectory computations fromthe four vertices of an arrival grid cell (unfilled circles in Fig. 3.1b) to their departure points (filled1note that CSLAM may also be cast in flux-form (Harris et al., 2011)41Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solvercircles in Fig. 3.1b). The departure cell area is then approximated using straight lines as cell edges2(dark grey region δA in Fig. 3.1b). To integrate the field φ n over δA, CSLAM implements aremapping algorithm that consists of a piecewise quasi-biparabolic subgrid-cell-reconstruction ofφ n in the two coordinates as in Nair and Machenhauer (2002) with an additional cross term asdescribed in Jablonowski (2004) that helps smooth subgrid distribution near sharp gradients,φ n(x,z) =〈φ n〉+axx+bx(112− x2)+azz+bz(112− z2)+12(cxz + czx)xz (3.10)where coefficients ax, bx, az, bz of the reconstructed parabolic function in the two coordinates areobtained as in Nair and Machenhauer (2002), and the cross-term coefficients cxz and czx are obtainedas in Jablonowski (2004). An average of the two coefficients of the cross term, cxz and czx, is takento avoid a directional bias (Jablonowski, 2004). The cell-average value over the Eulerian grid cellis denoted as〈φ n〉.The integration of the reconstruction function over the departure cell area is then computed.The area integration in CSLAM is transformed into a series of line integrals using the Gauss-Greentheorem, and involves solving for a set of weights w(i, j) that depends only on the departure cellboundary. As described in Lauritzen et al. (2010), the discrete conservative transport scheme fordeparture cell k is∫δA∗φ ndA =Lk∑l=1[∑i+ j≤2c(i, j)l w(i, j)kl](3.11)where c(0,0)l , c(1,0)l , c(0,1)l , c(2,0)l , c(0,2)l are the coefficients for the constant, x, z, x2, and z2 termsrespectively, c(1,1)l is the coefficient for the xz term in (3.10), and l is the index for the Eulerian gridcell(s) with which departure cell k overlaps (of a total of Lk overlapping Eulerian grid cells). Thepartitioning of the areal integration into computation of coefficients and weights greatly enhancesthe transport scheme’s computational efficiency for multi-tracer transport, as the weights can bereused for the remapping of all tracer species in the model. For full details on the basic CSLAMscheme, see Lauritzen et al. (2010); for high-resolution spherical implementations of CSLAM, thereader is referred to the modifications to the scheme documented in Erath et al. (2013). A rigorousassessment of the accuracy of linear transport using CSLAM [for the test case in Lauritzen et al.(2012)] and a comparison of CSLAM to a collection of state-of-the-art transport schemes can befound in Lauritzen et al. (2014).2higher-order edge approximations have been explored in Ullrich et al. (2012)42Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver3.3.2 Trajectory algorithmTo find the departure cell area, we trace the vertices of each arrival grid cell back one time step ∆tusing a trajectory algorithm described in Lauritzen et al. (2006). The trajectory is approximatedand split into two segments: departure grid point to trajectory midpoint, and trajectory midpoint toarrival grid point. The split-trajectory approximation facilitates the semi-implicit formulation of theflux-form conservation equation (section 3.3.4).The displacement in the two linear segments are determined using velocities at time-level n andvelocities extrapolated to time-level n+1, respectively. The first segment (from the departure pointposition rnD to midpoint trajectory rn+1/2D/2 ) is approximated asrn+1/2D/2 = rnD +∆t2vnD, (3.12)We iterate (3.12) three times to increase the accuracy of the computation of vnD. At each iteration, thevelocities are interpolated to the estimated departure location using bicubic Lagrange interpolation.The second segment (from midpoint trajectory rn+1/2D/2 to the arrival point rn+1) is approximatedusingrn+1/2D/2 = rn+1−∆t2v˜n+1, (3.13)where v˜n+1 is evaluated at the arrival grid point and denote velocities extrapolated to time-leveln+1 using a two-time-level extrapolationv˜n+1 = 2vn−vn−1. (3.14)To find rnD, we take the sum of the two half-trajectories [(3.12) and (3.13)],rnD = rn+1−∆t2(vnD + v˜n+1). (3.15)Higher-order approximations to the trajectory can be made by including an acceleration term asdescribed in McGregor (1993). Tests including an acceleration term (not shown) showed that sucha higher-order approximation made little difference to the solutions for this suite of tests.3.3.3 Discretization of the momentum equationsThe momentum equations are solved using the traditional semi-Lagrangian semi-implicit method,where material derivatives such as du/dt = ∂u/∂ t + u∂u/∂x + w∂u/∂ z and dw/dt = ∂w/∂ t +43Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solveru∂w/∂x+w∂w/∂ z [lhs of (3.1) and (3.2), respectively] are computed using a grid-point interpola-tion to the departure point. The two-time-level discretizations of the momentum equations areun+1A =[u−∆t(1−β2)(piρm)xγRdδxΘ′]nD+∆t(Fu)nD−∆t(1+β2)(pinρnm)xAγRdδxΘ′n+1A , (3.16)andwn+1A =[w−∆t(1−β2)(piρm)zγRdδzΘ′]nD+∆tρnmz[gρdpi ′pi −gρ′m]n+1/2D/2z−∆t(1+β2)(pinρnm)zAγRdδzΘ′n+1A+∆t(Fw)nD, (3.17)where the subscripts D, D/2 and A denote evaluation at the departure, midpoint trajectory, andarrival grid points, respectively, and the superscripts denote the time level. The spatial operators aredefined as(·)x=12((·)i,k +(·)i+1,k), (3.18)(·)z=12((·)i,k +(·)i,k+1), (3.19)δx(·) =(·)i+1,k− (·)i,k∆x, and (3.20)δz(·) =(·)i,k+1− (·)i,k∆z. (3.21)The gradient terms responsible for the fast-moving acoustic waves are solved implicitly with theoption of off-centering by setting β 6= 0. Numerical diffusion is represented in Fu and Fw in the44Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverform of second-order diffusion with physical viscosity ν ,F(·) = ν[δ 2x (·)+δ 2z (·)]. (3.22)The buoyancy terms in the vertical momentum equation are solved explicitly by extrapolating totime level n+1/2 using(·)n+1/2 =32(·)n−12(·)n−1, (3.23)and then interpolated to the midpoint trajectory. One way to evaluate the buoyancy term implicitlyis to concurrently update the density and pressure perturbation variables (ρ ′m and pi ′, respectively) atevery iteration of Θ˜′m in the Helmholtz solver. This implicit treatment of the buoyancy term involvesupdating the density perturbation using un+1 and wn+1 guesses at that iteration, and we have yet tofind a feasible way to incorporate this in the Helmholtz solver that ensures convergence at large timesteps. The implicit treatment of the buoyancy terms will be the scope of future work.3.3.4 Discretization of the thermodynamic equationIn our nonhydrostatic solver, we form and solve an implicit equation for Θn+1m . The implicit equationis formed in two steps. First, we compute the explicit solution of the flux-form thermodynamicequation using the conservative transport scheme CSLAM,Θˆn+1m =Θn+1m,exp +∆t2[∇eul · (Θnmvn)−∇lag · (Θnmvn)]δA∗∆A+∆t[FnΘm]δA∗∆A, (3.24)where the notation[·]denotes departure cell averages. The first term on the rhs of (3.24), Θn+1m,exp,is the explicit CSLAM update. The second term is a predictor-corrector term integrated over thedeparture cell to account for the discrepancy between the discrete Eulerian and Lagrangian fluxdivergences in the semi-implicit flux-form correction term. Similarly, in FΘm , second-order diffusion(with mixing coefficient given by ν times the Prandtl number) and the diabatic tendency from themicrophysical scheme are evaluated explicitly and integrated over the departure cell area. Since thepredictor-corrector and the forcing terms depend only on values at the previous time level, they canbe evaluated along with Θn+1m,exp in a single call to CSLAM, giving Θˆn+1m . Then, to allow for couplingto the momentum equations, a semi-implicit flux-form correction term is used to form the implicit45Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverequationΘ˜n+1m =Θˆn+1m −∆t2[∇eul · (Θˆn+1m vn+1)−∇lag · (Θˆn+1m v˜n+1)], (3.25)where Θ˜n+1m is the value of Θm at the new time level except for a final saturation adjustment that takesplace at the end of the time step to correct the diabatic tendency using the microphysics scheme.The new tendency is then carried over to the next time step to be used as an estimate of the diabaticterm in (3.24).The form of the semi-implicit correction term [square-bracketed terms in (3.25)] stems from thesplit-divergence approximation used in the trajectory computation. The semi-implicit discretizationfor Θn+1m is based on the flux-form scheme presented in Wong et al. (2013) for the height equationin their shallow-water equations solver. The flux-form scheme is based on the derivation of a sim-ilar semi-implicit discretization for the shallow-water model found in Lauritzen et al. (2006), butthe latter scheme uses a time-independent reference state, with which it becomes difficult to en-sure numerical consistency and maintain conservative properties (discussed in section 3.4). Insteadof using a time-independent reference state, we form the semi-implicit correction term using theexplicit solution Θˆn+1m from (3.24).The semi-implicit correction term is defined as the difference between an Eulerian flux diver-gence and a Lagrangian flux divergence. On an Arakawa C-grid, these would be defined as∇eul · (Θmv) =1∆x[(Θmxu)r− (Θmxu)l]+1∆z[(Θmzw)t − (Θmzw)b], (3.26)and∇lag · (Θmv) =1∆x∆z[ΘmxFr−ΘmxFl+ΘmzFt −ΘmzFb], (3.27)respectively, and F(·) are Lagrangian flux areas, where the subscripts r, l, t, b denote the right, left,top, and bottom cell faces of an Eulerian grid cell (Fig. 3.2). We use an exact computation of theLagrangian flux divergence in an Eulerian manner, where Lagrangian flux areas F(·) through eachcell face are defined asFr = urzz∆z− (uc2wc3−uc3wc2)∆t/2, (3.28)46Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverFigure 3.2: Geometric representation of the Lagrangian flux divergence, defined as the flux-area difference between the Eulerian arrival grid cell (solid square) and the departure cell(dashed polygon) in one time step. Velocities associated with the Eulerian grid cell at thecell faces (ul,ur,wt ,wb) and cell vertices (uc,wc)i for i = 1,2,3,4 are also shown. Whitearrows indicate the computed trajectories of each departure grid cell vertex.Fl = ulzz∆z− (uc1wc4−uc4wc1)∆t/2, (3.29)Ft = wtxx∆x− (uc3wc4−uc4wc3)∆t/2, (3.30)Fb = wbxx∆x− (uc2wc1−uc1wc2)∆t/2, (3.31)where the spatial operators are defined as(·)xx=14((·)i−1,k +2(·)i,k +(·)i+1,k), (3.32)(·)zz=14((·)i,k−1 +2(·)i,k +(·)i,k+1). (3.33)The terms proportional to ∆t/2 correct for the geometric differences between the Eulerian andLagrangian flux divergences (shaded areas in Fig. 3.2). [For full details on the derivation of F and∇lag · (Θmv), see Wong et al. (2013)].Using (3.26) and (3.27), the explicit equation for Θˆn+1m (3.24) and implicit equation for Θ˜n+1m47Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver(3.25) can be rewritten asΘˆn+1m =Θn+1m,exp +∆t2[∇eul · (Θnmv′n)]δA∗∆A+∆t[FnΘm]δA∗∆A, (3.34)andΘ˜n+1m =Θˆn+1m −∆t2[∇eul · (Θˆn+1m v′n+1)], (3.35)respectively, where v′ is a corrective velocity and∇eul · (Θmv′) =1∆x[Θxm(ur−Fr/∆z)−Θxm(ul−Fl/∆z)]+1∆z[Θzm(wt −Ft/∆x)−Θzm(wb−Fb/∆x)]. (3.36)3.3.5 Helmholtz equationThe Helmholtz equation with variable coefficients for the semi-implicit problem is solved using aconjugate-residual solver. Substitution of the momentum equations (3.16) and (3.17) into (3.35)forms the Helmholtz equation for Θ˜′n+1m ,−(∆t2)2γRd(1+β )[δx(Θˆn+1mpinρnmxδxΘ˜′n+1m )+δz(Θˆn+1mpinρnmzδzΘ˜′n+1m )]+ Θ˜′n+1m= RΘ−∆t2(1+β )[δx(Θˆn+1mxRu)+δz(Θˆn+1mzRw)]. (3.37)The terms Ru, Rw, and RΘ represent the known terms in (3.16), (3.17), and (3.35). The explicitsolution Θˆn+1m from CSLAM is computed prior to solving (3.37).Using the explicit solution Θˆn+1m allows for a straightforward and consistent formulation betweenthe thermodynamic and continuity equations, as long as the reconstruction of Θm is performed in aconsistent manner. To ensure this, in CSLAM we follow Nair and Lauritzen (2010) in separating thesubgrid-cell reconstructions for ρd and θm, and compute the second-order reconstruction function48Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverΘm(x,z) asΘm(x,z) =〈ρd〉θm(x,z)+〈θm〉(ρd(x,z)−〈ρd〉), (3.38)where〈ρd〉and〈θm〉are Eulerian grid cell values, and ρd(x,z) and θm(x,z) are reconstructionfunctions. To check for consistency, we substitute a field of constant θm, i.e. θm(x,z) =〈θm〉= 1,in (3.38) and see that the rhs of (3.38) properly reduces to ρd(x,z).In summary, the solution procedure for obtaining solutions for Θ˜′n+1m , un+1, and wn+1, is as fol-lows: (i) obtain solution for Θ˜′n+1m by solving the Helmholtz equation (4.26), (ii) substitute solutionfor Θ˜′n+1m into (3.16) and (3.17) to obtain solutions for un+1 and wn+1, respectively, and (iii) recalcu-late Θ˜′n+1m using un+1 and wn+1 to eliminate any roundoff errors. This solution procedure is similarto that used in Wong et al. (2013) for the shallow-water equations.3.3.6 Discretization of the continuity equationWe ensure that the flux-form thermodynamic equation is consistent with the continuity equation byusing the same numerical scheme, with the inclusion of the semi-implicit correction terms in thecontinuity equation. Again, we first use CSLAM to obtain the explicit solution ρˆn+1d in a similarmanner as in (3.34),ρˆn+1d =ρn+1d,exp +∆t2[∇eul · (ρnd v′n)]δA∗∆A. (3.39)Then, we add the semi-implicit correction term to (3.39) to be consistent with (3.35),ρn+1d =ρˆn+1d −∆t2[∇eul · (ρˆn+1d v′n+1)]. (3.40)The new time-level correction term is evaluated by back-substituting the solution of the velocityfield vn+1 into v′n+1.3.3.7 Discretization of moisture conservation equationsThe flux variables of mixing ratios of water vapor qv, cloud water qc, and rainwater qr are in-cluded as prognostic variables in the nonhydrostatic solver. Moist mass conservation equationsare integrated using CSLAM. To ensure moisture conservation, numerical consistency between thecontinuity equation and the moisture conservation equations needs to be ensured. Numerical in-consistency between the continuity equation and other scalar conservation equations can lead tospurious generation or removal of scalar mass, despite using inherently mass-conserving advection49Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverschemes.A consistent formulation of the moisture conservation equations using the scheme in Wong et al.(2013) for the flux variables Q j = ρdq j where q j = (qv,qc,qr) isQˆn+1j = Qn+1j,exp +∆t2[∇eul · (Qnjv′n)]δA∗∆A+∆t[Fnq j]δA∗∆A, (3.41)Q˜n+1j = Qˆn+1j −∆t2[∇eul · (Qˆn+1j v′n+1)](3.42)where v′n, v′n+1 and the computations for ∇eul · (·) are the same as in (3.40). The explicit CSLAMsolution Qˆn+1j (3.41) is computed using a consistent reconstruction as in (3.38). Fq j representssecond-order diffusion with a mixing coefficient same as that for Θm, and any diabatic tendenciesfrom the microphysics.3.3.8 Diabatic processesMicrophysical processes are modelled using the simple warm-rain Kessler parameterization, as de-scribed in Klemp and Wilhelmson (1978). In the evaluation of the thermodynamic and moistureconservation equations, the diabatic forcing is approximated in FΘm and FQ j [(3.34) and (3.41), re-spectively] using the most up-to-date values integrated over the departure cell. These values arethen removed from the solutions prior to calling the Kessler microphysics scheme. The includedmicrophysical processes are (1) condensation of water vapor into cloud-water, (2) autoconversionby diffusion and collection of cloud-water into rain-water, (3) evaporation of cloud-water and rain,and (4) precipitation of rain which is removed when it reaches the surface. These microphysicalprocesses are computed as a final adjustment at the end of the time step, advancing Θ˜n+1m and Q˜n+1jto Θn+1m and Qn+1j in a manner that is consistent with saturation conditions at the new time level.3.3.9 Diagnostic equation of statePressure is a diagnostic variable computed using the equation of state (3.6),p = p0(RdΘmp0)γ(3.43)where p0 is the reference surface pressure set to 100 kPa.50Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver3.3.10 Consistency and shape-preservationIn the CSLAM reconstruction step, we reconstruct Q j using (3.38) described in section 3.3.5 to en-sure consistency. To ensure shape preservation, we follow the two steps as described in Wong et al.(2013). First, we use the simple 2D filter by Barth and Jespersen (1989) that searches for new localminima and maxima in the reconstruction function of a scalar field such as moisture mixing ratio q j,and scales the function if these values exceed those in the neighbouring cell. For chemistry applica-tions, preservation of linear correlations in tracers is important, and it has been found that the limiterpreserves linear correlations between tracers, whereas typically linear correlation is only preservedwhen the limiter is not applied (Lauritzen et al., 2014). Second, to ensure shape-preservation in theflux-divergence terms, we compute the upwind moist species mixing ratio q∗j by first decoupling Q jfrom ρd . Then, the flux divergences are computed by centering density to each of the cell faces, i.e.∇eul ·ρdq jv′ =1∆x[(ρxdq∗ju′)r− (ρxdq∗ju′)l]+1∆z[(ρzdq∗jw′)t − (ρzdq∗jw′)b]. (3.44)The upwind q∗j values are determined using v′.3.4 Desirable properties of CSLAM-NHThe flux-form nonhydrostatic semi-implicit CISL solver CSLAM-NH has six main advantages anddesirable properties: (i) inherently mass-conserving using the conservative semi-Lagrangian trans-port scheme CSLAM, (ii) ensures numerically consistent transport, (iii) independent of a meanreference state, (iv) shape-preserving, and (v) like typical semi-implicit solvers, CSLAM-NH re-quires solving a single linear Helmholtz equation and (vi) a single application of CSLAM at eachtime step.CSLAM-NH uses a formulation of the discretized continuity equation that ensures numeri-cal consistency for a cell-integrated semi-Lagrangian (CISL) solver. In CSLAM-NH, a Helmholtzequation for the potential temperature perturbation is solved. Traditionally, to avoid solving a non-linear Helmholtz equation, the flux divergence term that is coupled to the momentum equations is51Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solveroften first linearized around a mean reference state Θm,ref, e.g.Θn+1m =Θn+1m,exp−∆t2[∇eul · (Θm,refv′n+1)]+∆t2[∇eul · (Θm,refv′n)]δA∗∆A+∆t[FnΘm]δA∗∆A, (3.45)where Θm,ref is a mean reference state that is often time-independent and varies with height. Achoice of reference state can be the hydrostatic background state ρdθ . The scheme (3.45) is anonhydrostatic extension to the SWE semi-implicit CISL continuity equation in Lauritzen et al.(2006).To ensure conservation of potential temperature, it is important for the discrete thermodynamicequation to be numerically consistent with the discrete continuity equation. One can include similarsemi-implicit correction terms as in (3.45) and discretize the continuity equation asρn+1d = ρn+1d,exp−∆t2[∇eul · (ρd,refv′n+1)]+∆t2[∇eul · (ρd,refv′n)]δA∗∆A. (3.46)However, note that the above scheme is only strictly consistent with (3.45) for Θm,ref = ρd,refθm,which is difficult to achieve as θm is unknown.Transported material, such as moisture and passive tracers with some mixing ratio q, are oftensolved explicitly using the CISL transport scheme, i.e.,φ n+1 = φ n+1exp +∆t[Fnφ]δA∗∆A, (3.47)where φ = ρdq is the scalar mass and[Fnφ]represents diffusion and any diabatic tendency evaluatedat time level n over the departure cell. Such explicit schemes would lead to numerical inconsistencybetween the discrete CISL continuity equation (3.46) and the discrete constituent mass conservationequations such as (3.47). If the discrete conservation equation is consistent with the discrete conti-nuity equation, the former should reduce to the latter when q is a constant, and an initially spatiallyuniform passive tracer field should remain so. The inconsistent flux-divergence correction term in(3.46) can spuriously generate or remove moisture or tracer mass in the model.Alternatively, one can formulate the discrete scalar conservation equation in a manner similar52Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverto (3.46) by including the flux-divergence correction terms,φ n+1 = φ n+1exp −∆t2[∇eul · (φrefv′n+1)]+∆t2[∇eul · (φrefv′n)]δA∗∆A+∆t[Fnφ]δA∗∆A. (3.48)However, similar to the case for Θm,ref, determining an appropriate choice for reference state φref isdifficult, making a numerically consistent formulation such as (3.48) hard to implement.The formulations we present for the thermodynamic, density, and moisture conservation equa-tions [(3.35), (3.40), and (3.42), respectively] are all numerically consistent with one another. Theseconsistent formulations are made possible by avoiding the use of a mean reference state. In ourformulation, we use the explicit CSLAM solution instead of a mean reference state in the flux-divergence correction terms. This approach eliminates the difficult choice of a mean reference stateφref for moisture or tracer mass.Even if an appropriate choice of φref can be found, using a time-independent mean referencestate in (3.48) can be problematic for regions with little moisture or tracer mass (φˆ n+1  1). De-pending on the magnitude of φref, the flux divergences are likely nonzero for a divergent flow andcan, therefore, generate or remove unphysical mass (Lauritzen et al., 2008). In the nonhydrostaticsolver presented here, by computing the flux divergences of the explicit solution φˆ n+1, the magni-tude of the flux divergences are better approximated for regions with little moisture or tracer mass.As scalar mass conservation is not guaranteed in an inconsistent solver, these solvers also gen-erally do not preserve the shape of scalar fields such as mixing ratios, even when shape-preservingfilters are applied to the transport scheme. The implications are that the scalar field may no longer bepositive-definite, and new unphysical minima and maxima can occur due to under- and overshoot-ing, respectively. The consistent and shape-preserving transport in the proposed solver ensures thatno new (unphysical) minimum and maximum (within machine roundoff) will occur.3.5 Idealized test casesTwo of the three idealized test cases presented, namely a density current simulation and a gravitywave simulation, are commonly used as benchmarks for testing nonhydrostatic solvers. The thirdidealized test case is a 2D squall line simulation, where the stability of the model is tested withlatent heating modeled by a simple warm-rain microphysics scheme. In addition to comparing withavailable solutions in the literature, comparisons with an Eulerian split-explicit model with 2nd-53Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverorder advection [2D version of the solver described in Klemp et al. (2007)] are also presented.3.5.1 Density currentThe nonlinear density current test case, described in Straka et al. (1992), is widely used as abenchmark test for nonhydrostatic solvers (e.g. Klemp et al., 2007; Xue et al., 2000). An ini-tial cold temperature perturbation is centered in the domain, and the negatively buoyant cold airdescends to the surface and forms symmetric density currents propagating in opposite directions.Straka et al. (1992) provides a well-documented comparison among various compressible and quasi-compressible solvers as well as a high-resolution benchmark solution.The numerical domain is centered at x = 0.0 km, with −25.6 km ≤ x ≤ 25.6 km and 0 ≤ z ≤6.4 km. As described in Straka et al. (1992), the initial condition is given by a temperature pertur-bation ∆T , where∆T =0.0◦C if L≥ 1.0−15.0◦C[cos(piL)+1.0]/2 if L < 1.0,(3.49)where L =√[(x− xc)/xr]2 +[(z− zc)/zr]2 where (xc,zc) = (0.0,3.0) km is the center of the per-turbation, and its width and depth are given by xr = 4.0 km and zr = 2.0 km. The surface temper-ature θ0 is at 300 K in a horizontally homogeneous and neutral environment. A constant physicalviscosity of 75 m2 s−1 is used. The domain is an x-periodic channel with reflective boundaryconditions along the upper and lower boundaries as specified by Straka et al. (1992) that require∂u/∂ z = w = ∂ρ/∂ z = ∂Θ/∂ z = 0.Following Straka et al. (1992), we simulate the density current test case using grid spacings∆x = ∆z = 400,200,100,50,25 m, with Eulerian time step sizes of ∆t = 4,2,1,0.5,and 0.25 s,respectively. Figure 3.3 shows the potential temperature perturbation (θ ′) from its mean state fromCSLAM-NH and the Eulerian split-explicit scheme with 2nd-order advection at the simulation endtime of 15 min using different model resolutions.The density current is clearly under-resolved using a 400 m-grid spacing, with only the mainrotor marginally resolved (7 km ≤ x ≤ 9 km). A grid-spacing of 200 m gives a better resolution ofthe main rotor as well as a second rotor (11 km≤ x≤ 12 km); however the leading third rotor is stillunder-resolved. For resolutions finer than ∆x = ∆z = 100 m, all three rotors are well-resolved withthe solutions converging and indistinguishable by eye between 50 m and 25 m grid spacings. Thedifferences among the model resolutions agree well with those documented in other nonhydrostaticsolvers such as in Straka et al. (1992), Xue et al. (2000), Skamarock and Klemp (2008), and Melvinet al. (2010).Positions of the density current front (specified to be at θ ′= -1 K), the minimum and maximum54Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver Δx = Δz = 400 mz (km)12345Δx = Δz = 400 mΔx = Δz = 200 mz (km)12345Δx = Δz = 200 mΔx = Δz = 100 mz (km)12345Δx = Δz = 100 mΔx = Δz = 50 mz (km)12345Δx = Δz = 50 mCSLAM-NH Eulerian Split-ExplicitΔx = Δz = 25 mx (km)z (km)0 2 4 6 8 10 12 14 1612345−10−8−6−4−2x (km)2 4 6 8 10 12 14 16Δx = Δz = 25 mperturbation θ (K)Figure 3.3: Potential temperature perturbation (K) after 15 min. Contour intervals are every 1K, starting at 0.5 K. Mean wind U = 0 m s−1.θ ′ values in the domain, and ∑θ ′sampled for all θ ′sampled and θ ′sampled > 0, and ∑θ ′2sampled are shown inTable 3.1. We also compare the results with those using SLICE [Table II of Melvin et al. (2010)]and REFC25, the 25 m reference solution in Table IV of Straka et al. (1992). In computing thesummation statistics, θ ′sampled from each of the CSLAM-NH runs (except for the 400 m grid-spacingrun) are sampled at 200 m resolution. This sampling is done so that we can make a direct compar-ison with the statistics of REFC25 in Straka et al. (1992) (where they sampled REFC25 at 200 mresolution). Statistics from the 25 m solution agree closely with the nonhydrostatic SLICE model,with a similar slight discrepancy in the density front location when compared to REFC25. BothCSLAM-NH and SLICE are semi-Lagrangian models with inherent dissipation and order of ac-curacy different from REFC25, an Eulerian compressible solver with 2nd-order advection; thesedifferences could lead to the slight discrepancy in the density front locations. In addition to modeldifferences, like the SLICE model, a different time step size is used to compute the 25 m solution(16 times larger than that used to compute REFC25). At coarser resolutions (400 m and 200 m),55Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverGrid Time θ ′min θ ′max Front ∑θ ′sampled ∑θ ′sampled ∑θ ′2sampledsize (m) step size (s) (K) (K) location (K) (for θ ′ > 0) (K2)(m) (K)400 4 -10.339 0.6804 14248 — — —200 2 -10.746 0.0846 14938 -1293.82 4.4398×10−1 5634.92100 1 -9.7694 0.0006 15234 -1361.41 1.8114×10−4 6127.90100 4 -9.6985 0.0053 15256 -1360.73 6.7741×10−3 6182.0350 0.5 -9.7078 0.0000 15360 -1394.93 2.0562×10−5 6395.6325 0.25 -9.7323 0.0000 15391 -1411.62 3.2974×10−8 6516.33SLICE400 4 -5.6608 0.3674 13572 — — —SLICE200 2 -8.0958 0.1226 14768 — — —SLICE100 1 -9.8574 0.0995 15182 — — —SLICE50 0.5 -9.4995 0.0626 15334 — — —SLICE25 0.25 -9.6548 0.0048 15390 — — —REFC25 1.5625×10−2 -9.7738 0.0000 15537 -1427.10 0.0000 6613.62Table 3.1: Statistics for the density current simulations at time 15 min using CSLAM-NH atvarious grid resolutions and time step sizes. Comparison values from the nonhydrostaticx-z solver using SLICE in Melvin et al. (2010) are also presented. REFC25 are valuestaken from Straka et al. (1992). θ ′sampled are solutions sampled at 200 m for comparisonwith values in Straka et al. (1992).the minimum θ ′ values are colder than those in SLICE; the front locations therefore also travelledfarther out from the centerline. Analytically, the maximum θ ′ should remain zero throughout thesimulation, as is the case in the higher resolution runs (50 m and 25 m). The contribution of positiveθ ′ values in ∑θ ′sampled is also small at these resolutions (in the order of 1×10−5 K and 1×10−8 Krespectively), increasing up to the order of 1×10−1 K at 200 m. [Straka et al. (1992) only reportedvalues up to 4 decimal points.]For the next simulation, mean background wind of U = 20 m s−1 is applied to the describedtest case, as is done in Skamarock and Klemp (2008) to examine phase errors. Solutions fromCSLAM-NH and the Eulerian split-explicit 2nd-order advection scheme of both the left- and right-moving currents at time 15 min using ∆x = ∆z = 200, 100, and 50 m are shown in Fig. 3.4. Timestep sizes are the same as in Fig. 3.3. The solutions from CSLAM-NH in general show propersymmetry about the translating centerline, although very subtle differences between the secondaryrotors in the left- and right-moving currents are noticeable at 200 m and 100 m grid spacing. As acomparison, the Eulerian split-explicit 2nd-order advection scheme shows noticeably larger errorsin the right-moving current as expected due to the right-moving current moving at a greater speed56Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverCSLAM-NH Eulerian Split-Explicit−10−8−6−4−2perturbation θ (K)Δx = Δz = 200 mz (km)12345Δx = Δz = 100 mz (km)12345Δx = Δz = 50 mx (km)z (km)−16 −14 −12 −10 −8 −6 −4 −2 0 2 4 6 8 10 12 14 1612345 Δx = Δz = 50 mΔx = Δz = 200 mΔx = Δz = 100 m x (km)−16 −14 −12 −10 −8 −6 −4 −2 0 2 4 6 8 10 12 14 16Figure 3.4: Potential temperature perturbation (K) after 15 min. Contoured as in Fig. 3.3.Solution is translated using a mean wind U = 20 m s−1.than the other (causing larger advective phase errors).For this test case, we found that the maximum stable time step size in CSLAM-NH is double ofthat of the Eulerian scheme. Fig. 3.5 shows solutions for tests where U = 0 m s−1 at ∆x = ∆z = 100m using a time step size of 3 s and 4 s, whereas the maximum stable Eulerian time step size is∆t = 2 s. The solution using a large time step of 4 s is almost indistinguishable by eye from the 25m high-resolution solutions (Fig. 3.3). With mean advection (U = 20 m s−1), the maximum stabletime step in CSLAM-NH is 3 s. As we increase the time step size to 4 s, the phase error was largeenough to form unphysically steep gradients at the leading edge of the right-moving current, whichthen caused the violation of the Lipschitz stability condition. The maximum stable time step in theEulerian model is 1.5 s. Using a time step size of 5 s, instability was observed in the vicinity of theleading edge of the subsiding cold air for both cases with and without the mean wind.3.5.2 Gravity waveA second test case of a gravity wave in a periodic channel with solid, free-slip upper- and lower-boundary conditions is used to evaluate the nonhydrostatic solver. The test case is described inSkamarock and Klemp (1994), where they presented results for a Boussinesq atmosphere. The test57Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver z (km)12345x (km)z (km)0 2 4 6 8 10 12 14 1612345Δt = 3 sΔt = 4 s−10−8−6−4−2perturbation θ (K)Figure 3.5: Potential temperature perturbation (K) from CSLAM-NH after 15 min for gridspacing ∆x = ∆z = 100 m using time step sizes ∆t = 3 and 4 s. Mean wind U = 0 m s−1.The Eulerian split-explicit scheme (not plotted) was numerically unstable for these timesteps, as it required ∆t ≤ 2 s for numerical stability of this gravity current. Contoured asin Fig. 3.3.case is characterized by an initial potential temperature perturbation of amplitude ∆θ0,θ(x,z, t = 0) = ∆θ0sin(piz/H)1+(x− xc)2/a2. (3.50)where ∆θ0 = 10−2 K, a = 5 km is the half-width of the initial perturbation, H = 10 km is thetotal depth of the domain, and xc = 0.25L, where L = 300 km is the length of the domain. Thebackground atmospheric stratification has a constant Brunt-Va¨isa¨la frequency N = 10−2 s−1. Forone simulation, no mean wind (U = 0) is prescribed. The other simulation uses a mean wind ofU = 20 m s−1, advecting the solution to the right while the two gravity wave modes propagate inopposite directions. Again, the mean advection of the solution accentuates any advective phasespeed errors in the scheme. The boundary condition is implemented by linear extrapolating u,Θ,and ρ values into the boundary, consistent with the free-slip boundary conditions, and setting w = 0.We run the gravity wave test case at grid spacings ∆x = ∆z = 1 km, 500 m, and 250 m usingEulerian time step sizes ∆t = 12 s, 6 s, and 3 s, respectively. Solutions from CSLAM-NH atthe three resolutions for U = 0 (not shown) are indistinguishable by eye from the 250 m and 500m solutions for U = 20 m s−1 in Fig. 3.6 and compare well with those using the WRF-ARWmodel (solutions using the 5th- and 6th-order advection scheme are available at http://www.mmm.ucar.edu/projects/srnwp tests/IG waves/ig wave.html), with the 2nd-order advection scheme of thesame Eulerian split-explicit scheme, and with the SLICE nonhydrostatic vertical model in Melvin58Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverΔx = Δz = 1 kmz (km)2468Δx = Δz = 1 kmΔx = Δz = 500 mz (km)2468Δx = Δz = 500 mΔx = Δz = 250 mx (km)z (km)−150 −100 −50 0 50 100 1502468x (km)Δx = Δz = 250 m−150 −100 −50 0 50 100 150CSLAM-NH Eulerian Split-Explicit−1 0 1 2 x 10-3perturbation θ (K) Figure 3.6: Potential temperature perturbation (K) after 50 min. Contour intervals are every5× 10−4 K (zero contour line not plotted). Solid lines indicate positive contours anddashed lines indicate negative contours. Solution is translated using a mean wind U =20 m s−1. Horizontal axis has also been translated with the mean wind so the line ofsymmetry remains at x = 0.et al. (2010). In Skamarock and Klemp (1994), the solution presented for this nonhydrostatic testcase uses a Boussinesq model, where the symmetry of the analytic Boussinesq solution in both xand z is maintained. The density variation in the full Euler equations results in solutions that areasymmetric in z, as observed in the CSLAM-NH solutions, the 2nd-order Eulerian solutions, the5th-order Eulerian solutions, as well as the SLICE nonhydrostatic vertical model solutions.Like in the density current test, we impose a mean advection wind U = 20 m s−1 to examinephase errors. These tests are made at the same grid spacings and time step sizes as in the no meanwind case. The right- and left-moving waves from CSLAM-NH exhibit nearly perfect symme-59Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver(a)Δt = 60 sz (km)2468 (d)Δt = 24 s(b)Δt = 75 sz (km)2468 (e)Δt = 30 s(c)Δt = 100 sx (km)z (km)−150 −100 −50 0 50 100 1502468 (f)Δt = 40 sx (km)−150 −100 −50 0 50 100 15043perturbation θ (K)−1 0 1 2 x 10-3U = 0 m s-1 U = 20 m s-1Figure 3.7: Potential temperature perturbation (K) solutions of the gravity wave case usingincreasingly large CSLAM-NH time steps (∆x = ∆z = 1 km) where (a)-(c) U = 0 m s−1and (d)-(f) U = 20 m s −1. Contoured as in Fig. 3.6.try, indicating there is minimal phase error in the solutions. The Eulerian split-explicit 2nd-orderadvection scheme shows more noticeable phase errors (Fig. 3.6).Testing of CSLAM-NH using larger time steps in this gravity wave test case reveals a numericalstability condition that is sensitive to the stratification N. (We note that CSLAM-NH is uncondi-tionally stable for N = 0, i.e. for a near-pure advection case of the initial warm perturbation.) Weevaluate the maximum stable CSLAM-NH time step size for the gravity wave case with a meanadvection wind speed of U = 20 m s−1 (∆x = ∆z = 1 km) over a range of N (0.01, 0.015, and 0.02s−1). Since the gravity wave phase speed varies with N, we increase/decrease the simulation time60Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverlength as appropriate such that the gravity wave solutions are similar to those shown in Fig. 3.6;for example, for N = 0.015 s−1, the simulation time is reduced to 2000 s. Test results showed thatthe maximum stable CSLAM-NH time step sizes are ∆tmax = 38,35, and 32 s for N = 0.01,0.015,and 0.02 s−1, respectively, whereas in the case of the Eulerian split-explicit scheme, the maximumstable large time steps are found to be ∆t = 60,55, and 50 s (with small time step size of 2.4 s),respectively, limited by the stability condition of the advection scheme. The buoyancy terms inthe vertical momentum equation are integrated explicitly in CSLAM-NH, and handled implicitly inthe Eulerian scheme. When we remove the buoyancy terms from the implicit step and solve themexplicitly in the Eulerian model, the time step sizes required to obtain solutions of similar accuracyas those from the vertically implicit model are reduced by 20–35%, and are closer to those foundin CSLAM-NH. The devising of an integration scheme that handles the buoyancy terms implicitlyin CSLAM-NH will require a robust and stable way of updating the density perturbation in theHelmholtz solver, and this will be addressed in future work.3.5.3 2D (x-z) squall lineWe perform a test case of a 2D squall line as described in Weisman and Klemp (1982) to evaluatemass conservation, consistency, and shape-preservation in the nonhydrostatic solver, in addition totesting for any small-scale computational instability in the model due to latent heating.The numerical domain is centered at x = 0.0 km, with −100 km ≤ x ≤ 100 km and 0 ≤ z ≤20 km. As in Weisman and Klemp (1982), a conditionally unstable thermodynamic profile is usedto initialize the horizontally homogeneous environment. Constant physical horizontal and verticaleddy viscosities of 250 m2 s−1 are used. A warm thermal perturbation near the surface is prescribedto initiate convection (Weisman et al., 1988). The initial thermal perturbation has a maximum of∆θ0 = 3 K, centered at zc = 1.5 km and along the centerline (xc = 0) of the domain, with a horizontalradius xr of 10 km and a vertical radius zr of 1.5 km. The shape of the perturbation is a cosine hillgiven asθ(x,z, t = 0) =∆θ0 cos2(piL/2) ,L < 1.0,0 ,L≥ 1.0,(3.51)where L =√(x/xr)2 +[(z− zc)/zr]2.A weak vertical wind shear within a 2.5 km-layer at the surface is used to promote the growthof the squall line. The initial wind profile is given asu(z, t = 0) =u · (z/zts)−us ,z < zts,u−us ,z≥ zts,(3.52)61Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverwhere u = 12 m s−1, us = 10 m s−1, and zts = 2.5 km. The environmental potential temperature andrelative humidity profiles at the initial time areθ(z, t = 0) =θ0 +(θtr−θ0)(z/ztr)54 ,z≤ ztr,θtr exp[ gcpTtr(z− ztr)],z > ztr,(3.53)andRH(z, t = 0) =1− 34(z/ztr)54 ,z≤ ztr,0.25 ,z > ztr,(3.54)where θtr = 343 K, ztr = 12.0 km, and Ttr = 213 K are the potential temperature, geometric height,and actual temperature at the tropopause. The maximum water mixing ratio is capped at 14 g kg−1.The surface potential temperature θ0 = 300 K. The skew-T log-p diagram for this sounding can befound in Fig. 1 of Weisman and Klemp (1982). Numerical simulations (unless otherwise stated) usea grid spacing ∆x = ∆z = 500 m, a time step ∆t = 5 s, and a time-off-centering parameter of β = 0.1to maintain numerical stability. Like the gravity-wave case, the boundary condition is implementedby linear extrapolating u,Θ, and ρ values into the boundary and setting w = 0, consistent with thefree-slip boundary conditions.A comparison of the squall line development among CSLAM-NH (with shape preservation),the 5th-order split-explicit, and the 2nd-order split-explicit Eulerian models is presented in Fig.3.8. Instantaneous and accumulated surface precipitation integrated across the model domain arepresented in Fig. 3.9; also shown is the rate of condensation over the entire domain. Maximumupdraft velocity is shown in Fig. 3.10. The series of updraft velocity peaks highlight the continuoustriggering of new convective systems along the advancing front.All three models (CSLAM-NH, Eulerian 5th-order advection, and Eulerian 2nd-order advec-tion) show similar development of the convective system (Fig. 3.8). At 0.6 h, all three modelsshow an initial downshear orientation of the system due to the ambient wind shear. As the stormcontinues to develop with the cold pool strengthening behind the system (not shown), convergenceand enhanced uplift lead to the storm tilting in a near-upright position (T = 0.8 h). At 1.3 h, a newcell is triggered near the edge of the cold pool, where uplift of the warm moist air in the boundarylayer is enhanced. At 1.7 h, the cold pool is strong enough to generate a circulation such that thesystem develops an upshear orientation, as described in Rotunno et al. (1988). Comparing to thesimulations from the Eulerian 2nd-order model, those from CSLAM-NH show closer resemblanceto those from the Eulerian 5th-order model. The better agreement is also evident in the moisturestatistics (Fig. 3.9), especially in the accumulated surface precipitation amounts and condensationrate in the domain.62Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver x (km)−10 0 10 20 30 40x (km) −10 0 10 20 30 40x (km) −10 0 10 20 30 40T = 0.6 h T = 0.6 h T = 0.8 h T = 0.8 h T = 1.3 h T = 1.3 h T = 1.7 h T = 1.7 h−20 −10 0 10 20 30vertical velocity (m s-1)CSLAM-NH Eulerian 5th order Eulerian 2nd order T = 0.6 hT = 0.8 h T = 1.3 h T = 1.7 hheight (km)2468101214height (km)2468101214height (km)2468101214height (km)2468101214Figure 3.8: Vertical cross-sections of vertical velocity (color shading in m s−1) and solid con-tour of the convective cloud structure (qc = 0.1 g kg−1) at times 0.6, 0.8, 1.3, 1.7 h of thesimulation for the 500 m grid-spacing runs with a time step of 5.0 s from (left) CSLAM-NH, (middle) 5th-order split-explicit Eulerian model, and (right) 2nd-order split-explicitEulerian model.Focussing on the two models that show more comparable results, the first maximum updraftvelocities from CSLAM-NH (34.1 m s−1) is slightly greater than that from Eulerian 5th-order ad-vection (31.6 m s−1) (Fig. 3.10). CSLAM-NH appears to show a weaker second peak updraftvelocity (21.9 m s−1) than the Eulerian 5th-order model (28.3 m s−1); however, the stronger firstpeak (∼ 34 m s−1) and weaker second peak (∼ 25 m s−1) are also observed in a higher-resolutionsimulation using the Eulerian 5th-order model at a grid spacing of 250 m and large time step sizeof 2.5 s (dashed black line in Fig. 3.10). For comparison, maximum updraft from CSLAM-NH at63Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver0100200300Precipitation rate(kg s−1 ) Eulerian 5th−orderEulerian 2nd−orderCSLAM−NH0246810Accumulated precip.(105 kg)0 60 1200200400600Time (min)Condensation rate(kg s−1 )Figure 3.9: Moisture statistics including surface precipitation rate (kg s−1), accumulated sur-face precipitation (kg), and condensation rate (kg s−1) from the microphysics usingCSLAM-NH, Eulerian 5th-order horizontal advection, and Eulerian 2nd-order horizon-tal advection at ∆x = ∆z = 500 m.∆x = 250 m and ∆t = 2.5 s (red dashed line in Fig. 3.10) is also shown, and at the higher resolution,the two models agree very well with each other.The maximum stable time step in the Eulerian split-explicit 5th-order advection scheme is alarge time step of 20 s and acoustic time step size of 1.25 s. The maximum CSLAM-NH stabletime step is limited to 15 s due to the violation of the Lipschitz stability condition in the vicinityof the updraft when a larger time step is used (the instability occurs when the storm reaches itsfirst maximum vertical updraft, which generates a strong horizontal wind shear between the updraft64Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver 0 10 20 30 40 50 600Time (min)5101520253035Maximum updraft (m s-1)EulerianCSLAM−NHEulerianCSLAM-NHΔx=500 m, Δt=5 sΔx=250 m, Δt=2.5 sFigure 3.10: Updraft intensities using CSLAM-NH (red) and Eulerian 5th-order horizontaladvection (black) at ∆x = ∆z = 500 m and ∆t = 5 s (solid), and ∆x = ∆z = 250 mand ∆t = 2.5 s (dashed).and the neighbouring downdraft). In Fig. 3.11, we see at larger time step sizes, maximum updraftvelocities remain close to the small time-step solutions.With the 2D squall-line test case, we examine the shape-preservation properties of CSLAM-NHusing the shape-preserving scheme by Barth and Jespersen (1989) in the CSLAM transport schemeand the upwind scheme for the flux-correction terms in the transport equations. An analogousimplementation of these schemes for a shallow-water model is described in Wong et al. (2013).To verify that consistency is achieved, an additional passive tracer with mixing ratio r is intro-duced into the model. The passive tracer initially has a constant mixing ratio of r0 = 1.0 g kg−1 andwe form the discretized conservation equation as in (3.42). The minimum and maximum values ofr are maintained at 1.0 g kg−1 (up to machine roundoff of order 10−14) throughout the simulationusing the consistent formulation in CSLAM-NH.For a passive tracer that uses an inconsistent discrete conservation equation such as (3.47),unphysical minima and maxima of the passive tracer mixing ratio are generated (Fig. 3.12). At theend of the squall line simulation at 2 h, the minimum and maximum mixing ratios r are 0.986 gkg−1 and 1.021 g kg−1, respectively (i.e. the error is on the order of 1 part in 100). We note thatthe shape-preserving limiter described in Barth and Jespersen (1989) was also applied in CSLAMin this test. Due to numerical inconsistency, however, the limiter becomes ineffective agreeing withthe results in Wong et al. (2013). This discrepancy from constancy highlights the importance ofensuring numerical consistency to properly maintain conservation of moisture and tracer mass in a65Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver 0 10 20 30 40 50 60Time (min)05101520253035Maximum updraft (m s-1)Δt = 2.5 sΔt = 5 sΔt = 10 sΔt = 15 sΔt = 2.5 sΔt = 5 sΔt = 10 sΔt = 15 sEulerianCSLAM-NHFigure 3.11: Timing and intensity of the maximum vertical updraft using ∆x = ∆z = 500 mat different CSLAM-NH time step sizes (solid lines), as compared to the Eulerian 5th-order horizontal advection vertical velocity (dashed lines). (Only first hour is plotted.)semi-implicit CISL solver.3.6 SummaryA new cell-integrated semi-Lagrangian (CISL) nonhydrostatic atmospheric solver, CSLAM-NH,for the moist Euler equations is introduced in this paper. The two-dimensional (x-z) Cartesiannonhydrostatic solver uses a CISL transport scheme, CSLAM, for conservative transport. It alsoincorporates a new approach to ensure numerical consistency among the CISL continuity equationand the conservation equations for potential temperature, moisture species, and passive tracers. Asemi-implicit time integration scheme is used to stably handle the fast-moving acoustic waves in thecompressible system.Based on a recently tested shallow-water equations solver, the extended nonhydrostatic atmo-spheric solver presented here, CSLAM-NH, possesses a number of features ideal for weather andclimate modelling purposes. The solver:1. is inherently mass-conserving through the use of a conservative transport scheme CSLAM,2. ensures numerical consistency between the continuity equation and other scalar mass conser-vation equations (the lack of which may lead to violation of scalar mass conservation),3. does not depend on a mean reference state,66Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solverheight (km)2468101214x (km)height (km)−10 0 10 20 30 40 502468101214x (km)−10 0 10 20 30 40 50−3 −2 −1 0 1 2 3mixing ratio (10-3 g kg-1) T = 0.6 h T = 0.8 hT = 1.3 h T = 1.7 hFigure 3.12: Mixing ratio errors (g kg −1) due to numerical inconsistency associated with(3.47). The passive tracer is initialized with a uniform mixing ratio field of 1.0 g kg−1.The consistent formulation in CSLAM-NH (which does not use (3.47)) ensures mixingratio constancy of the same passive tracer up to machine roundoff of order 10−14 (notshown).4. can be easily implemented with existing shape-preserving filters to ensure shape-preservationof scalar fields,5. requires a single linear Helmholtz equation solution (as in typical semi-implicit solvers) pertime step [for an explicit treatment of the gravity wave terms]3, and6. requires a single application of CSLAM per time step.Here, we tested the nonhydrostatic extension for three idealized test cases: a density current,a gravity wave, and a squall line. To represent microphysical processes in the squall line test,the Kessler warm-rain microphysics parameterization scheme is coupled to the dynamics. The 2Dsolver currently does not admit flow in the y-direction, and therefore, Coriolis terms are neglected;however, the tests we present allow for sufficient testing of typical meteorological flows. Resultscompare well with other existing Eulerian (such as WRF-ARW) and nonhydrostatic CISL solvers3not in manuscript67Chapter 3: Extension to a two-dimensional (x-z) nonhydrostatic atmospheric solver(such as the nonhydrostatic SLICE model). In the density current and gravity wave tests, we seethat CSLAM-NH allows for stable time steps two times larger than that in an Eulerian model. Inthe highly-nonhydrostatic flow of the squall line test case, the maximum stable time step size isof similar magnitude as the Eulerian split-explicit model. The strong wind shear across the stormupdraft imposes a time step limit in CSLAM-NH due to the Lipschitz stability condition (violationof which leads to the crossing of trajectories).Plans to extend the nonhydrostatic solver to include orographic influences are also underway.This work involves transformation of the nonhydrostatic equations into a terrain-following heightcoordinate. In traditional semi-Lagrangian semi-implicit solvers, flow over topography has beenfound to trigger spurious resonances and time off-centering in the implicit scheme has been foundto eliminate these noises. Thus far, without orography, we have found that our nonhydrostatic solveronly requires time off-centering (β = 0.1) in the squall line case to maintain numerical stability. Thenonhydrostatic solver with orography will allow us to test the conservative and consistent transportand stability of the new semi-implicit CISL discretization under the influence of a terrain-followingcoordinate.68Chapter 4Extension to a two-dimensional (x-z)nonhydrostatic solver with transformedvertical coordinates4.1 IntroductionSemi-Lagrangian semi-implicit (SLSI) schemes have been widely used in climate and numericalweather prediction (NWP) models since the pioneering work of Robert (1981) and Robert et al.(1985). The fully compressible nonhydrostatic equations permit fast-moving waves which limit themodel time step size. The combination of a semi-Lagrangian advection scheme with semi-implicittreatment of these waves allows for larger stable time steps, and therefore, increased computationalefficiency. Conservative semi-Lagrangian advection schemes, also known as cell-integrated semi-Lagrangian (CISL) transport schemes, are finite-volume methods that inherently conserve mass bytracking individual grid cells each time step (Rancic, 1992; Laprise and Plante, 1995; Machenhauerand Olk, 1997; Zerroukat et al., 2002; Nair and Machenhauer, 2002; Lauritzen et al., 2010). CISLtransport schemes allow for locally (and thus globally) conservative transport of total fluid mass(such as dry air in the atmosphere) and constituent (i.e., moisture and tracer) mass.However, some CISL schemes lack consistency between the numerical representation of thetotal dry air mass conservation, which we will refer to as the continuity equation, and constituentmass conservation equations (Jo¨ckel et al., 2001; Zhang et al., 2008; Wong et al., 2013). Numericalconsistency in the discrete tracer conservation equation requires the equation for a constant tracerfield to correspond numerically to the discrete mass continuity equation; this consistency ensuresthat an initially spatially uniform passive tracer field will remain so.To allow for large advection time steps, Lauritzen et al. (2010) developed a CISL transportscheme called the CSLAM (‘Conservative Semi-LAgrangian Multi-tracer’) transport scheme. TheCSLAM scheme has recently been implemented in the National Center for Atmospheric Research69Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates(NCAR) High-Order Methods Modeling Environment (HOMME) and was found to be an efficientand highly scalable transport scheme for atmospheric tracers (Erath et al., 2012). To ensure consis-tent numerical representations of the continuity equation and other scalar conservation equations,Wong et al. (2014) proposed a new discretization of the semi-implicit CISL continuity equationusing CSLAM. They showed that the new formulation can be straightforwardly extended to thescalar conservation equations in a fully consistent manner. Wong et al. (2013) also showed that anydiscrepancy between the numerical schemes can lead to spurious generation or removal of scalarmass. We refer to this nonhydrostatic atmospheric solver with conservative and consistent transportas CSLAM-NH.Idealized 2D benchmark test cases for a density current, gravity wave, as well as a squall line,using CSLAM-NH have been performed in Wong et al. (2014). These test cases used flat bottomboundary conditions for simplicity. In the real atmosphere, the bottom fluid boundary is often notflat. Mountains act as a stationary forcing and generate horizontally and vertically propagating in-ternal gravity waves in the atmosphere. Under certain atmospheric conditions they can also inducehighly nonlinear flows such as wave amplification and breaking. Numerical simulations of thesemountain waves have been extensively studied by many (e.g. Klemp and Lilly, 1978; Peltier andClark, 1982; Durran and Klemp, 1983; Durran, 1986) and several of these cases have become bench-mark tests in model development and intercomparison studies (e.g. Pinty et al., 1995; Bonaventura,2000; Xue et al., 2000; Doyle et al., 2000; Melvin et al., 2010). To further develop CSLAM-NH as aviable nonhydrostatic atmospheric solver, we have incorporated orography into the model and haveconducted a suite of these mountain-wave cases documented in the literature. The test suite includeslinear hydrostatic and nonhydrostatic dry mountain waves, a highly nonlinear dry mountain wavewith amplification and overturning of the waves, and a moist mountain flow with cloud and rainformation.The paper is organized as follows. A model description of CSLAM-NH is given in section 4.2.In section 4.3, simulations from the suite of idealized mountain wave tests are presented. Finally, asummary is given in section 4.4.4.2 Model description4.2.1 Governing equationsThe major modification to the model prognostic equations described in Wong et al. (2014) is thetransformation of the vertical coordinate from geometric height to a terrain-following height coor-dinate. In addition to this modification, we have also included the treatment of the gravity-wave70Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesterms in the implicit solver. The previous version of CSLAM-NH solves the buoyancy terms in thevertical momentum equation explicitly using a two time-level extrapolation scheme. For a gravity-wave test originally proposed in Skamarock and Klemp (1994), the time-step limit was found to berestricted by the explicit treatment of these buoyancy terms (Wong et al., 2014). To circumvent thistime-step restriction, an iterative approach is used to include these terms in the implicit solve. Wewill focus on the description of these two modifications and provide a basic description of the solver[readers are referred to Wong et al. (2014) for a more detailed description].A height-based coordinate is used to avoid the complication of a time-varying vertical coor-dinate system, as is the case with mass (pressure) coordinates. The use of terrain-following co-ordinates substantially simplifies the bottom boundary condition when topography is present. Forcell-integrated semi-Lagrangian advection, in a geometric height coordinate, approximated depar-ture cell boundaries may intersect the orography and create more complex cell configurations (e.g.more cell edges/vertices, which complicate the subgrid-cell reconstruction). On a computationalgrid defined by terrain-following vertical coordinates, however, the lowest cell boundaries will al-ways remain at the surface.Following Gal-Chen and Somerville (1975), the 2D terrain-following height coordinate ζ isexpressed using the linear transformationζ = ztz−h(x)zt −h(x), (4.1)where z(x,ζ ) is the physical height, h(x) is the terrain profile, and zt is the height of the model top[with the bottom defined as z = h(x)].The 2D governing equations expressed in (x,ζ )-coordinates are:∂u∂ t +(u∂u∂x)ζ+ω ∂u∂ζ =−1ρ˜mγRdpi[(∂Θ′m∂x)ζ+∂ (ζxΘ′m)∂ζ]+Fu, (4.2)∂w∂ t +(u∂w∂x)ζ+ω ∂w∂ζ =−1ρ˜m[γRdpi∂ (ζzΘ′m)∂ζ −gρ˜dpi ′pi +gρ˜′m]+Fw, (4.3)∂Θm∂ t +(∇ ·vΘm)ζ = FΘ, (4.4)∂ ρ˜d∂ t +(∇ ·vρ˜d)ζ = 0, (4.5)∂Q j∂ t +(∇ ·vQ j)ζ = FQ j (4.6)71Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesp = p0(RdΘmp0)γ, (4.7)where v = (u,w) is the horizontal and vertical wind components, pi = (p/p0)κ is the Exner function,p0 = 100 kPa is the reference pressure, Rd = 287 J kg−1 K−1 is the gas constant for dry air, cp =1003 J kg−1 K−1 is the specific heat for dry air at constant pressure, cv = 717 J kg−1 K−1 is thespecific heat of dry air at constant volume, and the ratios κ = Rd/cp ≈ 0.286 and γ = cp/cv ≈1.4. Perturbation variables from a time-independent hydrostatically balanced background state areused to reduce numerical errors in the calculations of the pressure-gradient terms (Klemp et al.,2007). The hydrostatically balanced background state is defined as dp(z)/dz =−ρd(z)g. Flux-formvariables are coupled to a scaled dry density adjusted to the transformed coordinate, ρ˜d = ρd/ζz, i.e.Θm = ρ˜dθm and Q j = ρ˜dq j. The notations (ζx,ζz) refer to the spatial derivatives of ζ . Perturbationvariables (primed) are defined via Θm = ρd(z)θ(z)+Θ′m, pi = pi+pi ′, ρd = ρd(z)+ρ ′d , and the moistdensity ρ˜m = ρ˜d(1+qv +qc +qr), where qv, qc, and qr are the mixing ratios for water vapor, cloud,and rainwater, respectively. The modified potential temperature θm is defined as θm = θ(1+ a′qv)where a′ ≡ Rv/Rd ' 1.61.Notations (·)ζ denote evaluation at constant ζ , and (∇ ·vb)ζ = δx(ub)+δζ (ωb) for any scalarvariable b. The variable ω = dζ/dt is the vertical motion perpendicular to the coordinate surface.For simplicity, we assume a non-rotating atmosphere. The terms Fu and Fw represent diffusion, andFΘ and FQ j represent diffusion as well as any diabatic source terms from parameterized physics.The governing equations used in CSLAM-NH include the advective form of the momentumequations [(4.2), (4.3)] so that we can use a traditional semi-Lagrangian discretization. The flux-form advection for potential temperature, density, and moisture/passive scalar variables [(4.4), (4.5),(4.6), respectively] are solved using the conservative semi-Lagrangian scheme CSLAM. Pressureis a diagnostic variable given by the equation of state (4.7). Following Klemp et al. (2007), thepressure-gradient terms are written in terms of potential temperature. The recasting allows for cou-pling of the implicit pressure-gradient terms with the flux divergence term in the potential temper-ature equation. The compressible nonhydrostatic equation set is still exact and no approximationshave been applied.4.2.2 CSLAM — a cell-integrated semi-Lagrangian transport schemeWhen advection terms are evaluated using an Eulerian scheme, the model time step sizes are re-stricted by the well-known Courant stability condition. To allow for larger advective time steps, thenonhydrostatic solver uses a cell-integrated semi-Lagrangian (CISL) transport scheme called theCSLAM (Conservative Semi-LAgrangian Multi-tracer) transport scheme developed by Lauritzenet al. (2010). This inherently conservative (both locally and globally) transport scheme is used to72Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates(a) Non-self-intersecting departure cell (b) Self-intersecting departure cellzxFigure 4.1: Discrete departure cells in CSLAM-NH are approximated using straight edges(shaded in grey). The departure cell vertices (black circles) are computed usingbackward-in-time trajectories (arrows) from the vertices (white circles) of the Eulerianarrival grid cell (white box). The CSLAM transport scheme is stable as long as the dis-crete departure grid cells are (a) non-self-intersecting, and becomes problematic if (b)the departure cell self-intersects since the scheme is no longer mass-conserving.solve the continuity and potential temperature equations, and for transport of any moist species orother tracers.The stability criterion for the CSLAM transport scheme is limited by the trajectory approxi-mations of the grid-cell vertices. To ensure stability in traditional semi-Lagrangian schemes, theLipschitz stability condition requires that, in 1D, no trajectories in the space-time domain shouldintersect one another (Smolarkiewicz and Pudykiewicz, 1992). In the CSLAM scheme, the stabilitycondition is slightly different in that the trajectories of neighbouring vertices may cross, as longas the discrete departure cells remain non-self-intersecting. In all test cases presented here, lineartrajectories as described in Wong et al. (2014) are assumed. Figure 4.1a shows a discrete arrivalgrid cell (white box) originating from a non-self-intersecting discrete departure cell (grey box) withstraight edges that are computed using the approximated displacement over one time step (arrows).As long as the departure cells remain non-self-intersecting, the scheme is stable and ensures globalmass conservation. In Fig. 4.1b, a more distorted flow causes the departure cell edges to inter-sect. This ‘twisting’ of the departure cell causes adjacent departure cells to overlap. In such a case,the scheme is no longer mass-conserving and becomes unstable. The stability and accuracy of the73Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesCSLAM scheme in highly deformed flows may be improved by using higher-order trajectory ap-proximations and/or higher-order approximations of departure cell boundaries. One such exampleis to use the parabolic (curved) departure cell edges that account for acceleration in the trajectoryapproximations developed by Ullrich et al. (2012). In the present study, we did not test any geo-metrical definitions other than quadrilateral departure cells, but the option could be explored in thefuture.4.2.3 Discretized momentum equationsThe momentum equations are solved in a traditional semi-Lagrangian semi-implicit manner, wherethe total derivatives du/dt and dw/dt are computed using a grid-point interpolation to the departurepoint. Bicubic Lagrange interpolation is used for all departure point evaluations. The two time-leveldiscretizations of the momentum equations areun+1A = unD +∆t(Fu)nD−∆t2{γRdpiρ˜mx[(δxΘ′m)ζ +δζ (ζxΘ′mxζ)]}nD−∆t2{γRdpinρ˜nmx[(δxΘ′m)ζ +δζ (ζxΘ′mxζ)]n+1A},(4.8)and(1+µ∆t)wn+1A = wnD +∆t(Fw)nD−∆t2{γRdpiρ˜mζ[δζ (ζzΘ′m)]−1ρ˜m[gρ˜dpi ′pi −gρ˜′m]ζ}nD−∆t2{γRdpinρ˜nmζ[δζ (ζzΘ′m)]n+1A−1ρ˜nm[gρ˜dpi ′pi −gρ˜′m]n+1Aζ},(4.9)where subscripts D and A denote evaluation at the departure and arrival grid points, respectively,and superscripts denote the time level. To reduce gravity wave reflection at the upper boundary,a Rayleigh damping term −µw is added to the vertical momentum equation, where the dampingcoefficient µ is a function of height z and applied in the top layers of the domain. The spatialaveraging operators are defined as(·)x=12((·)i,k +(·)i+1,k), and (4.10)74Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates(·)ζ=12((·)i,k +(·)i,k+1), (4.11)and gradient operators asδx(·) =(·)i+1,k− (·)i,k∆x, and (4.12)δζ (·) =(·)i,k+1− (·)i,k∆ζ .(4.13)The prognostic variable for vertical motion perpendicular to the terrain-following vertical coordinateζ isωn+1 =ζxun+1 +ζzwn+1. (4.14)We use the following notations to combine the known rhs terms in the momentum equations:RU ≡ unD +∆t(Fu)nD−∆t2{γRdpiρ˜mx[(δxΘ′m)ζ +δζ (ζxΘ′mxζ)]}nD,(4.15)RW ≡ wnD +∆t(Fw)nD−∆t2{γRdpiρ˜mζ[δζ (ζzΘ′m)]−1ρ˜m[gρ˜dpi ′pi −gρ˜′m]ζ}nD,(4.16)the implicit pressure-gradient terms:PGU =−∆t2{γRdpinρ˜nmx[(δxΘ′m)ζ +δζ (ζxΘ′mxζ)]n+1A}, (4.17)PGW =−∆t2{γRdpinρ˜nmζ[δζ (ζzΘ′m)]n+1A}, (4.18)and the implicit half of the buoyancy term:BW =∆t21ρ˜nm[gρ˜dpi ′pi −gρ˜′m]n+1Aζ. (4.19)75Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesThe subscripts in the notations denote the momentum equations to which the terms belong. Using(4.8), (4.9), and the notations above, the vertical momentum equation can be rewritten as:ωn+1 = ζxRU +ζz(1+µ∆t)−1(RW +BW)︸ ︷︷ ︸RΩ+ζxPGU +ζz(1+µ∆t)−1PGW.(4.20)The notation RΩ is used here to represent all the known terms plus the implicit buoyancy term (4.19).Often, off-centering of the time-averaged terms is needed in semi-Lagrangian semi-implicittime-stepping schemes to help eliminate computational noise, especially when orographic forcingis present and at large Courant numbers (e.g. Rivest et al., 1994). In CSLAM-NH, no off-centeringwas needed to attain the numerical stability in the solver for the test cases presented here.4.2.4 Conservative and consistent flux-form equationsAs noted by Lauritzen et al. (2006) and demonstrated in Wong et al. (2013) and Wong et al. (2014),when a numerical scheme different from the one used to evaluate the continuity equation is usedto transport scalar variables, conservation of scalar mass will no longer be guaranteed, despite theuse of a conservative transport scheme. The problem of numerical consistency in cell-integratedsemi-Lagrangian schemes is resolved through the use of a new flux-form CISL continuity equationintroduced in Wong et al. (2013) for the shallow-water equations and tested for a 2D nonhydrostaticatmosphere without topography (Wong et al., 2014). The new flux-form CISL continuity equationallows for a straightforward implementation of a CISL scalar transport scheme that ensures numer-ical consistency. Here, we further test the proposed formulation based on the CSLAM transportscheme for 2D idealized cases over mountains.The potential temperature, continuity, and scalar-mass conservation equations are all solvedconsistently using the same numerical scheme presented in Wong et al. (2014):Θˆn+1m =Θn+1m,exp +∆t2[∇eul · (v′nΘˆnm]δA∗∆A+∆tFnΘmδA∗∆A(4.21)andΘ˜n+1m = Θˆn+1m −∆t2[∇eul · (v′n+1Θˆn+1m )]. (4.22)The flux divergence in terms of a corrective velocity v′ in the semi-implicit correction term is76Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesdefined as:∇eul · (Θmv′) =1∆x[Θxm(ur−Fr/∆ζ )−Θxm(ul−Fl/∆ζ )]+1∆ζ [Θζm(ωt −Fl/∆x)−Θzm(ωb−F/∆x)],(4.23)where F =F (u,ω) are Lagrangian flux areas, computed as in Wong et al. (2014). The velocitiesur, ul , ωt , and ωb are staggered velocities at the cell faces (Arakawa C grid).In the semi-implicit flux-form equation, instead of linearizing around a mean reference state,we utilize Θˆn+1m using the CSLAM transport scheme to ensure consistency of the semi-implicitcorrection term among all the scalar flux-form equations. Included in this CSLAM computationare all the terms to be integrated over the departure cell: the explicit conservative CSLAM solution(Θn+1m,exp), a predictor-corrector term (the flux divergence term at time level n), explicit diffusion, anddiabatic tendency (the latter two are combined in FΘm). The diabatic tendencies are approximatedusing values at the previous time level. The resulting approximation Θˆn+1m (4.21) is then used in(4.22). The solution from (4.22) is the solution from the dynamics, and prior to any adjustment tosaturation by a moist microphysics scheme.Consistent formulations of the continuity equation and scalar mass conservation equations arestraightforwardly discretized as:φˆ n+1 = φ n+1exp +∆t2[∇eul · (v′nφˆ n]δA∗∆A+∆tFnφδA∗∆A,(4.24)andφ˜ n+1 = φˆ n+1− ∆t2[∇eul · (v′n+1φˆ n+1)], (4.25)where φ = ρ˜d or Q j. Similar to (4.21), (4.24) combines the advected quantities in the explicitsolution using the CSLAM scheme, the predictor-corrector flux divergence term, diffusion, anddiabatic tendencies from the previous time level. The solution of the velocities at the new time levelare used to compute the flux divergence in (4.25).77Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates4.2.5 Helmholtz equationBy eliminating v′n+1 in the potential temperature equation (4.22) using the horizontal and verticalmomentum equations, we form the Helmholtz equation:Θ˜′n+1m +∆t2{δx(PGU Θˆn+1mx)+δζ[(ζxPGUxζ+ζz(1+µ∆t)−1PGW)Θˆn+1mζ]}=Θˆ′n+1m −∆t2[δx(RU Θˆn+1mx)+δζ (RΩΘˆn+1mζ)], (4.26)where PGU and PGW are, as defined earlier, the implicit pressure-gradient terms expressed as func-tions of Θ˜n+1m . All the terms on the rhs of (4.26) are pre-computed at the beginning of each timestep, and the implicit buoyancy term in RΩ is updated at each iteration of the Helmholtz equationsolver (described next).4.2.6 Iterative centered-implicit time-stepping schemeThe compressible Euler equations permit fast horizontally- and vertically-propagating acoustic andgravity waves. To alleviate the time-step limit due to acoustic waves, in the previous version ofCSLAM-NH (Wong et al., 2014), an implicit time-stepping scheme was used to solve the pressure-gradient and mass-divergence terms. The remaining buoyancy terms were evaluated explicitly usinga two time-level extrapolation scheme. The semi-implicit time integration scheme allowed the useof time steps much larger than those allowed in a classical explicit scheme, which would otherwisehave been restricted by the speed of sound. The buoyancy terms responsible for gravity waves,however, imposed a restriction to the maximum stable time step.Instead of evaluating the gravity wave terms explicitly using time extrapolation, we use an itera-tive approach for a more accurate and implicit treatment of these terms. The solution procedure canbe summarized in two main components as follows. First, the departure cell areas are approximatedusing backward trajectories from the arrival grid cell vertices. The forcing terms (RU , RΩ) and theexplicit departure cell-averaged potential temperature, Θˆn+1m [using (4.21)] are evaluated and formthe RHS of the Helmholtz equation (4.26). The implicit buoyancy term BW in RΩ is evaluated at timelevel n as an initial estimate. The explicit departure cell-averaged density ρˆn+1d [using (4.24)] is alsopre-computed. The second component involves solving the linear Helmholtz equation for Θ˜′n+1m ;here we use a conjugate gradient residual solver (Skamarock et al., 1997). The solution Θ˜′n+1m isthen back-substituted into the momentum equations (4.8) and (4.20) to get un+1 and ωn+1, respec-tively. Finally, the implicit buoyancy term BW (4.19) in RΩ is updated using (i) ρ˜′n+1m by evaluating(4.25) and (ii) pi ′n+1 directly from Θ˜′n+1m using the equation of state pi ≡ (RdΘm/p0)R/cv . At the endof the second component, the trajectories and forcing terms (first component of the procedure) are78Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesrecomputed using the latest solution of un+1 and ωn+1.Depending on the test case, two to four iterations of each component are performed. For thenonlinear flow tests, iterating more than twice did not further improve the maximum stable timestep size. For the linear cases, the maximum time step can be further increased by performing moreiterations (iterating more than four times does not further improve stability). At each iteration, theHelmholtz solver converges progressively faster (since the latest estimate of Θ˜′n+1m is used as thestarting point). The iterative scheme is used for advancing the dry dynamics; after which, tracersare advected using (4.25) and the moist physics are called (once at each time step).The use of an iterative centered-implicit scheme is found to substantially increase the stabletime step size in CSLAM-NH at the expense of solving the Helmholtz equation more than onceper time step. To demonstrate this behaviour, we conduct the gravity-wave test originally proposedin Skamarock and Klemp (1994), using CSLAM-NH as was done in Wong et al. (2014) with agrid spacing of ∆x = ∆z = 1 km and an imposed mean wind U = 20 m s−1. Wong et al. (2014)used an explicit treatment of the buoyancy terms and found that the maximum stable time step wasrestricted to ∆t = 38 s [at a nominal advective Courant number (Cr) of 0.76]. In the current versionof the iterative centered-implicit CSLAM-NH, we have found that for the same simulation, themaximum stable time step increased to 100 s, roughly by a factor of 2.6 (Cr = 2). For comparison,the maximum stable time step for an Eulerian split-explicit third-order Runge-Kutta time steppingscheme was 60 s (Cr = 1.2) (Wong et al., 2014).Similar iterative approaches were found to improve numerical stability in other semi-Lagrangiansolvers. In the Canadian Global Environmental Multiscale (GEM) model, Coˆte´ et al. (1998) dis-cretized the governing equations in a fully implicit manner and used an iterative procedure to avoidsolving a nonlinear Helmholtz equation. This procedure is also implemented in Melvin et al. (2010)for the vertical-slice nonhydrostatic solver using the ‘Semi-Lagrangian Inherently Conserving andEfficient’ (SLICE) transport scheme. An alternative predictor-corrector (thus, also iterative) ap-proach was tested in the European Centre for Medium-Range Weather Forecasts (ECMWF) Inte-grated Forecast System (IFS) model by Cullen (2001). In that study, a positive improvement in ac-curacy was noticeable only when the advective velocities, in addition to the buoyancy terms, wereiterated. Using an idealized analysis of acoustic modes in a 1D nonhydrostatic vertical column,Cordero et al. (2005) demonstrated the impact of using time-extrapolated and -interpolated trajec-tory computations on the numerical stability of a semi-Lagrangian centered-semi-implicit scheme.When extrapolation and large time steps were used for the trajectories, the vertical structure of theacoustic modes were found to be distorted (with spurious zeros forming with time). The time in-terpolation scheme on the other hand was found to be stable in all cases. The idealized analysisby Cordero et al. (2005) supports the findings in Cullen (2001) and the method used in Coˆte´ et al.79Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates(1998), with a recommendation for time-interpolated trajectory computations (e.g. by repeating thefirst component of the CSLAM-NH solution procedure).The disadvantage of this approach is that the linear Helmholtz equation for potential temperatureΘm is solved a number of times with the buoyancy terms updated at the end of each iteration.However, the increased stability will allow a larger time step to be used and can help offset theadded computational expense of solving the dry dynamics (calculated once at each time step). Afterthe dry dynamics, the solver then advects passive tracers (once at each time step). With a largertime-step size, the total number of times tracers are advected during the entire simulation is reduced.Therefore, the overall execution time spent on scalar transport is also reduced. This reduction mayhave a significant impact on computational time, especially when the number of tracers used inchemistry applications is large.4.2.7 Boundary conditionsPeriodic-in-x and free-slip top and bottom boundary conditions are applied in all our tests. Thevertical velocities at the top and bottom boundaries are set to ω = 0, and ensures no normal fluxthrough them. The boundary conditions are implemented by extrapolating Θm, ρd , and u into theboundary.4.2.8 Implicit Rayleigh dampingTo prevent the reflection of vertically-propagating gravity waves along the rigid model top, a damp-ing term, −µw, is added in the vertical momentum equation based on the scheme proposed inKlemp et al. (2008) and implemented in Melvin et al. (2010). The damping profile µ(z) proposedby Klemp and Lilly (1978) is used:µ(z) =µmax sin2(pi2z−zdzt−zd)if z > zd ,0 if z≤ zd .(4.27)The profile is characterized by a gradual increase of viscosity with height, which is desirable toprevent any reflections that would otherwise occur from a sharp increase in viscosity. The dampinglayer starts from a user-specified height zd and extends to the top of the domain zt . The valuesof µmax are specified for each test case. Klemp et al. (2008) analyzed the reflection properties ofthis implicit Rayleigh damping layer. The scheme they proposed is slightly different from the oneapplied here; in particular, they proposed implementing the damping term as an adjustment step.The adjustment step approach includes an extra damping term that resembles vertical diffusion, inaddition to the effect of a damping term −µw added directly in the vertical momentum equation.80Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesFor smaller (nonhydrostatic) horizontal scales, however, the effect of the damping term dominatesand there is little difference between the two approaches. Klemp et al. (2008) also provided someguidance in selecting a suitable µmax value based on the reflection properties of the scheme. Basedon their experimentation, smaller damping coefficients µmax should be used for smaller dominanthorizontal wavelengths.4.3 Idealized test cases: results4.3.1 Linear mountain waves over bell-shaped mountainTo test the response of the nonhydrostatic solver to orographic forcing, two adiabatic linear mountain-wave simulations are conducted first. Both cases assume a simple hill profile h(x) of a witch-of-Agnesi curve, defined ash(x) =hma2x2 +a2, (4.28)with a small amplitude hm = 1 m but different half-widths, a. Gravity waves generated by flowmoving over a wide hill under conditions where U/Na 1 are approximately hydrostatic and arevertically propagating (Smith, 1979). We simulate flow with a constant upstream wind speed U = 20m s−1, in an atmosphere initially in hydrostatic balance and isothermal at T = 250 K (equivalentto a stratification of N2 ≈ 3.83× 10−4 s−2, where N2 = g2/cpT for an isothermal atmosphere).The mountain half-width is set at a = 10 km to give U/Na = 0.1, such that nonhydrostatic effectsare small for this broad low hill. The analytical steady-state solution of vertical velocity usingthe linear theory and fast-Fourier transform algorithm described in Smith (1980) is plotted for thisflow in Fig. 4.2. The physical domain is 120 km wide and 25 km deep. The simulation is run for anondimensional time Ut/a= 30 (equivalent to t = 15000 s) to ensure the solution has reached steadystate. The numerical domain has dimensions 120 × 100 (∆x = 1 km and ∆z = 250 m). A Rayleighdamping layer (µmax = 0.1 s−1) is implemented in the top 10 km of the domain (approximately 1.5times the vertical wavelength, λz = 2piU/N = 6.4 km).Results from simulations using a small Courant number Cr= 0.2 and large Cr= 1.5 are shown inFig. 4.3. An upstream tilt of the phase lines is observed, corresponding to energy originating fromthe ground (the mountain) and propagating upwards. As expected, the amplitude of the verticalvelocity also increases with height (∝ ρ−1/2), corresponding to the effect of wave amplification dueto decreasing density at higher altitudes. The slight downstream tilt of the wave pattern with heightis due to weak nonhydrostatic influences and is also observed in other nonhydrostatic models forthe same test case (e.g. Melvin et al., 2010). The modelled vertical velocities of Fig. 4.3 comparevery well with the exact linear solution shown in Fig. 4.2.81Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinateshorizontal distance (km)height (km) −40 −20 0 20 400510−4−2024x 10−3w (m s-1)Figure 4.2: Analytical steady-state solution for the linear hydrostatic mountain waves obtainedfollowing Smith (1980).height (km) 0510−4−2024x 10horizontal distance (km)height (km) −40 −20 0 20 400510-3w (m s-1)w (m s-1)Cr = 1.5Cr = 0.2Figure 4.3: CSLAM-NH simulation of a linear hydrostatic wave (U/Na≈ 0.1) for a low widemountain showing vertical velocity w (in m s−1) after T = 15000 s using (top) ∆t = 10 s(Cr = 0.2) and (bottom) ∆t = 75 s (Cr = 1.5). Contour interval is 5×10−4 m s−1. Meanwind is from left to right.For a narrow mountain (U/Na ≈ 1), the mountain waves are strongly nonhydrostatic. Thesewaves are highly dispersive, with shorter horizontal scales propagating farther downstream withheight, and scales less than 2piU/N becoming evanescent. To simulate such a flow, the half-width aof the mountain is reduced to 1 km, and the initial background state has a constant stratification ofN2 = 1×10−4 s−2. The impinging flow has an upstream wind speed of U = 10 m s−1 (U/Na = 1).The analytical steady-state solution of vertical velocity obtained following Smith (1980) is plottedfor this flow in Fig. 4.4. The domain is 144 km wide and 25 km deep. The numerical domain hasdimensions 360 × 100 grid cells (∆x = 400 m and ∆z = 250 m). Even though the mountain is only82Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinateshorizontal distance (km)height (km) −10 0 10 20 300510−505x 10-3w (m s-1)Figure 4.4: Analytical steady-state solution for the linear nonhydrostatic mountain waves ob-tained following Smith (1980).height (km) 0510−5−2.502.55 x10−3horizontal distance (km)height (km) −10 0 10 20 300510w (m s-1) w (m s-1)Cr = 0.125Cr = 1.5Figure 4.5: CSLAM-NH simulation of a linear nonhydrostatic wave (U/Na = 1) for a narrowmountain showing vertical velocity w (in m s−1) after T = 9000 s using (top) ∆t = 5 s(Cr = 0.125) and (bottom) ∆t = 60 s (Cr=1.5). Contour interval is 6×10−4 m s−1. Meanwind is from left to right.marginally resolved at this model resolution, the configuration is kept for comparison purposes withMelvin et al. (2010). The Rayleigh damping layer is applied to the top 13 km of the domain (twicethe length of λz) with µmax = 0.05 s−1.Results from CSLAM-NH are shown in Fig. 4.5 for two different time step sizes (Cr = 0.125 andCr = 1.5). As expected, far downstream of the mountain, the solutions exhibit a more pronounced(relative to the linear hydrostatic case of Figs. 4.3 and 4.2) downstream tilt of the phase due to thenonhydrostatic component of the waves, and the waves also decay with height (note: only part ofthe domain is shown in the figure). This solution also compares well with the analytical solution83Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesshown in Fig. 4.4.4.3.2 Downslope windstormTo test the nonhydrostatic solver in a highly nonlinear flow, a simulation of the famous downs-lope windstorm that occurred on 11 January 1972 in Boulder, Colorado (Lilly, 1978) is conducted.Strong surface winds, gusting to 55 m s−1, were observed in Boulder, Colorado on that day. Thewindstorm has been a long-standing case for theory development and numerical model verification,e.g., Klemp and Lilly (1978); Peltier and Clark (1979); Durran (1986). More recently, Doyle et al.(2000) carried out a model intercomparison study of 11 different high-resolution models to assesstheir ability in numerically simulating the wave breaking process of this windstorm. Prior to Doyleet al. (2000), smoothed soundings were used to initialize the models; here, we use the same 12 UTC11 January 1972 Grand Junction, Colorado sounding as in Doyle et al. (2000), where they showedthat a more realistic simulation of the windstorm was generated.The numerical set-up is based on Doyle et al. (2000). The mountain half-width is 10 km witha height of 2 km. The domain is 240 km wide and 25 km deep. The numerical domain dimensionsare 240 × 125 grid cells (∆x = 1 km and ∆z = 200 m). The time step sizes used in Doyle et al.(2000) ranged from 1 s to 12 s [the largest time step size of 12 s was run using the Durran andKlemp (1983) model that uses a time-splitting scheme with a small time step of 3 s]. To comparethe results with these models, a similar time step size of 10 s is used in CSLAM-NH. A Rayleighdamping layer (µmax = 0.05 s−1) is applied only in the top 7 km of the domain (18≤ z≤ 25 km) toprevent the damping of the physically significant wave breaking in the lower stratosphere. A fourth-order horizontal smoothing filter is applied with a coefficient KD = 1×109 m4 s−1, which smoothsout any small-scale variations and helps maintain numerical stability in the model. Unlike most ofthe models in that study, no turbulence parameterization or any other explicit diffusion was used;turbulent dissipation is solely dependent on the hyperviscosity applied and any inherent numericaldissipation associated with the model discretization.In the results presented within Doyle et al. (2000), all models produced significant strengtheningof the winds on the lee of the mountain and wave breaking in the upper troposphere and stratosphereat time 3 h. Despite using identical initial conditions, however, significant differences were foundamong the model results due to differences in the model formulations (e.g. spatial and temporaldiscretizations, type of explicit diffusion used), as well as the nonlinearity of the flow.The CSLAM-NH results at 3 h are presented in Fig. 4.6. The wave breaking regions in CSLAM-NH can be identified as the adiabatic (well-mixed) regions (Fig. 4.6a) and highly turbulent areas(Richardson number, Ri < 0.25) (Fig. 4.6c, where to be consistent with Doyle et al. (2000), theRichardson number sgn(Ri)|Ri|0.5 is plotted). The Richardson number used in Fig. 4.6c is the bulk84Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates0 0 0000000000000000000000000000088888888888 8888888888888888161616161616161616166161616 1616161616161624242424242424242424242424323232 32323232323232324040 404040484848565664−16−8−8−8−8−8−8−8288 288296296296296304304304312312320320320 320328 328336336344344344352352352352360360360360368 368368376 376376384 384384384392 392392392400 400408408408416 416416424 424 424432 432 432440440440448 448448456 456456464 464464472472472480480480480488488488488496 496 496504504 504512512512520 520520528 528 528536 536 536544 544 544552552 552560 560 560568 568 568576 576 576584 584 584592 592600 600608 608616 616624 624632 632640 640648 648656 656664 664672 672Horizontal distance (km)0510152025Height (km)Height (km)Height (km)−100 −50 0 50 100(c) sgn(Ri) |Ri|0.505101520250510152025(a) θ(b) uFigure 4.6: CSLAM-NH simulation for the 1972 Boulder windstorm case (a) potential tem-perature θ (K) (with contour interval of 8 K), (b) horizontal velocity U (m s−1) (withcontour interval of 8 m s−1), and (c) Richardson number sgn(Ri)|Ri|0.5 (−5≤ Ri ≤ 1 areplotted with contour interval of 0.5, following Doyle et al. (2000); grey contour showsnegative values) at T= 3 h using a time step ∆t = 10 s.85Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesRichardson number (dry),Ri =g/θ(∆θ/∆z)(∆u/∆z)2. (4.29)For locally statically stable air (Ri > 0), the critical Richardson number at which wind shear is strongenough to sustain turbulence and overcome the damping by negative buoyancy is 0.25. The wavebreaking regions appear to be in the vicinity of 12 ≤ z ≤ 16 km and 17 ≤ z ≤ 20 km, comparableto the results in the intercomparison study. An initial critical level at z = 21 km (where U = 0) isalso found to be damping in the CSLAM-NH model simulation, and traps the vertically-propagatinggravity waves. The damping effect is evident in the smooth isentropes and lack of turbulence (largeRi, not contoured) above that height.The lateral position of the hydraulic jumps at time 3 h varied among the models given in Doyleet al. (2000), with several occurring over the leeslope and others farther downstream. The associatedmaximum leeslope winds from the 11 models were found to range from 43 m s−1 to 86 m s−1. InCSLAM-NH, the hydraulic jump feature is found on the leeslope, and the simulated maximumdownslope wind speed at the surface (lowest model level) is located at 10.5 km downstream fromthe mountain crest at 56.6 m s−1 (Fig. 4.6b). Flow features aloft such as the flow reversal at5≤ z≤ 10 km that were present in many of the models in Doyle et al. (2000), is also present in theCSLAM-NH results. This weakening of the winds above the hydraulic jump was also observed inthe aircraft flight data analysis (see, e.g., Fig. 2b in Doyle et al., 2000).The hyperviscosity coefficients used by the models in the model intercomparison study rangedfrom 1.1× 108 to 5.0× 109 m4 s−1. Time series of simulated maximum downslope wind speedsusing different diffusion coefficients in CSLAM-NH are given in Fig. 4.7. Results from varyingthe horizontal smoothing coefficient from 1.5× 108 to 1× 109 m4 s−1 show slight variation in thesimulated maximum downslope wind speed, with values at time 3 h ranging from 53.1 m s−1 to 56.6m s−1. The impact of using different magnitudes of horizontal smoothing is apparent once the wavesbegin to break, giving a maximum range of predicted downslope wind speeds of approximately12 m s−1. The general trend of the downslope windstorm development, however, is similar withmaximum surface winds of the simulation occurring at around 3 h, with weakening thereafter dueto the limited horizontal extent of the domain and periodic lateral boundaries.The maximum stable time step in CSLAM-NH for this wave breaking case is 20 s (when KD =5.0× 109 m4 s−1 is applied). With a time step larger than 20 s, the errors of the linear trajectoryapproximations become large enough that the departure cells self-intersect as illustrated in Fig. 4.8.In this case, the flow is characterized by a strong horizontal shear of the vertical wind speeds as wellas large vertical Courant number, Crz ≈ 3 (not shown). Higher-order cell-edge approximations havebeen explored by Ullrich et al. (2012), and may help alleviate the time step limit by increasing the86Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates0 0.5 1 1.5 2 2.5 3 3.5 410203040506070Time (mins)Maximum downslope wind speed(m s-1) KD = 1.5 × 108 = 2.5 × 108 = 5.0 × 108 = 10 × 108KDKDKDFigure 4.7: Time series of the simulated maximum CSLAM-NH downslope wind speeds (ms−1) for the 1972 Boulder windstorm case using different horizontal smoothing coeffi-cients, KD (m4 s−1).accuracy of the area integration. Overall, CSLAM-NH is able to generate comparable results to theother models, and at a maximum time step size that is roughly double of those used in Doyle et al.(2000).4.3.3 Moist flow over a mountain in a nearly neutral environmentThe nonhydrostatic solver is tested for another nonlinear flow, but in this case, we also include theeffects from moist processes. A simulation of saturated flow over a mountain in an initially nearlyneutral environment is conducted. This test case also demonstrates the ability of the solver in pro-ducing realistic orographic precipitation. The simulation is based on the test cases presented inMiglietta and Rotunno (2005). Moisture in the atmosphere is an important factor in modifying flowover topography. Durran and Klemp (1983) studied the influence of moisture on mountain wavesusing numerical simulations. In both a linear mountain-wave test and a downslope-windstorm test,they found that the inclusion of upstream moisture can greatly reduce the amplitude of these wavesrelative to their dry analogs. As the mountain enhances lifting of the moist flow over the windwardside, condensation commonly occurs, leading to clouds and precipitation. The downstream evap-oration of these clouds and precipitation can reduce the static stability at these altitudes, and theair can become desaturated on the lee side of the mountain due to rainout processes and diabaticwarming in the descent.For a nearly neutral flow, Miglietta and Rotunno (2005) simulated the transition of saturated airupstream to unsaturated air downstream due to diabatic warming in the downward motion on thelee. The inverse Froude number Nmhm/U is near zero, indicating that the resistance due to gravity87Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinateszxFigure 4.8: A self-intersecting departure cell (highlighted in red with vertices marked by blackcircles) in CSLAM-NH when a large time step size of 25 s is used for the strongly shearedflow in the 1972 Boulder downslope windstorm case. Black circles indicate departuregrid cell vertices, and white circles the Eulerian arrival grid cell vertices. Arrows sym-bolize the computed backward-in-time trajectories. Trajectories and the arrival grid cellassociated with the self-intersecting departure cell are highlighted in red.is minimal and the flow can freely translate over the mountain.To include moisture effects when determining local static stability, Lalas and Einaudi (1974),and later verified by Durran and Klemp (1982), derived an expression for the moist Brunt-Va¨isa¨la¨frequency:N2m = gΓ(d lnθdz+LvcpTdqsdz)−g1+qwdqwdz, (4.30)where T is the absolute temperature, qs is the saturated water vapour mixing ratio, qw is the totalwater mixing ratio, andΓ=1+ 1qs+ε∂qs∂ lnθ∣∣pi1+ LvcpT∂qs∂ lnθ∣∣pi, (4.31)is the ratio of the moist to dry adiabatic lapse rates. (All other variables are as defined previously.)88Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesA few more specifics regarding the generation of the moist neutral sounding that supplementsthe derivation presented in Miglietta and Rotunno (2005) are given. Following the procedure inMiglietta and Rotunno (2005), to generate the initial sounding of a specific Nm, a first-order or-dinary differential equation is solved. To create the initial nearly neutral sounding, the first-orderordinary differential equation for potential temperature is solved iteratively based on a specifiedsurface temperature (15 ◦C) and reference pressure (p0 = 100 kPa). To be consistent with Migliettaand Rotunno (2005), the Wexler’s formula for saturated vapour pressure (in mb) is used:es(T ) = 6.11exp(17.67T −273.15KT −29.65K). (4.32)The definition for qs = εes/(p−es), where ε = Rd/Rv is used to derive ∂qs/∂ lnθ |pi in (4.31). First,differentiating qs with respect to lnθ at constant pi gives∂qs∂ lnθ∣∣∣∣pi=qsespp− es∂es∂ lnθ∣∣∣∣pi, (4.33)and differentiating (4.32) (using T = piθ ) with respect to lnθ at constant pi gives the expression∂es∂ lnθ∣∣∣∣pi= esT17.67(243.15)(T −29.65)2. (4.34)To find θ(z) for a specific Nm, (4.30) must be iterated to convergence (10−12) at each pressure level(or height) since qs = qs(pi,θ) and Γ are also functions of the unknown. To get the pressure at eachheight, the hydrostatic equation is used:pi j+1−pi j∆z=−gcp1+qwzθ zm. (4.35)For each model level j, the discrete form of the ODE solving for θ at height z islnθ j+1 +LvcpTz(qs, j+1−qs, j)− (qw, j+1−qw, j)(1(1+qw)Γ)z= lnθ j +N2mgΓz(z j+1− z j). (4.36)[Note: Miglietta and Rotunno (2005) expresses this equation in terms of (T, p)).] Other aspectsof the Kessler microphysics scheme also require modification. Following Miglietta and Rotunno(2005), no autoconversion from cloud water to rain is permitted in the first five hours to allow forinitial adjustment of the flow to the impulsive introduction of terrain. For a consistent definition of89Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesqs throughout the model, the production of cloud water due to saturation is also modified:∣∣∣∣dqsdt∣∣∣∣=qv−qs1+ Lvcp( ∂qs∂T)|p, (4.37)where, based on the Wexler’s equation for es:∂qs∂T∣∣∣∣p=qs pp− es17.67(243.5)(T −29.65)2. (4.38)Miglietta and Rotunno (2005) used a small N2m = 3× 10−6 s−2 to represent a nearly neutraltroposphere due to the limitations of the single machine precision accuracy of their model. Theyfound that using any smaller Nm led to solutions that were apparently convectively unstable. TheCSLAM-NH solver has machine double precision accuracy, so for this moist neutral flow case, aninitial Nm = 0 in the troposphere is applied. Simulations using a range of ‘small’ N2m ∼ 10−11 in theinitialization step show solutions similar to applying Nm = 0 and resemble those in Miglietta andRotunno (2005) more than using their N2m = 3×10−6 s−2. An initial N2m = 4.84×10−4 s−2 is usedfor the isothermal stratosphere.Two mountain cases with different heights, hm = 700 m and 2 km, were chosen from Migliettaand Rotunno (2005) for their distinct differences in orographic distribution of moisture. Both testcases are run using the Witch-of-Agnesi curve with a half-width of 10 km. The same numericaldomain that is 800 km wide and 20 km deep is used, and the grid dimensions are 400 × 80 gridcells (∆x = 2 km and ∆z = 250 m). In both cases, a mean wind U = 10 m s−1 is applied. Theatmosphere is initially saturated (qv ≡ qs) with constant cloud-water mixing ratio qc = 0.05 g kg−1set everywhere in the domain to prevent the atmosphere from becoming subsaturated due to theimpulsive introduction of the mountain at initial time. The Rayleigh damping layer (µmax = 0.1s−1) is applied in the top 5 km of the domain. Second-order filters in the horizontal and verticaldirections are applied with coefficients 3000 and 3 m2 s−1, respectively. The Prandtl number is 3.This configuration is the same as that in Miglietta and Rotunno (2005).Both cases suggest a desaturation of the air downstream of the mountain with time. Migliettaand Rotunno (2005) noticed in their simulations that for intermediate mountain heights (500 m≤ hm ≤ 1500 m) the unsaturated region downstream of the hill unexpectedly extends upstream aswell. A later study by Keller et al. (2012) showed that this upstream extent of the subsaturated airis due to local adiabatic descent and warming caused by a transient upstream-propagating gravitywave, a fundamental feature of a two-layer two-dimensional atmosphere with topography intro-duced impulsively. The purpose of performing the test case as prescribed in Miglietta and Rotunno90Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates−100 −50 0 50 100051015Horizontal distance (km)Height (km) 00.10.5qc (g kg-1)Figure 4.9: CSLAM-NH cloud-water mixing ratio (g kg−1) at time 5 h from an initially sat-urated nearly neutral flow (with an initial qc = 0.05 g kg−1) over a 700 m hill. Whiteregion above ground indicates subsaturated air (qc ≡ 0).(2005) with a mountain 700 m high is to ensure that CSLAM-NH can generate comparable resultsto models used in the literature, such as that in Miglietta and Rotunno (2005) who used the Weatherand Research Forecasting (WRF) model v1.3.Fig. 4.9 shows the solution from CSLAM-NH [c.f. Fig. 5d of Miglietta and Rotunno (2005)]using a time step size of 20 s. The white region indicates subsaturated air, as described previously.Although the upstream region of the subsaturated air in Miglietta and Rotunno (2005) extends far-ther upstream (x = −100 km) than that found using CSLAM-NH, the solution from CSLAM-NHcompares very well with that obtained in Miglietta and Rotunno (2005). The maximum stableCSLAM-NH time step size is 50 s with two iterations of both components in the iterative centered-implicit scheme.Fig. 4.10a shows the CSLAM-NH cloud-water mixing ratio at time 5 h 10 mins (10 minsafter autoconversion of rain is permitted) of a simulation using ∆t = 20 s for the large-amplitudemountain (hm = 2 km) case [c.f. black contours in Fig. 8a in Miglietta and Rotunno (2005)]. Similarto the results presented in Miglietta and Rotunno (2005), no upstream region of the subsaturated airis found. In addition, the formation of convective cells due to the reduction of local static stabilitydownstream of the mountain is also detected in the CSLAM-NH simulation. The instability is foundto be primarily associated with a hydraulic jump feature downwind. Fig. 4.10b shows the rainwater91Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesmixing ratio for the same simulation time as in Fig. 4.10a.Compared to the results in Miglietta and Rotunno (2005), CSLAM-NH indicates more rainspillover to the lee of the mountain (c.f. grey contours in Fig. 8a in Miglietta and Rotunno (2005)).Simulation of our case using an Eulerian split-explicit model similar to the one used in Miglietta andRotunno (2005) shows virtually the same distributions of cloud- and rainwater as in the CSLAM-NH simulation (Fig. 4.10). The similarity of the CSLAM-NH solution to that of the second Eulerianmodel seems to suggest that the discrepancy is not specific to CSLAM-NH and may be related tocertain aspects of the initialization procedure. Miglietta and Rotunno (2005) suggested that theirsimulations were sensitive to small changes to Nm ≈ 0. If there is a slight discrepancy between thevalue of Nm specified in the initialization than that used in Miglietta and Rotunno (2005), the flowdynamics may be altered. In this situation, the greater amount of waterloading to the lee of themountain could imply a lower effective terrain height (Nmhm/U), i.e. lower stability and/or strongerwinds, such that the advection of the hydrometeors happens at a faster time scale than the fallout ofprecipitation.−40 −20 0 20 40 60 800123456789Horizontal distance (km)Height (km) 00.20.61−40 −20 0 20 40 60 800123456789Horizontal distance (km)Height (km) 00.20.40.60.81(a) qc (g kg-1) (b) qr (g kg-1)Figure 4.10: (a) CSLAM-NH cloud-water mixing ratio qc (g kg−1) at time 5 h 10 mins froman initially saturated nearly neutral flow (with an initial qc = 0.05 g kg−1) over a 2 kmmountain. White region above ground indicates subsaturated air (qc ≡ 0 g kg−1) (b)CSLAM-NH rain-water mixing ratio qr (g kg−1) at the same simulation time.92Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates4.3.4 Computational performance (This section is not in the manuscript.)The computational expense associated with the iterative centered-implicit scheme in CSLAM-NHis examined. The iterative scheme allows a larger CSLAM-NH time-step size, at the expense ofsolving the Helmholtz equation more than once. Taking a larger time step increases computationalefficiency by requiring fewer time steps; however, it may also slow down the convergence rate inthe Helmholtz solver. The moist neutral flow described in section 4.3.3 has the greatest complexityof topographical and latent heating effects and will be used to examine some of the computationalaspects in CSLAM-NH.To examine whether the larger time steps allowable by the iterative scheme actually save com-putational time, simulations using CSLAM-NH are made for a range of time-step sizes allowableby the iterative solution procedure. Fig. 4.11 shows the computational times required for each time-step size as compared to using ∆t = 20 s (which is the time-step size used in section 4.3.3). Withoutusing the iterative centered-implicit scheme, the maximum stable time-step size is 5 s. The smallertime step required a total computational time roughly 23% more (dashed line in Fig. 4.11) than thatin section 4.3.3. The solid line shows the computational times using two iterations for time-stepsizes between 10 and 50 s.Using the iterative scheme to double the time-step size to ∆t = 10 s increases the computationaltime by 9% as compared to not using the iterative scheme. The added computational expense stemsfrom the iterations and the slower convergence in the Helmholtz solver. Comparing to ∆t = 20 s,the time step is reduced by half. However, a reduction of time-step size by half does not double thecomputational time, which indicates that some components of the solver are running more efficientlyat smaller time steps, such as the elliptic solver attaining faster convergence.When the time-step size is increased to 30 s, an improvement in the solver efficiency is ob-served. The larger time-step size is able to reduce the computational time by 13%. When thetime-step size is increased further to ∆t = 40 s, the computational time is reduced by another 4%,and by another 3% when increased to ∆t = 50 s. The negative impact of the slower Helmholtz solverconvergence rate on the model performance is evident at the larger time-step sizes. Although thereis still improvement in the net efficiency when using large time steps, the benefit is greatly hinderedby the elliptic solver. In comparison to the non-iterative stable time-step size, the iterative schemeimproves the overall efficiency by about 43%. In the iterative scheme, the need to recompute thetrajectories (briefly discussed in section 4.2.6) was found to help stabilize CSLAM-NH in the nu-merical experiments. The stability properties of the trajectory recomputations in CSLAM-NH maybe worth further investigation.The distribution of the computational expense of the current CSLAM-NH solver can be sep-93Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates10 20 30 40 500.80.911.11.21.31.4Time step size (s)Normalized computational timeFigure 4.11: Computational times for the moist neutral flow over a 2 km mountain using dif-ferent CSLAM-NH time-step sizes (∆t). Times are normalized to that using ∆t = 20 sfor comparison purposes. The solid line indicates computational times using the itera-tive centered-implicit scheme. Dashed line indicates using the non-iterative CSLAM-NH with ∆t = 5 s.arated into two main components: the dry dynamics and the moist physics. The computationalexpense is heavily skewed towards the elliptic solver. The total computational time for the moistneutral flow case using a time-step size of 50 s (372 time steps) is distributed as follows. (Note: thefollowing discussion accounts for 99.7% of the total computational time, while the neglected 0.3%is from the initialization procedure and a limited number of I/O processes.) The computationalexpense is dominated by the solution procedure of the dry dynamics (96.5%) and in particular, bythe elliptic solver, which accounts for 86.1% of the total computational time. The physics compo-nent is defined as the components unrelated to the dry dynamics, i.e. all tracer transport processesand subgrid-scale physics, such as the microphysics parameterization scheme. In CSLAM-NH, thiscomponent accounts for the remaining 3.2% of the total computational time.In the current development of CSLAM-NH, the conjugate-residual method in Skamarock et al.(1997) is used to solve the Helmholtz equation that arises from the semi-implicit time-steppingscheme. The conjugate-residual solver is flexible and allows the Helmholtz equation to have avariable coefficient structure. However, without a good preconditioner, the scheme has a slow con-vergence rate. A possible way of improving the efficiency of the Helmholtz solver in CSLAM-NHis by relaxing the stopping criterion. Analogous to the results found in Smolarkiewicz et al. (1997)for the anelastic equations, Skamarock et al. (1997) noted that for elastic systems (such as that inCSLAM-NH), the stopping criteria in conjugate residual solvers need not converge to machine pre-94Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatescision, which is the current implementation in CSLAM-NH. Rather, the stopping criterion can bedetermined through a more physical justification of a small change in the divergence solution rela-tive to the mean flow. An examination of whether CSLAM-NH model behaviour changes with morerelaxed levels of convergence will help determine if the stopping criterion can be eased off to im-prove efficiency. In CSLAM-NH, each call to the Helmholtz solver in the iterative centered-implicitscheme may also require a different level of convergence.Other elliptic solvers that could be explored to increase the efficiency include the use of pre-conditioners in residual solvers (Thomas et al., 1997) and the multigrid method (Fulton et al., 1986;Buckeridge, 2010). In a multi-processor environment, an efficient parallelizable geometric multi-grid solver has recently been explored by Mu¨ller and Scheichl (2014). The authors demonstratedthat comparing to other elliptic solvers (algebraic multigrid methods and preconditioned conjugateresidual solvers), the geometric multigrid method is roughly 5 to 10 times faster. The use of sucha solver can likely help improve the efficiency of CSLAM-NH by an order of magnitude and for itto become on a par with that of the Eulerian solver, closing the gap between the two computationaltimes needed when the elliptic solve dominates the run time in CSLAM-NH (e.g., Fig. 4.12a). Asthe purpose of this dissertation is to develop a nonhydrostatic semi-implicit cell-integrated semi-Lagrangian solver, and not Helmholtz solvers (as these are an entire field in itself), the developmentof efficient Helmholtz solvers is not explored here. However, when CSLAM-NH is to be imple-mented in a more comprehensive weather and climate model, advances in the field such as those inMu¨ller and Scheichl (2014) can be considered.Passive tracers are not physically coupled to the dry dynamics, but are passively advected usingthe solution of the flow at the end of each dry dynamics step. Using a single processor, the introduc-tion of more passive tracers, in addition to the existing moisture species (water vapour, cloudwater,and rainwater), therefore should increase the computational time in a linear fashion. Shown in Fig.4.12a are the CPU times for CSLAM-NH (∆t = 50 s) and an Eulerian solver. This Eulerian solveris the same as the second comparison Eulerian solver described in section 4.3.3 and is run with thesame number of tracers with ∆t = 20 s (and a small time step ∆τ = 2 s). As expected, the timesincrease linearly with each additional tracer. The maximum number of 16 passive tracers is limitedby the available memory on the workstation.As mentioned earlier, for a small number of tracers, the computational time in CSLAM-NHis dominated by the elliptic solver. As the number of tracers in the model increases, the com-putational expense shifts towards expenses related to the tracer transport scheme. In each call tothe CSLAM-NH tracer transport (advection) subroutine, computations using the explicit CSLAMtransport scheme and of the corrective fluxes for the “semi-implicit correction term” are made. Thecomputations of the corrective fluxes increases the computational expense by about 26% per time95Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinates0 2 4 6 8 10 12 140 2 4 6 8 10 12 14 16 18Number of passive tracersCPU time (s) x102160 2 4 6 8 10CPU time (s) 0 200 400 600 80024681012Number of passive tracersT CSLAM−NH/T Eulerianx103(a) (b)(c)CSLAM−NHEulerianCSLAM−NHEulerianFigure 4.12: Comparison of the computational times using CSLAM-NH and an Eulerian split-explicit solver. (a) CPU times (in seconds) for the two solvers for a range of tracercounts. (b) Projected computational times using best linear fits. (c) Ratio of CPU timesof CSLAM-NH to the Eulerian solver.step (in return for maintaining consistent transport with air density) as compared to using only theexplicit CSLAM transport scheme.Because the total required computational time is almost directly proportional to the numberof tracers, one can approximate the time required (T ) for a number of tracers (#tracers) using thefollowing best linear fits:TCSLAM−NH = 10.017 ·#tracers +1497.75 (r2 = 0.9895), (4.39)TEulerian = 5.060 ·#tracers +132.63 (r2 = 0.9987). (4.40)The best linear fits are plotted as dashed lines in Fig. 4.12b. The y-intercept represents the compu-tational base time of the respective solver (mostly dependent on the dry dynamics; in CSLAM-NH,mostly on the Helmholtz solver), and the slope indicates the cost associated with advecting addi-tional tracers. For this case (i.e., the test flow and the corresponding maximum stable time-stepsize), the cost of the cell-integrated transport scheme in CSLAM-NH is roughly double of that in96Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesthe Eulerian solver (estimated using the ratio of the two slopes). The cost of the tracer transportschemes can also be visualized by projecting the computational time for a very large number oftracers, at which point, tracer advection dominates the computational expense. For example, Fig.4.12c shows that, when projected for 800 tracers, the ratio of the two solvers’ computational timesapproaches two.For the current version of CSLAM-NH to be competitive with the Eulerian split-explicit solverfor a large number of tracers, the maximum CSLAM-NH stable time-step size (based on this moistneutral flow case) will need to be roughly five times of that of the Eulerian solver. To demonstratethis, the Eulerian time step is reduced to ∆t = 10 s. The linear relationship between CPU time andnumber of tracers for the Eulerian solver [as in (4.39) and (4.40)] isTEulerian = 10.653 ·#tracers +218.41 (r2 = 0.9977). (4.41)For the nonlinear test cases, the maximum stable CSLAM-NH time-step size is found to be lim-ited by self-intersecting departure cells (as in Fig. 4.8). Application of artificial dissipation inthese regions may help prevent the occurrence of these departure cells, and increase the stability ofCSLAM-NH.The results in this section indicate that the cost effectiveness of scalar transport in CSLAM-NH is limited in a single-processor setting. The Eulerian split-explicit solver uses small time stepsto handle fast-moving acoustic waves and a Runge-Kutta 3rd-order time-integration scheme, thatrequires solving the dry dynamics three times. Despite having to transport passive scalars threetimes in the Runge-Kutta time-stepping scheme, the Eulerian solver is still more efficient. The scalartransport code in CSLAM-NH is likely under-optimized in minor aspects. For example, for subgrid-cell reconstruction purposes, ghost cells are filled at each time step in the advection subroutine inCSLAM-NH. The Eulerian solver does not include filling of any ghost cells (whereas in a limitedarea model, the lateral boundary conditions will likely require some ghost cell specifications).The performance of CSLAM-NH will likely compare very differently on multiple processors,based on the results related to the recent effort of Erath et al. (2012). The authors implemented andoptimized the CSLAM transport scheme to run in HOMME, one of the dynamical core options inCAM (see the introduction of this chapter for definitions of the acronyms). CAM with the dynamicalcore HOMME is collectively referred as ‘CAM-SE’ [since HOMME uses a spectral element (SE)method] (Dennis et al., 2012). Comparing to the existing CAM-SE tracer transport scheme, whichis a three-stage second-order Runge-Kutta time-stepping scheme, the CSLAM transport scheme canuse a larger time step and requires only one call to CSLAM per time step (as opposed to three callsin the existing scheme). The computational performance of using the CSLAM transport scheme was97Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesfound to be much more efficient and more scalable than the CAM-SE advection scheme. The authorsfound that the CSLAM scheme was able to outperform the existing scheme, not because of fewerfloating-point operations, but largely because the CAM-SE advection scheme required much more(approximately ten times) communication time than CSLAM when a large number of processorsare used. Therefore, in a parallel environment and for a large number of tracers, CSLAM-NH maystill have an advantage over an Eulerian solver that uses a multi-stage Runge-Kutta time-steppingtransport scheme.Finally, another factor related to the computational efficiency in the solvers is the number oftimes the physics parameterization schemes are called. In both solvers, a straightforward time-splitting approach is used to couple the dynamics to the physics (in this case, to the cloud micro-physics subroutine). The time-splitting approach is defined as the sequential calling of the dynamicsand physics, so that the solution from the dry dynamics is used as input to the physical parameteri-zation schemes. In both solvers, the microphysics scheme is called once at each time step. Becauselarger time steps are used in CSLAM-NH, fewer calls are made to the microphysics scheme, whichhelp reduce the computational expense.4.4 SummaryA nonhydrostatic atmospheric solver (CSLAM-NH) that uses a new discrete formulation of thesemi-implicit continuity equation for cell-integrated semi-Lagrangian transport schemes is furthertested for flows over idealized orography. Here, the solver using the Conservative Semi-LAgrangianMulti-tracer (CSLAM) transport scheme is tested against various idealized mountain-wave casesand exhibits accurate and stable behaviour under the influence of a terrain-following height coordi-nate. An implicit Rayleigh damping layer is also implemented in this extended version of CSLAM-NH to help prevent unphysical reflection of vertically-propagating gravity waves at the model top.The new discrete semi-implicit continuity equation used in CSLAM-NH allows for a straight-forward implementation of consistent flux-form equations for scalars in the model. This aspect ofthe model is important in ensuring inherent mass conservation of these scalars, such as moistureand chemical species, and may prove to be important in longer NWP and climate simulations. Thetime integration of both the gravity and acoustic waves are handled implicitly in the solver using aniterative centered-implicit scheme. The iterative scheme allows for larger maximum stable time stepsizes at the expense of solving the linear Helmholtz problem in the dynamics more than once. Asthe same time step size is used for tracer transport in CSLAM-NH, these larger time steps will con-sequently reduce the number of tracer advection steps per simulation, which may help compensatefor the added expense. Moreover, at each time step, the solution procedure for the dry dynamics98Chapter 4: Extension to a two-dimensional (x-z) nonhydrostatic solver with transformed vertical coordinatesis carried out only once, and then followed by the transport of any tracers in the solver. In largeclimate and chemistry models, the computational cost associated with the parameterized physicsand the transport of the many (order of 102) tracer species is likely to outweigh that associated withthe dynamics. The computational ‘burden’ of tracer transport makes the selection of a multi-tracerefficient transport scheme, such as the CSLAM transport scheme, important.Three idealized test cases available from the literature were used to verify the stability and accu-racy of the proposed solver over topography. Simulations of linear hydrostatic and nonhydrostaticmountain waves compared well with numerical solutions from the literature. The simulation of ahighly nonlinear wave-breaking case of the 11 January 1972 Boulder windstorm highlighted theability of the solver to handle highly-sheared flow at large time steps. Due to the strong nonlinear-ity of the flow, the simulations from the models used in the intercomparison study of Doyle et al.(2000) varied in their fine-scale features. Although there is limited predictability of the precisionof these features, all models, including CSLAM-NH (the simulation of which is presented here),showed similar main features of the windstorm, such as the locations of the wave breaking regionsand hydraulic jump downstream of the mountain. Finally, moist nearly neutral orographic flowsbased on Miglietta and Rotunno (2005) are tested. Two mountain profiles were used: a lower 700-m tall mountain and a much higher 2 km mountain. For the lower mountain case, CSLAM-NHshows comparable results with those in Miglietta and Rotunno (2005), including downstream andupstream regions of subsaturated air. For the higher mountain case, there is more rain spillover tothe leeside of the mountain as compared to the results presented in Miglietta and Rotunno (2005).However, similar solutions are found using another comparison Eulerian split-explicit model, whichsuggests that certain aspects (e.g., initialization) of the model other than model formulation may becausing the discrepancy, and that the discrepancy is not specific to CSLAM-NH.In its current state of development, CSLAM-NH is a two-dimensional prototypical nonhydro-static atmospheric solver in Cartesian geometry that has shown promising potential for weatherand climate applications. Attractive features of this solver include the consistent formulation ofthe semi-implicit cell-integrated semi-Lagrangian continuity and scalar conservation equations, inconjunction with the inherently conservative multi-tracer CSLAM transport scheme. Further devel-opment work (e.g. implementation of three-dimensional CSLAM transport, extension of the schemeto a sphere) remains for the solver to be implemented as a dynamical core in a full NWP and climatemodel.99Chapter 5Conclusion5.1 Summary of contributionsIn this dissertation, an alternative cell-integrated semi-Lagrangian (CISL) semi-implicit nonhydro-static solver for the atmosphere is presented. Existing CISL semi-implicit solvers either do not useconsistent numerical schemes in transporting air and tracer mass and require iterations to improveconsistency between the discrete equations (Thuburn et al., 2010), or solve a nonlinear Helmholtzequation, as proposed by Thuburn (2008).A novel, easy-to-implement flux-form semi-implicit time-stepping scheme for CISL solvers isproposed here. This scheme features a new discrete semi-implicit continuity equation that can bestraightforwardly extended to other constituent-conservation equations. The approach works byeliminating the dependence on a mean reference state. The implicit system still forms a linearHelmholtz equation that can be solved using existing elliptic solvers. The proposed solver is testedin two dimensions using nine idealized test cases. Most of these cases are used by other modeldevelopers as benchmark tests for dynamical cores in atmospheric models. Results from these testcases show that the new nonhydrostatic solver produces comparable results, and shows potential forfurther development and implementation in a more comprehensive weather model.The development of the two-dimensional nonhydrostatic prototype solver for the atmospherecan be summarized as follows.• The initial testbed for the prototype solver is the shallow-water equations (SWE) system.Chapter 2 provides a detailed description of the proposed discrete semi-implicit height equa-tion. The SWE solver utilizes the conservative semi-Lagrangian transport scheme CSLAMto ensure local mass conservation. The SWE solver, named CSLAM-SW (after the transportscheme), uses a new semi-implicit formulation to ensure that air mass and tracer constituentsare transported in a consistent manner. CSLAM-SW is shown to be stable and accurate at aCourant number of 2.5 for a highly nonlinear barotropic jet. Results from a different solverwith inconsistent transport schemes showed that numerical errors in the passive advection of100Chapter 5: Conclusionthe specific concentration of a constituent grows with time step size and nonlinearity of theflow. Moreover, the application of shape-preserving schemes to the tracer mixing ratio be-comes ineffective. Using CSLAM-SW, these errors are successfully eliminated (reduced tomachine roundoff).• Nonhydrostatic modelling using the CSLAM transport scheme and a semi-implicit time in-tegration scheme is performed for the first time (Chapter 3). The proposed discrete semi-implicit height equation for the SWE solver in Chapter 2 is extended and tested for the two-dimensional nonhydrostatic moist Euler equations in the vertical plane. The nonhydrostaticsolver is coupled to an existing microphysics scheme to represent subgrid-scale diabatic ef-fects due to phase changes. CSLAM-NH solutions for a density-current test case and gravity-wave test compare well with those presented in the literature. The gravity-wave test highlightsthe need to integrate the buoyancy terms in the vertical momentum equation in an implicitmanner (addressed in Chapter 4). A third idealized test case of a 2D squall line shows com-parable results to an existing Eulerian nonhydrostatic solver [similar to the one used in theWeather Research and Forecasting Advanced Research and Weather (WRF-ARW) model].The mixing ratio error induced by the lack of numerical consistency at the end of a two-hoursquall line simulation is shown to be on the order of 1 part in 100.• Test cases with the inclusion of topography are presented in Chapter 4. The test cases pre-sented in this chapter include linear mountain waves, a nonlinear downslope windstorm, andorographic cloud formation. Results from CSLAM-NH for the linear mountain waves and thenonlinear downslope windstorm are comparable to those found in the literature. In the downs-lope windstorm case, the maximum vertical Courant number was approximately three, beforeCSLAM-NH eventually becomes unstable due to self-intersecting departure cells. These self-intersecting cells appear in regions with strong wind gradient, such as in the wave breakingregion in the downslope windstorm simulation. In one of the moist orographic flow cases, thespatial distribution of rain differs from that presented in the literature, but a second Euleriansolver shows virtually the same results as CSLAM-NH. This second set of results indicatesthat factors other than the solver component may be in play, e.g., sensitivity to the initial moiststatic stability.In summary, the nonhydrostatic solver CSLAM-NH possesses the following desirable features:1. It inherently conserves mass using the ‘Conservative Semi-LAgrangian Multi-tracer’ (CSLAM)transport scheme.101Chapter 5: Conclusion2. It consistently transports scalar variables using a new semi-implicit formulation of the conti-nuity equation.3. It does not depend on a mean reference state.4. It prevents spurious generation/removal of tracer mass through the use of shape-preservingfilters.5. It requires an elliptic solver that resembles those in traditional semi-Lagrangian semi-implicitsolvers.6. It requires a single call to the conservative semi-Lagrangian transport scheme at each timestep for each tracer.5.2 Impact on the broader communityThere is growing research interest in the influence of trace gases and aerosol concentrations on theatmosphere. Trace gases, such as carbon dioxide, methane, and ozone, make up only a very smallportion of atmospheric constituents. However, their limited concentrations can still have significantcontributions to the radiation budget in the atmosphere, as some are highly radiatively-absorbingspecies (Petty, 2006). Other aerosols such as dust particles and sea spray can be advected by thejet stream and be carried far from their source locations. These aerosols, aside from having a directimpact on the radiation budget, also have an indirect effect through cloud microphysical processes.Dust particles and sea spray can serve as cloud condensation nuclei or ice nuclei, and enhanceprecipitation formation. For example, Creamean et al. (2013), through winter observations from afield campaign called Calwater in the Sierra Nevada region, showed that the long-range transportof dust and biological aerosols from Asia and the Sahara desert likely contributed to the positivefeedback mechanism in orographic precipitation formation in the region.The impact of these aerosols and trace gases on the global radiation and water budget is of-ten studied through atmospheric chemistry and global climate models (GCMs). The representationof physical and chemical processes related to these trace gases in current GCMs and atmosphericchemistry models is becoming more sophisticated. The number of prognostic aerosols and chemicalspecies in current comprehensive atmospheric models far exceeds the number of scalars (roughly10) required by the model dynamics alone. For example, in the fully-coupled online chemistrymodel in the Weather Research and Forecasting – Chemistry model (WRF-CHEM) (Grell et al.,2005), the basic configuration involves 39 basic prognostic chemical species, 34 aerosols, and about10 other variables associated with organic aerosols. (An online chemistry model has the advan-tage of allowing a feedback mechanism for the aerosol radiative impact back to the meteorological102Chapter 5: Conclusionmodel.) In the chemistry version of NCAR’s Community Atmospheric Model (CAM1), the numberof prognostic tracers is on the order of 100.Cloud microphysical schemes in NWP models are also becoming more sophisticated. Cloudmicrophysical schemes are explicit representations of cloud processes in a model (as opposed toimplicit representations using convective parameterizations). In bulk microphysical schemes, hy-drometeor interactions are represented by pre-defined (typically empirical) functional relations forthe particle-size distribution. The simplest microphysical scheme only predicts hydrometeor mix-ing ratios (these are referred to as single-moment schemes), and only predicts for warm-rain speciessuch as cloud and rain droplets. More advanced microphysical schemes now typically includesix classes of hydrometeors, with the addition of ice, snow, graupel, and hail. In double-momentschemes, the concentrations are also predicted to help give a better estimate of particle size distri-butions (Stensrud, 2007). The total number of prognostic variables then increases from two to 12.The growing demand for representation of these tracers and cloud processes in our models calls formore accurate and efficient solvers.Another desirable property when modelling multiple tracers is the ability to maintain relativetracer concentration ratios of long-lived chemical tracers. The inability to preserve the functional re-lations among tracers can be interpreted as erroneous numerical mixing of these tracers (Lauritzenet al., 2011). The relationships among tracers can be linear or nonlinear. Only fully Lagrangiantransport schemes (such as those used in direct numerical simulations) are likely able to handleboth, whereas Eulerian and semi-Lagrangian schemes should ideally preserve linear relations. Re-call that shape-preservation is a desirable property in the case of tracer transport. Ovtchinnikov andEaster (2009) identified a problem when monotonocity constraints are applied in the transport ofindividual tracers. The shape-preserving filters studied were found to be nonlinear algorithms thatcould not preserve relative tracer concentrations. By transporting three inter-related aerosol andcloud particles in an idealized cloud simulation, they found that the percentage error in the (sup-posedly constant) total sum of the three mixing ratios can reach 30%. Based on the same notion,inconsistent transport can also lead to the destruction of the (ideally conserved) total mixing ratiosof inter-related tracers. Testing of functional relation preservation in CSLAM-NH is outside thescope of the dissertation. However, the shape-preserving filter described in Barth and Jespersen(1989) and implemented in the CSLAM component of this solver has been shown to preserve tracerrelations fairly well in a recent comparison study of transport schemes by Lauritzen et al. (2014).Recent testing using the CSLAM transport scheme for advection of tracers in CAM showed that1The Community Atmosphere Model (CAM) is developed at the National Center for Atmospheric Research (NCAR)for the weather and climate research communities. CAM is the atmospheric component of the broader Community EarthSystem Model (CESM) [one of the participating climate models used in the Intergovernmental Panel on Climate Changereport (AR5)].103Chapter 5: Conclusionthe transport scheme is computationally efficient and highly scalable in a parallelized environment(Erath et al., 2012). The adaptability of CSLAM for a large number of tracers is shown through theirsimulations with up to 800 tracers. However, the numerical scheme used to solve the air densitycontinuity equation is not consistent with that for the tracer continuity equation. The continuityequation in Erath et al. (2012) is solved using the spectral element dynamical core (HOMME). Theissue of numerical consistency therefore remains. The method proposed in this dissertation can beused to provide consistent transport and ensure tracer mass conservation in the model.5.3 Further workSeveral aspects of the CSLAM-NH dynamical solver developed in this work will warrant moreresearch if it is to be deployed in a more comprehensive weather and climate model. The maintopics include (i) dimensionality, (ii) potential ways to enhance CSLAM-NH stability, and (iii)cost-effectiveness in a parallelized environment.5.3.1 DimensionalityAt its current state of development, CSLAM-NH remains a two-dimensional prototype solver in thevertical plane. Despite being a two-dimensional solver, the idealized test cases chosen in this disser-tation are representative of typical meteorological flows. The governing nonhydrostatic equationspermit vertical acceleration; the solver demonstrates that it is able to accurately and stably handlestrong vertical motion. The nonhydrostatic prototype solver in the vertical plane (Chapters 3 and4) is, therefore, a logical step forward from a shallow-water model (Chapter 2), in the developmentof a full atmospheric solver. Alas, the world is not flat; the atmosphere is three-dimensional (orfour, including temporal evolution) and the solver will need to include this third spatial dimension.The dimensionality of CSLAM-NH is limited by the inherently conservative transport scheme. Toextend a transport scheme to three dimensions, various paths are given in Lauritzen et al. (2011)(their Chapter 8), and a brief description is given here.The first option is a very rigorous approach of extending the CSLAM transport scheme to carryout three-dimensional remapping. The fully semi-Lagrangian version of the CSLAM transportscheme was developed as a two-dimensional transport scheme. Thus far, the CSLAM transportin the vertical plane tested here is a straightforward implementation of the existing infrastructureusing a terrain-following vertical coordinate. The work required to extend the transport schemeto three dimensions will require new reconstruction and remapping methods. In two dimensions,the transport scheme requires a 2D reconstruction function and a remapping scheme that relies onknowledge of the departure cell boundary. In three dimensions, using a similar notion of applying104Chapter 5: Conclusionthe Gauss-Green theorem will require some surface integral of the departure cell (assuming planarand quadrilateral faces in three dimensions, the cell will be a hexahedron). The surface integral canbe approximated as the sum of the area integral over each face, which can each be decomposedinto separate line integrals as is done in the 2D version of CSLAM. Creation of the search algo-rithm for the cell boundaries (i.e., each face’s boundaries) in three dimensions will be a challengingtask. Garimella et al. (2007) has, however, demonstrated that such a remapping algorithm works inCartesian geometry using piecewise linear reconstruction functions.The second option is to employ a Lagrangian vertical coordinate as is done in several fullysemi-Lagrangian solvers, and originally proposed by Lin (2004). The use of a Lagrangian verticalcoordinate assumes that transport is done along a horizontal surface that moves with the flow. Thisassumption eliminates the vertical advection terms in the discretized equations of motion, and theremaining horizontal advection is calculated by the two-dimensional CSLAM transport scheme. Forthe hydrostatic approximation, as is used in Lin (2004), the vertical coordinate (pressure) position ofeach model level can be determined using the hydrostatic equation, with the model top specified ata fixed pressure. The distribution of the prognostic variables in the Lagrangian column will becomeless resolved in regions that are more vertically divergent. Periodic remapping from the verticalLagrangian column back to an Eulerian reference grid is therefore needed. For a nonhydrostaticsolver that uses a terrain-following height coordinate, the vertical position can be determined usingthe prognostic contravariant vertical velocity. In nonhydrostatic flows, vertical remapping will needto be performed more frequently than for hydrostatic flows because strong vertical accelerations willlikely lead to more deformed flows. The advantage of this approach over the fully three-dimensionalCSLAM transport is that this implementation would require very little modification to the methodproposed in this dissertation.The third option is a ‘cascade approach’, similar to those employed in Nair et al. (2002) and Zer-roukat et al. (2002). The cascade approach is similar to using a Lagrangian vertical coordinate whereremapping is not done in a fully three-dimensional manner. In Nair et al. (2002) and Zerroukat et al.(2002), sequential application of one-dimensional conservative remapping along each coordinate isperformed to achieve conservative remapping in two- to three-dimensions. The cascade approachutilizes an intermediate grid as a “stepping stone” for remapping along each coordinate from theEulerian grid to the Lagrangian grid. In CSLAM-NH, one could apply two-dimensional remappingonto an intermediate grid in x-z, and then apply one-dimensional conservative remapping in the y-direction from the intermediate grid to the Lagrangian grid (or intermediate grid in x-y space, andthen, in the z-direction). As mentioned in Lauritzen (2007), directional biases may occur in split-ting schemes, but can be reduced by alternating the sequence of splitting (directions) every time-stepand/or using smaller (such as halved) time-steps per sweep to reduce the errors. A comparative study105Chapter 5: Conclusionof different splitting methods, such as that in Schneider and Bott (2014), for extending CSLAM-NHto three dimensions would be useful.5.3.2 Enhancement of CSLAM-NH stabilityModel time steps are limited by two main factors: accuracy and stability. The main advantage ofusing semi-Lagrangian schemes is their larger stable time steps, such that accuracy becomes theonly limiting factor when choosing time step sizes. In other words, the model time step size ischosen to acquire an accurate representation of the physical phenomena of interest.For the more complex and nonlinear test cases presented, the stable CSLAM-NH time steps areoften limited by self-intersecting departure cells in CSLAM, which is a stability issue. An importantcriterion for any numerical scheme is its robustness in terms of numerical stability, and this criterionis especially important when model output is used as a part of an operational forecasting system. Inan operational forecast model, tuning of various diffusion options often allows enhanced numericalstability by damping unwanted waves. In the current test cases, minimal explicit diffusion is usedin CSLAM-NH.Possible improvements to the stability of CSLAM-NH are described next. As mentioned inChapter 4, higher-order departure cell-edge approximations may help alleviate the time step limit.The recommendation stems from the work of Ullrich et al. (2012). They introduced quadratic fits(as opposed to straight line approximations) of the cell edges and showed that using quadratic-fitcell boundaries led to more accurate transport of a tracer in a strongly nonlinear sheared flow. Onthe contrary, tracing back more points along the cell edge and connecting them with line segmentsshowed less of an improvement. In addition to using curved cell edges, they also used a fourth-orderRunge-Kutta trajectory algorithm to find the departure point of each cell vertex. Their approach wasimplemented for a simplified flux-form version of the CSLAM transport scheme, but the methodcan also be used in the fully semi-Lagrangian version of CSLAM (as used here).Explicit diffusion is often used in numerical weather prediction models to control spurious noiseand to enhance the stability of a numerical scheme. The introduction of dissipation mechanisms inCSLAM-NH to control the appearance of these self-intersecting departure cells might also be ex-plored. For example, in the downslope windstorm case, unlike 90% of the models in Doyle et al.(2000), CSLAM-NH did not use a turbulence parameterization scheme. Self-intersection of de-parture cells is mostly due to strong wind gradients, which often generate turbulence in the realatmosphere (analogous to the Lipschitz stability condition), and some parameterization of turbu-lence will likely help increase the stability of CSLAM-NH.In the case of the detection of a self-intersecting departure cell, the two problematic cell verticescan be swapped to “untwist” the self-intersecting departure cell (Lauritzen, 2013, personal commu-106Chapter 5: Conclusionnication). Like other filters, the application of such a “smoother” will need to be carefully examinedto determine its damping effect. The application will likely be analogous to applying diffusion tosmooth out local regions with strong wind gradients. When such an algorithm is used, it may benecessary to apply a correction to the velocities of the associated grid points to maintain consistencybetween the trajectory computations and the divergence terms.5.3.3 Cost-effectiveness in a parallel environmentAs shown in Chapter 4, the efficiency of the current prototype solver can be further improved.Computational expenses in CSLAM-NH are dominated by the elliptic solver. The elliptic equationin CSLAM-NH is constructed to be similar to those in traditional semi-Lagrangian semi-implicitsolvers. Recently, Mu¨ller and Scheichl (2014) performed a comprehensive comparison of a range ofalgorithms typically used to solve the Helmholtz equation and tested the elliptic solvers in a parallelenvironment with up to 65536 cores, in anticipation of global semi-Lagrangian semi-implicit solversbeing applied for high-resolution simulations. The authors tested and optimized five algorithms,which include two variants of algebraic multigrid solvers, two variants of preconditioned conjugate-gradient solvers, and a geometric multigrid solver. The authors performed a series of weak scalingassumption tests using all five algorithms and strong scaling tests with the geometric multigridsolver.In the weak scaling tests, the number of grid cells (219) per processor remained unchanged.The problem size was then increased (effectively doubling the resolution of each subsequent test).The number of processors was also increased until all 65536 cores were in use. The time stepsize was decreased accordingly so that the Courant number of the problem stayed the same. Theauthors found that three of the five tested solvers showed generally good weak scalability, with thegeometric multigrid solver performing the best. The latter took only 0.14 to 0.17 s per iteration, andsix iterations to converge, using a total time of 0.86 to 1.06 s. Compared to the other four methods,the elliptic solve of their geometric multigrid method is roughly 5 to 10 times faster.The elliptic solver in operational traditional semi-Lagrangian semi-implicit NWP models islikely optimized for operational use. Since the Helmholtz equation in CSLAM-NH is formulated tobe similar to those used in these traditional solvers, existing optimized techniques can be applied.Newer developments, such as those described above from Mu¨ller and Scheichl (2014), are alsoapplicable when updating and optimizing the Helmholtz solver in CSLAM-NH to more moderntechniques.As mentioned in section 5.2, the CSLAM transport scheme has been implemented in the CAM-SE model for tracer advection. Erath et al. (2012) re-designed the CSLAM scheme to work ef-ficiently in a parallel environment by minimizing the amount of inter-processor communication,107Chapter 5: Conclusionwhich was mostly needed in the subgrid-cell reconstruction calculations. They found that the com-putational time needed to advect a second tracer was only 1/23 of the time needed for advecting thefirst, valid up to 1014 processors. With more processors (tested for up to 4056 processors) and alarge number of tracers, the communication time limits the scalability slightly. In the latter case, thecomputational time needed to advect a second tracer was 1/14 of the time for the first tracer.5.4 SummaryCell-integrated semi-Lagrangian transport schemes are relatively recent developments within thefield of atmospheric solvers. Traditional semi-Lagrangian transport schemes are much more widelyused in operational forecasting. At the time of writing this dissertation, many operational NWPcentres are aware of the benefits of using an inherently conservative semi-Lagrangian transportscheme, and therefore, have developed and/or implemented some version of a cell-integrated semi-Lagrangian solver. However, these solvers have only been applied in research mode, and the devel-opment of CSLAM-NH in this dissertation fits into this category. Currently, very few operationalglobal models are nonhydrostatic, but this will likely change: current global model horizontal reso-lution is already approaching 25 km.As research interest in the impact of trace gases and aerosols continues to grow, accurate tracertransport will become more important. Not only may the neglect of numerically consistent transportlead to spurious generation and removal of tracer mass, but also shape-preservation may no longer beensured. CSLAM-NH uses a consistent and conservative transport scheme for all scalars. Consistenttransport schemes will likely be important in ensuring other (stricter) properties such as functionalrelation preservation of long-lived chemical constituents.The consistent semi-implicit continuity equation formulation has only been tested in CSLAM-NH thus far. The implementation of the consistent formulation in other CISL solvers will determinethe generalizability of this scheme.108BibliographyArakawa, A. and V. R. Lamb, 1977: Computational design of the basic dynamical processes of theUCLA general circulation model. Methods in Computational Physics, 17, Academic Press,173–265.Barth, T. J. and D. C. Jespersen, 1989: The design and application of upwind schemes onunstructured meshes. 27th Aerospace Sciences Meeting, 89 (89-0366).Bermejo, R. and A. Staniforth, 1991: The Conversion of Semi-Lagrangian Advection Schemes toQuasi-Monotone Schemes. Mon. Wea. Rev., 120, 2622–2632.Bonaventura, L., 2000: A Semi-implicit Semi-Lagrangian Scheme Using the Height Coordinate fora Nonhydrostatic and Fully Elastic Model of Atmospheric Flows. J. Comp. Phys., 158, 186–213.Bubnova´, R., G. Hello, P. Be´nard, and J.-F. Geleyn, 1995: Integration of the Fully ElasticEquations Cast in the Hydrostatic Pressure Terrain-Following Coordinate in the Framework ofthe ARPEGE/Aladin NWP System. Mon. Wea. Rev., 123, 515–535.Buckeridge, S., 2010: Numerical solution of weather and climate systems. Ph.D. thesis, Universityof Bath, 228 pp.Cordero, E., N. Wood, and A. Staniforth, 2005: Impact of semi-Lagrangian trajectories on thediscrete normal modes of a non-hydrostatic vertical-column model. Q. J. R. Meteorol. Soc., 131,93–108.Coˆte´, J., S. Gravel, A. Me´thot, A. Patoine, M. Roch, and A. Staniforth, 1998: The operationalCMC-MRB global environmental multiscale (GEM) model. Part I: Design considerations andformulation. Mon. Wea. Rev., 126, 1373–1395.Creamean, J. M., et al., 2013: Dust and biological aerosols from the Sahara and Asia influenceprecipitation in the Western US. Science, 339, 1572–1578.Cullen, M. J. P., 2001: Alternative implementations of the semi-Lagrangian semi-implicit schemesin the ECMWF model. Q. J. R. Meteorol. Soc., 127, 2787–2802.Dennis, J. M., et al., 2012: CAM-SE: A scalable spectral element dynamical core for theCommunity Atmosphere Model. Int. J. High Perform. C., 26, 74–89.Doyle, J. D., et al., 2000: An intercomparison of model-predicted wave breaking for the 11January 1972 Boulder windstorm. Mon. Wea. Rev., 128, 901–914.Dukowicz, J. K. and J. W. Kodis, 1986: Accurate conservative remapping (rezoning) for arbitraryLagrangian-Eulerian computations. SIAM J. Sci. Stat. Comp., 8, 305.109Durran, D. R., 1986: Another Look at Downslope Windstorms. Part I: The Development ofAnalogs to Supercritical Flow in an Infinitely Deep, Continuously Stratified Fluid. J. Atmos.Sci., 43, 2527–2543.Durran, D. R., 2010: Numerical Methods for Fluid Dy.ics – With Applications to Geophysics. 2ded., Springer, 516 pp.Durran, D. R. and J. B. Klemp, 1982: On the effects of moisture on the Brunt-Va¨isa¨la¨ frequency. J.Atmos. Sci., 39, 2152–2158.Durran, D. R. and J. B. Klemp, 1983: A Compressible Model for the Simulation of MoistMountain Waves. Mon. Wea. Rev., 111, 2341–2361.Erath, C., P. H. Lauritzen, J. H. Garcia, and H. M. Tufo, 2012: Integrating a scalable and effcientsemi-Lagrangian multi-tracer transport scheme in HOMME. Procedia Computer Science, 9,994–1003.Erath, C., P. H. Lauritzen, and H. M. Tufo, 2013: On mass-conservation in high-orderhigh-resolution rigorous remapping schemes on the sphere. Mon. Wea. Rev., 141, 2128–2133.Fulton, S. R., P. E. Ciesielski, and W. H. Schubert, 1986: Multigrid methods for elliptic problems:A review. Mon. Wea. Rev., 114, 943–959.Gal-Chen, T. and R. Somerville, 1975: On the use of a coordinate transformation for the solutionof the Navier-Stokes equations. J. Comput. Phys.Garimella, R., M. Kucharik, and M. Shashkov, 2007: An efficient linearity and bound preservingconservative interpolation (remapping) on polyhedral meshes. Computers & fluids, 36, 224–237.Grell, G. A., S. E. Peckham, R. Schmitz, S. A. McKeen, G. Frost, W. C. Skamarock, and B. Eder,2005: Fully coupled “online” chemistry within the WRF model. Atmos. Environ., 39,6957–6975.Harris, L. M., P. H. Lauritzen, and R. Mittal, 2011: A flux-form version of the conservativesemi-Lagrangian multi-tracer transport scheme (CSLAM) on the cubed sphere grid. J. Comput.Phys., 230, 1215–1237.Jablonowski, C., 2004: Adaptive grids in weather and climate modeling. Ph.D. thesis, Universityof Michigan, 292 pp.JMA, 2013: Outline of the operational numerical weather prediction at the Japan MeteorologicalAgency. Appendix to the WMO Technical Progress Report on GDPFS and NWP. Tech. rep.,JMA. URL http://www.jma.go.jp/jma/jma-eng/jma-center/nwp/outline2013-nwp/index.htm.Jo¨ckel, P., R. von Kuhlmann, M. Lawrence, B. Steil, C. Brenninkmeijer, P. Crutzen, P. Rasch, andB. Eaton, 2001: On a fundamental problem in implementing flux-form advection schemes fortracer transport in 3-dimensional general circulation and chemistry transport models. Q. J. R.Meteorol. Soc., 127, 1035–1052.110Keller, T. L., R. Rotunno, M. Steiner, and R. D. Sharman, 2012: Upstream-Propagating WaveModes in Moist and Dry Flow over Topography. J. Atmos. Sci., 69, 3060–3076.Klemp, J. B., J. Dudhia, and A. D. Hassiotis, 2008: An Upper Gravity-Wave Absorbing Layer forNWP Applications. Mon. Wea. Rev., 136, 3987–4004.Klemp, J. B. and D. K. Lilly, 1978: Numerical Simulation of Hydrostatic Mountain Waves. J.Atmos. Sci., 35, 78–107.Klemp, J. B., W. C. Skamarock, and J. Dudhia, 2007: Conservative split-explicit time integrationmethods for the compressible nonhydrostatic equations. Mon. Wea. Rev., 135, 2897–2913.Klemp, J. B. and R. B. Wilhelmson, 1978: The Simulation of Three-Dimensional ConvectiveStorm Dynamics. J. Atmos. Sci., 35, 1070–1096.Kwizak, M. and A. Robert, 1971: A semi-implicit scheme for grid point atmospheric models of theprimitive equations. Mon. Wea. Rev., 99, 32–36.Lalas, D. P. and F. Einaudi, 1974: On the correct use of the wet adiabatic lapse rate in stabilitycriteria of a saturated atmosphere. J. Appl. Meteor., 13, 318–324.Laprise, J. and A. Plante, 1995: A class of semi-Lagrangian integrated-mass (SLIM) numericaltransport algorithms. Mon. Wea. Rev., 123, 553–565.Lauritzen, P., 2007: A Stability Analysis of Finite-Volume Advection Schemes Permitting LongTime Steps. Mon. Wea. Rev., 135, 2658–2673.Lauritzen, P. H., 2005: An Inherently Mass-conservative Semi-implicit Semi-lagrangian Model.Ph.D. thesis, University of Copenhagen, Denmark, 283 pp., [Available fromhttp://www.cgd.ucar.edu/cms/pel/papers/phd.pdf.].Lauritzen, P. H., E. Kaas, and B. Machenhauer, 2006: A mass-conservative semi-implicitsemi-Lagrangian limited-area shallow-water model on the sphere. Mon. Wea. Rev., 134,1205–1221.Lauritzen, P. H., E. Kaas, B. Machenhauer, and K. Lindberg, 2008: A mass-conservative version ofthe semi-implicit semi-Lagrangian HIRLAM. Q. J. R. Meteorol. Soc., 134, 1583–1595.Lauritzen, P. H., R. D. Nair, and P. A. Ullrich, 2010: A conservative semi-Lagrangian multi-tracertransport scheme (CSLAM) on the cubed-sphere grid. J. Comput. Phys., 229, 1401–1424.Lauritzen, P. H., W. C. Skamarock, M. J. Prather, and M. A. Taylor, 2012: A standard test casesuite for two-dimensional linear transport on the sphere. Geosci. Model Dev., 5, 887–901.Lauritzen, P. H., P. A. Ullrich, and R. D. Nair, 2011: Atmospheric Transport Schemes: DesirableProperties and a Semi-Lagrangian View on Finite-Volume Discretizations. NumericalTechniques for Global Atmospheric Models, P. Lauritzen, C. Jablonowski, M. Taylor, R. Nair,111T. J. Barth, M. Griebel, D. E. Keyes, R. M. Nieminen, D. Roose, and T. Schlick, Eds., SpringerBerlin Heidelberg, Lecture Notes in Computational Science and Engineering, Vol. 80, 185–250,URL http://dx.doi.org/10.1007/978-3-642-11640-7 8.Lauritzen, P. H., et al., 2014: A standard test case suite for two-dimensional linear transport on thesphere: results from a collection of state-of-the-art schemes. Geosci. Model Dev., 7, 105–145.Lilly, D. K., 1978: A severe downslope windstorm and aircraft turbulence event induced by amountain wave. J. Atmos. Sci., 35, 59–77.Lin, S.-J., 2004: A “vertically Lagrangian” finite-volume dynamical core for global models. Mon.Wea. Rev., 132, 2293–2307.Machenhauer, B., E. Kaas, and P. H. Lauritzen, 2009: Finite volume meteorology. ComputationalMethods for the Atmosphere and the Oceans: Special Volume, R. Temam and J. Tribbia, Eds.,Elsevier, Handb. Numer. Anal., 3–120.Machenhauer, B. and M. Olk, 1997: The implementation of the semi-implicit scheme incell-integrated semi-Lagrangian models. Atmos.-Ocean, 35 (special issue), 103–126.McGregor, J. L., 1993: Economical determination of departure points for semi-Lagrangian models.Mon. Wea. Rev., 121, 221–230.Melvin, T., M. Dubal, N. Wood, A. Staniforth, and M. Zerroukat, 2010: An inherentlymass-conserving iterative semi-implicit semi-Lagrangian discretization of the non-hydrostaticvertical-slice equations. Q. J. R. Meteorol. Soc., 136, 799–814.Miglietta, M. and R. Rotunno, 2005: Simulations of Moist Nearly Neutral Flow over a Ridge. J.Atmos. Sci., 62, 1410–1427.Moorthi, S., R. W. Higgins, and J. R. Bates, 1995: A global multilevel atmospheric model using avector semi-Lagrangian finite-difference scheme. Part II: Version with physics. Mon. Wea. Rev.,123, 1523–1541.Mu¨ller, E. and R. Scheichl, 2014: Massively parallel solvers for elliptic PDEs in numericalweather and climate prediction. Q. J. R. Meteorol. Soc., doi:doi:10.1002/qj.2327.Nair, R. and B. Machenhauer, 2002: The mass-conservative cell-integrated semi-Lagrangianadvection scheme on the sphere. Mon. Wea. Rev., 130, 649–667.Nair, R. D. and P. H. Lauritzen, 2010: A class of deformational flow test cases for linear transportproblems on the sphere. J. Comput. Phys., 229, 8868–8887.Nair, R. D., J. S. Scroggs, and F. H. Semazzi, 2002: Efficient conservative global transportschemes for climate and atmospheric chemistry models. Mon. Wea. Rev., 130, 2059–2073.112Ovtchinnikov, M. and R. Easter, 2009: Nonlinear advection algorithms applied to interrelatedtracers: errors and implications for modeling aerosol-cloud interactions. Mon. Wea. Rev., 137,632–644.Peltier, W. and T. Clark, 1982: Nonlinear mountain waves in two and three spatial dimensions. Q.J. R. Meteorol. Soc., 109, 527–548.Peltier, W. R. and T. L. Clark, 1979: The Evolution and Stability of Finite-Amplitude MountainWaves. Part II: Surface Wave Drag and Severe Downslope Windstorms. J. Atmos. Sci., 36,1498–1529.Petty, G. W., 2006: A first course in atmospheric radiation. Sundog Publishing, Wisconsin, USA.Pinty, J. P., R. Benoit, E. Richard, and R. Laprise, 1995: Simple tests of a semi-implicitsemi-Lagrangian model on 2D mountain wave problems. Mon. Wea. Rev., 123, 3042–3058.Poulin, F. and G. Flierl, 2003: The nonlinear evolution of barotropically unstable jets. J. Phys.Oceanogr., 33, 2173–2192.Priestley, A., 1993: A quasi-conservative version of the semi-Lagrangian advection scheme. Mon.Wea. Rev., 121, 621–629.Rancic, M., 1992: Semi-Lagrangian piecewise biparabolic scheme for two-dimensional horizontaladvection of a passive scalar. Mon. Wea. Rev., 120, 1394–1406.Randall, D., 1994: Geostrophic adjustment and the finite-difference shallow-water equations. Mon.Wea. Rev., 122, 1371–1377.Rasch, P. and D. Williamson, 1990a: Computational aspects of moisture transport in global-modelsof the atmosphere. Q. J. R. Meteorol. Soc., 116, 1071–1090.Rasch, P. J. and D. L. Williamson, 1990b: On shape-preserving interpolation and semi-Lagrangiantransport. SIAM J. Sci. Stat. Comp., 11, 656–687.Rivest, C., A. Staniforth, and A. Robert, 1994: Spurious resonant response of semi-Lagrangiandiscretizations to orographic forcing: Diagnosis and solution. Mon. Wea. Rev., 122, 366–376.Robert, A., 1981: A stable numerical integration scheme for the primitive meteorologicalequations. Atmos.-Ocean, 19, 35–46.Robert, A., T. Yee, and H. Ritchie, 1985: A semi-Lagrangian and semi-implicitnumerical-integration scheme for multilevel atmospheric models. Mon. Wea. Rev., 113, 388–394.RosHydromet, 2011: Annual Joint WMO Technical Progress Report on the Global DataProcessing and Forecasting System (GDPFS) and related numerical weather prediction researchactivities for 2011. Tech. rep., Federal Service for Hydrometeorology and EnvironmentalMonitoring of the Russian Federation. URL http://www.wmo.int/pages/prog/www/Documents/PublicWeb/www/gdpfs/GDPFS-NWP Annualreports11/2011 RussianFed EN.doc.113Rotunno, R., J. B. Klemp, and M. L. Weisman, 1988: A theory for strong, long-lived squall lines.J. Atmos. Sci., 45, 463–485.Schneider, W. and A. Bott, 2014: On the time-splitting errors of one-dimensional advectionschemes in numerical weather prediction models; a comparative study. Q. J. R. Meteorol. Soc.,in press, doi:doi:10.1002/qj.2301.Skamarock, W. C. and J. B. Klemp, 1994: Efficiency and accuracy of the Klemp-Wilhelmsontime-splitting technique. Mon. Wea. Rev., 122, 2623–2630.Skamarock, W. C. and J. B. Klemp, 2008: A time-split nonhydrostatic atmospheric model forweather research and forecasting applications. J. Comput. Phys., 227, 3465–3485.Skamarock, W. C., J. B. Klemp, and M. G. Duda, 2012: A multiscale nonhydrostatic atmosphericmodel using centroidal Voronoi tesselations and C-grid staggering. Mon. Wea. Rev., 140,3090–3105.Skamarock, W. C., P. K. Smolarkiewicz, and J. B. Klemp, 1997: Preconditioned conjugate-residualsolvers for Helmholtz equations in nonhydrostatic models. Mon. Wea. Rev., 125, 587–599.Smith, R. B., 1979: The influence of mountains on the atmosphere. Adv. Geophys., 21, 87–230.Smith, R. B., 1980: Linear theory of stratified hydrostatic flow past an isolated mountain. Tellus,32, 348–364.Smolarkiewicz, P. and J. Pudykiewicz, 1992: A class of Semi-Langragian approximations forfluids. J. Atmos. Sci., 49, 2082–2096.Smolarkiewicz, P. K., V. Grubisˇic´, and L. G. Margolin, 1997: On forward-in-time differencing forfluids: Stopping criteria for iterative solutions of anelastic pressure equations. Mon. Wea. Rev.,125, 647–654.Staniforth, A. and J. Coˆte´, 1991: Semi-Lagrangian integration schemes for atmospheric models – areview. Mon. Wea. Rev., 119, 2206–2223.Stensrud, D. J., 2007: Parameterization schemes: keys to understanding numerical weatherprediction models. Cambridge University Press, Cambridge, UK.Straka, J., R. B. Wilhelmson, L. J. Wicker, J. R. Anderson, and K. K. Droegemeier, 1992:Numerical solutions of a non-linear density current: A benchmark solution and comparisons.Int. J. Numer. Methods Fluids, 17, 1–22.Thomas, S. J., A. V. Malevsky, M. Desgagne´, R. Benoit, P. Pellerin, and M. Valin, 1997: Massivelyparallel implementation of the mesoscale compressible community model. Parallel Comput., 23,2143–2160.Thuburn, J., 2008: A fully implicit, mass-conserving, semi-Lagrangian scheme for the f-planeshallow-water equations. Int. J. Numer. Methods Fluids, 56, 1047–1059.114Thuburn, J., M. Zerroukat, N. Wood, and A. Staniforth, 2010: Coupling a mass-conservingsemi-Lagrangian scheme (SLICE) to a semi-implicit discretization of the shallow-waterequations: Minimizing the dependence on a reference atmosphere. Q. J. R. Meteorol. Soc., 136,146–154.Ullrich, P. A., P. H. Lauritzen, and C. Jablonowski, 2012: Some considerations for high-order‘incremental remap’-based transport schemes: edges, reconstructions, and area integration. Int.J. Numer. Methods Fluids, 71, 1131–1151.Weisman, M. L. and J. B. Klemp, 1982: The dependence of numerically simulated convectivestorms on vertical wind shear and buoyancy. Mon. Wea. Rev., 110, 504–520.Weisman, M. L., J. B. Klemp, and R. Rotunno, 1988: Structure and evolution of numericallysimulated squall lines. J. Atmos. Sci., 45, 1990–2013.Williamson, D. and P. Rasch, 1994: Water vapor transport in the NCAR CCM2. Tellus A, 46,34–51.Williamson, D. L., J. T. Kiehl, V. Ramanathan, R. E. Dickinson, and J. J. Hack, 1987: Descriptionof NCAR community climate model (CCM1). Tech. rep., National Center for AtmosphericResearch.Wong, M., W. C. Skamarock, P. H. Lauritzen, J. B. Klemp, and R. B. Stull, 2014: A compressiblenonhydrostatic cell-integrated semi-Lagrangian semi-implicit solver (CSLAM-NH) withconsistent and conservative transport. Mon. Wea. Rev. (in press).Wong, M., W. C. Skamarock, P. H. Lauritzen, and R. B. Stull, 2013: A cell-integratedsemi-Lagrangian semi-implicit shallow-water model (CSLAM-SW) with conservative andconsistent transport. Mon. Wea. Rev., 141, 2545–2560.Xue, M., K. K. Droegemeier, and V. Wong, 2000: The Advanced Regional Prediction System(ARPS) - A multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I:Model dynamics and verification. Meteorol. Atmos. Phys., 75, 161–193.Yang, X., J. Hu, D. Chen, H. Zhang, X. Shen, J. Chen, and L. Ji, 2008: Verification of GRAPESunified global and regional numerical weather prediction model dynamic core. Chinese Sci.Bull., 53, 3458–3464.Zerroukat, M., N. Wood, and A. Staniforth, 2002: SLICE: A semi-Lagrangian inherentlyconserving and efficient scheme for transport problems. Q. J. R. Meteorol. Soc., 128, 2801–2820.Zerroukat, M., N. Wood, A. Staniforth, A. White, and J. Thuburn, 2009: An inherentlymass-conserving semi-implicit semi-Lagrangian discretisation of the shallow-water equations onthe sphere. Q. J. R. Meteorol. Soc., 135, 1104–1116.Zhang, K., H. Wan, B. Wang, and M. Zhang, 2008: Consistency problem with tracer advection inthe atmospheric model GAMIL. Adv. Atmos. Sci., 25, 306–318.115Appendix ASupporting MaterialsA.1 Numerical schemes for comparisonA.1.1 A two-time-level traditional semi-Lagrangian semi-implicit modelA traditional grid-point semi-implicit semi-Lagrangian model on a staggered C-grid is constructedfor comparison purposes. The scheme uses a forward-in-time off-centering parameter β for numer-ical stability purposes. The discretized system is given byun+1A = ∆t(1+β2)[f vxy−g′δxh]n+1A+Rnu, (A.1)vn+1A = ∆t(1+β2)[− f uxy−g′δyh]n+1A+Rnv , (A.2)hn+1A =−∆t(1+β2)H0(δxu+δyv)n+1A+Rnh +Rn+ 12h , (A.3)whereRnu = und +∆t(1−β2)[f vxy−g′δxh]nd, (A.4)Rnv = vnd +∆t(1−β2)[− f uxy−g′δyh]nd, (A.5)Rnh = hnd−∆t(1−β2)H0(δxu+δyv)nd, (A.6)Rn+ 12h =−∆t(h′δxu+h′δyv)n+ 12d/2, (A.7)and h′ = h−H0. The operators are defined asδxφ =φi, j−φi−1, j∆x; δyφ =φi, j−φi, j−1∆y, (A.8)116φ x = 12(φi, j +φi+1, j), (A.9)φ xy = φ xy= φ yx=14(φi, j +φi, j+1 +φi+1, j +φi+1, j+1). (A.10)The Rn terms define the known terms that are evaluated at time level n and interpolated to thedeparture point. The Rn+12 term is the nonlinear term evaluated by extrapolating values from timelevel n and n− 1 to time level n+ 12 , and interpolated to the estimated mid-point trajectory. Thetime-off-centering parameter β is set to 0.1 for all runs.A.1.2 An Eulerian leap-frog semi-implicit advective modelThe Eulerian C-grid staggering model uses the semi-implicit leap-frog time-stepping scheme andmomentum equations in the advective form. The model has an Asselin time-filter and a time-off-centering parameter (β = 0.1) to eliminate spurious oscillations. Numerical viscosity is alsoapplied for certain test cases (see section 4b). Using the same notations as for the traditional semi-Lagrangian model, the discretized system is given byun+1 = ∆t(1+β)(f vxy−gδxh)n+1+Ru, (A.11)vn+1 = ∆t(1+β)(− f uxy−gδyh)n+1+Rv, (A.12)hn+1 =−∆t(1+β)H0(δxu+δyv)n+1+Rh, (A.13)whereRu = un−1−2∆t(uδxu+ vδyu)n+∆t(1−β )(f vxy−gδxh)n−1, (A.14)Rv = vn−1−2∆t(uδxv+ vδyv)n+∆t(1−β)(− f uxy−gδyh)n−1, (A.15)Rh = hn−1−∆t(1−β)H0(δxu+δyv)n−1−2∆t(h′δxu+h′δyv)n+ 12. (A.16)117"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2014-05"@en ; edm:isShownAt "10.14288/1.0167310"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Atmospheric Science"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivs 2.5 Canada"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/2.5/ca/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "A fully-compressible nonhydrostatic cell-integrated semi-Lagrangian atmospheric solver with conservative and consistent transport"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/46519"@en .