## Abstract

The development and verification of a new version of the DieCAST ocean circulation model to be referred to as CANDIE (Canadian Diecast) are considered. Both CANDIE and DieCAST have many features in common with the well-known Modular Ocean Model (MOM) of the Geophysical Fluid Dynamics Laboratory. Of particular relevance to the present study are the rigid-lid approximation and the use of standard Cartesian coordinates. The DieCAST formulation in terms of the surface pressure, rather than the volume transport streamfunction, is also used in CANDIE to reduce numerical sensitivity to ocean depth variations. The major difference between MOM and DieCAST is the use of a mixed C and A grid formulation in DieCAST rather than the B grid formulation used in MOM. CANDIE differs from DieCAST in the use of a standard C grid formulation and a reduction in the magnitude of the time truncation error associated with the implicit treatment of the Coriolis force. The implementation of the rigid-lid approximation is retained.

To assess the reliability of both DieCAST and CANDIE, the authors have applied these models to a problem used by Haidvogel and Beckmann in a comparison of several different model formulations. The tests include the influence of a steep-sided coastal canyon that represents a significant challenge for the step topography of Cartesian coordinate models. Haidvogel and Beckmann’s tests show general agreement between models based on topography-following coordinates, but significantly different results were obtained with MOM. The results of DieCAST for the homogeneous test case also differ substantially from those of the *σ*-coordinate models, largely due to dissipation associated with low-order interpolations used adjacent to solid boundaries in DieCAST. However, the results of CANDIE are in good agreement with those of the *σ*-coordinate models for both homogeneous and stratified coastal canyon experiments. These results clearly demonstrate that the differences found for both MOM and DieCAST are not due to intrinsic limitations associated with the use of Cartesian coordinates.

## 1. Introduction

The choice of a vertical coordinate system for use in a numerical model has been a controversial issue in the ocean modeling community. The geopotential or *z*-level coordinate is simplest and has been widely used in the past. However, *z*-level models must represent the real bathymetry by a series of steps, which may lead to large truncation errors over steep topography (e.g., Gerdes 1993; Adcroft et al. 1997; Gnanadesikan and Pacanowski 1997). An alternative approach is to use terrain-following coordinates with the lowest surface conforming to the bottom. These coordinates are specifically designed to give improved representation of variable bottom topography, but they are subject to their own forms of systematic error, most notably the errors in the horizontal pressure gradients associated with the combination of stratification and bottom topography (Haney 1991; Mellor et al. 1994).

Recognizing the need for systematic comparison of coastal ocean models, Haidvogel and Beckmann (1998) (henceforth referred to as HB) use an idealized coastal canyon problem to compare the results of various models. Their results show that most of the models produce similar overall circulation patterns. However, the well-known Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model (MOM) (e.g., Bryan 1969; Cox 1984; Semtner 1986) gave results substantially different from those produced by other models. Since the GFDL model was the only *z*-coordinate model that contributed to this intercomparison, the question naturally arises whether the undesirable features associated with the GFDL model are common to other *z*-level models.

