Abstract

Cut cells use regular or nearly regular polygonal cells to describe fields. For a given orography, some cells may be completely under the mountain, some completely above the mountain, and some are partially filled with air. While there are reports indicating considerably improved simulations with cut cells, inaccuracies may arise with some approximations, producing noise in fields near the surface. This behavior may depend strongly on the approximations made for the advection terms near the surface. This paper investigates the accuracy of advection for numerical schemes for a nondivergent flow near a mountain surface. The schemes use C-grid staggering with densities located at cell centers or on the corners of cells. Also, a nonconserving scheme is considered, which was used in the past with real-data cut-cell simulations. Since the cut cells near the surface create an irregular resolution, the accuracy and order of some approximations may break down near the surface. The objective of this paper is to find schemes having the same accuracy for advection near the surface as in the interior of the domain. As a test problem, uniform advection by a nondivergent velocity field is used with a 45° slope mountain (represented as a straight line) on a rectangular grid. Along the surface a sequence of triangular and pentagonal cells of quite different sizes are generated. Some schemes being discussed for cut cells lead to inaccurate and noisy solutions for this perfectly smooth mountain. A scheme using piecewise linear basis functions in a C grid with density points at the cell corners avoids these inaccuracies.

1. Introduction

Terrain-following (TF) coordinates provide a good approximation to mountain-induced flow if the approximation condition dh ≪ dz is met, where dh is the maximum change of orographic height between neighboring grid cells and dz is the vertical grid length. Real data models rarely meet this condition (Damrath et al. 2000), as forecast scores depend very much on having mountains in the model as realistic as possible (Wallace et al. 1983; Tibaldi 1986). This means that orographic filtering is kept to a minimum (Davies and Brown 2001). The resolution (after filtering) of the orography should roughly correspond to the effective resolution, as for example proposed by Skamarock (2004). The errors associated with TF coordinates have been known for a long time (Sundquist 1967). They tend to be larger with increasing resolution and can even lead to instabilities. The rather strong impact of various boundary representations was investigated by Webster et al. (2003). Within the TF coordinate approach, work has been done to reduce such errors (e.g., see the discussion in Shaw and Weller 2016). However, errors coming from TF coordinates when the condition dh ≪ dz is violated persist, and Steppeler et al. (2006) showed in a series of real data forecasts, that cut cells are able to improve forecasts considerably and produce good results without using any orographic filtering.

Adcroft et al. (1997) discussed a number of simplifications of the cut-cell approach, including the step approach. This approach was shown by Gallus and Klemp (2000) to be nonconvergent. This error results from the orography being represented as a piecewise constant function. When using cut cells with piecewise linear functions, Steppeler et al. (2002) showed that the solution converged and represented well the known linear solution for a shallow smooth mountain. A number of other errors are encountered, depending on approximations with the cut-cell approach. Shaw and Weller (2016) used horizontal cell boundaries being not quite horizontal and showed that the atmosphere at rest was not realized exactly. By comparing these results with Steppeler et al. (2006), it appears that these deviations from the atmosphere at rest produce velocities larger than the mountain/valley winds caused by radiation during day and night. Therefore, even small deviations from horizontal cell boundaries should be avoided with cut cells. The method outlined in section 4, using the corner grid, has horizontal lower, and upper, cell boundaries and therefore in a full application may be expected to realize the atmosphere at rest exactly. Because of the irregular resolution of cut cells near the mountain surface, errors and noise may appear in solutions for smooth mountains. The measures used to avoid such errors and noise are described in Steppeler et al. (2006) and Steppeler et al. (2011, 2013). The amount of noise error seems to depend critically on the approximation of the advection terms.

Good et al. (2014) investigated the advection errors near, but not immediately on the mountain surface and found much reduced error for the cut cells, compared to the TF coordinate. Shaw and Weller (2016) showed that a considerable reduction of error is achieved when the direction of advection coincides with coordinate lines. It may be difficult to find coordinates being suitable for all kinds of meteorological flow. However, the orographic surface is a line, which may be followed often by the flow. It, therefore, seems important to have a good accuracy of advection for this situation. This may also be at the heart of accuracy problems caused by the lower surface, as the advection term is also important when computing gravity waves.

In this paper we test the accuracy of 2D advection for flow at the orographic surface. A simple test is used, with a mountain slope of 45°, dx = dz, and a homogeneous wind field parallel to the terrain. The linear advection is conserving if div(U, W) = 0. The schemes are for a C grid with densities at cell centers or at corners. The methods proposed are conserving and derived as continuous spectral elements, describing fields within cells as polynomials. The C-grid schemes are low-order versions of these methods.

As we are interested in the effect of spatial discretization, all time integration is done by the highly accurate Runge–Kutta (RK4) method.

2. The L-Galerkin method

The approach to cut-cell modeling here is based on the local Galerkin method. This uses the representation of fields by piecewise polynomial and continuous functions and neglects approximations on discontinuous function spaces. Widely used local Galerkin methods using continuous or discontinuous basis functions are known under the names of continuous Galerkin (CG) or spectral element methods (see Taylor et al. 1997). To describe the C grid by basis functions, some discontinuities are admitted. However, the fluxes used will always be continuous splines. The CG method allows a high order of approximation in combination with conservation properties. For our simple advection problem we consider the conservation of mass. The following description includes the high-order field representations, even though in sections 3 and 4 only polynomials of degree 0 or 1 will be used.

