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.
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.
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
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.
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
where the vertical decay functions are given by
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.
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.
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 x–z 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:
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
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
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
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.
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
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
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:
A tracer with density ϕ is positioned upstream above the height of the terrain. It has the following shape:
with radius r given by
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
where ϕ is the numerical tracer value, is the analytic value, and is the cell volume.
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
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
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:
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.
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
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
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.
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
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
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.
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.
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.
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:
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
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:
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.
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.
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.
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.
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:
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:
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:
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:
where subscript F denotes midpoint interpolation from two surrounding cell values onto face f:
where denotes the two cells sharing face f. Here is the consistent cell center reconstruction of from surrounding values of :
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:
The divergence is discretized using Gauss’s divergence theorem so that
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:
Equation (A6) can now be approximated by
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:
where the diagonal matrix and the nonorthogonal correction is . The orthogonal part, , can be treated implicitly in the Helmholtz equation:
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.