One goal of the present study is to compare the results presented by HB with at least one additional model based on *z* coordinates. For this purpose we chose to start with the C grid version of the DieCAST model (Dietrich/Center for Air–Sea Technology; see Dietrich et al. 1987). DieCAST has been successfully applied to the wind-driven circulation in Lake Neuchatel, Switzerland (Zuur and Dietrich 1990), and loop current eddies in the Gulf of Mexico (Dietrich and Lin 1994) and is also being applied to study the buoyancy and wind-driven circulation of the North Atlantic (D. Dietrich 1997, personal communication). DieCAST is a primitive equation, *z*-level model that uses the hydrostatic, Boussinesq, and rigid-lid approximations. In these respects, it is similar to the GFDL model. However, DieCAST differs from MOM in that it is formulated in terms of the surface pressure rather than the volume transport streamfunction. It should be noted that the latest version of MOM includes an optional surface pressure formulation (e.g., Pacanowski 1996). Therefore, unless otherwise stated, MOM in this paper refers to the traditional version used in the HB test cases, which is formulated in terms of the streamfunction. DieCAST also differs from MOM in that it uses a mixed C and A grid discretization rather than the B grid discretization used in MOM. Note that the rigid-lid approximation used in DieCAST and MOM eliminates fast-surface gravity waves, which limits the application of both models to frequencies such that *ω*/*f* ≪ (*r*_{e}/*L*)^{2}, where *ω* is the frequency, *f* is the Coriolis parameter, *r*_{e} is the external Rossby radius, and *L* is the length scale of the motion (Bryan 1969).

Initial results obtained with the DieCAST model for the coastal canyon test cases were not very satisfactory. This has led to some straightforward but important modifications of the original DieCAST code. The resulting circulation model will be referred to as CANDIE: the Canadian version of the DieCAST model. The second purpose of this paper is to introduce the new CANDIE model and discuss its similarities to and differences from DieCAST. The final purpose is to demonstrate the good agreement between the results obtained with CANDIE and those obtained with models utilizing topography-following coordinates.

The organization of this paper is as follows. The basic governing equations are reviewed in section 2. Time and space discretizations are discussed in sections 3 and 4. A brief discussion of the differences between CANDIE and DieCAST is given in section 5. We apply CANDIE to HB’s coastal canyon test problems and compare results with those produced by other coastal ocean models in section 6. An example of DieCAST results is also shown in section 6. A summary and discussion of results are presented in section 7.

## 2. The governing equations

We consider the three-dimensional (3D) primitive equations for an incompressible, stratified fluid using the rigid-lid, Boussinesq, and hydrostatic approximations. They are essentially the same as those considered by Dietrich (1992), except that spherical polar coordinates are used here:

where *u, υ, w* are the east (*λ*), north (*ϕ*) and vertical (*z*) components of the velocity; *p* is pressure; *ρ* is density; *T* and *S* are potential temperature and salinity; *K*_{m} and *K*_{h} are vertical eddy viscosity and diffusivity coefficients; *f* is the Coriolis parameter; *ρ*_{o} is a reference density; *R* and *g* are the earth’s radius and gravitational acceleration; ℒ is an advection operator defined as

and 𝒟_{m} and 𝒟_{h} are diffusion operators defined as

where *A*_{m} and *A*_{h} are horizontal eddy viscosity and diffusivity coefficients, respectively.

Under the hydrostatic approximation (3), the pressure field at depth *z* can be expressed in terms of the pressure at the rigid upper surface, *p*_{s}, plus that due to the fluid between the rigid lid and the depth *z*:

Lateral boundary conditions are required to close the above governing equations. No normal flow across solid boundaries is imposed, and the component of horizontal velocity tangent to solid walls satisfies either free-slip or no-slip boundary conditions. In the case of free-slip boundary conditions the tangential stress at vertical boundaries is set to zero, and for the no-slip boundary condition, the tangential velocity at the boundary is set to zero. Note that since the real bathymetry is represented by a series of steps in the *z*-level model, lateral boundary conditions are required not only at the coast but also at all submerged vertical boundaries.

The boundary conditions for the vertical velocity component *w* are based on the rigid-lid approximation and the condition of no normal flow at the bottom. For general topography these conditions can be expressed as

where the subscripts 0 and −*h* indicate evaluation at the sea surface *z* = 0 and at the sea bottom *z* = −*h,* respectively.

Integrating the continuity equation (4) and using (11) gives

Equation (14) holds for general topography and is valid in both *z*- and *σ*-coordinate formulations. Note that for the steplike topography of *z*-coordinate models, *h* is uniform over each individual grid cell and there is no flow through solid boundaries, so that in the model *w* is set to zero at the bottom.

## 3. Time discretization

As discussed in the introduction, both DieCAST and CANDIE are similar to the GFDL model, most notably in the use of Cartesian coordinates and the rigid-lid approximation. However, unlike the GFDL model, which uses a volume transport streamfunction, both DieCAST and CANDIE are formulated in terms of the surface pressure. The main advantages of using a pressure formulation are the increased ease with which islands are handled and the reduced numerical sensitivity to depth variations. The latter attribute is due to the appearance of *h*(*λ, ϕ*) [rather than 1/*h*(*λ, ϕ*) as in the streamfunction formulation] in the differential equation determining the surface pressure and reduces the need to smooth the model bottom topography to maintain numerical stability (e.g., Dukowicz et al. 1993). Again, we emphasize that in all references to the GFDL model, we are referring to the particular version used in the test cases presented by HB, which use the streamfunction formulation. Madala and Piacsek (1977), Killworth et al. (1991), Dukowicz et al. (1993), Dukowicz and Smith (1994), and Pacanowski (1996) have each reformulated the GFDL model using different pressure formulations, and it would be of interest to explore how these models perform on Haidvogel and Beckmann’s test cases. However, our interest here is in the approach used by Dietrich, which is distinct from each of those mentioned above. This approach is presented below.

Using forward time stepping for horizontal diffusion terms, centered time stepping for advection and pressure variations, backward (or fully implicit) time stepping for vertical diffusion, and either centered or backward time stepping for the Coriolis term, the time discretization of (1), (2), (6), and (7) can be written as

where Δ*t* is the time step and superscripts *n* − 1, *n,* and *n* + 1 represent the values of quantities at the corresponding time steps. The quantity *γ* is a measure of implicitness for the Coriolis term: *γ* = 0 gives centered (leapfrog) time differencing and *γ* = 1 gives backward (fully implicit) treatment of the Coriolis term. Dietrich et al. (1987) noted that, when the Ekman layer is resolved, an extremely short time step would be required if an explicit scheme were used for the vertical diffusion terms. To avoid this undesirable restriction, while maintaining an accurate representation of the Ekman layer dynamics, Dietrich et al. (1990) suggested that both the Coriolis and the vertical diffusion terms should be treated implicitly. The equations derived below allow for either implicit or explicit treatment of the Coriolis terms, but they treat the vertical mixing in both the momentum and tracer equations implicitly.

Clearly, since the Coriolis term may include a fully implicit contribution and the vertical diffusion is always treated implicitly, the state variables cannot be updated in the straightforward manner implied by the above equations. In practice, the updating is done through a sequence of steps equivalent to solving the above equations in combination with the continuity and hydrostatic equations at appropriate time levels.

Let us rewrite *u*^{n+1}, *υ*^{n+1}, *T*^{n+1}, and *S*^{n+1} in terms of a set of trial variables:

where the initial trial fields *ũ*^{n+1}, *υ̃*^{n+1}, *T̃*^{n+1}, and *S̃*^{n+1} are updated to include the effects of horizontal diffusion (forward), nonlinearity (centered), the pressure gradient force (centered baroclinic and forward surface contributions), and the contribution to the Coriolis force based on centered differencing. The fields *û*^{n+1}, *υ̂*^{n+1}, *T̂*^{n+1}, and *Ŝ*^{n+1} include the additional effects of vertical diffusion (backward) and the contribution to the Coriolis force based on backward differencing. Finally, *δu*^{n+1} and *δυ*^{n+1} represent the corrections to the velocity components required to allow for a time-centered surface pressure.

The initial trial fields *ũ*^{n+1}, *υ̃*^{n+1}, *T̃*^{n+1}, and *S̃*^{n+1} are given by

Each of the terms on the right sides of (23)–(26) are expressed in terms of previously determined state variables and hence are trivially calculated. Note the appearance of *p*^{n−1}_{s} rather than *p*^{n}_{s} on the right-hand sides of (23) and (24). The quantity *p*^{n}_{s} is not yet known, but it will be determined and properly accounted for in the final step of the solution procedure.

The advanced trial fields *û*^{n+1}, *υ̂*^{n+1}, *T̂*^{n+1}, and *Ŝ*^{n+1} are then determined by updating *ũ*^{n+1}, *υ̃*^{n+1}, *T̃*^{n+1}, and *S̃*^{n+1} to give the best current estimates of the effects of vertical diffusion and Coriolis:

where *F* = 2*f*Δ*t.* Note that we have used ∂*u*^{n+1}/∂*z* ≡ ∂*û*^{n+1}/∂*z* and ∂*υ*^{n+1}/∂*z* ≡ ∂*û*^{n+1}/∂*z* in the vertical diffusion terms on the right sides of (27) and (28) since the pressure correction terms *δu*^{n+1} and *δυ*^{n+1} are depth independent [see (33) and (34)]. Also note that *T*^{n+1} ≡ *T̂*^{n+1} and *S*^{n+1} ≡ *Ŝ*^{n+1} [see (21) and (22)] have been used in the vertical diffusion terms on the right-hand sides of (29) and (30).

At this point, *T̂* and *Ŝ* are completely updated, but *ũ* and *υ̃* are based on the use of *p*^{n−1}_{s} rather than *p*^{n}_{s} on the right sides of (23) and (24), and the effect spills over into (27) and (28) through the use of (*ũ*^{n+1}, *υ̃*^{n+1}) on the right sides of these equations. It remains to adjust the pressure terms to the time level *n* and update the velocities appropriately in order to achieve a fully time-centered pressure formulation.

Subtracting (27) from (15), then using (23) to eliminate *ũ*^{n+1}, and performing the corresponding operations on (16), (24), and (28) gives

in which *δ*(*u, υ*)^{n+1} = (*u, υ*)^{n+1} − (*û, **υ̂*)^{n+1} and *δp*^{n}_{s} = *p*^{n}_{s} − *p*^{n−1}_{s}. It follows that

Note that *δp*^{n}_{s}, *δu*^{n+1}, and *δυ*^{n+1} are all independent of *z.* Hence, they make no contribution to the vertical diffusion terms that are thus fully accounted by (27) and (28).

where the final equality is a convenient shorthand notation: the penultimate form is used in all numerical calculations. Note that (*û*^{n+1}, *υ̂*^{n+1}) does not in general satisfy the depth-integrated continuity equation. Therefore, *ŵ*^{n+1}_{−h} is not always zero, which leads to nonzero corrections *δ*(*u, υ*)^{n+1} and *δp*^{n}_{s}.

Equation (35) is the general equation determining the surface pressure correction *δp*^{n}_{s}. For the special case *γ* = 0 (centered Coriolis), this equation simplifies to

Equations (35) and (36) are analogous to Eq. (22) in Dukowicz et al. (1993).

Equations (35) and (36) are forced 2D elliptic equations for the surface pressure correction *δp*^{n}_{s}, with the forcing determined by the net divergence (or convergence) implied by the advanced trial horizontal velocity field (*û*^{n+1}, *υ̂*^{n+1}) over each vertical water column. Equation (35) or (36) determines the change in the surface pressure required to adjust the vertical velocity at the bottom to zero, as required by the steplike bottom topography, and is solved subject to the boundary condition determined from (33) and (34) with zero normal velocity through the boundaries of the model domain. Once *δp*^{n}_{s} has been determined from (35) or (36), the barotropic horizontal velocity corrections *δu*^{n+1} and *δυ*^{n+1} are determined by (33) and (34). Note that while the pressure is fully updated only to the time level *n,* the velocity, temperature, and salinity are each fully updated to the time level *n* + 1, as required by the above scheme.

It is of significance that the above equation for the surface pressure is linear. All nonlinear effects are included through the forcing term, which is fully determined at the central time level *n* so that no iterations are required. As a consequence, the surface pressure can be determined very efficiently using standard routines for elliptic equations. Also note that both DieCAST and CANDIE use the modified version of the stabilized error vector propagation technique described by Madala (1978) to determine the surface pressure correction *δp*^{n}_{s}. Roache (1995) discusses this computationally efficient and numerically stable elliptic marching method in some detail.

In summary, each time step proceeds as follows. We first determine the baroclinic contribution *p*_{b} to the pressure field at time level *n* and add it to the surface pressure from time level *n* − 1 to give an initial estimate of the pressure field at time level *n.* The initial trial fields (*ũ*^{n+1}, *υ̃*^{n+1}), *T̃*^{n+1}, and *S̃*^{n+1} are determined to allow for the pressure gradient, diffusive horizontal fluxes, advection, and an estimate of the Coriolis force based on the state variables at the previous two time steps. We then calculate the advanced trial fields (*û*^{n+1}, *υ̂*^{n+1}), *T̂*^{n+1}, and *Ŝ*^{n+1} to include vertical diffusion and any implicit contribution to the Coriolis terms. The surface pressure correction *δp*^{n}_{s} is then determined to eliminate the depth integrated horizontal divergence of (*û*^{n+1}, *υ̂*^{n+1}) over each water column. After numerically solving the governing 2D elliptic equation for *δp*^{n}_{s}, the barotropic velocity corrections (*δu*^{n+1}, *δυ*^{n+1}) are readily calculated through (33) and (34). Adding these velocity corrections to (*û*^{n+1}, *υ̂*^{n+1}) yields the total horizontal velocity components (*u*^{n+1}, *υ*^{n+1}) at the end of the time step. The vertical velocity component *w*^{n+1} is then determined at each *z* level by the horizontal divergence of (*u*^{n+1}, *υ*^{n+1}) through the continuity equation (4).

## 4. Spatial discretization and treatment of the Coriolis force

DieCAST and CANDIE both use a finite difference scheme based on Cartesian coordinates with unevenly spaced *z* levels in the vertical. The Arakawa C grid is used for the spatial discretization with state variables *u, υ, w,* and *p* defined on the staggered grid illustrated in Fig. 1. Note that temperature, salinity, and density are defined at *p* points and vertical variations in *p* are estimated from the overlying density field. The advection operator is conservative and is standard for the C grid (see Dietrich 1992).

The C grid is attractive and widely used in the community (e.g., MICOM of Bleck et al. 1992 and the MIT model of Marshall et al. 1997) because of the ease with which the control volume approach can be implemented. The main weakness of the C grid, on the other hand, is that the horizontal staggering of the *u* and *υ* points causes difficulties in estimating the Coriolis terms. DieCAST and CANDIE use different approaches to deal with this difficulty.

In DieCAST, the Coriolis force is treated implicitly (*γ* = 1). The approach suggested by Dietrich et al. (1987) and Dietrich et al. (1990) uses a blend of A and C grids, thus avoiding the computational difficulty mentioned above. The key to this approach is to interpolate the trial velocity components *ũ*^{n+1} and *υ̃*^{n+1} to the *p* points and update the advanced trial velocity components at the *p* points by integrating equations that include both vertical diffusion and Coriolis terms (the Ekman spiral equations). The updated velocity is then interpolated back to the staggered *u* and *υ* points. Note, however, that the blend of A and C grids may introduce significant numerical dissipation, as we now discuss.

Let (*ũ*^{n+1}_{c}, *υ̃*^{n+1}_{c}), calculated using (23)–(26), represent the trial velocity components at the *p* points. For example, they could be calculated using two-point averaging:

The indices *i, j,* and *k* denote the east, north, and vertical coordinates, respectively. The advanced trial velocity components at the *p* points (*û*^{n+1}_{c}, *υ̂*^{n+1}_{c}) are determined by the implicit equations [see (27) and (28)]:

These velocity components are then interpolated back to the *u* and *υ* points again to give *û*^{n+1} and *υ̂*^{n+1} at the *u* and *υ* points.

Results show that the smoothing involved in interpolating to the *p* points and back again introduces numerical dissipation. For example, if *γF* = 0 and *K*_{m} = 0, differences between *û* and *ũ* are due solely to the interpolation scheme. It is easily shown that linear interpolation to the A grid and then back to the C grid is equivalent to replacing *ũ*_{i,j,k} by 0.25*ũ*_{i−1,j,k} + 0.5*ũ*_{i,j,k} + 0.25*ũ*_{i+1,j,k} [i.e., numerical filter weights of (1, 2, 1)/4]. The smoothing effect is greater for waves with shorter wavelengths. Figure 2 shows the numerical dissipation of wave amplitudes *Â/Ã* as a function of the normalized wavenumber *k*Δ*x,* where Δ*x* is the grid spacing. The amplitude reduction is about 50% for waves with wavelengths of 4Δ*x* (i.e., *k*Δ*x* = *π*/2). Further, this dissipation is clearly magnified by reducing Δ*t* (Dietrich et al. 1990) since the interpolations are carried out more times over a given time interval. The smoothing effect of the interpolation scheme discussed above is nearly equivalent to harmonic lateral dissipation (∂*u*/∂*t* = *A*_{h}∂^{2}*u*/∂*x*^{2}) with the nondimensionalized diffusion coefficient *γ*_{h} = *A*_{h}Δ*t*/(Δ*x*)^{2} (see Fig. 2).