A general outline of the L-Galerkin approach is given in the  appendix. The  appendix describes the L-Galerkin approach in combination with cut cells and the formalism presented there allows schemes other than the CG scheme. With the cut-cell method the dynamic amplitudes and field representation by polynomials even for cut cells is done within the associated uncut cells. This means that dynamic amplitudes under the mountain are possible, when the interpolation provided by this amplitude has an impact for the field above the mountain. The equations of motion are, however, approximated only above the mountain, which means, that the Galerkin integrals are taken only over the area of a cell above the mountain. For low-order basis function this will be explained in detail in sections 3 and 4. For the cut-cell method, amplitudes are given on the corners or centers of an uncut cell, even when the cell is cut. An alternative possibility is to just use small cells near the boundary. This is an L-Galerkin procedure on irregular cells and is not considered here.

As piecewise polynomials can always be described by a set of grid points, any L-Galerkin scheme can be considered as a finite-difference scheme and vice versa. If irregular cells occur and a finite-difference scheme is given for the regular cell, the L-Galerkin interpretation allows a natural generalization of this finite-difference scheme to the irregular cell. The latter property is used in this paper. Even for uncut cells many more schemes are possible than known under the name CG. The CG or spectral element scheme used by Taylor et al. (1997) is the only continuous L-Galerkin method widely applied and will be referred to as the standard L-Galerkin method. For fields represented by piecewise quadratic basis functions the  appendix defines two families of schemes under the names o2o2 and o2o3.

In sections 3 and 4 it will be shown how the standard L-Galerkin method using piecewise linear or constant basis functions is equivalent to centered differences for uncut cells and provides naturally the finite-difference equations for the cut cells. L-Galerkin cut-cell methods are local approximations. Thus, the matrix inversions involved in older finite-element models are avoided, which require nonlocal approximations that are problematic for use on massively parallel computers. The aim of such simple models is to work toward realistic models, recognizing that considerable future work would be required in order to generalize the scheme for use in the full model.

A realistic model application will probably use an existing forecast model with TF coordinates as a starting point, to provide features like physics treatment and data handling, which may be called the background model. Only a small percentage of the code in the background model needs to be altered to make the change to cut cells. Any of the L-Galerkin methods conceivable according to the  appendix will fit into any realistic finite-difference model. The only modeling concept not easily compatible with L-Galerkin cut cells is the global spectral method. Thus, it appears possible, in principle, to create a C-grid cut-cell scheme for a local finite-difference or spectral-element model and vice versa. However, for practical reasons it may be necessary to choose the cut-cell dynamics similar to that of the background TF model. So, for example, a C-grid background model should be developed into a C-grid cut-cell model.

The L-Galerkin concept as presented in the  appendix is flexible enough to emulate for uncut cells most of the spatial discretizations in practical use, such as centered differences using different kinds of staggering or spectral elements of a different order. An L-Galerkin method can be chosen to be identical to that of a given TF model on flat terrain.

Before describing specific numerical experiments, it is worth emphasizing that the stability constraints for irregular grids cannot be estimated by evaluating a Courant number based on the smallest occurring dx. As an example consider the 1D advection equation:

 
formula

Define the grid determining the cell intervals as xi+1/2 = xi−1/2 +dxi. The field ρi is defined at the interval center grid xi = (xi−1/2 + xi+1/2)/2. The scheme tested is defining the field ρi+1/2 at the interval ends by linear interpolation:

 
formula

Choose for even i and for odd i. Such grids occur by computing a plane wave on a hexagonal grid (with and ), where the grid points are defined on the corners of hexagonal cells. This example allows an analytic stability analysis. The finite-difference equations can be transformed into

 
formula

From this arithmetic transformation it follows that the scheme is stable for CFL < 1 with leapfrog time integration or CFL < 2.8 with RK4, when CFL is computed using . It follows also, that the scheme is second order on this grid. In particular, for with fixed , the time step does not have to go to 0 for stable integrations. Therefore, for a general cut-cell grid a correct stability analysis for irregular grids is necessary.

3. 2D advection in C grid, density at cell centers

This section investigates advection in the x direction and along the diagonal using schemes with densities at cell centers, as summarized in Table 1. We are interested in the accuracy of advection near orography. We use the C grid, with U defined to the right and left of density points and W at the top and bottom. There are two choices for cells with the C grid: cells with density defined at cell centers (cell-centered grid) and cells with density defined at the cell corners (corner grid). In both cases the velocities U and W are defined on cell boundaries. This section deals with the cell-centered definition of densities, which naturally allows the use of the finite-volume method for discretization in irregular cells. Discretization for the corner grid will be treated in section 4.

Table 1.

Error information for numerical experiments.

Error information for numerical experiments.
Error information for numerical experiments.

The 2D advection equation is conserving if the velocity field is divergence free, div(U, W) = 0:

 
formula

At lateral boundaries the field is put to 0. For these pure advection tests, a uniform 300 × 300 points grid is used with dx = dz = 1. A Gaussian density function is initially centered at x = z = 100 and transported by a constant velocity field either with U = 1, W = 0 (such that at t = 100 the advected field should be centered at x = 200, z = 100) or with U = W = 1 (such that at t = 100 the advected field should be centered at x = z = 200). The initial density field is

 
formula

