Abstract

Terrain-following coordinates are widely used in operational models but the cut-cell method has been proposed as an alternative that can more accurately represent atmospheric dynamics over steep orography. Because the type of grid is usually chosen during model implementation, it becomes necessary to use different models to compare the accuracy of different grids. In contrast, here a C-grid finite-volume model enables a like-for-like comparison of terrain-following and cut-cell grids. A series of standard two-dimensional tests using idealized terrain are performed: tracer advection in a prescribed horizontal velocity field, a test starting from resting initial conditions, and orographically induced gravity waves described by nonhydrostatic dynamics. In addition, three new tests are formulated: a more challenging resting atmosphere case, and two new advection tests having a velocity field that is everywhere tangential to the terrain-following coordinate surfaces. These new tests present a challenge on cut-cell grids. The results of the advection tests demonstrate that accuracy depends primarily upon alignment of the flow with the grid rather than grid orthogonality. A resting atmosphere is well maintained on all grids. In the gravity waves test, results on all grids are in good agreement with existing results from the literature, although terrain-following velocity fields lead to errors on cut-cell grids. Because of semi-implicit time stepping and an upwind-biased, explicit advection scheme, there are no time step restrictions associated with small cut cells. In contradiction to other studies, no significant advantages of cut cells or smoothed coordinates are found.

1. Introduction

Representing orography accurately in numerical weather prediction systems is necessary to model downslope winds and local precipitation. Orography also exerts strong nonlocal influences: from the latent heat release due to convection, by directly forcing gravity waves and planetary waves, and by the atmospheric response to form drag and gravity wave drag. There are two main approaches to representing orography on a grid: terrain-following layers and cut cells, with the immersed (or embedded) boundary method (Simon et al. 2012) being similar to a cut-cell approach. All methods align cells in vertical columns. Most models are designed for a particular type of grid, and the study by Good et al. (2014) compared cut-cell results with terrain-following solutions implemented within different models. Instead, this study uses a single model to enable a like-for-like comparison between solutions using terrain-following and cut-cell grids.

With increasing horizontal model resolution, the discrete representation of terrain can become steeper, making accurate calculation of the horizontal pressure gradient more difficult when using terrain-following layers (Gary 1973; Steppeler et al. 2002). Numerical errors in this calculation result in spurious winds and can cause numerical instability (Fast 2003; Webster et al. 2003). Cut-cell methods seek to reduce the error that is associated with steep orography.

With terrain-following (TF) layers the terrain’s influence decays with height so that the bottommost layers follow the underlying surface closely while the uppermost layers are flat. There are two main approaches to minimizing errors associated with TF layers. First, by smoothing the effects of terrain with height, the influence of the terrain is reduced, hence errors in the calculated horizontal pressure gradient are also reduced aloft (Schär et al. 2002; Leuenberger et al. 2010; Klemp 2011). However, the error is not reduced at the ground where steep terrain remains unmodified.

Second, numerical errors can also be reduced by improving the accuracy in calculating the horizontal pressure gradient itself. Terrain-following layers are usually implemented using a coordinate transformation onto a rectangular computational domain, which introduces metric terms into the equations of motion. The techniques proposed by Klemp (2011) and Zängl (2012) both involve interpolation onto z levels in order to calculate the horizontal pressure gradient. This gave them the flexibility to design more accurate horizontal pressure gradient discretizations using more appropriate stencils. The technique proposed by Weller and Shahrokhi (2014) involved calculating pressure gradients in the direction aligned with the grid, thus ensuring curl-free pressure gradients and improving accuracy.

Despite their associated numerical errors, TF layers are in widespread operational use (Steppeler et al. 2003). They are attractive because their rectangular structure is simple to process by computer and link with parameterizations, and boundary layer resolution can be increased with variable spacing of vertical layers (Schär et al. 2002).

Cut cells is an alternative method in which the grid does not follow the terrain but, instead, cells that lie entirely below the terrain are removed, and those that intersect the surface are modified in shape so that they more closely fit the terrain. The resulting grid is orthogonal everywhere except near cells that have been cut. Hence, errors are still introduced when calculating the horizontal pressure gradient between cut and uncut cells.

The cut-cell method can create some very small cells that reduce computational efficiency (Klein et al. 2009), and several approaches have been tried to alleviate the problem. Yamazaki and Satomura (2010) combine small cells with horizontally or vertically adjacent cells. Steppeler et al. (2002) employ a thin-wall approximation to increase the computational volume of small cells without altering the terrain. Jebens et al. (2011) avoid the time step restriction associated with explicit schemes by using an implicit method for cut cells and a semiexplicit method elsewhere.

Some studies have shown examples where cut cells produce more accurate results when compared to TF coordinates. Spurious winds seen in TF coordinates are not present with cut cells and errors do not increase with steeper terrain (Good et al. 2014). A comparison of TF and cut cells using real initial data by Steppeler et al. (2013) found that 5-day forecasts of precipitation and wind over Asia in January 1989 were more accurate in the cut-cell model, although this result was dependent on using an old version of a model.

Another alternative method is the eta coordinate, described by Mesinger et al. (1988). This transformation, expressed in pressure coordinates, quantizes the surface pressure at each grid box using prescribed geometric heights. This results in terrain profiles having a staircase pattern which is known as “step” orography. The eta coordinate improves the accuracy of the horizontal pressure gradient calculation compared to the sigma coordinate (Mesinger et al. 1988).

In an experiment of orographically induced gravity waves, Gallus and Klemp (2000) found that horizontal flow along the lee slope was artificially weak in the eta model. Mesinger et al. (2012) offer an explanation for this behavior: air flowing along the lee slope cannot travel diagonally downward but must first travel horizontally, then vertically downward. However, lee slope winds are weakened because some of the air continues to be transported horizontally aloft.

Mesinger et al. (2012) refined the formulation to allow diagonal transport of momentum and temperature immediately above sloping terrain. This arrangement is similar to the finite-volume cut-cell method. The new method improved test results, increasing lee slope winds by 4–5 m s−1 (Mesinger et al. 2012).

This study uses a modified version of the fully compressible model from Weller and Shahrokhi (2014) to enable a like-for-like comparison between terrain-following and cut-cell grids for idealized, two-dimensional test cases from the literature. Section 2 presents the formulation of the terrain-following and cut-cell grids used in the experiments that follow. In section 3 we give the governing equations and outline the model from Weller and Shahrokhi (2014). Section 4 analyzes the results from three advection tests, a test of a stably stratified atmosphere initially at rest, and orographically induced gravity waves. Concluding remarks are made in section 5.