Two approaches have been suggested to reduce this numerical dissipation. The standard version of DieCAST uses a fourth-order interpolation scheme to reduce the dissipation associated with each pair of interpolations. Dietrich (1993) shows test cases for a doubly periodic domain that demonstrate the utility of this approach for regions that are removed from horizontal boundaries. The effective filter weights for the fourth-order interpolation scheme are (1, −18, 63, 164, 63, −18, 1)/256, corresponding to reduced but still significant numerical dissipation. Figure 2 shows that the amplitude reduction using the fourth-order interpolation is reduced to about 20% for waves with wavelengths of 4Δ*x.* For comparison, we also plot the amplitude reduction corresponding to biharmonic lateral dissipation (∂*u*/∂*t* = −*A*_{b}∂^{4}*u*/∂*x*^{4}) with the nondimensionalized diffusion coefficient *γ*_{b} = *A*_{b}Δ*t*/(Δ*x*)^{4} (see Fig. 2). Note, however, that near horizontal boundaries the original version of DieCAST uses a second-order scheme, and the associated dissipation remains substantial (see the discussion of the canyon test problem in the next section). Also, this approach still results in dissipation, which increases as 1/Δ*t.*

A second modification, which further reduces the dissipation, is to interpolate only the *changes* in velocity at *p* points back to the staggered *u* and *υ* points (e.g., Dietrich et al. 1990). Let