with xi = i × dx, zk = k × dz, dx = dz = 1, and dxmid = 8 for all experiments.

For the finite-volume integration of (4), we need amplitudes of ρ at cell edges, xi+1/2 = (i+1/2)dx, zk+1/2 = (k+1/2)dz, which must be interpolated using values from neighboring cells. A simple assumption used by Steppeler et al. (2002), Yamazaki and Satomura (2008), Lock et al. (2012), and others is the cell-centered approach, where the locations of the scalar variables are at the centers of the uncut cell, such that

 
formula

The assumption that the position is at the cell centers in cut cells is not first-order accurate near surfaces. We consider advection over a distance of x = 100, meaning 100 time steps for dt = 1. Figure 1a shows the result for advection in the x direction (U = 1, W = 0) in the absence of any terrain.

Fig. 1.

2D advection in the absence of terrain for t = 0, t = 50, and t = 100. The arrows indicate the wind direction (a) velocity in x direction: (left) the initial field and results at (middle) t = 50 and (right) t = 100. (b) As in (a), but with the lower half of the initial field removed. (c) As in (a), but for advection along the diagonal. (d) As in (b), but for advection along the diagonal and the rotated initial condition with the lower half removed. The contour interval is 0.1.

Fig. 1.

2D advection in the absence of terrain for t = 0, t = 50, and t = 100. The arrows indicate the wind direction (a) velocity in x direction: (left) the initial field and results at (middle) t = 50 and (right) t = 100. (b) As in (a), but with the lower half of the initial field removed. (c) As in (a), but for advection along the diagonal. (d) As in (b), but for advection along the diagonal and the rotated initial condition with the lower half removed. The contour interval is 0.1.

As shown by Shaw and Weller (2016), for advection along coordinate lines, accurate solutions can be obtained for fields that are not smooth in the direction orthogonal to the coordinate lines. To demonstrate this, Fig. 1b shows the result for the initial values, for otherwise, where are the values used in (5). This confirms the notion of Shaw and Weller (2016) that for flow along coordinate lines strong gradients normal to the coordinate lines may be present, and the advection may still be accurate. For example, when the coordinate lines are horizontal, clouds in a model may have sharp vertical structures and still may be advected correctly in the horizontal. With TF coordinates over nonflat terrain this is not possible. Because of this effect, Good et al. (2014) found improved advection in higher layers of their cut-cell model, as compared to the TF coordinate.

When the flow does not follow coordinate lines, sharp gradients normal to the flow lines may lead to significant numerical error. Figures 1c and 1d correspond to Figs. 1a and 1b for a flow along the diagonal U = W = 1. The initial density field in Fig. 2d is rotated 45° from that shown in Fig. 1b. It is given by (5) for k > i or k = i and 0 for k < i. With the advection now oblique to the coordinate lines, the structure loses its sharp gradient and is elongated in the direction of the flow. Objective error information at t = 100 for the experiments shown in Fig. 1, named 1a–d, is given in Table 1 in terms of the maximum error , the minimum error , and the root-mean-square error , where represents the exact solution.

Fig. 2.

Orographically restricted grid. (a) δ = 0. (b) As in (a), but for δ = ½. The point A is a point (a) on the boundary or (b) outside the boundary, whose associated basis function reaches into the model area and therefore A carries a dynamic amplitude. The point B is a point inside the model area involving a cut cell. The finite-difference equation for this point is different from that of interior points. The point C is an interior point involving uncut cells only.

Fig. 2.

Orographically restricted grid. (a) δ = 0. (b) As in (a), but for δ = ½. The point A is a point (a) on the boundary or (b) outside the boundary, whose associated basis function reaches into the model area and therefore A carries a dynamic amplitude. The point B is a point inside the model area involving a cut cell. The finite-difference equation for this point is different from that of interior points. The point C is an interior point involving uncut cells only.

Exploiting the accuracy of flow along coordinate lines is typically difficult except for along horizontal coordinate lines. Near an orographic surface, the flow may be nearly parallel to the surface and it is of interest to investigate the transport along such surfaces, in particular for cut-cell grids. We choose a 45° sloped mountain surface and a homogeneous diagonal velocity field with U = W = 1, which also satisfies the zero normal velocity boundary condition. As this is a perfectly smooth mountain, it is important that a scheme performs well for this test. Because of the simplicity of the test, however, this may not be sufficient. Tests using a curved surface will need to be investigated as the next step. As cut cells create a rather irregular grid at the surface, not all schemes may be expected to perform well with this test. The aim of this investigation is to find schemes performing well for transport at the surface.

Figures 2a and 2b show two examples of the grid used. The parameter δ is used to shift the orographic surface by a distance δdx and allows the creation of triangular cells of arbitrarily small size. Figure 2b shows also that the boundary cells may consist of triangles and pentagons of rather different size. Even in the special case δ = 0 shown in Fig. 2a the area differences between the triangles and quadrilaterals are 1:2.