2. Grids

Here we describe the formulation of the terrain-following grids and the method of cut-cell grid construction. The techniques presented are used to define the grids for the experiments in section 4.

Gal-Chen and Somerville (1975) proposed a basic terrain-following (BTF) coordinate defined as

 
formula

where, in two dimensions, is the physical height of the Cartesian coordinate surface at the model level with transformed height , H is the height of the domain, and is the height of the terrain surface. In this formulation z varies between h and H and ranges from 0 to H. Using this coordinate, the terrain’s influence decays linearly with height but disappears only at the top of the domain. An example is shown in Fig. 1a.

Fig. 1.

Examples of (a) BTF, (b) SLEVE, and (c) a cut-cell grid, showing cell edges in the lowest four layers. The full two-dimensional grids are 20 km wide and 20 km high. SLEVE parameters are specified in the resting atmosphere test in section 4c. The cut-cell grid was created by intersecting the terrain surface with a regular grid as described in section 2. Axes are in units of m.

Fig. 1.

Examples of (a) BTF, (b) SLEVE, and (c) a cut-cell grid, showing cell edges in the lowest four layers. The full two-dimensional grids are 20 km wide and 20 km high. SLEVE parameters are specified in the resting atmosphere test in section 4c. The cut-cell grid was created by intersecting the terrain surface with a regular grid as described in section 2. Axes are in units of m.

The smooth level vertical (SLEVE) coordinate proposed by Schär et al. (2002) achieves a more regular TF grid in the middle and top of the domain than the BTF coordinate. The terrain height is split into large-scale and small-scale components, and , such that , with each component having a different exponential decay. The transformation is defined as

 
formula

where the vertical decay functions are given by

 
formula

and and are the scale heights of large-scale and small-scale terrain, respectively. The exponent n was introduced by Leuenberger et al. (2010) in order to increase cell thickness in the layers nearest the ground, allowing longer time steps. Leuenberger et al. (2010) found the exponent has an optimal value of . Choosing gives the decay functions used by Schär et al. (2002). An example of the SLEVE grid can be seen in Fig. 1b.

Most implementations of terrain-following layers use a coordinate system that makes the computational domain rectangular, but introduces metric terms into the equations of motion. Instead, the model employed in this study uses Cartesian coordinates and nonorthogonal grids. By doing so, results from the same model can be compared between terrain-following and cut-cell grids without modifying the equation set or discretization.

Cut-cell grids are generated in a different way to the typical shaving technique described by Adcroft et al. (1997). Starting from a uniform grid, all cell vertices that lie beneath the orography are moved up to the surface. Additionally, to avoid creating very thin cells, all vertices up to above the orography are moved down to the surface. Where all four of a cell’s vertices are moved, the cell has zero volume and so it is removed. Where two vertices at the same horizontal location are moved up to the surface they will occupy the same point; this results in a zero-length edge that is removed to create a triangular cell. Figure 2 shows how a 2 × 3 cell, uniform grid is transformed into a cut-cell grid. Cells and are removed because they have zero volume, and the zero-length edge at point q is removed to create a triangular cell . Point p is moved down because it is within of the surface, avoiding the creation of a very thin cell.

Fig. 2.

Illustration of a cut-cell grid (a) before and (b) after construction. The terrain surface, denoted by a heavy dotted line, intersects a uniform rectangular grid comprising six cells, . The cell vertices, marked by open circles, are moved to the points at which the terrain intersects vertical cell edges, marked by filled circles. Cells that have no volume are removed. Where a cell has two vertices occupying the same point, the zero-length edge that joins those vertices is removed. In this illustration, cells and are removed because they have no volume, and the zero-length edge at point q is removed to create a triangular cell, . Point p is moved down because it is within of the surface, avoiding the creation of a thin cell.

Fig. 2.

Illustration of a cut-cell grid (a) before and (b) after construction. The terrain surface, denoted by a heavy dotted line, intersects a uniform rectangular grid comprising six cells, . The cell vertices, marked by open circles, are moved to the points at which the terrain intersects vertical cell edges, marked by filled circles. Cells that have no volume are removed. Where a cell has two vertices occupying the same point, the zero-length edge that joins those vertices is removed. In this illustration, cells and are removed because they have no volume, and the zero-length edge at point q is removed to create a triangular cell, . Point p is moved down because it is within of the surface, avoiding the creation of a thin cell.

Some small cells are generated but, unlike most cut-cell grids, cells are typically made smaller in height but their width is unaltered. (A grid that has these thin cells can be seen in Fig. 5c.) This technique has the advantage that cells are not shortened in the direction of flow and so there should be no additional constraints on the advective Courant number.

3. Models

Three models are used for the test cases in this study: two linear advection models and a model of the fully compressible Euler equations. All are operated in a two-dimensional xz slice configuration.

The two finite-volume models make use of the upwind-biased multidimensional cubic advection scheme from Weller and Shahrokhi (2014), which is nonmonotonic and not flux corrected. The scheme uses a least squares approach to fit a multidimensional polynomial over an upwind-biased stencil that contains more cells than the number of polynomial coefficients. This fit is used to interpolate cell values onto face values for the discretization of the advection term using Guass’s divergence theorem. Following Lashley (2002) and Weller et al. (2009), the two cells on either side of the face we are interpolating onto are weighted in the least squares fit so that the fit goes nearly exactly through these cell centers but does not go exactly through the other points. This method worked well when used for terrain-following meshes by Weller and Shahrokhi (2014) but can be unstable in the presence of very small cut cells. This is because the least squares fit can generate a larger interpolation weight for the downwind cell than the upwind cell. To overcome this problem, wherever a large downwind cell interpolation weight is calculated by the least squares fit, the weighting of the upwind cell is increased for the least squares fitting and the fit is recalculated. This procedure is repeated until the interpolation weight of the upwind cell is greater than the interpolation weight of the downwind cell. More details of this approach and a study of its behavior is the subject of a future publication.

a. Finite-volume linear advection model

The first model discretizes the linear advection equation in flux form:

 
formula

where is a prescribed velocity field and the tracer density ϕ is interpolated onto cell faces using one of two schemes: first, the centered linear scheme, which takes the average of the two neighboring cell values; second, the upwind-biased cubic scheme. The time derivative is solved using a three-stage, second order Runge–Kutta scheme defined as

 
formula
 