represent the changes in the trial velocity components at *p* points due to the Coriolis and vertical diffusion terms. Interpolating (*δu*_{c}, *δυ*_{c}) back to the *u* and *υ* points, the advanced trial velocity components on the C grid are determined by

If this approach is used, then the dissipation associated with the interpolations is reduced to zero for the special case *γF* ≡ 0 and *K*_{m} ≡ 0 considered above. Implementation of this approach in the standard DieCAST model can reduce numerical dissipation significantly, particularly when small time steps are required, but some dissipation remains. Further work is required to improve the accuracy adjacent to boundaries. This issue warrants further investigation as it may be critical for problems that are strongly influenced by the boundary conditions (e.g., Haidvogel et al. 1992). An example from the DieCAST model that includes both of the above modifications will be discussed in section 6.

One further point should be mentioned regarding the treatment of the Coriolis term in DieCAST. Although *û*^{n+1}_{c,i,j,k} and *υ̂*^{n+1}_{c,i,j,k} are updated using (39) and (40) with *γ* = 1, the surface pressure is computed using (36) and the associated velocity corrections are computed using (33) and (34) with *γ* = 0, respectively, that is, as if the Coriolis term were being treated explicitly. As noted by Dietrich et al. (1987), this leads to an error of order *F* = 2*f*Δ*t.* J. Sheng et al. (1997, personal communication) give an example where this error has significant effect even for *F* substantially less than 1.

In CANDIE we choose to treat the Coriolis force explicitly (*γ* = 0) and use the standard four-point averaging of *u* (*υ*) to determine appropriate estimates at the *υ* (*u*) locations for use in (23) and (24). This method is computationally inefficient if the Coriolis force is treated implicitly since it would require *û*^{n+1} and *υ̂*^{n+1} to be updated at all grid points simultaneously [see Xu (1994) for an example where this is done].

## 5. The differences between CANDIE and DieCAST

In the next section we will show results obtained from both DieCAST (e.g., Dietrich et al. 1990) and CANDIE for the test problems considered by HB. As mentioned in the introduction, CANDIE is a derivative of DieCAST obtained by implementing some straightforward modifications. Many of the changes have been cosmetic, but a few have important effects on the comparisons to be shown. The modifications listed below are the differences between DieCAST and CANDIE that have consequences for the results of these models. The first two points are important for the comparisons shown in the next section. The third point deals with improved accuracy of the sidewall boundary conditions, and the fourth point is a matter of numerical efficiency for a specific parameter choice that we use in the next section. The final two points do not affect results presented in this paper. They are important, however, for other applications of CANDIE.

CANDIE uses a standard C grid treatment of the model variables, whereas DieCAST uses a blend of the A and C grids. In particular, the vertical diffusion equations in CANDIE are integrated at