For the first tests we use density amplitudes with index i, k being located only at cell centers, even when the cells are cut. Amplitudes at cell boundaries have half indices i + 1/2, k for flux in x direction and i, k + 1/2 for the flux in z direction. Since we use constant velocities, staggered representation does not make a difference for the velocities. However, for the interpolated fluxes and in the x and z directions, we use basis functions that are piecewise constant normal to the direction of the flux and piecewise linear in the flux direction. These are the standard assumptions (e.g., see Steppeler 1989). For cell-centered schemes it is necessary to approximate densities at flux locations. In the present paper we consider equally weighted interpolations only:

 
formula

Using , (7) allows fluxes to be computed at cell edges. For cell edges that are cut by the terrain, we define or as the part of a cell edge being in the atmosphere, and we define weights ; . The time derivative is computed by the following finite-volume equation:

 
formula

where Fi,k is the part of the cell area being in the atmosphere above the surface.

The first experiment uses the “thin wall” approximation from Steppeler et al. (2002) together with the orography shown in Fig. 2a (orography going through grid points, δ = 0). The results reported in Steppeler et al. (2002) use the thin-wall approximation for all terms except advection. For advection the fix described below is used. Without this fix for the advection term noisy results were obtained. The thin-wall approximation consists of extending the triangles in Fig. 2a to quadrilaterals and using the interior areas of the uncut cells in the finite-volume approximation. The result is shown in Fig. 3a for an advection over 50dx and 100dx. As compared to free advection (Fig. 1c) there are substantial errors, which in a nonlinear setting could lead to noise not limited to the surface. These errors can be understood as the cut cells have an effective velocity reduced by ½ as compared to the cut cells. In initial experiments leading to Steppeler et al. (2002), the thin-wall approximation led to noisy gravity waves generated by a shallow mountain. Some details of this approach were only reported later, in Steppeler et al. (2006) and Steppeler et al. (2011, 2013). As this cut-cell model did not use the conserving equations, the thin-wall approximation was still used for the pressure gradient in the momentum equations and the divergence term in the pressure equation, but the advection terms were represented by the approximation described in the following.

Fig. 3.

Advection for C-grid schemes with density defined at cell centers: (a) thin-wall approximation, initial field, and after transport over 50dx and 100dx. (b) As in (a), but using the method of posing boundary values for advection tendencies. This result is independent of δ between 0 and 1. (c) As in (a), but for correct cell volume dx × dz/2 and flux interpolation based on equal distance cell centers. (d) As in (c), but for δ = ½ (the grid shown in Fig. 2b).

Fig. 3.

Advection for C-grid schemes with density defined at cell centers: (a) thin-wall approximation, initial field, and after transport over 50dx and 100dx. (b) As in (a), but using the method of posing boundary values for advection tendencies. This result is independent of δ between 0 and 1. (c) As in (a), but for correct cell volume dx × dz/2 and flux interpolation based on equal distance cell centers. (d) As in (c), but for δ = ½ (the grid shown in Fig. 2b).

To fix the noise problems with the thin walls, the method of posing boundary values for advection tendencies was used. This can be described as follows:

 
formula

For the special mountain used here, i = k indicates points on the mountain surface. In the following, k being a function of i will indicate a boundary point.

The result is shown in Fig. 3b. The solution is slightly broadened at the surface, but the result may be acceptable since the correction is very easy to apply in a 3D model. The method does not depend on details of the cut cells between the grids shown in Fig. 2a or 2b. The results are identical and not shown for δ = ½.

In the present test of a straight mountain, the posing of boundary values is unique, as the atmosphere is only on one side of the mountain. In cases more general than treated here, it should be noted that the value of ρ inside the mountain may be double valued. Consider as an example a very thin mountain positioned at i = i0. In this case, which is not further investigated here, (9) has to be replaced by

 
formula

In a 3D application, such as Steppeler et al. (2006) up to four values may appear at thin mountain peaks. As these values are diagnostic, there is no long-term storage necessary for this and, therefore, no complex programming is necessary.

When using the cells with the correct volume, rather than the thin-wall approximation, the interpolation of the fluxes at the cell boundaries being used in (2) must be defined for cells cut by the surface, which in our case are those with k = i. A possible approximation is that the densities are defined at the centers of the uncut cells, as was done by Steppeler et al. (2002), Lock et al. (2012), and J. Klemp (2015, unpublished manuscript). All tests presented here are based on this approximation. This means that (7) is used also for k = i. With this approximation (equal distance cell centers) for the fluxes and the correct cell volumes dx × dz/2, the results for the advection are given in Fig. 3c, again for the grid shown in Fig. 2a, with orography going through the density points. For this grid the errors appear small.

Again with the equal distance cell centers, the advection is repeated with δ = ½ for the grid shown in Fig. 2b. The result is shown in Fig. 3d and exhibits large errors, so this scheme cannot be considered suitable. The large errors are caused by a problem of consistency. Similar results were obtained after halving the time step. The error information for the results shown in Fig. 3 is given in Table 1 for the experiments 3a–d.

It is unknown if other assumptions concerning the position of density amplitudes, together with different interpolation procedures to the cell surfaces, could produce better results. In this paper the improvement of accuracy will be achieved by the transition to the grid where the densities are piecewise linear splines.

4. Advection in C grid with density at cell corners

