## Abstract

Cut-cell meshes present an attractive alternative to terrain-following coordinates for the representation of topography within atmospheric flow simulations, particularly in regions of steep topographic gradients. In this paper, we present an explicit two-dimensional method for the numerical solution on such meshes of atmospheric flow equations including gravitational sources. This method is fully conservative and allows for time steps determined by the regular grid spacing, avoiding potential stability issues due to arbitrarily small boundary cells. We believe that the scheme is unique in that it is developed within a dimensionally split framework, in which each coordinate direction in the flow is solved independently at each time step. Other notable features of the scheme are: (i) its conceptual and practical simplicity, (ii) its flexibility with regard to the one-dimensional flux approximation scheme employed, and (iii) the well-balancing of the gravitational sources allowing for stable simulation of near-hydrostatic flows. The presented method is applied to a selection of test problems including buoyant bubble rise interacting with geometry and lee-wave generation due to topography.

## 1. Introduction

For many years, the method of choice for the representation of topography in atmospheric simulations has been the terrain-following coordinate approach (e.g. Philips 1957; Gal-Chen & Somerville 1975; Simmons & Burridge 1981; Zhu *et al.* 1992; Webster *et al*. 1999). However, improvements in computing power and algorithms in recent years have made ever-increasing horizontal resolutions accessible within such simulations. At these high resolutions, relatively steep gradients in orography are evident and, as a result, the use of the terrain-following approach typically introduces significant truncation errors in these regions, as discussed in Sundqvist (1976) and Janjic (1989).

Regular Cartesian meshes present many advantages for numerical simulation due to their underlying structure, including efficiency of both storage and computation, and trivial partitioning for parallel processing. Such meshes are also attractive for their simple compatibility with adaptive mesh refinement and dimensional splitting techniques. These grids are, however, not naturally suitable for the representation of geometrically complex boundary geometries. A first approximation of the effects of irregular boundaries within regular grid solvers may be achieved by the use of *step boundaries* (see Egger (1974) and Mesinger *et al.* (1988) for applications in the context of atmospheric flows). Here the boundary is approximated by the nearest regular grid cell edges, hence reducing the problem to the solution purely on regular grid cells. However, these methods result in a particularly coarse representation of the true geometry and, for example, have been demonstrated to predict incorrectly the generation of lee waves in the wake of mountains in Gallus & Klemp (2000).

*Cut-cell* techniques (also referred to as *Cartesian-grid methods*, *embedded boundaries* or *shaved cells*) offer a more accurate method for representing this irregular boundary. Here, a regular grid is intersected with the irregular geometry of the physical simulation domain, resulting in a collection of control volumes that are entirely regular away from the domain boundaries, while accurately representing the boundary through a layer of irregular ‘cut cells’ adjacent to the boundary (figure 1). The grid generation procedure is particularly simple and automated, allowing for efficient application on arbitrary geometries in multiple dimensions. First considered for the solution of the potential flow equations in Purvis & Burkhalter (1979), the approach has been successfully applied to varied and complex geometries for a range of different flow types (e.g. Quirk 1994; Pember *et al.* 1995; Almgren *et al.* 1997; Aftosmis *et al.* 1998; Causon *et al.* 2000; Colella *et al.* 2006), though it has only relatively recently been incorporated into atmospheric simulations (Steppeler *et al.* 2002; Rosatti *et al.* 2005; Gatti-Bono & Colella 2006; Yamazaki & Satomura 2008). The methodology is also attractive for the related field of oceanographic simulation as demonstrated in Adcroft *et al.* (1997) and Barad (2006) where even more significant gradients in bathymetry are evident.

The primary disadvantage associated with the cut-cell grid is the generation of arbitrarily small and highly distorted volumes adjacent to the boundaries of the domain, frequently referred to as the *small-cell problem*. Such cells may lead to stability and efficiency issues for numerical methods and therefore require special attention. Various numerical approaches have been introduced to resolve this issue for compressible flow simulation, including the merging of these small irregular cells with neighbouring regular cells in, for example, Clarke *et al.* (1986), Quirk (1994) and Yamazaki & Satomura (2008), the use of implicit methods stable for large time steps (LeVeque & Shyue 1996; Murman *et al.* 2003), or extending the influence of the cut cells and redistributing the resultant conservation error as in Pember *et al.* (1995). A novel approach based on constructing boundary-aligned boxes of equivalent length to regular cells, combined with rotated grid methods, has also been proposed in Helzel *et al.* (2005).