*u*and*υ*points separately, whereas in DieCAST they are calculated at the pressure points. It follows that in CANDIE interpolations between A and C grids are avoided.Coriolis terms in CANDIE are treated explicitly (

*γ*= 0) and computed using standard four-point averaging. As pointed out at the end of the last section, the surface pressure in the present version of DieCAST is determined by integrating (36) even though the Coriolis terms are treated implicitly. This approximation may result in significant error if implicit time stepping is used to increase Δ*t*such that*F*is not much less than unity. The tests with DieCAST reported below have*F*= 0.168 for the homogeneous case and 0.017 for the stratified case.At boundaries, the normal velocity is set to zero in both DieCAST and CANDIE. The normal flux of momentum for both the tangential and normal velocities is also set to zero in the original DieCAST code. CANDIE calculates the horizontal fluxes using a corrected representation of the coastal boundary conditions, in which either free-slip or no-slip boundary conditions are allowed for. The results using DieCAST reported in the next section also use the correct representation for the free-slip and no-slip boundary conditions.

In CANDIE, the equations dealing with the implicit treatment of vertical diffusion are solved through the inversion of a simple tridiagonal matrix equation. An iterative approach used in DieCAST was found to converge very slowly when the vertical eddy viscosity was given a very large value.

In contrast with the original version of DieCAST, CANDIE does not use the “swamp layer” numeric approach, in which land areas are replaced by shallow water regions. Currents over land areas are always set to zero in CANDIE. An iterative solution of the pressure equation on a regular domain is used to eliminate divergence errors in the depth-integrated volume transport adjacent to land boundaries. The most recent versions of DieCAST have also been modified to use this approach (Dietrich et al. 1996).

DieCAST does not include an explicit convection scheme. Tests with CANDIE gave poor results in areas where static instability develops. The vertical convection scheme discussed in the appendix of Wright and Stocker (1992) has been implemented in CANDIE to efficiently remove static instabilities.

## 6. Coastal canyon test problems

To provide a standard test case for the intercomparison of coastal ocean models, HB recently formulated a problem involving the rectification of oscillatory wind-forced flow over a coastal canyon. The results produced by a representative selection of models used by the international community of ocean modelers are compared in their paper. The models considered include the GFDL Modular Ocean Model (GFDLM; Bryan 1969; Cox 1984), the GeoHydrodynamics and Environmental Research Model (GHERM; Beckers 1991), the Miami Isopycnic Coordinate Ocean Model (MICOM; Bleck et al. 1992), the S-Coordinate Rutgers University Model (SCRUM; Song and Haidvogel 1994), the Spectral Element Model (SEOM; Iskandarani et al. 1995), the Semi-Lagrangian Shallow Water Model (SLSWM; Sanderson 1997); the Semi-Spectral Primitive Equation Model (SPEM; Haidvogel et al. 1986), and the Princeton Ocean Model (POM; Blumberg and Mellor 1987). A brief review of the dynamic and computational attributes of the above models can be found in HB.

The geometry used in the test problem consists of a periodic channel with a shelflike depth profile in the cross-channel direction, an isolated cross-shelf canyon, and inshore and offshore solid walls (see Fig. 3). The dimensions of the channel are 128 km in the along-channel direction and 96 km in the cross-channel direction. Results are determined for a local *f*-plane approximation. The coordinates *x* and *y* are aligned with the along-channel and cross-channel directions, respectively.

The circulation in the coastal channel is driven by a periodic along-channel wind stress, *τ*^{x}, given by