In this section we consider schemes where the density is defined at cell corners. As our test problem uses constant U and W, the special C-grid basis functions for velocities need not be used. Here U and W are given at any point of the model area. The density as a function of x and z is given as bilinear interpolation in each quadrilateral cell. As seen from Fig. 2, triangular cells may also occur, for which triangular interpolation is used. This is the first-order interpolation defined by the corner amplitudes. For the pentagons in Fig. 2b we interpolate within quadrilaterals using bilinear functions defined by the four corner amplitudes. Note that for pentagons, amplitudes and interpolations are always defined for the corresponding uncut cells.

As seen from Figs. 2a and 2b, in our particular example only triangular cells with the right angle on top will occur. Figure 4 indicates the numbering of corners for a triangular and a quadrilateral cell. This numbering avoids double indexing in the equations. For example in Fig. 2a, a triangular cell at the lower boundary (defined by indices i, i) is defined by its corners (i, k), (i, k + 1), and (i + 1, k + 1) with k = i. An example is indicated in Fig. 2 by the letter A. As shown in Fig. 4, these corners can be referred to as 1, 2, and 3. The index of a cell is defined as that of the lower-left corner (index 1 in Fig. 4). The above-mentioned cell therefore has index j = (i, k) with k = i. It should be kept in mind that the associated uncut cells are used for carrying the amplitudes. The smaller areas of the cut cells are used only when computing the Galerkin integrals.

Fig. 4.

The numbering of corners belonging to (left) square and (right) triangular cells.

Fig. 4.

The numbering of corners belonging to (left) square and (right) triangular cells.

To perform the discretization it is necessary to know the corners for each cell and the neighboring cells that share each corner. Because of the simplicity of our example, these neighborhood relations are directly apparent in Fig. 2 and a formalism for unstructured grids is not necessary. As seen in Fig. 2, there are four neighboring cells for interior grid points and three for boundary points. For example the boundary point with index (i, i), indicated by A in Fig. 2a has three neighboring cells, two triangles, and one quadrilateral. The indices representing cell corners are (i − 1, i − 1), (i − 1, i), and (i, i), with the index numbering as defined in Fig. 4.

Consider as an example a triangular cell i, k. We indicate the cell corner indices by li such that l1 = (i, k), l2 = (i, k + 1), and l3 = (i + 1, k + 1). For each corner of each cell a basis function fj,l (x, z) is defined, being different from 0 only in the cell j. Inside quadrilateral cells fj,l is the bilinear function being 1 at corner l and 0 at all other corners. The bilinear polynomial is of the form p = p0 + p1x + p2z + p3xz. For triangular cells linear functions in x, z are fitted. The basis functions for the triangular cells Cj, j = (i, i) (see Figs. 2a and 4) are defined as

 
formula

The fj,l are 1 at the selected vertex and 0 at the others. Outside the selected triangle j, the fj,l are defined to be 0. We define the characteristic function of cell Cj as follows:

 
formula

A discontinuous function is given by a set of amplitudes Aj,l, j = (i, k):

 
formula

Another representation, obtained by rearranging the sum in (13), is obtained by considering the index j = i, k to be the index of a point. Let be the amplitudes belonging to the point j = i, k. Here is defined by = Aj, where j′ is one of the four or three neighboring cells of point j = i, k and i′ is the corner of cell j′ being identical to j. The numbering of the corners is indicated in Fig. 4. For a quadrilateral cell j = (i, k) in the interior of the model area, the neighboring cells are, numbered clockwise: (i − 1, k − 1), (i − 1, k), (i, k), and (i, k − 1). Combining this with Fig. 4, we obtain the following:

 
formula

The basis functions can be renumbered in a similar way:

 
formula

Using (13), (14), and (15), we obtain the alternative representation of (13), which is just a rearrangement of the sum in (13):

 
formula

Here the superscript d is used to indicate that we are dealing with discontinuous functions. Continuous functions are obtained, when the A in the second sum does not depend on l:

 
formula

It is also possible to define a function for any point j = i, k as follows:

 
formula

From (18), we obtain the representation for a continuous function:

 
formula

For points j = (i, k) with triangular and quadrilateral neighboring cells, (13) to (19) can be written in an analogous way.

The application of the right-hand side of (4) leads to discontinuous functions, represented according to (16). To obtain a time derivative to be used with the Runge–Kutta method, must be approximated by . This can be done in a variety of ways. Here we use the standard version of the spectral element scheme (Taylor et al. 1997). The derivation of the method in terms of spectral elements results from the observation that the Gauss–Lobatto points for linear functions are just the end points of the interval in question. We obtain the amplitude as a weighted average of the discontinuous amplitudes , using weights :

 
formula

The values of the weights follow from the requirement of mass conservation. We require that, for a field where A is different from 0 at only one point j, the mass of the discontinuous function is identical to the mass of the approximation . The weights then become

 
formula

For dx = dz = 1 the weights for the quadrilaterals and triangles are as follows:

 
formula
 
formula

For quadrilateral and triangular cells the derivative in x or z at the center of a boundary cell is given by a centered difference with the target point at the middle of the boundary cell:

 
formula

The weights to obtain x derivatives at corner points from these directional derivatives are given in (19) and (20).

Let us first consider the grid in Fig. 2a and define the fluxes . For interior points j = i, k, surrounded by quadrilaterals only:

 
formula

which is the centered difference scheme for uncut cells.

For grid points (i − 1, i) in the second diagonal from the row of surface points (i, i) we obtain the following:

 
formula