formula
 
formula

where at time level n. This time-stepping scheme is used for consistency with the trapezoidal implicit scheme used for the fully compressible model, described in section 3c. To ensure that the discrete velocity field is nondivergent, velocities are prescribed at cell faces by differencing the streamfunction along the edges from stored at cell vertices.

b. Finite-difference linear advection model

The second model is a modified version of the linear advection model first used by Schär et al. (2002). It uses terrain-following coordinates and it is configured with leapfrog time stepping and either second-order centered differences, or a fourth-order centered difference scheme given by

 
formula
 
formula

and similarly for .

Once again, velocity fields are prescribed using a streamfunction defined at cell vertices [referred to as double staggered grid points by Schär et al. (2002)]. The original version of the code effectively smoothed the orography, interpolating the geometric height z at doubly staggered points from values at adjacent half levels in order to calculate the streamfunction. The modified version used here directly calculates the height at vertices to enable comparisons with the finite-volume model solutions.

c. Finite-volume fully compressible model

The third model is taken from Weller and Shahrokhi (2014) that details a discretization of the fully compressible Euler equations, given by

 
formula
 
formula
 
formula
 
formula

where ρ is the density, is the velocity field, is the gravitational acceleration, is the heat capacity at constant pressure, is the potential temperature, T is the temperature, p is the pressure, is a reference pressure, is the Exner function of pressure, and is the gas constant to heat capacity ratio. The variable μ is a damping function used for the sponge layer in the gravity waves test in section 4d.

The fully compressible model uses the C-grid staggering in the horizontal and the Lorenz staggering in the vertical such that θ, ρ, and are stored at cell centroids and the covariant component of velocity at cell faces. The model is configured without Coriolis forces.

Acoustic and gravity waves are treated implicitly and advection is treated explicitly. The trapezoidal implicit treatment of fast waves and the Hodge operator suitable for nonorthogonal grids are described in the  appendix. To avoid time-splitting errors between the advection and the fast waves, the advection is time stepped using a three-stage, second-order Runge–Kutta scheme. The advection terms of the momentum and θ equations, (7a) and (7c), are discretized in flux form using the upwind-biased cubic scheme.

4. Results