with *τ*_{o} = 10^{−4} sin(2*πt*/*T*_{w}) Pa kg^{−1} m^{3}, *L*_{y} = 96 km, *L*_{w} = 20 km, and *T*_{w} = 10 days. The along-channel wind stress is relatively uniform near the coast but decreases significantly with *y* over the shelf break region and is very weak over the deep water near the offshore wall (see Fig. 3). Detailed discussions of the mean circulation over a coastal canyon driven by a periodic, along-channel wind forcing can be found in Haidvogel and Brink (1986).

### a. A homogeneous coastal canyon experiment

We first consider a homogeneous coastal canyon experiment by setting the density to be uniform everywhere. We use Δ*x* = Δ*y* = 2 km; constant Coriolis parameter *f* = 10^{−4} s^{−1}; a linear bottom stress coefficient, *r*_{b} = 3 × 10^{−4} m s^{−1}; *g* = 9.81 m s^{−2}; and *A*_{m} = 5 m^{2} s^{−1}. The time step is Δ*t* = 1728 s or 50 time steps per day, which is comparable to that used by MOM and SPEM. Twenty unevenly spaced *z* levels were used in the vertical with greater resolution near the upper surface. Cell boundaries are at depths of 0, 10, 20, 30, 40, 60, 79, 107, 149, 209, 295, 417, 585, 807, 1090, 1430, 1812, 2208, 3104, and 4000 m. The solid lines in Figs. 3c and 3d are the bottom profiles corresponding to *x* = 63 km and *y* = 19 km, respectively, and the dotted lines correspond to the model representations of the bottom topography for these sections.

Note that the homogeneous model intercomparisons made by HB considered depth-independent (2D) circulation by applying the surface wind and bottom stresses as uniform body forces over the water column. An alternative approach, which we use here, is to run a 3D model with very strong vertical mixing and determine the depth-mean flow from the 3D model results. To ensure vertically uniform horizontal currents we use a vertical eddy viscosity coefficient of 10^{4} m^{2} s^{−1}. With this eddy viscosity, vertical shear is eliminated over a depth of 1000 m on a timescale of order 100 s, so that the entire water column is effectively homogenized every time step.

CANDIE was integrated from a state of rest for a period of 120 days. The model results over the last 30 days (equivalent to three oscillation periods) were averaged to yield time mean (residual) fields as suggested by HB. The residual surface currents calculated using CANDIE (see Fig. 4) are dominated by a net flow in the direction of shelf wave propagation (the prograde direction) and an anticyclonic gyre over the shallow water region near the canyon. The currents over the deep water are relatively weak. The overall circulation patterns produced by CANDIE are consistent with those estimated by *σ*-coordinate models such as SPEM and MICOM (compare Fig. 4 with Fig. 7 of HB).

To quantitatively assess the skill of the CANDIE model, we calculated the along-channel, depth- and time-averaged residual *u*^{xzt} and compared it with the results obtained from the other models considered by HB (Fig. 5; the data points for the coastal models other than CANDIE were read from Fig. 6 of HB). In the case of free-slip boundary conditions, the cross-channel profile of *u*^{xzt} produced by CANDIE agrees reasonably well with those produced by all other models except SLSWM (see Fig. 5a). Unfortunately, there were no GFDL model results available in this case since the free-slip boundary conditions are difficult to implement for the B grid GFDL model. The CANDIE results in Fig. 5a are, therefore, particularly important since they represent the only results for this test case obtained from a *z*-coordinate model. The maximum surface residual |**v**(*z* = 0)^{t}|_{max} produced by CANDIE is about 12.9 cm s^{−1}; the maximum along-channel, depth- and time-averaged current, *u*^{xzt}_{max}, is about 3.7 cm s^{−1}; and the net transport through the channel *hu*^{xzt} is about 0.325 Sv (1 Sv = 10^{6} m^{3} s^{−1}). These values again compare well with the results produced by all other coastal models except SLSWM (see Table 1).

In contrast with the results for the free-slip case, large differences occur among the cross-channel profiles of *u*^{xzt} produced by the various models in the case of no-slip boundary conditions (cf. Figs. 5a,b). Nevertheless, the results produced by CANDIE are relatively smooth and most consistent with the SPEM results. The results produced by the GFDL model (which only allows no-slip boundary conditions) are less smooth in spite of the use of higher horizontal eddy viscosity (15 m^{2} s^{−1}), which was required for model stability. The maximum surface residual |**v**(*z* = 0)^{t}|_{max} predicted by CANDIE is about 11.8 cm s^{−1}, the maximum value of *u*^{xzt} is about 2.8 cm s^{−1}, and the net transport *hu*^{xzt} is about 0.21 Sv in the no-slip case. These values are again very close to those produced by SPEM (see Table 1).

We have also used CANDIE to investigate the sensitivity of model results to the value of horizontal eddy viscosity. Figure 6 shows the cross-channel profiles of *u*^{xzt} using *A*_{m} = 5 and 15 m^{2} s^{−1}. Clearly, *u*^{xzt} is substantially weaker for higher horizontal eddy viscosity, particularly in the no-slip case. When *A*_{m} is increased from 5 to 15 m^{2} s^{−1}, the maximum value of *u*^{xzt}_{max} is reduced by about 20% in the free-slip case and by about 55% in the no-slip case. This result supports the speculation by Haidvogel and Beckmann that the weaker currents determined by the GFDL model may be a consequence of the larger value of *A*_{m} required for model stability.

An additional question raised by the results of the GFDL model shown by HB is the role of vertical resolution in determining the model results. The models based on topography-following coordinates have an advantage in this respect since they can reasonably represent rapid, but smooth, depth variations even with a single model grid cell. On the other hand, *z*-coordinate models require numerous levels to represent the bottom topography even though there may be no vertical structure in the current field. For example, the GFDL model requires more than 50 vertical levels to achieve convergence for the topography used in the present test case. Thus, it is clear that the response of a homogeneous fluid to depth-independent forcing is most efficiently studied with the use of topography-following coordinates. However, more realistic problems will require increased vertical resolution for either type of model formulation, and it is of interest to consider whether or not a *z*-coordinate model can adequately represent the abrupt topography of this problem with vertical resolution typical of that used in realistic model simulations.

To test the sensitivity of the CANDIE results to the number of *z* levels used, we redid the test problem using 10 and 40 *z* levels, respectively, and compared the results with those obtained with 20 *z* levels (a typical value used in realistic simulations) in Fig. 7. The cell boundaries are at depths of 0, 20, 40, 60, 100, 200, 500, 1000, 2000, and 4000 m for the 10-level case; and at 0, 2.5, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 69, 79, 93, 107, 128, 149, 179, 209, 252, 295, 356, 417, 501, 585, 696, 807, 949, 1090, 1260, 1430, 1620, 1812, 2010, 2208, 2656, 3104, 3552, and 4000 m for the 40-level case.

First, consider the comparison between results obtained with 20 and 40 vertical levels. Although there are some differences in *u*^{xzt}, particularly at the coast where the relative depth changes between adjacent grid cells are largest, the differences are much smaller than the differences between the results of the different *σ*-coordinate models considered by HB. Even with 10 levels in the vertical the results are within the range of uncertainty represented by the differences between the results of the different *σ*-coordinate models. It does not appear that the resolution of bottom topography by *z*-coordinate models is a major limiting factor in determining accurate model results for this particular test problem. Since this test is designed to represent the effects of abrupt topography, this is clearly an encouraging result for *z*-coordinate models.

Finally, we consider the significance of the differences between DieCAST and CANDIE. Results obtained from DieCAST, using the explicit Coriolis term and a second-order scheme to interpolate the *changes* in velocity at *p* points back to the *u* and *υ* points discussed in section 4, are presented in Fig. 6. These results are not significantly changed if the second-order interpolations are replaced with fourth-order interpolations away from the boundaries. Clearly, DieCAST underestimates the residual *u*^{xzt} in both the free-slip and no-slip cases, and it produces weak retrograde flows near the inshore wall that are not present in the results from CANDIE, indicating that numerical dissipation is still large near vertical boundaries in DieCAST. The differences are even more substantial if the second-order interpolations are retained and the total velocity at *p* points is interpolated back to the *u* and *υ* points.

### b. A stratified coastal canyon experiment

We now consider the stratified coastal canyon experiment discussed by HB. Note that the combined effects of steep topography and strong stratification have been a great challenge to all the coastal ocean models (Haidvogel and Beckmann 1998). To begin the process of identifying key model sensitivities, HB simply reformatted the homogeneous form stress problem discussed in the last section by relaxing the assumption of uniform density and body forcing.

The water density in this experiment is initially horizontally uniform but vertically stratified according to

in units of kg m^{−3}, with *H*_{ρ} = 800 m and *ρ*_{o} = 1000 kg m^{−3}. This gives a first baroclinic Rossby radius of deformation of 40 km in the deepest part of the channel.

The vertical eddy viscosity is given by

in units of m^{2} s^{−1} with *H*_{k} = 50 m. Therefore, the vertical eddy viscosity is maximum and approximately equal to 10.5 × 10^{−3} m^{2} s^{−1} at the top and bottom, and it decays into the interior with an *e*-folding scale of *H*_{k} (Haidvogel and Beckmann 1998).

As in the homogeneous problem, we use 20 unevenly spaced *z* levels and 2-km resolution in the horizontal directions, but with a reduced time step of Δ*t* = 172.8 s, or 500 timesteps per day. The horizontal eddy viscosity and diffusivity coefficients are both set to 20 m^{2} s^{−1} and the vertical diffusivity coefficient is set to 10^{−4} m^{2} s^{−1} (see Table 2). The other model parameters are all set to the same values as in the homogeneous coastal canyon experiment.

Figure 8 shows four snapshots of the evolving horizontal currents at a depth of 100 m for the case of free-slip boundary conditions. These results illustrate the evolution over a single 10-day cycle of the periodic wind forcing near the end of the model run. The wind switches from westward to eastward at day 100 and returns to being westward after day 105. The currents at day 102 are relatively weak everywhere except in the mouth of the canyon. At days 104 and 106, the currents remain weak in the far field but are relatively strong over the canyon and the shelf break. By day 108, the currents are again relatively weak everywhere except in the mouth of the canyon.

Figure 9 shows the corresponding evolution of the density anomaly (Δ*ρ,* the difference between the density field at the given time and the initial density field) at the depth of 100 m. The far-field density anomalies are relatively small. At day 102, the near-field density anomaly is positive close to the downstream canyon wall but negative over the other areas. At days 104 and 106, the negative density anomalies are of similar magnitude but have developed many small-scale features in the mouth of the canyon. The density anomaly at day 108 is similar to that at day 106, except for a positive density anomaly in the mouth of the canyon. Note that the winds are downwelling (upwelling) favorable on the southern (northern) side of the channel between days 101 and 105, and the situation is reversed from day 106 to day 110. The presence of positive density anomalies in the northern region on days 104 and 106, and in the southern region on day 108, can each be attributed primarily to preceding periods of upwelling favorable winds. However, the picture is clearly complicated by the persistent advection by the residual flow field and the variable bottom topography.

The model results between days 90 and 120 were used to estimate the rectified flow for this period, as in the homogeneous experiment. However, as pointed out by HB, secular trends associated with the continuing evolution of the density field remain present at day 120, so these results do not correspond to a system in equilibrium.

The time mean currents at the depth of 100 m are weak in the far field but relatively strong over the shelf break and in the mouth of the canyon (see Fig. 10). Note that the results in both free-slip and no-slip cases are very similar, except that the currents are weaker in the no-slip case. The overall features of the mean flow at 100 m produced by CANDIE are quantitatively comparable with those produced by POM and SPEM (compare our Fig. 10 with Fig. 9 of HB).

The time mean density anomaly at 100 m is small in the far field but large over the shelf break, particularly on the upstream canyon wall, which is very similar to the results produced by SPEM and GHERM. Over the downstream canyon wall and in the mouth of the canyon, however, the time mean density anomaly produced by CANDIE is negative and near zero, quite different from the positive density anomaly produced by GHERM, SPEM, and POM. It seems likely that the differences in the representation of bottom topography make a significant contribution to these model discrepancies. Given that the *z*-coordinate model uses a rather crude representation of the bottom, it is tempting to attribute the observed differences to the inadequacies of this representation. However, the possibility of significant differences associated with pressure gradient errors in the *σ*-coordinate models cannot be ruled out. Overall, it is encouraging that the results agree as well as they do.

The differences between model results may also be partially a consequence of the different friction parameterizations used in the various models [SPEM used biharmonic friction parallel to *σ* surfaces and POM used the approach of Smagorinsky (1963) for the lateral diffusivity with zero vertical diffusivity, while CANDIE was run with Laplacian lateral diffusion (see Table 2)]. To demonstrate the model sensitivity to these subgrid-scale parameterizations, we consider the results obtained with CANDIE using different parameterizations. The results are discussed in terms of three bulk quantities: the maximum residual at 100 m ; the maximum along-channel, depth- and time-averaged current *u*^{xzt}_{max}; and the net transport through the channel *hu*^{xzt}.

We first ran CANDIE using different vertical eddy diffusivity coefficients, *K*_{h} = 10^{−5}, 10^{−4}, and 10^{−3} m^{2} s^{−1}, with the horizontal eddy viscosity and diffusivity coefficients both fixed at 20 m^{2} s^{−1}. The maximum residual and the net transport *hu*^{xzt} both increase with the increased vertical eddy diffusivity coefficient (see Table 3). The maximum residual *u*^{xzt}_{max}, on the other hand, is smallest for *K*_{h} = 10^{−4} m^{2} s^{−1} and largest for *K*_{h} = 10^{−3} m^{2} s^{−1}.

Figure 11 shows the cross-channel profiles of *u*^{xzt} for the three different values of *K*_{h}. In contrast to the results of the homogeneous experiment, *u*^{xzt} is quite large near the coast in the stratified experiment. It decreases significantly over the inner coastal region and then increases again to a maximum value at about 18 km offshore. The cross-channel structures are similar for free-slip and no-slip boundary conditions, but the residual is prograde everywhere in the case of free-slip boundary conditions, whereas there is a band of weak retrograde flow centered at about 8 km offshore for the no-slip case. The differences in *u*^{xzt} are not significant between *K*_{h} = 10^{−5} and 10^{−4} m^{2} s^{−1}, suggesting that diffusion of this magnitude does not play a dominant role in determining the evolution over a timescale of order 100 days. For *K*_{h} = 10^{−3} m^{2} s^{−1}, the solution shows a shift toward the results obtained for a homogeneous fluid, indicating that the homogenization over a period of 100 days is significant for this value of *K*_{h}. Note that for *K*_{h} = 10^{−3} m^{2} s^{−1} stratification will be diffusively eliminated over a 100-m depth interval on a timescale of order 100 days.

We also ran CANDIE using horizontal eddy viscosity and diffusivity coefficients of 10, 20, and 40 m^{2} s^{−1} with the vertical eddy diffusivity coefficient fixed at 10^{−4} m^{2} s^{−1}. As expected, the net transport *hu*^{xzt} and the maximum residual *u*^{xzt}_{max} both decrease with the increased lateral mixing coefficient (Table 4). With no-slip boundary conditions, the maximum residual also decreases with the increased horizontal mixing. However, with free-slip boundary conditions, the maximum residual is smallest using *A*_{m} = *A*_{h} = 10 m^{2} s^{−1}. Figure 12 shows the cross-channel profiles of *u*^{xzt} for three different values of *A*_{m} and *A*_{h}. It is apparent that the cross-channel variations of *u*^{xzt} are very sensitive to variations in the horizontal eddy viscosity and diffusivity coefficients over this range.

The small-scale density features shown in Fig. 9 are also sensitive to the variations in the horizontal eddy viscosity and diffusivity coefficients. Figure 13 shows the evolution of the density anomaly at the depth of 100 m with *K*_{h} = 10^{−4} m^{2} s^{−1} but using *A*_{m} = *A*_{h} = 40 m^{2} s^{−1}. The large-scale patterns in Figs. 9 and 13 are quantitatively similar. However, the small-scale features present in Fig. 9 are substantially reduced in Fig. 13. Note that the only differences between these two figures are that the horizontal eddy viscosity and diffusivity corresponding to Fig. 13 are larger by a factor of 2 than those used in Fig. 9.

While the above results are broadly consistent with expectations, there is a disconcerting aspect of this test problem that deserves mention. In both Figs. 11 and 12, it is apparent that the coastal boundary layer in which the tangential velocity is reduced to zero in the no-slip case has not been resolved by the 2-km resolution used in the tests. An estimate of the distance to which the coastal boundary effect diffuses over the forcing timescale gives some idea of the resolution that would be required very near the wall to properly represent its effects. Taking a timescale of *τ*_{F} = 10/*π* days and *A*_{m} = 20 m^{2} s^{−1}, we estimate this distance to be of order (*A*_{m}*τ*_{F})^{1/2} ≈ 2 km. Resolving this boundary layer would require several grid cells within this distance of the southern boundary. Unfortunately, this would introduce a difference between our results and those presented by HB that would confuse the present intercomparison, but it would clearly be of interest to examine the effect of increasing the resolution near the southern boundary.

## 7. Discussion and conclusions

The primary purposes of this manuscript have been to present a revised discussion of the basic modeling approach used in the DieCAST model developed by Dietrich (e.g., Dietrich et al. 1990) and then to verify the model results against the test cases presented by HB. One of the distinct features of the DieCAST model is its usage of a surface pressure formulation rather than a volume transport streamfunction. The main advantage of the pressure formulation is its elimination of computational errors associated with ocean-depth fluctuations that appear in the streamfunction formulation (e.g., Dukowicz et al. 1993). Note, however, that the latest version of MOM can use the pressure formulation.

The original version of DieCAST, however, did not perform well on HB’s test cases. The dominant source of differences between DieCAST and the other model results has been traced to excess momentum dissipation near vertical boundaries associated with the interpolations from the C grid to the A grid and back again. Away from boundaries, this problem is substantially reduced (although not eliminated) by using higher-order interpolations and/or interpolating only *changes* in velocity back to the C grid. However, the higher-order interpolations require significant modifications if they are to be used effectively near vertical boundaries, and it is found that simply interpolating velocity changes is not sufficient to eliminate the excessive dissipation experienced for the test cases considered. We have also found that the implicit treatment of the Coriolis force used in the original version of DieCAST suffers from an unnecessarily large time truncation error, which may have significant consequences even in the absence of boundaries.

The problems noted above have led to the development of a modified version of DieCAST that we refer to as CANDIE: the Canadian version of DieCAST. CANDIE retains the important formulation of the rigid lid in terms of the surface pressure, but it replaces the mixed A and C grid formulation for the implicit treatment of the Coriolis force with a standard explicit treatment on the C grid. This simple modification eliminates the excess dissipation associated with coastal boundaries and avoids the error in DieCAST’s implicit treatment of the Coriolis force. We have also shown how the pressure equation must be modified in order to accurately implement an implicit treatment of Coriolis effects within the present model formulation. Several other modifications of the original DieCAST code are discussed in section 5, but these are of lesser importance for the present application.

There has been much debate about advantages and disadvantages of using either a *z*-coordinate model or a *σ*-coordinate model for studies of ocean circulation. The primary concern about *z*-coordinate models is that the representation of irregular topography by a series of steps may introduce large truncation errors in the steep topography regions and poorly resolved bottom in the weak slope regions [e.g., see Gnanadesikan and Pacanowski 1997). For *σ*-coordinate models, the concern is that truncation errors associated with the calculation of the horizontal component of the baroclinic pressure gradients in a model with coordinates inclined to the horizontal may cause significant errors.

In assessing predictive skills of widely used coastal ocean models, Haidvogel and Beckmann (1998) recently proposed an idealized coastal canyon test case and compared the results of these models. While terrain-following models in these test cases are in qualitative structural agreement, the well-known *z*-coordinate GFDL Modular Ocean Model has great difficulty in reproducing similar circulation patterns. Particularly in the stratified coastal canyon test case, MOM significantly underestimated the residual flow over the canyon. These discrepancies emphasize the question of whether or not *z*-coordinate models accurately represent the flow near steep topography with the vertical resolution typically used in ocean models.

The generally good agreement shown here between the results obtained with the *z*-coordinate CANDIE model and the *σ*-coordinate models in Haidvogel and Beckmann’s coastal canyon test cases should be encouraging for users of either class of models. The fact that results produced by CANDIE are very similar to those produced by terrain-following models such as SPEM, SCRUM, and POM in the coastal canyon problems, clearly indicates that the use of *z* coordinates does not have such serious negative effects on the model results as sometimes thought. Similarly, the agreement between the residual circulation patterns obtained for the baroclinic test problem with the two different classes of models indicates that the pressure gradient errors associated with the *σ*-coordinate models may not be a cause for great concern. Unfortunately, the lack of analytical solutions for these test problems does not permit any definitive statements regarding the specific deficiencies of either class of models.

## Acknowledgments

We are grateful to Keith R. Thompson and Tony Song for making useful comments and suggestions during the course of this work. We received financial support from the Natural Sciences and Engineering Research Council of Canada, the Canadian Institute for Climate Studies, and the Environmental Research Program sponsored by IBM.

## REFERENCES

## Footnotes

*Corresponding author address:* Dr. Jinyu Sheng, Dalhousie University, Dept. of Oceanography, Halifax, NS B3H 4J1 Canada.

Email: sheng@phy.ocean.dal.ca