For points (i, i) on the boundary the difference equations are as follows:

 
formula

The coefficients are computed according to (21)(23). Note that the derivative in each triangular cell is constant, with (12) giving the amplitudes to be applied in each corner. For quadrilateral cells the derivative is linear.

The equations for the fluxes are analogous to (25)(27). The results using these equations for advection on the grid in Fig. 2a (orography going through corners) are shown in Fig. 5a. These results are similar in quality to those shown in Fig. 3b.

Fig. 5.

Advection in the orographically restricted grid for the C grid defining density at the corners of cells: (a) using the grid shown in Fig. 2a and (b) using the grid shown in Fig. 2b.

Fig. 5.

Advection in the orographically restricted grid for the C grid defining density at the corners of cells: (a) using the grid shown in Fig. 2a and (b) using the grid shown in Fig. 2b.

An alternative choice for the basis function for quadrilaterals can be obtained by dividing each quadrilateral into two triangles, where here we choose the dividing line to go from i, k to i + 1, k + 1. This will lead to six neighboring points. The only differences are that in (23) the integral is ⅞ for corners 1 and 3 and ⅙ for corners 2 and 4. In this case, for example, the finite-difference equation in (25) changes to

 
formula

which again recovers the centered difference approximation for uncut cells. For the cases involving cut cells, the weights given in (26) and (27) are slightly changed. The results for this 2D advection test (not shown) are very similar to those shown in Fig. 5a and are of the same quality. For the purpose of this paper the division of quadrilaterals into triangles has the advantage that the analytic computation of the weights is easier than using bilinear basis functions.

There is an arithmetic transformation of the difference equations given above, which allows one to see immediately that the combinations of the one-sided differences given above are stable. The equations are averages of centered differences along the flow direction. For example for the triangular cell (i, i) on the surface, whose lower-left corner is indicated by A in Fig. 2a, the divergence or time derivative of density is given by

 
formula

With the definition , (26) for j = i − 1, i can be rewritten as

 
formula

and (27) for j = (i, i) becomes

 
formula

For the grid shown in Fig. 2b, small triangles and pentagons appear as cells in addition to quadrilaterals. For the purposes of computing the Galerkin integrals in order to compute time derivatives, we subdivide each pentagon into three quadrilaterals and one triangle, such that we consider the same cell structure is used as in the grid shown in Fig. 2a. Dynamic amplitudes are defined only for the uncut cells of size dx × dz. For the nodes of the small areas it is possible to compute amplitudes diagnostically using the interpolation in the associated large uncut cell. These diagnostic amplitudes can be used to compute the integrals analytically over the whole cut cell.

For purposes of interpolation and defining basis functions, the cut cells are extended to their uncut area. No division into smaller cells is used for this purpose. The interpolation by basis functions is the same as for the grid in Fig. 2a. The amplitudes are defined in the uncut cells, meaning, that the divergence inside each cell is obtained by the same formula as for the grid in Fig. 2a. However when the weights in (20) are computed by the Galerkin integrals in (22) and (23), the integration extends only over the cut cells, not the corresponding uncut cells, resulting in different weights in (30) and (31). For the grid in Fig. 2b, the finite-difference equations corresponding to (30) and (31) become

 
formula
 
formula

Figure 5b shows the advection results using the grid in Fig. 2b. Comparing with Fig. 5a, the accuracy is not degraded by the transition from the grid shown in Fig. 2a to the grid shown in Fig. 2b. Objective error information for the results shown in Fig. 5 is given in Table 1 for experiments 5a and 5b.

So far only an analytic computation of Galerkin integrals was used, which is possible for the special cases shown in Fig. 2. General Galerkin integrals may have to be computed numerically.

Equations (30)(33) and the results in Fig. 5 indicate in a heuristic way that the time step for stability does not go to 0 for small cut cells. Another case involving rather small triangles occurs with δ = 0.9. The resulting dynamic equations are as follows:

 
formula
 
formula

The plotted advection results (not shown) look very similar to those of Fig. 5b. As a stability analysis for the irregular cut-cell resolution is beyond the scope of this paper, we can only give heuristic arguments for the stability of this method. In this simulation experiment, (34) and (35) seem to provide a stable integration with dt = 1 in spite of the small triangular cells involved. This is consistent with the stability analysis given in section 2 for a 1D case with alternating large and small grid intervals.

The limiting case of vanishing 1 − δ is obtained by considering 1 − δ to be a very small number and neglecting all terms proportional to 1 − δ:

 
formula
 
formula

The time derivative amplitude in (37) defines the field for an area proportional to 1 − δ, which is approaching 0. Integrations could be done with dt = 1 and the results (not shown) are comparable to those in Fig. 5b.

To achieve second-order differencing for irregular cells, it is possible to define the fluxes as piecewise second-order polynomials. Because of the simplicity of our test case, (32) and (33) show that the piecewise linear fluxes already lead to second-order differences.

5. Conclusions

Using a very simple advection test along a straight-line 45° mountain, a number of numerical schemes used on cut-cell grids exhibit large errors and generate noise for advection along the surface. Two schemes investigated here appear to be free of these inaccuracies. One of these is the method of posing boundary values for the advection tendencies. This method can be used with the thin-wall approximation for the pressure gradient. It is nonconserving. The other method to avoid inaccuracies and noise at the lower boundary uses a C-grid representation with the densities defined at the corners of cells. The finite differences for irregular cut cells derive from the spectral element approximation being applied to piecewise linear basis functions.