We present here an alternative numerical scheme for dealing with this small-cell problem, which we believe to be relatively simple both conceptually and in implementation. The method is developed within a dimensionally split framework, allowing for its coupling with any of a wide range of existing dimensionally split regular-grid explicit solvers. At its heart lies a one-dimensional flux stabilization constructed to ensure that the resultant update is stable for a time step independent of cell volume, combined with an auxiliary ‘splitting’ flux determined by considering a boundary-tangential flow problem, the effects of which cancel out over the splitting cycle. The resultant scheme is fully conservative and first-order accurate near boundaries. Here we apply the approach to the solution of fully compressible atmospheric flow problems and combine it with the well-balancing strategy of Botta *et al.* (2004) to allow for the accurate simulation of near-hydrostatic flows. Results are presented for (i) static and dynamic flows for a sample ‘zeppelin’ geometry embedded into the atmosphere and (ii) the oft-employed examples of lee-wave generation of Durran (1986).

## 2. Mesh generation

In two dimensions a cut-cell grid is generated by intersecting a regular grid of *N*_{i}×*N*_{j} regular control volumes of size Δ*x*×Δ*y* with the physical domain boundary (figure 1). Here we assume that this boundary may be represented by a level set specified on the nodes of the regular grid and construct a piecewise linear representation of the interface within each regular cell. For the method described here, the required geometric details for the irregular cells are as follows:

(i)

*α*_{i,j}is the volume of the irregular cell relative to that of a regular cell,(ii)

*β*_{i±1/2,j}and*β*_{i,j±1/2}are the area of the cell interface in each coordinate direction relative to a regular cell interface,(iii) is the cell centre of mass,

(iv) is the centre of the boundary interface within the cell, and

(v) is the normal to the boundary interface within the cell.

These quantities are all computed simply from the level-set representation of the topography, and the details are omitted here for the sake of brevity.

## 3. Model description

We assume a dry non-hydrostatic atmosphere with no rotation and consider the fully compressible system of Euler equations with gravity as a model for flow scales of 10–100 km horizontally and 5–50 km vertically. Those are the characteristic scales of flows considered in this paper. While this fully compressible model is computationally inefficient for long-term simulation of atmospheric flows because the time step is limited by the sound speed, it provides a valid starting point for the assessment of the cut-cell approach. Hence, the system we wish to solve is
3.1
3.2
3.3where the flow variables *ρ*, **V**=(*u*,*v*), *p* and *θ* are density, velocity, pressure and potential temperature, respectively, and *Φ* is a time-independent gravitational field, with in the present study. An ideal-gas equation of state is assumed, giving the pressure
3.4where *γ*=1.4 is the ratio of specific heats for air.

## 4. Numerical method

To solve the system of equations (3.1)–(3.3) numerically on the chosen finite-volume discretization, we wish to employ dimensional splitting to obtain a series of conservative one-dimensional problems. Therefore, we consider first the numerical solution of the one-dimensional problem in §4*a*. The extension of this approach to multiple dimensions within the dimensional splitting is then introduced in §4*b*.

In the context of atmospheric flows, care is also needed in computing the pressure divergence and gravitational source terms in equation (3.2) in order to ensure that the important hydrostatic balance is maintained. We therefore implement the well-balancing approach of Botta *et al.* (2004) within the cut-cell mesh framework, as detailed in §4*c*.

### (a) One-dimensional method

#### (i) Regular cells