A series of two-dimensional tests are performed over idealized orography. For each test, results on terrain-following and cut-cell grids are compared. The first test from Schär et al. (2002) advects a tracer in a horizontal velocity field. Second, a new tracer advection test is formulated employing a terrain-following velocity field to challenge the advection scheme on cut-cell grids. The third test solves the Euler equations for a stably stratified atmosphere initially at rest, following Klemp (2011). Fourth, as specified by Schär et al. (2002), a test of orographically induced gravity waves is performed. Finally, another advection test is formulated that transports a stably stratified thermal profile in a terrain-following velocity field. No explicit diffusion is used in any of the tests. (The OpenFOAM implementation of the numerical model, grid generation utilities and test cases are available at https://github.com/hertzsprung/tf-cutcell-comparison/tree/shaw-weller-2015-mwr.)

a. Horizontal advection

Following Schär et al. (2002), a tracer is transported above wave-shaped terrain by solving the advection equation for a prescribed horizontal wind. This test challenges the accuracy of the advection scheme in the presence of grid distortions.

The domain width is 301 km, taken as the horizontal distance between the inlet and outlet boundaries. The domain is 25 km high, discretized onto a grid with and . Note that Schär et al. (2002) measured the domain width as 300 km between the outermost cell centers where tracer values are specified. Both formulations create a cell center (or mass point) rather than a cell face (or horizontal velocity point) over the top of the highest peak, which is crucial for the accuracy of the centered advection schemes.

The terrain is wave shaped, specified by the surface height h, such that

 
formula

where

 
formula

where is the mountain envelope half-width, is the maximum mountain height, is the wavelength, , and . On the SLEVE grid, the large-scale component is given by and is the large-scale height, and is the small-scale height. The optimization of SLEVE by Leuenberger et al. (2010) is not used, so the exponent .

The wind is entirely horizontal and is prescribed as

 
formula

where u0 = 10 m s−1, and . This results in a constant wind above , and zero flow at 4 km and below.

The discrete velocity field is defined using a streamfunction . Given that , the streamfunction is found by vertical integration of the velocity profile:

 
formula

A tracer with density ϕ is positioned upstream above the height of the terrain. It has the following shape:

 
formula

with radius r given by

 
formula

where and are the horizontal and vertical half-widths, respectively, and ϕ0 = 1 kg m−3 is the maximum density of the tracer. At , the tracer is centered at so that the tracer is upwind of the mountain and well above the maximum terrain height of 3 km. Analytic solutions can be found by setting the tracer center such that . Tests are integrated forward in time for 10 000 s with a time step of .

The test was executed on the BTF, SLEVE, and cut-cell grids using a centered linear scheme and the upwind-biased cubic scheme. Results were also obtained on BTF and SLEVE grids with the fourth-order scheme from Schär et al. (2002) using the modified version of their code.

Minimum and maximum tracer values and error norms on the BTF, SLEVE, cut-cell, and regular grids are summarized in Table 1, where the error norm is defined as

 
formula

where ϕ is the numerical tracer value, is the analytic value, and is the cell volume.

Table 1.

Minimum and maximum tracer densities (kg m−3) and error norms, defined by Eq. (13), at t = 10 000 s in the horizontal and terrain-following tracer advection tests using centered linear and cubic upwind-biased schemes. For the horizontal advection test, error norms, minimum and maximum values are given for the fourth-order scheme using the modified code from Schär et al. (2002).

Minimum and maximum tracer densities (kg m−3) and  error norms, defined by Eq. (13), at t = 10 000 s in the horizontal and terrain-following tracer advection tests using centered linear and cubic upwind-biased schemes. For the horizontal advection test,  error norms, minimum and maximum values are given for the fourth-order scheme using the modified code from Schär et al. (2002).
Minimum and maximum tracer densities (kg m−3) and  error norms, defined by Eq. (13), at t = 10 000 s in the horizontal and terrain-following tracer advection tests using centered linear and cubic upwind-biased schemes. For the horizontal advection test,  error norms, minimum and maximum values are given for the fourth-order scheme using the modified code from Schär et al. (2002).

The results of the cubic upwind-biased scheme on TF and regular grids are comparable with those for the fourth-order centered scheme from Schär et al. (2002). The error is largest on the BTF grid with , but is significantly reduced on the SLEVE grid with . Advection is most accurate on the cut-cell grid, with approximately half of that on the SLEVE grid. Tracer minima and maxima for the centered linear and fourth-order schemes are lower than those presented by Schär et al. (2002) because no interpolation is used to calculate the streamfunction.

The results of the horizontal advection test show that numerical errors are due either to misalignment of the flow with the grid, or to grid distortions. In the following section, we propose a new test in order to identify the cause of the errors.

b. Terrain-following advection

In the horizontal advection test, results were least accurate on the BTF grid, where the grid was most nonorthogonal and flow was misaligned with the grid layers. Here, we formulate a new tracer advection test in which the velocity field is everywhere tangential to the basic terrain-following coordinate surfaces. On the BTF grid, the flow is then aligned with the grid, but the data in the multidimensional advection stencil are not uniformly distributed because the BTF grid is nonorthogonal. Conversely, on the cut-cell grid, the flow is misaligned with the grid but, except in the lowest layers, the grid is orthogonal. This test determines whether the primary source of numerical error is due to nonorthogonality or misalignment of the flow with grid layers.

The spatial domain, mountain profile, initial tracer profile, and discretization are the same as those in the horizontal tracer advection test, except for the time step . The velocity field is defined using a streamfunction so that the discrete velocity field is nondivergent and follows the BTF coordinate surfaces given by Eq. (1) such that

 
formula

where u0 = 10 m s−1, which is the horizontal wind speed where . The horizontal and vertical components of velocity, u and w, are then given by

 
formula
 
formula

Unlike the horizontal advection test, flow extends from the top of the domain all the way to the ground. The discrete velocity field is calculated using the streamfunction in the same way as the horizontal advection test.

At t = 10 000 s the tracer has passed over the mountain. The horizontal position of the tracer center can be calculated by integrating along the trajectory to find t, the time taken to pass from one side of the mountain to the other:

 
formula
 
formula
 
formula

Hence, we find that x(t = 10 000 s) = 51 577.4 m. Because the velocity field is nondivergent, the flow accelerates over mountain ridges and the tracer travels 1577.4 m farther compared to advection in the purely horizontal velocity field. Tracer height is unchanged downwind of the mountains because advection is parallel to the coordinate surfaces.

Tracer contours at t = 0, 5000, and 10 000 s are shown in Fig. 3 using the centered linear scheme on the BTF grid and cut-cell grid (Figs. 3a and 3b, respectively). At , the tracer is distorted by the terrain-following velocity field. On the BTF grid, the tracer correctly returns to its original shape having cleared the mountain by t = 10 000 s, but this is not the case with centered linear scheme on the cut-cell grid. Here, the tracer has spread vertically due to increased numerical errors when the tracer is transported between layers. Dispersion errors are apparent with grid-scale oscillations that travel in the opposite direction to the wind (Fig. 3d) and some artifacts remain above the mountain peak.

Fig. 3.

Tracer contours advected in a terrain-following velocity field at t = 0, 5000, and 10 000 s using the centered linear scheme on (a) the BTF grid and (b) the cut-cell grid with contour intervals every 0.1. Errors at t = 10 000 s are shown for (c) the centered linear scheme on the BTF grid, (d) the centered linear scheme on the cut-cell grid, (e) the upwind-biased cubic scheme on the BTF grid, and (f) the upwind-biased cubic scheme on the cut-cell grid with contour intervals every 0.01. Negative contours are denoted by dotted lines. The terrain profile is also shown immediately above the x axis.

Fig. 3.

Tracer contours advected in a terrain-following velocity field at t = 0, 5000, and 10 000 s using the centered linear scheme on (a) the BTF grid and (b) the cut-cell grid with contour intervals every 0.1. Errors at t = 10 000 s are shown for (c) the centered linear scheme on the BTF grid, (d) the centered linear scheme on the cut-cell grid, (e) the upwind-biased cubic scheme on the BTF grid, and (f) the upwind-biased cubic scheme on the cut-cell grid with contour intervals every 0.01. Negative contours are denoted by dotted lines. The terrain profile is also shown immediately above the x axis.

A small improvement is obtained on the BTF grid by using the upwind-biased cubic scheme: as seen in Fig. 3e, errors are less than 0.02 in magnitude and errors are confined to the expected region of the tracer. However, results are substantially improved by using the upwind-biased cubic scheme on the cut-cell grid (Fig. 3f). Results on the SLEVE grid are comparable to those on the cut-cell grid except that the artifacts above the mountain peak with the centered linear scheme on the cut-cell grid are not present on the SLEVE grid (not shown).

The errors and tracer extrema for this test are compared with the horizontal advection results in Table 1. In the terrain-following velocity field, tracer accuracy is greatest on the BTF grid. Errors are about 10 times larger on the SLEVE and cut-cell grids compared to the BTF grid.

We conclude from this test that accuracy depends upon alignment of the flow with the grid, and accuracy is not significantly reduced by grid distortions. Error on the BTF grid in the terrain-following advection test is comparable with the error on the SLEVE grid in the horizontal advection test.

c. Stratified atmosphere initially at rest

An idealized terrain profile is defined along with a stably stratified atmosphere at rest in hydrostatic balance. The analytic solution is time invariant, but numerical errors in calculating the pressure gradient can give rise to spurious velocities that become more severe over steeper terrain (Klemp 2011). Cut-cell grids are often suggested as a technique for reducing these spurious circulations (Yamazaki and Satomura 2010; Jebens et al. 2011; Good et al. 2014).

The test setup follows the specification by Klemp (2011). The domain is 200 km wide and 20 km high, and the grid resolution is . All boundary conditions are no normal flow.

The wave-shaped mountain profile has a surface height h given by

 
formula

where is the mountain half-width, is the maximum mountain height, and is the wavelength. For the optimized SLEVE grid, the large-scale component is specified as

 
formula

and, following Leuenberger et al. (2010), is the large-scale height, is the small-scale height, and the optimal exponent value of is used.

Tests were performed with two different stability profiles, both having an initial potential temperature field has and a constant static stability with Brunt–Väisälä frequency N = 0.01 s−1 everywhere, except for a more stable layer of N = 0.02 s−1. Figure 4a shows where this inversion layer is positioned in the two tests: the “high inversion” test follows Klemp (2011), placing the layer between ; the “low inversion” test is designed to challenge the pressure gradient calculations on the cut-cell grid by placing the inversion layer between so that it intersects the terrain.

Fig. 4.

Setup and results of a stratified atmosphere initially at rest. Tests are performed on four grids for two different stability profiles: (a) the placement of the inversion layer in the two profiles. The low inversion is positioned so that it intersects the terrain, shown immediately above the x axis. In each test, the inversion layer has a Brunt–Väisälä frequency N = 0.02 s−1, and N = 0.01 s−1 elsewhere. (b) The maximum magnitude of spurious vertical velocity w (m s−1) with results on BTF, SLEVE, and cut-cell grids using the model from Weller and Shahrokhi (2014), which includes a curl-free pressure gradient formulation. Results for the high inversion test are shown with solid lines and results for the low inversion test are shown with dashed lines.

Fig. 4.

Setup and results of a stratified atmosphere initially at rest. Tests are performed on four grids for two different stability profiles: (a) the placement of the inversion layer in the two profiles. The low inversion is positioned so that it intersects the terrain, shown immediately above the x axis. In each test, the inversion layer has a Brunt–Väisälä frequency N = 0.02 s−1, and N = 0.01 s−1 elsewhere. (b) The maximum magnitude of spurious vertical velocity w (m s−1) with results on BTF, SLEVE, and cut-cell grids using the model from Weller and Shahrokhi (2014), which includes a curl-free pressure gradient formulation. Results for the high inversion test are shown with solid lines and results for the low inversion test are shown with dashed lines.

The Exner function of pressure is calculated so that it is in discrete hydrostatic balance in the vertical direction (Weller and Shahrokhi 2014). The damping function μ is set to 0 s−1. Unlike Klemp (2011), there is no eddy diffusion in the equation set.

The test was integrated forward by 5 h using a time step on the BTF, SLEVE, and cut-cell grids. Maximum vertical velocities are presented in Fig. 4b and are similar on the BTF, SLEVE, and cut-cell grids. For the high inversion test, the largest vertical velocity of 0.37 m s−1 was found on the BTF grid after 400 s, compared with a maximum of ~7 m s−1 found by Klemp (2011) using their improved horizontal pressure gradient formulation. Errors are two orders of magnitude smaller on the cut-cell grid with vertical velocities of ~1 × 10−4 m s−1, but this advantage is lost when the inversion layer is lowered to intersect the terrain. Unlike the result from Klemp (2011), the SLEVE grid does not further reduce vertical velocities compared to the BTF grid. This implies that the numerics we are using are less sensitive to grid distortions.

Good et al. (2014) found the maximum vertical velocity in their cut-cell model was 1 × 10−12 m s−1, which is better than any result obtained here. It is worth noting that our model stores values at the geometric center of cut cells, whereas the model used by Good et al. (2014) has cell centers at the center of the uncut cell, resulting in the center of some cut cells being below the ground (S.-J. Lock 2014, personal communication). This means that the grid is effectively regular when calculating horizontal and vertical gradients. This would account for the very small velocities found by Good et al. (2014).

The results in Fig. 4b have maximum errors that are comparable with Weller and Shahrokhi (2014) but, due to the more stable split into implicitly and explicitly treated terms (described in the  appendix), the errors decay over time due to the dissipative nature of the advection scheme.

In summary, we reproduce the result found by Good et al. (2014) that cut cells can reduce spurious velocities over orography. However, in addition, we find that, with the right numerics, these errors can be very small on a BTF grid. We also find that, if changes in stratification intersect cut cells, spurious velocities on cut-cell grids are comparable with those on TF grids.

d. Gravity waves

The test originally specified by Schär et al. (2002) prescribes flow over terrain with small-scale and large-scale undulations, which induces propagating and evanescent gravity waves.

Following Melvin et al. (2010), the domain is 300 km wide and 30 km high. The mountain profile has the same form as Eq. (20), but the gravity waves tests have a mountain height of . As in the resting atmosphere test, is the mountain half-width and is the wavelength.

A uniform horizontal wind (u, w) = (10, 0) m s−1 is prescribed in the interior domain and at the inlet boundary. No normal flow is imposed at the top and bottom boundaries and the velocity field has a zero gradient outlet boundary condition.

The initial thermodynamic conditions have constant static stability with N = 0.01 s−1 everywhere, such that

 
formula

where the temperature at is . Potential temperature values are prescribed at the inlet and upper boundary using Eq. (22), and a zero gradient boundary condition is applied at the outlet. At the ground, fixed gradients are imposed by calculating the component of normal to each face using the vertical derivative of Eq. (22). For the Exner function of pressure, hydrostatic balance is prescribed on top and bottom boundaries and the inlet and outlet are zero normal gradient.

Sponge layers are added to the upper 10 km and leftmost 10 km at the inlet boundary to damp the reflection of waves. The damping function μ is adapted from Melvin et al. (2010) such that

 
formula
 
formula
 
formula

where is the damping coefficient, is the bottom of the sponge layer, is the top of the domain, x0 = −150 km is the leftmost limit of the domain, and xI = −140 km is the rightmost extent of the inlet sponge layer. The sponge layer is only active on faces whose normal is vertical so that it damps vertical momentum only.

Note that, while the domain itself is 30 km in height, for the purposes of generating BTF grids, the domain height is set to 20 km because the sponge layer occupies the uppermost 10 km.

The simulation is integrated forward by 5 h and the time step, , is chosen so that it scales linearly with spatial resolution and, following the original test specified by Schär et al. (2002), when . Test results are compared between the BTF and cut-cell grids at several resolutions. The spatial and temporal resolutions tested are shown in Table 2. The lowest resolution is the same as that used by Schär et al. (2002), and higher resolutions preserve the same aspect ratio. The vertical resolution is chosen to test various configurations of cut-cell grid. At , the mountain lies entirely within the lowest layer of cells, while at and the mountain peak is aligned with the interface between layers. With increasing resolutions up to , the orography intersects more layers and becomes better resolved. Three of the cut-cell grids are shown in Fig. 5 at , , and . Small cells are visible on the 150-m grid but, on the 200-m grid, the small cells are merged with those in the layer above.

Table 2.

Spatial and temporal resolutions used in the gravity waves test. The resolution of has the same parameters as the original test case specified by Schär et al. (2002). At other resolutions, the vertical resolution is prescribed, and horizontal and temporal resolutions are calculated to preserve the same aspect ratios as the original test case.

Spatial and temporal resolutions used in the gravity waves test. The resolution of  has the same parameters as the original test case specified by Schär et al. (2002). At other resolutions, the vertical resolution is prescribed, and horizontal and temporal resolutions are calculated to preserve the same aspect ratios as the original test case.
Spatial and temporal resolutions used in the gravity waves test. The resolution of  has the same parameters as the original test case specified by Schär et al. (2002). At other resolutions, the vertical resolution is prescribed, and horizontal and temporal resolutions are calculated to preserve the same aspect ratios as the original test case.
Fig. 5.

Cut-cell grids used for the gravity waves and thermal advection tests at (a) , (b) , and (c) . The mountain peak . At and , the grid creation process has merged small cells with the cells in the layer above but, at , small cells have been retained. The full two-dimensional grids are 300 km wide and 30 km high. Axes are in units of m.

Fig. 5.

Cut-cell grids used for the gravity waves and thermal advection tests at (a) , (b) , and (c) . The mountain peak . At and , the grid creation process has merged small cells with the cells in the layer above but, at , small cells have been retained. The full two-dimensional grids are 300 km wide and 30 km high. Axes are in units of m.

The ratio of minimum and maximum cell areas in the various grids is shown in Table 3, providing an indication of size of the smallest cut cells. As expected, there is almost no variation in cell sizes on the BTF grids. Small cells are generated on cut-cell grids at resolutions finer than in which the terrain intersects grid layers.

Table 3.

Cell area ratios of BTF and cut-cell grids used in the gravity waves and thermal advection tests. Cell sizes are almost uniform on BTF grids, but for the cut-cell grids the cell area ratio gives an indication of the smallest cell sizes.

Cell area ratios of BTF and cut-cell grids used in the gravity waves and thermal advection tests. Cell sizes are almost uniform on BTF grids, but for the cut-cell grids the cell area ratio gives an indication of the smallest cell sizes.
Cell area ratios of BTF and cut-cell grids used in the gravity waves and thermal advection tests. Cell sizes are almost uniform on BTF grids, but for the cut-cell grids the cell area ratio gives an indication of the smallest cell sizes.

At , vertical velocities on the BTF and cut-cell grids are visually indistinguishable (not shown). They agree with the high-resolution mass-conserving semi-implicit semi-Lagrangian solution from Melvin et al. (2010). The initial thermal profile is subtracted from the potential temperature field at the end of the integration to reveal the structure of thermal anomalies. The anomalies on the BTF grid with are shown in Fig. 6. A vertical profile is taken at , marked by the dashed line in Fig. 6, with results shown for the BTF grids in Fig. 7a and on the cut-cell grids in Fig. 7b. The position is chosen to be far away from the mountain where the gravity wave amplitude is small in order to better reveal numerical errors. On all grids, potential temperature differences increase with height in the lowest 1200 m at , in agreement with the results seen in Fig. 6. Results are seen to converge on all grids, with the exception of small errors in the lowest layers on the cut-cell grids.

Fig. 6.

Differences in potential temperature between the start and end of the gravity waves test on the BTF grid with . The dashed line at marks the position of the vertical profile in Fig. 7. Axes are in units of m.

Fig. 6.

Differences in potential temperature between the start and end of the gravity waves test on the BTF grid with . The dashed line at marks the position of the vertical profile in Fig. 7. Axes are in units of m.

Fig. 7.

Vertical profiles of potential temperature differences between the start and end of the gravity waves test on (a) the BTF grid and (b) the cut-cell grid. Results are compared with thermal advection tests results, showing differences in potential temperature between the numeric and analytic solutions at t = 18 000 s on (c) the BTF grid and (d) the cut-cell grid. The results are convergent, except for errors found in the lowest layers on the cut-cell grids.

Fig. 7.

Vertical profiles of potential temperature differences between the start and end of the gravity waves test on (a) the BTF grid and (b) the cut-cell grid. Results are compared with thermal advection tests results, showing differences in potential temperature between the numeric and analytic solutions at t = 18 000 s on (c) the BTF grid and (d) the cut-cell grid. The results are convergent, except for errors found in the lowest layers on the cut-cell grids.

To summarize, results of the gravity waves test on all grids are in good agreement with the reference solution from Melvin et al. (2010). The potential temperature field converges, though errors are found in the lowest layers on the cut-cell grids. The source of the errors in the cut-cell grids will be investigated further with an advection test in section 4e.

e. Terrain-following advection of thermal profile

The potential temperature anomalies in the gravity waves test do not converge with resolution when using the cut-cell grids. This may be due to differences in the wind fields between grids, or errors in the advection of potential temperature, among other possible causes. To help establish the primary source of error, a new advection test is formulated in which the initial potential temperature field from the gravity waves test is used. To eliminate any differences in wind fields, the field is advected in a fixed, terrain-following velocity field that mimics the flow in the gravity waves test.

The spatial domain, mountain profile, grid resolutions, and time steps are the same as those in the gravity waves test in section 4d. The terrain-following velocity field is defined by the streamfunction as follows:

 
formula

where is the level at which the terrain-following layers become flat; the domain height is . For , the u and w components of velocity are given by Eq. (15), but has the same form as Eq. (20), hence the derivative is

 
formula

for , , and .

The potential temperature field θ, and its boundary conditions, are the same as those of the initial potential temperature field in the gravity waves test. Following the gravity waves test, the simulation is integrated forward by 18 000 s, by which time the potential temperature initially upwind of the mountain will have cleared the mountain range. Hence, the analytic solution can be found by considering the vertical displacement of the thermal profile by the terrain-following velocity field:

 
formula

where the potential temperature at , , and the transform is given by rearranging Eq. (1).

Enlargements of the error field near the mountain are shown in Fig. 8 at with contours of potential temperature overlaid. Errors are only just visible on the BTF grid with an error of 1.12 × 10−7. However, on the cut-cell grid, the error is about 10 times larger. Advection errors are apparent around mountainous terrain, with small cells having some of the largest errors. These errors are advected horizontally along the lee slope, forming stripes. The same error structure is present on all cut-cell grids.

Fig. 8.

Error in potential temperature (measured in K) in the thermal advection test at a resolution of on (a) the BTF grid and (b) the cut-cell grid. Errors are negligible on the BTF grid, but on the cut-cell grid errors are generated near mountainous terrain and are advected horizontally on the lee side. Contours of the potential temperature field at t = 18 000 s are overlaid. Axes are in units of m.

Fig. 8.

Error in potential temperature (measured in K) in the thermal advection test at a resolution of on (a) the BTF grid and (b) the cut-cell grid. Errors are negligible on the BTF grid, but on the cut-cell grid errors are generated near mountainous terrain and are advected horizontally on the lee side. Contours of the potential temperature field at t = 18 000 s are overlaid. Axes are in units of m.

For comparison with the potential temperature anomalies in the gravity waves test, vertical profiles of potential temperature error are taken at . As seen in Fig. 7c, errors are negligible on the BTF grids, but Fig. 7d reveals significant errors in the lowest layers of the cut-cell grids that were advected from the mountain peaks.

While the magnitude and structure of error on the cut-cell grids in this test differs from potential temperature anomalies in the gravity waves test, results on the BTF grids are in close agreement in both tests but not on the cut-cell grids. Therefore, it is likely that anomalies on the cut-cell grids in the gravity waves test are primarily due to errors in the advection of potential temperature through cut cells.

5. Conclusions

We have presented a like-for-like comparison between terrain-following and cut-cell grids using a single model. Accuracy on the BTF, SLEVE, and cut-cell grids was evaluated in a series of two-dimensional tests.

Across all tests, a high degree of accuracy was achieved for all grids. Even on the highly distorted BTF grid errors were often small in the tests presented here. In the first two tests, tracers were advected by horizontal and terrain-following velocity fields. We found that the accuracy of the upwind-biased cubic advection scheme depended upon alignment of the flow with the grid rather than on grid distortions. Spurious vertical velocities in the resting atmosphere tests were similar on terrain-following and cut-cell grids. In the gravity waves test, vertical velocities were in good agreement with the reference solution from Melvin et al. (2010) across all grids.

Cut-cell grids reduced errors in the horizontal advection test. Conversely, in the terrain-following tracer advection test, errors were large on the SLEVE and cut-cell grids where velocities were misaligned with the grids. Errors were also large on the cut-cell grids in the terrain-following thermal advection test. This result suggests that, in the gravity waves test, potential temperature errors in the cut-cell grids are primarily due to advection errors.

The cubic upwind-biased advection scheme takes an approach for treating small cut cells that differs from other existing approaches by adjusting weightings to ensure that advection remains upwind-biased near small cells. (The implementation of this technique in OpenFOAM is available at https://github.com/hertzsprung/AtmosFOAM/tree/shaw-weller-2015-mwr and will be described in greater detail a future publication.) Combined with semi-implicit time stepping and a new cut-cell generation technique that preserves cell length in the direction of the flow, small cells did not impose additional time step constraints. By using a suitable multidimensional advection scheme and a curl-free pressure gradient formulation, we did not find significant advantages of cut cells or smoothed coordinate systems unlike Good et al. (2014), Klemp (2011), and Schär et al. (2002). In contrast, errors that do not reduce with resolution are on cut-cell grids. No significant problems were found when using basic terrain-following grids.

Acknowledgments

I am grateful to my cosupervisors John Methven and Terry Davies for their valuable input, and to Christoph Schär (ETH) for his assistance in reproducing his advection test results. The authors thank the two anonymous reviewers whose comments greatly improved this manuscript. I am thankful for the NERC studentship, which helped fund this work. Weller is funded by NERC Grant NE/H015698/1.

APPENDIX

Semi-Implicit Treatment of the Hodge Operator

To ensure curl-free pressure gradients, following Weller and Shahrokhi (2014), the covariant momentum component, that is the momentum at the cell face in the direction between cell centers, is used as the prognostic variable for velocity:

 
formula

where is the vector between cell centers and subscript f means “at face f.” The contravariant momentum component, that is the flux across faces, is a diagnostic variable:

 
formula

where is the outward-pointing normal vector to face f with magnitude equal to the area of the face. If U is the vector of all values of and V is the vector of all values of then we can define the Hodge operator as a matrix that transforms V to U:

 
formula

For energy conservation, Thuburn and Cotter (2012) showed that the Hodge operator must be symmetric and positive definite. We define a symmetric suitable for arbitrary 3D meshes:

 
formula

where subscript F denotes midpoint interpolation from two surrounding cell values onto face f:

 
formula

where denotes the two cells sharing face f. Here is the consistent cell center reconstruction of from surrounding values of :

 
formula

where is a tensor and so the inversion of the tensor sum is a local operation that can be calculated once for every cell in the grid before time stepping begins. The implied by this reconstruction of U is likely to be positive definite for meshes with sufficiently low nonorthogonality, although this has not been proved.

The semi-implicit technique involves combining the momentum [Eq. (7a)], continuity [Eq. (7b)], and θ [Eq. (7c)] equations and the equation of state [Eq. (7d)] to form a Helmholtz equation to be solved implicitly, as described by Weller and Shahrokhi (2014). The semi-implicit solution technique with a Hodge operator can be defined by considering only a discretized form of the continuity equation:

 
formula

The divergence is discretized using Gauss’s divergence theorem so that

 
formula

where is the volume of cell c, denotes the faces of cell c, and if points outward from the cell and −1 otherwise. Equation (A7) is now a sum over a sum since is one element of a matrix–vector multiply. To simplify the construction of the matrix for the Helmholtz problem, only the diagonal terms of are treated implicitly. Therefore, is separated into a diagonal and off-diagonal matrix:

 
formula

Equation (A6) can now be approximated by

 
formula

where superscript denotes lagged values taken from a previous iteration or from a previous stage of a Runge–Kutta scheme. This was the approach taken by Weller and Shahrokhi (2014). However, the numerical solution of Eq. (A9) turns out to be unstable when using a large time step on highly nonorthogonal meshes associated with terrain-following layers over steep orography. Improved stability and energy conservation can be achieved by splitting into a diagonal component, which would be correct on an orthogonal grid, and a nonorthogonal correction:

 
formula

where the diagonal matrix and the nonorthogonal correction is . The orthogonal part, , can be treated implicitly in the Helmholtz equation:

 
formula

This form is used for the solutions of the Euler equations in this paper and is stable, with good energy conservation for all of the tests presented.

REFERENCES

REFERENCES
Adcroft
,
A.
,
C.
Hill
, and
J.
Marshall
,
1997
:
Representation of topography by shaved cells in a height coordinate ocean model
.
Mon. Wea. Rev.
,
125
,
2293
2315
, doi:.
Fast
,
J. D.
,
2003
:
Forecasts of valley circulations using the terrain-following and step-mountain vertical coordinates in the Meso-Eta model
.
Wea. Forecasting
,
18
,
1192
1206
, doi:.
Gal-Chen
,
T.
, and
R. C.
Somerville
,
1975
:
On the use of a coordinate transformation for the solution of the Navier-Stokes equations
.
J. Comput. Phys.
,
17
,
209
228
, doi:.
Gallus
,
W. A.
, and
J. B.
Klemp
,
2000
:
Behavior of flow over step orography
.
Mon. Wea. Rev.
,
128
,
1153
1164
, doi:.
Gary
,
J. M.
,
1973
:
Estimate of truncation error in transformed coordinate, primitive equation atmospheric models
.
J. Atmos. Sci.
,
30
,
223
233
, doi:.
Good
,
B.
,
A.
Gadian
,
S.-J.
Lock
, and
A.
Ross
,
2014
:
Performance of the cut-cell method of representing orography in idealized simulations
.
Atmos. Sci. Lett.
,
15
,
44
49
, doi:.
Jebens
,
S.
,
O.
Knoth
, and
R.
Weiner
,
2011
:
Partially implicit peer methods for the compressible Euler equations
.
J. Comput. Phys.
,
230
,
4955
4974
, doi:.
Klein
,
R.
,
K.
Bates
, and
N.
Nikiforakis
,
2009
:
Well-balanced compressible cut-cell simulation of atmospheric flow
.
Philos. Trans. Roy. Soc. London
,
367
,
4559
4575
, doi:.
Klemp
,
J. B.
,
2011
:
A terrain-following coordinate with smoothed coordinate surfaces
.
Mon. Wea. Rev.
,
139
,
2163
2169
, doi:.
Lashley
,
R. K.
,
2002
: Automatic generation of accurate advection schemes on unstructured grids and their application to meteorological problems. Ph.D. thesis, University of Reading, 247 pp. [Available online at http://ethos.bl.uk/OrderDetails.do?uin=uk.bl.ethos.394176.]
Leuenberger
,
D.
,
M.
Koller
,
O.
Fuhrer
, and
C.
Schär
,
2010
:
A generalization of the SLEVE vertical coordinate
.
Mon. Wea. Rev.
,
138
,
3683
3689
, doi:.
Melvin
,
T.
,
M.
Dubal
,
N.
Wood
,
A.
Staniforth
, and
M.
Zerroukat
,
2010
:
An inherently mass-conserving iterative semi-implicit semi-Lagrangian discretization of the non-hydrostatic vertical-slice equations
.
Quart. J. Roy. Meteor. Soc.
,
136
,
799
814
, doi:.
Mesinger
,
F.
,
Z. I.
Janjić
,
S.
Ničkovic
,
D.
Gavrilov
, and
D. G.
Deaven
,
1988
:
The step-mountain coordinate: Model description and performance for cases of alpine lee cyclogenesis and for a case of an Appalachian redevelopment
.
Mon. Wea. Rev.
,
116
,
1493
1518
, doi:.
Mesinger
,
F.
, and Coauthors
,
2012
:
An upgraded version of the Eta model
.
Meteor. Atmos. Phys.
,
116
,
63
79
, doi:.
Schär
,
C.
,
D.
Leuenberger
,
O.
Fuhrer
,
D.
Lüthi
, and
C.
Girard
,
2002
:
A new terrain-following vertical coordinate formulation for atmospheric prediction models
.
Mon. Wea. Rev.
,
130
,
2459
2480
, doi:.
Simon
,
J. S.
,
K. A.
Lundquist
, and
F. K.
Chow
,
2012
: Application of the immersed boundary method to simulations of flow over steep, mountainous terrain. 15th Conf. on Mountain Meteorology, Steamboat, CO, Amer. Meteor. Soc., P45. [Available online at https://ams.confex.com/ams/15MountMet/webprogram/Paper210300.html.]
Steppeler
,
J.
,
H.-W.
Bitzer
,
M.
Minotte
, and
L.
Bonaventura
,
2002
:
Nonhydrostatic atmospheric modeling using a z-coordinate representation
.
Mon. Wea. Rev.
,
130
,
2143
2149
, doi:.
Steppeler
,
J.
,
R.
Hess
,
U.
Schättler
, and
L.
Bonaventura
,
2003
:
Review of numerical methods for nonhydrostatic weather prediction models
.
Meteor. Atmos. Phys.
,
82
,
287
301
, doi:.
Steppeler
,
J.
,
S.-H.
Park
, and
A.
Dobler
,
2013
:
Forecasts covering one month using a cut-cell model
.
Geosci. Model Dev.
,
6
,
875
882
, doi:.
Thuburn
,
J.
, and
C.
Cotter
,
2012
:
A framework for mimetic discretization of the rotating shallow-water equations on arbitrary polygonal grids
.
SIAM J. Sci. Comput.
,
34
,
B203
B225
, doi:.
Webster
,
S.
,
A.
Brown
,
D.
Cameron
, and
C.
Jones
,
2003
:
Improvements to the representation of orography in the Met Office Unified Model
.
Quart. J. Roy. Meteor. Soc.
,
129
,
1989
2010
, doi:.
Weller
,
H.
, and
A.
Shahrokhi
,
2014
:
Curl-free pressure gradients over orography in a solution of the fully compressible Euler equations with implicit treatment of acoustic and gravity waves
.
Mon. Wea. Rev.
,
142
,
4439
4457
, doi:.
Weller
,
H.
,
H. G.
Weller
, and
A.
Fournier
,
2009
:
Voronoi, Delaunay, and block-structured mesh refinement for solution of the shallow-water equations on the sphere
.
Mon. Wea. Rev.
,
137
,
4208
4224
, doi:.
Yamazaki
,
H.
, and
T.
Satomura
,
2010
:
Nonhydrostatic atmospheric modeling using a combined Cartesian grid
.
Mon. Wea. Rev.
,
138
,
3932
3945
, doi:.
Zängl
,
G.
,
2012
:
Extending the numerical stability limit of terrain-following coordinate models over steep slopes
.
Mon. Wea. Rev.
,
140
,
3722
3733
, doi:.
This article is licensed under a Creative Commons Attribution 4.0 license.