Acknowledgments

Funding for this research was provided by the National Center for Atmospheric Research through support from the National Science Foundation under Cooperative Agreement AGS-0856145.

APPENDIX

General Outline of the L-Galerkin Method with Cut Cells

The L-Galerkin procedure is given in its general form for the case of cut cells. We use the name L-Galerkin scheme for any numerical methods based on field representations by piecewise polynomial representations on polygonal grids. The standard CG method used by Taylor et al. (1997) or the third-order method presented by Steppeler (1976) are examples. This general outline of the 2D case will include the cut cells and it becomes clear that a great variety of approaches even for uncut cells is possible, which for the most part are unexplored.

In sections 2 and 3 examples of low-order basis functions, leading to a C-grid approximation are given and for example in Taylor et al. (1997) high-order schemes are described using only uncut cells. The formalism described in the following admits many more possibilities than have been investigated so far. For simplicity of notation the general outline of the 2D L-Galerkin approach is given for an unstructured grid, but our example uses the special case of a structured grid. With cut cells there is an orography profile, which may be any smooth curve. In our example this is a piecewise linear spline O(x) defined on the horizontal coordinate. Only points above O(x) are considered to be in the atmosphere. Most horizontal discretization schemes in practical use can be considered as L-Galerkin schemes. A short general outline of the L-Galerkin method will be given for the 2D case.

Let Cj be a division of the model area M into polygonal cells: For (x, z) points in the model area, approximation basis functions bj,l (x, z) are defined by bj,l (x, z) = 0 for , and bj,l (x, z) = pl (x, z) for , with pl (x, z) being a polynomial of the form pl (x, z) = . For a field f(x, z) we have the piecewise polynomial approximation f(x, z) = . While CG methods have been applied for rather high orders, such as , we consider low-order schemes with being 0, 1, 2, and 3 to be sufficiently accurate. The amplitudes Aj,l are the dynamic variables of the L-Galerkin scheme. For an L-Galerkin scheme there is always an equivalent grid point representation by choosing (xi,k, zi,k) for each cell collocation point. For a field f(x, z), the gridpoint values are fi,k = . When choosing identical ranges for the indices k and l, this equation can be inverted and, therefore, the gridpoint representation of a field is equivalent to the polynomial representation.

A popular example is the regular square spectral element grid. For the order two representation this grid is given by equidistant grid (xi, zk), with xi = idx, zk = kdz. The cell for order two is defined by using even i and k as corners. Each cell contains 3 × 3 = 9 points, where some involving an even index are shared with neighboring cells. The polynomial basis for the full grid is p0 = 1, p1 = x, p2 = x2, p3 = z, p4 = xz, p5 = xz2, p6 = z2, p7 = x2z, and p8 = x2z2. The equidistant structured collocation grid defined above has nine points per cell and the gridpoint values give an equivalent representation as compared to the polynomial amplitudes Aj,l. Note that index j = (i′, k′) points to a cell, l points to a polynomial, and (i, k) points to a collocation grid point. If we define the index of a cell to be that of the lower-left corner, cell indices for our example are given by even i′ and j′.

Corner points (i and k even) are shared by four cells and edge points (i even and k odd or vice versa) are shared by two cells. By choosing the cell or corner amplitudes independently in each cell, discontinuous functions are defined. A continuous field is obtained by allowing only one amplitude at the corner and edge points that are shared by several cells. For unstructured mesh programming, the administration of cells and their interior, edge, and corner points is one of the technical challenges. With a structured grid, as used in the example of a second-order field representation, the administration of corner and edge points is much easier.

The fluxes F(x, z) may be approximated in the same function space as the fields f(x, z) and this is required in the standard L-Galerkin (CG) scheme of Taylor et al. (1997). Forming the derivative of f(x, z) will create discontinuities at cell corners and edges. The standard L-Galerkin operator approximates this in the space of continuous functions. If o2 is the space of piecewise second-order polynomials defined above, we may call this an o2o2 method, as both the flux F(x, z) and the field f(x, z) are approximated in the o2 space. As an alternative, it is also possible to approximate the flux in the higher order o3 space. It is possible to define F(x, z) such that it is differentiable at corner points. The flux derivative is then in the same space as the original field and can, therefore, be used in time marching schemes without further approximation.

Not all schemes possible within the definitions above are stable. The standard L-Galerkin procedure can be unstable when the collocation points are not chosen as Gauss–Lobatto points. However, in addition to the standard L-Galerkin procedure, we found a number of stable schemes of o3o3 and o2o3 type, as well as the low-order approach described in this paper. While a systematic description of many L-Galerkin schemes is beyond the scope of this paper, a few technicalities are mentioned, which have potential for increasing numerical efficiency and accuracy.

As gridpoint representations and function space representations by amplitudes are equivalent, all known operators, such as spatial averaging can immediately be applied in function space. With L-Galerkin methods the collocation operator is often used. If the formulation of the flux as a function of the fields involves a product, this product is taken only at certain collocation points and within a cell an interpolation using the function space, such as the o2 space, is done. Applying such operators in the same way as with Arakawa finite-difference schemes makes a large number of cut-cell L-Galerkin schemes possible, which could be investigated for their efficiency and accuracy.