We begin our discussion by describing the method employed here for the one-dimensional problem on regular cells. We denote by the vector of conserved state variables at time *t*^{n} for the cell *V*_{i}. The standard explicit conservative update for the numerical solution in the cell is then given by
4.1where Δ*t* is the chosen time step, *t*^{n+1}=*t*^{n}+Δ*t* and **f**_{i±1/2} are numerical approximations of the flux at the left and right boundaries of the cell, fixed over this time step. This explicit update is then stable for a time step Δ*t*≤Δ*x*/*λ*_{i}, where *λ*_{i} is the speed of the fastest wave in the cell (equivalent to the largest eigenvalue of the Jacobian of the flux function).

Various methods are available for the approximation of the fluxes **f**_{i±1/2}. Here, the state differences are decomposed in terms of characteristic variables *q*_{k} as , where **r**_{k} are the right eigenvectors of the flux Jacobian. The local double logarithmic reconstruction (LDLR) of Artebrant & Schroll (2006) is then applied to produce a spatial reconstruction of these characteristic variables across each cell. From this reconstruction, half-time-step interface states are computed for each cell. The left and right states from neighbouring cells then define a Riemann problem at each cell interface and the Harten–Lax–van Leer–Einfeldt–Munz (HLLEM) approximate Riemann solver of Einfeldt (1988) is employed to approximate the solution to this problem and obtain the required interface flux. The LDLR reconstruction is chosen here because of its superior performance for low Mach number flows, though it may be substituted by other schemes such as the monotone upstream-centred scheme for conservation laws (MUSCL) reconstruction.

#### (ii) Cut-cell flux stabilization

As illustrated in figure 2, the same approach may be considered for cells that are cut by the domain boundary. Consider a cut cell *V*_{i} of length fraction *α*_{i} positioned adjacent to regular cell *V*_{i−1}. Then the conservative update for the cut cell is given by
4.2where **f**_{B} is the boundary flux determined by the relevant boundary condition and is the effective flux between the cut cell and neighbouring regular cell.

By extending the ‘influence’ of the cut cell to the entire regular cell length, we obtain an update that is stable for the same time step as that for the regular cells,
4.3where **f**_{i−1/2} is computed in the same manner as for the regular cells, by obtaining left and right interface states and from reconstructions on the regular and cut cell, respectively, and solving approximately the resulting Riemann problem. This reconstruction is limited, so no overshoots over the adjacent cell data can occur, and the numerical flux **f**_{i−1/2} remains well defined and bounded as long as the neighbouring cell averages remain controlled.

To ensure that conservation is maintained, this update must match the conservative update of equation (4.2). Equating the two yields an expression for the stabilized flux , 4.4

By applying this stabilized flux to the update for both the cut cell and its neighbouring regular cell, the scheme as a whole is hence both conservative and stable for the time step determined by the regular cell length.

We note that this stabilized flux also, importantly, respects the natural limits as the cut cell approaches a full or empty cell, i.e. when , , and when , .

### (b) Multi-dimensional extension

#### (i) Interface shielding

By applying dimensional splitting, we may reduce the two-dimensional problem to the solution of a series of one-dimensional strips of cells in each coordinate direction. In the case of the cut-cell grid, this is still not, however, a truly one-dimensional problem, because cells within this strip may be irregular and we must therefore account for this multi-dimensional nature in the calculation of the interface fluxes. This is achieved here by the decomposition of each cell interface neighbouring a cut cell into regions that are ‘shielded’ by the boundary and those that are ‘unshielded’, and computing individual fluxes and for each of these regions, as displayed in figure 3. Here we assume that we are solving in the *x*-direction and that the normal of the boundary in the current mixed cell is in the negative *i*-direction (i.e. that *β*_{i−1/2}>*β*_{i+1/2}); other directions are implemented by the natural change of indices. Also we restrict ourselves presently to the two-dimensional problem, though these ideas also extend naturally to three dimensions.

The total flux at the interface in question is then assembled as an area-weighted sum: 4.5

We note that, in the case of concave boundaries, we may encounter situations in which the interface is shielded by boundaries in both neighbouring cells. This situation may be dealt with by further subdividing the shielded area into singly shielded and doubly shielded areas and treating each of these individually. However, as the test problems investigated here involve only convex boundaries, we may ignore this case here in order to simplify the explanation.

The unshielded portion of the flux may then be taken to be the standard one-dimensional flux for regular cells since it does not ‘see’ the boundary, i.e.
The shielded portion of the flux is computed using the one-dimensional stabilization approach introduced in §4*a*, with the effective one-dimensional irregular cell length taken as the average distance to the interface in the current coordinate direction,
where is the shielded volume fraction. Hence

#### (ii) Splitting flux

In using directional operator splitting in conjunction with cut cells in a straightforward fashion, one faces a major obstacle that, if not addressed, renders the entire approach useless. We explain this on the basis of an example.

Consider a parallel, constant-velocity flow along an infinite slanted planar wall. Then the exact solution would simply maintain this steady-state flow. Given the initial data for this flow, we know the exact fluxes across any of the grid cell interfaces of the regular and cut cells, as well as the pressure on the wall, which matches the constant pressure in the rest of the flow field. The first ‘*x*-step’ of the operator splitting cycle would then yield a mass update for the cut cell in figure 3*a* according to
4.6and the constant state would be severely distorted by order *O*(1) changes immediately.

As far as we know, this issue has thus far prohibited the use of directional operator splitting in conjunction with cut-cell approaches. We overcome this difficulty here by re-interpreting the notion of operator splitting in the context of finite-volume methods as follows: in, say, an *x*-step, all grid cells receive updates that accumulate fluxes *across all of their interfaces*. In computing these updates, interfaces facing the *x*-direction are assigned fluxes obtained through the procedures detailed in the previous section, whereas those interfaces facing the *y*-direction are assigned an auxiliary flux , i.e. the exact Euler fluxes computed on the basis of the data from the previous time step.

The introduction of these auxiliary fluxes first of all does not introduce any inconsistencies. On regular cells, the auxiliary fluxes even cancel out in each step. For cut cells, the sum of all auxiliary fluxes accumulated over one operator splitting cycle also cancels out. As a proof of concept, the auxiliary fluxes guarantee that a constant parallel flow will be maintained up to machine accuracy in each split step.

It is a straightforward matter to show that the above is equivalent to applying the standard operator splitting technique, considering only those cell faces whose normal points in the *x*-direction, but using the *flux differences* instead of the unmodified numerical fluxes **f**_{i+1/2} alone. This is actually the version implemented to run the examples shown below.

In particular, prior to the first one-dimensional sweep of each time step, a reference state, is computed for each cut cell. This state is constructed as a tangential flow at the boundary within the cell, where the reference pressure is computed from the Riemann problem normal to the wall, and the velocity is obtained by retaining only the component of the cell velocity tangential to the boundary, i.e. . We may then define the auxiliary flux to be the flux function evaluated at this reference state, e.g.
4.7The update within the cut cell is then given by
4.8
4.9
4.10where *V* =Δ*x* Δ*y* is the regular cell volume. As discussed above, we note that the total contribution from the auxiliary flux cancels out over the full splitting cycle.

### (c) Gravity term well-balancing

As discussed in Botta *et al*. (2004), Greenberg & Le Roux (1996) and LeVeque (1998), nearly hydrostatic flows require careful treatment of the hydrostatic balance, particularly near steep topography. Numerical imbalances between the gravitational source and pressure divergence term may upset the balance, producing spurious vertical motions.

We follow the technique introduced in Botta *et al*. (2004), in which a local, time-varying hydrostatic reconstruction is employed within each cell in the computational domain. This hydrostatic reconstruction is dependent upon a local reconstruction of the entropy field within the cell, though, for many situations, a constant-entropy approximation is seen to be sufficient. Making this assumption, the resultant reconstruction may be computed from the cell average states **w**_{i,j}=(*ρ*,*p*,*u*,*v*)_{i,j} as
4.11where is a piecewise linear approximation of the potential temperature between grid cell centres and *ξ* is the offset from the cell centre in the direction of the gravity vector (here simply the *y* offset). Transforming back to conservative variables, , we can thus compute modified neighbouring states in which only the non-hydrostatic portion of the variation is retained as
4.12The chosen reconstruction method, here the LDLR, is then applied to produce a (suitably limited, if necessary) reconstruction of these modified states, (rather than a reconstruction of the full conservative states), and the cell interface states are obtained by summing the hydrostatic and deviation reconstruction at the relevant position; for example, the right interface state for the cell is given by
4.13

For cut cells, the above process remains entirely unchanged, and boundary interface states are reconstructed by the evaluation of equation (4.13) at the centre of the boundary interface .

The evaluation of the gravitational source term is performed using the same hydrostatic reconstruction as that employed for the flux approximation above, with the general form for the source on cell *V*_{i,j} given by
4.14

For regular cells, this simplifies to 4.15

For cut cells, we employ the same division into shielded and unshielded cell portions as described in the previous section. We are then able to compute individually the source term in these two portions. In the shielded portion, 4.16where the signs ± are determined by the sign of the boundary normal to the direction under consideration. In the unshielded portion, the source is equivalent to the regular cell source scaled by the unshielded interface area. This gives a total source term for the cut cell of 4.17

### (d) Summary

The following list provides a step-by-step summary of the method.

Compute and store the reference state

*w*^{0}for all cut cells.For each directional sweep of a Strang splitting:

(a) compute the hydrostatic reconstruction on all cells and evaluate at cell interfaces facing the current coordinate direction,

(b) divide each mixed cell into shielded and unshielded portions in the current coordinate direction, as illustrated in figure 3

*a*,(c) obtain cell interface states for both regular and boundary interfaces by applying the chosen spatial reconstruction to neighbour cell states modified to discount hydrostatic variations,

(d) solve the collection of Riemann problems determined by these cell interface states to approximate the regular fluxes at all cell interfaces normal to the current coordinate direction, including those for mixed cells where it is potentially unstable,

(e) for cut cells, compute the stabilized fluxes for all shielded portions of grid cell interfaces,

(f) assemble a total flux for these interfaces by an area-weighted sum of the regular flux and the stabilized flux for the unshielded and shielded portions, respectively, and

(g) update all cells using the explicit update formula incorporating the auxiliary flux computed at the reference state.

## 5. Results

### (a) Zeppelin test

As a test of the gravitational balancing within the cut-cell framework and a demonstration of the flexibility of the meshing technique, we first consider an artificial two-dimensional configuration comprising a stationary elliptic zeppelin suspended within a vertically stratified atmosphere as depicted in figure 4. The computational domain has dimensions *L*=2, *H*=1, corresponding to a dimensional domain size of 20 km×10 km, with solid wall boundary conditions imposed at the top and bottom of the domain. The zeppelin has dimensions *l*_{0}=0.466 and height *h*_{0}=0.333, also with solid wall boundary conditions at the surface. A periodic boundary condition is applied at the horizontal domain edges. The initial condition is that of a homentropic atmosphere at rest and therefore the state should theoretically remain constant for all future time. However, as demonstrated in Botta *et al*. (2004), the numerical simulation of such a situation using a method that is not correctly well-balanced results in the generation of relatively large vertical velocities. The stratification used for the initial condition is given by
5.1

Running this simulation for 500 computational time steps results in a maximal vertical velocity that remains below 4×10^{−15} throughout. These errors are at the order of the machine precision and are caused by numerical round-off errors in the finite-precision arithmetic, demonstrating that the balanced state is maintained over this period.

To introduce a dynamic component to this zeppelin configuration, we now consider a modification that introduces a circular region of increased potential temperature below the zeppelin surface. The potential temperature profile is given by
5.2where , with *R*_{B}=0.15, *x*_{B}=−0.23, *y*_{B}=0.2 the radius and position of the perturbation and *θ*_{B}=0.0333 the maximum temperature increase.

Figure 5 plots the computed potential temperature evolution at five representative time points in the simulation. The initial condition is seen in (figure 5*a*), with the increased temperature ‘bubble’ located beneath the left edge of the zeppelin. Buoyancy effects result in this bubble immediately beginning to rise in the atmosphere. In figure 5*b*, the bubble first interacts with the zeppelin surface. On impacting the surface, the high-temperature region is divided into two pieces (figure 5*c*), with the left piece continuing to travel clockwise around the boundary. Eventually, this packet of hotter fluid becomes detached from the zeppelin again (figure 5*d*) and continues to rise into the atmosphere, while being stretched horizontally and generating vortices at its extrema. Meanwhile, the second, smaller, portion of the bubble remains trapped beneath the zeppelin for longer, eventually leaving the surface in a similar manner on the opposite side of the obstruction. While it is not possible here to assess these results quantitatively, the qualitative behaviour is as expected.

### (b) Lee waves

To test the presented algorithm in the context of atmospheric flows, we consider the familiar example of the generation of gravity waves in the downstream flow of a mountain (Durran 1986), also previously considered as a test problem for cut-cell discretization in Gatti-Bono & Colella (2006) (hereafter GBC). The configuration for the problem is shown diagrammatically in figure 6. All parameters are taken to be the non-dimensional equivalents of those applied in GBC.

The lower boundary implemented in the simulations is that of a ‘witch of Agnesi’, with a height profile described by the function
5.3where *a* is the characteristic width and *h*_{0} is the maximum height of the hill, here taken as *a*=1 and *h*_{0}=0.06. The initial condition is that of a two-layer atmosphere in hydrostatic balance, with an interface between the two layers at *y*=*H*_{d}. The thermal profile in each of these layers is chosen such that the Brunt–Väisälä frequency is constant within the layer, with *N*=*N*_{l} in the lower layer and *N*=*N*_{u} in the upper. This results in the profile
5.4where *θ*_{0}=1 here. The pressure and density may then be computed from this temperature profile using the hydrostatic relations. Equivalent to GBC, we consider three different layer combinations denoted as cases (a)–(c), as detailed in table 1. The initial flow velocity is purely horizontal, given by *V*_{0}=(*u*_{0},0) with *u*_{0}=0.06816 for all the results plotted here.

The simulation domain has dimensions *x*∈[−7.2,12.8] and *y*∈[−0.04,1.28] and is discretized with a regular grid of 300×100 cells, giving Δ*x*=0.06 and Δ*y*=0.0132. Note that these cells are anisotropic, and, as a result, the computational time step is limited by the more restrictive vertical Courant–Friedrichs–Lewy (CFL) condition. Periodic boundary conditions are applied at the left and right edges of the domain and a solid wall condition used at the top boundary. To minimize gravity wave reflection from the top of the domain and transmission across the periodic boundary, an absorption layer is implemented on the top, left and right edges of the domain. Within this zone, source terms are added to the system to relax the flow towards the ambient atmospheric state as defined by the initial conditions. This absorption zone has width 3.0 on the left and right of the domain and height 0.63 on the top of the domain, and the relaxation source employed within the domain is given by
5.5where and are the relative distance of the cell position within the horizontal and vertical layers, respectively, and *η*_{x}=0.83 and *η*_{y}=0.17 are the relaxation rates in each of these.

The results of the simulations are plotted in figure 7 at an output time of *t*=293.4, equivalent to a dimensional time of *t*′=10 000 s as used in GBC. This time is reached after approximately 35 000 time steps at a CFL of 0.8. The results compare qualitatively well with those presented in both GBC and Durran (1986), with equivalent wavelengths and comparable amplitude. Case (b) shows decay of the amplitude of the gravity waves similar to that displayed in GBC, where computations were performed at a lower resolution of 128×64 cells. The likely cause of the more significant damping seen here is the shorter time step limited by the acoustic wave speed in the present study resulting in larger numerical diffusion; the simulations of GBC incorporated an implicit treatment of acoustics, allowing for a larger time step and hence lower diffusion.

## 6. Conclusion

We have presented a relatively simple approach for the simulation of atmospheric flows on two-dimensional cut-cell meshes. The method operates stably in conjunction with directional operator splitting, and incorporates the well-balancing strategy of Botta *et al*. (2004) for nearly hydrostatic atmospheric flows. The resultant solver has demonstrated its ability to solve simple atmospheric flow problems, with results that compare well with previous publications.

The extension of the approach to three dimensions and concave boundaries proceeds along the same conceptual lines introduced here for the two-dimensional version, though the notation and exposition are necessarily more involved. These extensions are a requirement to solve real-world atmospheric problems, and we anticipate a further publication detailing the process and demonstrating its application for idealized and real-world simulations.

As noted previously, the fully explicit solver implemented here is computationally inefficient for the practical solution of atmospheric flows due to a time step limited by the acoustic wave propagation speed. Of interest therefore is the development of anelastic and pseudo-incompressible number model solver versions of this scheme by combining ideas from Klein (2009) (sound-proof model solver technique) and Oevermann *et al*. (2009) (cut-cell solver for Poisson-type projection equations). Also, an extension to higher than first order is currently under development.

## Footnotes

One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.

- © 2009 The Royal Society