A large saving of computational effort is possible, when sparse grids in combination with serendipity function spaces are used. The above example of the o2 space can use a sparse grid, when using only corner and edge points in the gridpoint representation. For this example this means that the grid points j = (i, k) with odd i and k are not used. In the polynomial representation the basis function x2z2 is not used. This again leads to a one-to-one correspondence between polynomial amplitudes and the sparse grid. With finite elements serendipity sparse grids are standard (e.g., see Brenner and Scott 2002). For 3D and order three, the savings of computer time obtained by using sparse grids are substantial. Current applications of CG methods do not use sparse grids, indicating a large potential for increasing the efficiency for the higher-order methods. Basis functions of order up to 1, as used in this paper do not allow sparse grids.

REFERENCES

REFERENCES
Adcroft
,
A.
,
C.
Hill
, and
J.
Marshall
,
1997
:
Representation of topography by cells in a height coordinate model
.
Mon. Wea. Rev.
,
125
,
2293
2315
, doi:.
Brenner
,
S. C.
, and
L. R.
Scott
,
2002
: The Mathematical Theory of Finite Element Methods. 2nd ed. Springer, 392 pp.
Damrath
,
U.
,
G.
Doms
,
D.
Frühwald
,
E.
Heise
, and
B.
Richter
,
2000
:
Operational quantitative precipitation forecasting at the German Weather Service
.
J. Hydrol.
,
239
,
260
285
, doi:.
Davies
,
L. A.
, and
A. R.
Brown
,
2001
:
Assessment of which scales of orography can be credibly resolved in a model
.
Quart. J. Roy. Meteor. Soc.
,
127
,
1225
1236
, doi:.
Gallus
,
W. A.
, and
J. B.
Klemp
,
2000
:
Behavior of flow over step orography
.
Mon. Wea. Rev.
,
128
,
1153
1164
, doi:.
Good
,
B.
,
A.
Gadian
,
S.-J.
Lock
, and
A. N.
Ross
,
2014
:
Performance of representing orography in idealized simulations
.
Atmos. Sci. Lett.
,
15
,
44
49
, doi:.
Lock
,
S.-J.
,
H.
Bitzer
,
A.
Coals
,
A.
Gadian
, and
S.
Mobbs
,
2012
:
Demonstration of a cut cell representation of 3D orography for studies of atmospheric flows over very steep hills
.
Mon. Wea. Rev.
,
140
,
411
424
, doi:.
Shaw
,
J.
, and
H.
Weller
,
2016
:
Comparison of terrain following and cut cell grids using a nonhydrostatic model
.
Mon. Wea. Rev.
,
144
,
2085
2099
, doi:.
Skamarock
,
W. C.
,
2004
:
Evaluation of mesoscale NWP models using kinetic energy spectra
.
Mon. Wea. Rev.
,
132
,
3019
3032
, doi:.
Steppeler
,
J.
,
1976
:
The application of the second and third degree methods
.
J. Comput. Phys.
,
22
,
295
318
, doi:.
Steppeler
,
J.
,
1989
:
Analysis of group velocities of various finite element schemes
.
Beitr. Phys. Atmos.
,
62
,
151
161
.
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.
,
H. W.
Bitzer
,
Z.
Janjic
,
U.
Schättler
,
P.
Prohl
,
U.
Gjertsen
,
J.
Parfinievicz
, and
U.
Damrath
,
2006
:
Prediction of clouds and rain using a z-coordinate nonhydrostatic model
.
Mon. Wea. Rev.
,
134
,
3625
3643
, doi:.
Steppeler
,
J.
,
S.-H.
Park
, and
A.
Dobler
,
2011
:
A five day hindcast experiment using a cut cell z-coordinate model
.
Atmos. Sci. Lett.
,
12
,
340
344
, 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:.
Sundquist
,
H.
,
1967
:
On vertical interpolation and truncation in connection with the use of sigma system models
.
Atmosphere
,
14
,
37
52
.
Taylor
,
M.
,
J.
Tribbia
, and
M.
Iskandarani
,
1997
:
The spectral element method for the shallow water equations on the sphere
.
J. Comput. Phys.
,
130
,
92
108
, doi:.
Tibaldi
,
S.
,
1986
: Envelope orography and maintenance of quasi-stationary waves in the ECMWF model. Advances in Geophysics, Vol. 29, Academic Press, 339–374, doi:.
Wallace
,
J. M.
,
S.
Tibaldi
, and
A. J.
Simmons
,
1983
:
Reduction of systematic errors on the ECMWF model through the introduction of an envelope orography
.
Quart. J. Roy. Meteor. Soc.
,
109
,
683
717
, doi:.
Webster
,
S.
,
A. R.
Brown
,
D. R.
Cameron
, and
C. P.
Jones
,
2003
:
Improvements to the representation of Orography in the Met Office Unified Model
.
Quart. J. Roy. Meteor. Soc.
,
129
,
1989
2010
, doi:.
Yamazaki
,
H.
, and
T.
Satomura
,
2008
:
Vertically combined shaved cell model in a z-coordinate non-hydrostatic model
.
Atmos. Sci. Lett.
,
9
,
171
175
, doi:.

Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).