## Abstract

In a recent paper by Mellor et al., it was found that, in two-dimensional (*x, z*) applications with finite horizontal viscosity and zero diffusivity, the velocity error, associated with the evaluation of horizontal density or pressure gradients on a sigma coordinate grid, prognostically disappeared, leaving behind a small and physically insignificant distortion in the density field. The initial error is numerically consistent in that it decreases as the square of the grid increment size. In this paper, we label this error as a sigma error of the first kind.

In three-dimensional applications, the authors have encountered an error that did not disappear and that has not been understood by us or, apparently, others. This is a vorticity error that is labeled a sigma error of the second kind and is a subject of this paper. Although it does not prognostically disappear, it seems to be tolerably small. To evaluate these numerical errors, the authors have adopted the seamount problem initiated by Beckman and Haidvogel. It represents a stringent test case, as evidenced by their paper, wherein the model is initialized with horizontal isopycnals, zero velocity, and no forcing; then, any velocities that develop must be considered errors.

Two appendices are important adjuncts to the paper, the first providing theoretical confirmation and understanding of the numerical results, and the second delving into additional errors related to horizontal or isosigma diffusion. It is, however, shown that satisfactory numerical solutions are obtained with zero diffusivity.

## 1. Introduction

Numerical models using sigma coordinates have the capability of dealing with ocean applications with large topographic variability—ocean basins, shelf breaks, continental shelves, and estuaries and bays—and they can model bottom boundary layers. The latter are deemed to be important to the maintenance of baroclinicity in the ocean (Mellor and Wang 1996) and to the formation of deep water (Jungclaus and Backhaus 1994; Zavatarelli and Mellor 1995; Baringer and Price 1997a,b). An apparent disadvantage to the use of sigma models is the existence of a numerical, baroclinic pressure gradient error. The error can produce velocity errors that can be detected in the case of an initially horizontally homogeneous density field that, in theory, should produce zero velocities. A recent paper by Mellor et al. (1994, henceforth MEO) derived an expression for the error in determining the horizontal gradient of the density. Thus

Throughout the paper, Cartesian coordinates are denoted by (*x**, *y**, *z*) and sigma coordinates by (*x, y, σ*). In the above formula, *σ* = *z*/*H*(*x, y*), *H* is the bottom depth, and *b*′ ≡ *ρ*′*g*/*ρ*_{o}, where *ρ*′ is the density after subtraction of the horizontally averaged density, a strategy that reduces error (Gary 1973; Mellor et al. 1994). The grid-cell variation of *H*(*x, y*) is *δ*_{x}*H* due to the horizontal grid spacing *δx,* and *δσ* is the vertical grid spacing. The error in the *y* component of the density gradient is the same as (1) after appropriate substitution of *δy* and *δ*_{y}*H.* To find the error in the horizontal pressure gradient, the expression on the right of Eq. (1) is integrated with respect to *z.*

A concept of “hydrostatic inconsistency” has been advanced in several previous papers (Rousseau and Pham 1971; see also references in Haney 1991 and MEO) and defined by |(*σδ*_{x}*H*)/*Hδσ*)| > 1; we do not believe the concept is particularly meaningful. The error in (1) is actually proportional to (*δσ*)^{2} − *σ*^{2}(*δ*_{x}*H/H*)^{2} so that |(*σδ*_{x}*H*)/*Hδσ*)| = 1 represents a locus of minimum error. However, the important point is that the error limits to zero as (*δσ*)^{2} and (*δH*)^{2} → 0 on an error field of (*δσ*)^{2} versus (*δH*)^{2}. For typical applications, (*δσ*)^{2} is relatively small so that the gridpoint change in *H* dominates the error.

The model used to explore these matters in the first part of MEO was a *two-dimensional* version of the Princeton Ocean Model (POM) (Blumberg and Mellor 1987; Mellor 1996). Even without subtracting the area-averaged density from the initial density field, it was found that the initial velocity error, incurred after the model was spun up from rest diagnostically (density held constant), decayed to zero prognostically. The density field perturbed advectively to a new state that just canceled out the original numerical density error. Note again that this finding was a *two-dimensional* result, as is the error given by Eq. (1). Thus, the error that has been the focus of much previous discussion and journal exposure is, we believe, quite benign. We will henceforth refer to this as a sigma error of the first kind (SEFK).

In *three-dimensional* models, we have, in the past, detected an error that does not vanish prognostically. In this paper, we determine that it is a sigma coordinate error; it is a vorticity error; it is smaller than SEFK before decay; and it is nil for two-dimensional flows but not for three-dimensional flows. In the theoretical continuum limit, any initial error should decay advectively, as shown in appendix A. In a finite difference calculation, SEFK also decays; however, the vorticity error, which we will call a sigma error of the second kind (SESK), does not. We illustrate this error for the test seamount problem initiated by Beckman and Haidvogel (1993, henceforth BH).

In sections 2, 3, and 4, we specialize to a numerical model with zero horizontal diffusivity as in BH. Appendix A contains theory that aids qualitative understanding of the numerical results. In appendix B, an analysis and discussion of errors associated with horizontal and isosigma diffusion are provided.

## 2. The seamount problem with no forcing and zero diffusivity; SPEM results

The BH paper applied the Semi Spectral Primitive Equation Model (SPEM; see Haidvogel et al. 1991; Hedstrom 1990) to the case of the idealized seamount illustrated in Fig. 1. The SPEM model is a sigma coordinate model that has some similarities to the POM model. The major difference, relative to the topic of this paper, is that SPEM represents vertical variations of properties by a series expansion in Chebyshev polynomials, whereas the POM model discretizes vertical variations using finite difference algebra.

Beckman and Haidvogel (1993) applied their model to the case of the idealized seamount illustrated in Fig. 1. The seamount was defined by

A reference stratification was imposed according to

where Δ_{z}*ρ* is the density difference between the maximum water depth and the surface.

Their model was then linearized so as to preclude the need for horizontal diffusivity. A biharmonic viscosity was used. Beckman and Haidvogel (1993) next imposed an initial condition with a small perturbation density such that

and ran a series of numerical experiments with various Burger numbers and “slope parameters.” The Burger number, the ratio of the internal deformation radius to the horizontal advective scale *L* was defined by

where *N*^{2}_{o} ≡ *gH*^{−1}_{o}Δ_{z}*ρ*/*ρ*_{o} is the Brunt–Väisälä frequency. In (2) and (5), the depth is fixed at the value, *H*_{o} = 4500 m, *L* = 25 km, and the Coriolis parameter *f* = 10^{−4} s^{−1}. The slope parameter is a numerical resolution parameter defined by

where 0 < *r* < 1. Here *δH* is the difference in adjacent cell depths and *H* is the mean of the depths.

Beckman and Haidvogel (1993) ran the SPEM model with no external forcing so that, theoretically, there should be no motion. Motion, attributed to the sigma coordinate error, was obtained, however; the maximum velocity for the case *r*_{max} = 0.21 is plotted in Fig. 2 (BH’s Fig. 6) and the 10-day value is tabulated in Table 1 (their Table 2a with intermediate Burger numbers deleted). One sees that, for *S* ≥ 1, the error is linearly increasing at 10 days and the cases with *S* ≥ 6 have blown up or are in the process of blowing up. One is tempted to speculate that the other cases would suffer the same fate for long enough integration times. More recently, the error has been reduced by higher-order differencing of the horizontal pressure gradients (McCalpin 1994), and we have learned (D. Haidvogel 1996, personal communication) that SPEM has converted to vertical finite differencing. One supposes that the present paper will apply to the revised SPEM with lowest-order differencing.

## 3. The linearized POM model

We decided that BH had initiated a challenging test case and so configured the POM model for the seamount case. Beckman and Haidvogel (1993) linearized the model so that the horizontal diffusivity could be set to zero. Therefore, the governing equations in rectangular coordinates are

where *b*_{r}(*z*) ≡ *ρ*_{r}(*z*)*g*/*ρ*_{o} is the background density, *b*′*(*z*) ≡ *ρ*′*(*z*)*g*/*ρ*_{o} is the perturbation density, and *p** ≡ pressure/*ρ*_{o} is the kinematic perturbation pressure. Note that −∂*b*_{r}/∂*z* = *N*^{2} is the square of the Brunt–Väisälä frequency. The surface and bottom boundary conditions are *w** = ∂*η*/∂*t* and *w** = −*u**∂*H*/∂*x** − *v**∂*H*/∂*y**, respectively, and *η* is the surface elevation. The vertical viscosity and diffusivity were zero. The variables (*F*_{x}, *F*_{y}) denote horizontal viscous terms and, for this paper, they are similar but not identical to a conventional Laplacian form, as discussed below.

We now transform (7), (8a)–(8c), and (9) from the Cartesian coordinate system (*x**, *y**, *z*) to the sigma system (*x, y, σ*) according to

POM is a free surface model so that the surface elevation, *η*(*x, y, t*), is included in (10c). Otherwise, for the purposes of this paper, horizontal variations in *η* are insignificant relative to variations in *H* and can be excluded from discussion.

After some linearization to account for the fact that *η* ≪ *H,* the conversion relations for dependent variables are

Using (11a)–(11d), the transformed equations are

where *ω* ≡ *w* − *σu*∂*H*/∂*x* − *σv*∂*H*/∂*y* − (*σ* + 1)∂*η*/∂*t*;at the surface *ω*(*x, y,* 0) = 0 and at the bottom *ω*(*x, y,* −1) = 0.

In POM, the first two terms on the right side of (13a), using (13c), are written −*gH*∂*η*/∂*x* − *H *∫^{0}_{σ} [*H*∂*b*′/∂*x* − *σ*′(∂*H*/∂*x*)(∂*b*′/∂*σ*′)] *dσ*′ and similiarly in (13b). This is useful to specifically invoke the elevation *η* and to specifically identify the SEFK that occurs because the two terms of the integrand should exactly cancel for the case, *b*′ = *b*′(*σH*). For finite *δσ* and *δH,* they do not cancel and Eq. (1) is the error. On the other hand, the form in (13a), (13b) does simplify algebra and aids intuitive comprehension in section 4.

In accordance with Mellor and Blumberg (1985), the horizontal viscous terms are given by

The model grid replicated the BH grid. The number of grid points is 64 (*x* direction) × 48 (*y* direction) × 10 (sigma direction). The reference calculation used a grid size, *δx* and *δy,* that varied from 8 km far from the center of the seamount to 4 km close to the seamount, which yields a value of *r*_{max} = 0.21; the grid sizes may be varied in order to vary *r.* In this paper, POM uses 10 equally spaced sigma layers. The differences in the SPEM and POM were 1) the aforementioned vertical discretization; 2) SPEM has a rigid lid, whereas POM has a free surface; and 3) POM uses a harmonic viscosity, *A*_{M} = 2000 m^{2} s^{−1}, instead of a biharmonic viscosity, scaled here to conform approximately to the biharmonic coefficient used by BH at the grid size of 6 km. (Ordinarily, POM uses a nonlinear Smagorinsky diffusivity.) This *A*_{M} value is, however, quite large for the present grid size and we will shortly consider smaller values.

The result of the POM calculations is shown in Fig. 3 and Table 2. The behavior of the POM model differs significantly from that of Fig. 2 and Table 1. The POM errors are smaller and, more importantly, cease to grow after a short time so long as *S* ⩽ 8. Also, the POM model is able to deal with large values of *r*_{max}. Note that the maximum possible value with the present seamount geometry is 0.82; it is the case where the seamount is represented by a single grid point, where *H* = 450 m surrounded by grid points, and where *H* = 4500 m.

The linear system is numerically unstable for *S* ≥ 8;we have not identified the cause, but, in section 4, we will find that lowering the horizontal viscosity remedies this problem at a cost of increased numerical error. For *S* < 8, the error in Table 2 is very nearly a function of *r*_{max} and independent of *S* except for the fact that the *S* = 0 case seems to be peculiarly anomalous. However, this behavior can be justified theoretically, as shown in appendix A. Appendix A may be summarized as follows.

For

S= 0 and an error in the horizontal pressure gradient, the associated velocity field oscillates and decays in time and finally asymptotes to a constant value; the density field is unchanged. ForS> 0 and an error in the horizontal pressure gradient, the associated velocity field oscillates and decays in time and asymptotes to zero; a small compensating density change is created by the model.

In other words, there is a discontinuity in behavior of the system with respect to *S.* This partially explains the behavior of the solutions in Fig. 3. Of course, the *S* = 0 case is an artifact of the division of density into a background density and a disturbance density wherein the background density is nil and, relatively, the “disturbance density” is not small.

## 4. A sigma error of the second kind (SESK)

Although the velocity error exhibited in Fig. 3 and Table 2 are generally small, we find that they differ from that detailed in MEO, which disappeared prognostically in two-dimensional problems and which, according to appendix A, should decay. In Fig. 3, the errors do decay but asymptote to a nonzero value.

We now believe that there is another kind of error at play, a vorticity error. Thus, in forming vorticity from (8a), (8b), one uses ∂^{2}*p**/∂*x**∂*y** − ∂^{2}*p**/∂*y**∂*x** = 0; the same cancellation may be obtained analytically from the first two terms on the right sides of (13a), (13b). However, in the finite difference version, we find that cross-differentiation of (13a), (13b) results in a nonzero vorticity error.

Referring to Fig. 4, obtain the vorticity (actually the circulation after which division by *δx δy* yields vorticity) according to

(The first *i, j* location spans *u*_{i,j} in Fig. 4; the second spans *v*_{i,j}.) The right side is centrally differenced. The sigma coordinate, pressure gradient terms cancel, but the terms containing *b*′ and the depth gradients do not. Thus,

The vertical integral of (18), ∫^{0}_{−1} Γ_{p }*dσ,* for *S* = 4, *r*_{max} = 0.211, and *t* = 0 is plotted in Fig. 5. This error is labeled an SESK.

A vertically integrated diagnostic evaluation of (17) shows that the SESK is balanced by the viscous term; the Coriolis term, related to water column stretching, is small. This vorticity error is the *t* → ∞ asymptote displayed in Fig. 3. Except for the *S* = 0 case, it is independent of *S* simply because the initial disturbance *b*′ is independent of *S.* In Table 3, we give the results of calculations for a horizontal viscosity reduced by a factor of 10. For the smaller values of *r*_{max} one obtains a maximum velocity increased by a factor of 8 over that reported in Table 2, so there is near-linearity of error and inverse viscosity. Note that, somewhat surprisingly, the case *S* = 10 is now stable. Since this is a larger stratification than encountered in the ocean and a viscosity of 200 m^{2} s^{−1} is somewhat larger than typical for this resolution, we have not probed further to find the stabilty limits as a function of *S* and *A*_{M}.

Equation (18) is *the* error but difficult to comprehend. Therefore, expand the variables in a Taylor series around the center of the dashed stencil in Fig. 4b such that *H*_{i,j} − *H*_{i,j−1} = (*δ*_{y}*H*)_{o} + *δ*_{x}(*δ*_{y}*H*)/2, etc.; then

Now specialize for *b*′ = *b*′(*z*) = *b*′[*σH*(*x, y*)]. The first term on the right-hand side of (19) is nil, but the second term yields

In cylindrical coordinates *b*′ = *b*′(*r*) and *H* = *H*(*r*) [where *r* is the radial coordinate, not to be confused with (6)], so the problem is axisymmetric. However, transform to *H*(*x, y*), determine *δ*_{xy}*H* = (∂^{2}*H*/∂*x*∂*y*) *δx δy,* etc., and from (20), Γ_{p} = *F*(*r*) cos*θ* sin*θ* (cos^{2}*θ* − sin^{2}*θ*). Thus, the eight-lobed pattern in Fig. 5 is obtained; the temporally asymptotic streamfunction solution (not shown) also displays a similar pattern. One also sees that for two-dimensional problems *δ*_{xy}*H* is nil and so is the vorticity error. This explains why the SESK did not emerge in the two-dimensional calculations of Mellor et al. (1994). Note that, for fixed *H*(*x, y*) and *δx*/*δy* = constant, Γ_{p} ∼ *δx*^{4} as *δx* → 0.

## 5. Return to the nonlinear model

For easy reference, we repeat the full, nonlinear equations except that they are specialized to solve for density directly instead of temperature and salinity. Thus, if *D* ≡ *H* + η and we define the operator

then the equations of motion are

Note that, consistent with the previous linear system, *ρ*′ = *ρ* − *ρ*_{r} is used in the baroclinic integrals of (23a), (23b). Here, *F*_{x} and *F*_{y} are defined as in (15a), (15b), and (16). We have formally added isosigma diffusion to (24).

A primary finding is that, with the same initial conditions [(3) and (4)], this system of equations, when integrated numerically with *A*_{H} = 0, behaves almost exactly as did the linear system. The exception is the *S* = 0 case in Fig. 3, which now also decays to the asymptotic SESK value after about 5 days. Thus, in retrospect, there may have been no compelling need to deal with the linear equations. However, there are the virtues of connecting with the BH paper, the theoretical development of appendix A, and the very fact that we have established the equivalence of the linear and nonlinear sets of equations for the problem of this paper.

In appendix B, we find that horizontal and isosigma diffusion are additional sources of error. Therefore, we end this paper with a demonstration that this diffusion can be set to zero in a nonlinear calculation of the seamount problem with throughflow. We set a uniform velocity of 20 cm s^{−1} at inflow and radiation conditions for both the external and internal velocities at outflow. The initial density and density inflow, *ρ* = *ρ*_{r} + *ρ*′, is given by (3) and (4). The initial density field, after transforming to *σ* space, is subtracted prior to calculating the integrals in (23a), (23b) since its horizontal gradients in *z* space are nil. The run time was 30 days and there were no numerical instabilities. In Fig. 6, we show the differences between the instantaneous density at day = 2 and the initial, horizontally homogeneous, density at a depth of 2000 m. The SEFK velocity error is expected to vanish. For this calculation we have used the Smagorinsky diffusivity; in the vicinity of the seamount, it gave values in the range 50–120 m^{2} s^{−1}. By comparing the difference in calculated density and initial density (the difference being due to throughflow dynamics) with the results of Tables 2 and 3 (or Fig. 7 with *A*_{H} = 0 in appendix B), and since the error is linearly related to density perturbation and inverse viscosity, it is possible to make a rough estimate of the maximum SESK velocity error; it is less than 0.5 cm s^{−1}.

## 6. Summary

Two kinds of sigma coordinate errors, which we have labeled SEFK and SESK, have been identified. For the sigma coordinate error frequently discussed in the literature and evaluated in Eq. (1), the prognostic response from quiescent initial conditions is, first creation of a velocity error field, then creation of a compensating density perturbation field, and finally complete decay of the velocity error field. The SEFK velocity error vanishes due to advection of the density field, which causes a cancellation of the velocity error (Mellor et al. 1994), as demonstrated theoretically in appendix A. This SEFK behavior is operative for two- and three-dimensional problems. On the other hand, the SESK velocity error does not vanish but is small and can be made smaller by subtracting a horizontally averaged initial density before computing the baroclinic integrals. Smoothing the topography to reduce the largest value of |*δH*|/*H* can reduce error in both SEFK and SESK. Using a curvilinear grid that approximately follows bathymetric contours also reduces SESK.

As discussed in appendix B, addition of horizontal or isosigma diffusion adds other kinds of error. However, POM seems to execute smoothly with zero diffusion.

## Acknowledgments

This research was suported by Office of Naval Research Grant N0014-93-0037 and by the computational support and facillities of NOAA’s Geophysical Fluid Dynamics Laboratory. The constructive comments of reviewers resulted in an improved paper.

## REFERENCES

### APPENDIX A

#### The Oscillatory Decay of SEFK

Here we create an analytical model of the behavior of our numerical model. It will be justified only insofar as it conforms qualitatively to the behavior of our numerical solutions and helps explain their behavior. We will explicitly model the pressure (or density) gradient error due to bottom slope but, to simplify the mathematics, we will deal with a flat bottom. To justify this step, we have performed numerical calculations with the pressure gradient terms and error calculated using the seamount topography but otherwise setting *H* = constant = 4500 m in Eqs. (12), (13a)–(13c), and (14), thus mimicking the simplified analytical model. As expected, the solutions behave qualitatively the same as before, and, in fact, there are only 20%–30% quantitative differences.

##### The case of nonzero viscosity and zero diffusivity

We first deal with a linear system with nonzero viscosity and zero diffusivity. Thus we write equations (7), (8a)–(8c), and (9) according to

where we have eliminated pressure and at the same time introduced a “numerical error,” ɛ(*x, y, z*), which is spatially variable but constant temporally. The error is irrotational and thus models only SEFK.

For an initial *b*′(*x, y, z,* 0) = *b*^{′}_{o}(*z*), ∂*b*^{′}_{o}/∂*x* = ∂*b*^{′}_{o}/∂*y* = ∂*b*^{′}_{o}/∂*t* = 0 so that *b*^{′}_{o}(*z*) can be subtracted from *b*′ and added back after the time-dependent solution has been obtained. Thus the initial conditions are

For *S* = *N*^{2} = 0, it can be seen from (A3) and (A4d) that *b*′(*x, y, z, t*) = 0 and according to (A2a), (A2b), there must be a nonzero velocity error for all *t* > 0.

To determine the behavior for *S* > 0 is complicated. Nevertheless, solutions can be obtained as a particular solution, *u*_{p}, *v*_{p}, *w*_{p}, *b*^{′}_{p} = 0, 0, 0, −ɛ, plus a complicated Fourier manifold of homogeneous solutions, which can be obtained as a three-dimensional sum of modes of the form

Furthermore, the particular solution can also be projected on to a steady-state mode in the form of

where *E*_{n}(*z*) will subsequently be related to ɛ and *W*_{n}(*z*).

Essentially the problem is to find *b*^{′}_{h} so that *b*′ = *b*^{′}_{p} + *b*^{′}_{h} = 0 at *t* = 0 and then determine the behavior of *u, v, w, b*′ for *t* > 0. Now, substitution of (A5a)–(A5d) into the homogeneous form (ɛ = 0) of Eqs. (A2a), (A2b), and (A3) yields

where *υ* ≡ (ℓ^{2} + *m*^{2})*υ.* Substitution of the above into Eq. (A1) yields

where the dependencies *σ* = *σ*_{n}(ℓ, *m*), are understood. Subject to the boundary conditions, *W*_{n}(0) = *W*_{n}(−*H*) = 0, *W*_{n}(*z*), and *λ*^{2}_{n} may be determined as an eigenvalue problem (for *N*^{2} = constant = *N*^{2}_{o}, *W*_{n} ∝ sin*λ*_{n}*z,* and *λ*_{n} = *nπ*/*H* for *n* = 1, 2, 3, . . .).

Now the dispersion relation (A9) may be written

where

Finding the roots of the cubic in (A10) is complicated but for small *A,* we find that the three roots are

The imaginary parts of *σ* = *fω,* when inserted into (A5a)–(A5d), yields a zero or negative exponent for *e*^{iσt}.

The initial condition for *W*_{n} is nil, so the solutions corresponding to the first root can be discarded. The second and third roots can be written *σ*_{r} + *iσ*_{i} and −*σ*_{r} + *iσ*_{i} so that the *w* and *b*′ solutions are

where we have set *a*_{2} = −*a*_{1} = *i*/2. Since the homogeneous solutions may be multiplied by any constant, to satisfy (A4d), we have

and *E*_{n}(*z*) = *N*^{2}(*z*)*W*_{n}(*z*). If we invert

to determine *A*_{n}(ℓ, *m*), then *b*′ and *w* can be determined as well as *u* and *v* from (A7a), (A7b). In particular,

The double integrals over ℓ and *m* can be converted to discrete summations for a finite horizontal domain. For an infinite domain and, say, a horizontal, Gaussian error distribution, Fourier transforms are especially simple.

However, without expending further effort, we learn from Eq. (A15) that *b*′ asymptotes to a nonzero field as *t* → ∞ whereas, from (A16), *w,* and therefore *u* and *v,* first increase and then finally decay to zero. The motions are oscillatory. This behavior corresponds to the SEFK portion of the *S* > 0 plots of Fig. 3.

##### The case of nonzero diffussivity

In steady state it may be seen directly, for *N*^{2} ≠ 0, that *w* ≠ 0 so long as the perturbation density has a nonzero, steady-state field, as it does in the previous case for *α* = 0. In fact, the particular solutions are all nonzero, unlike the case for *α* = 0. A steady particular solution can be obtained from (A1), (A2a), (A2b), and (A3′). For example, we obtain the tranform of the particular solution for vertical velocity in the form

*U*_{pn}, *V*_{pn}, and *B*_{pn} may be similarly related to *W*_{pn} or *E*_{n}. The result is that the steady velocity components are nonzero unless *α* = 0. The dispersion relation (A10) for small *A* does not change significantly except that *A* itself becomes *A* = (ℓ^{2} + *m*^{2})(*υ* + *α*/2)/*f.* Thus, all homogeneous solutions decay and the particular solution becomes the finite, *t* → ∞, solution.

### APPENDIX B

#### The Linear Equations Including Isosigma and Horizontal Diffusivity

In this paper, we have studied the errors labeled SEFK and SESK. Aditional errors due to horizontal or isosigma diffusion exist and for completeness are discussed here.

We now amend Eq. (14) to include horizontal diffusivity so that

and

where

Note that the diffusion is based on *b*′; we return to this point later.

We have rerun the numerical model for the standard case, *S* = 4, *r*_{max} = 0.21, and, now, *A*_{H} = *A*_{M} = 2000 m^{2} s^{−1}. For large *t,* we now obtain the SEFK error, a finite difference numerical error, plus additional maximum velocity that we consider an error but that is built into the governing differential equation. It is due to the fact that gradients are reckoned along sigma surfaces so that, in a situation where isopycnals in Cartesian coordinates are horizontal, then there is effective vertical diffusion. This is accordance with the theory at the end of appendix A. We find that the total error has increased from 0.08 to 3.7 cm s^{−1}. Figure 7 presents results in the form of the maximum velocity error as a function of *A*_{H} and the inverse Prandtl number *A*_{H}/*A*_{M}. A lesson is that one might profitably set *A*_{H} < *A*_{M}.

How about rotating the diffusion so that it is horizontal? Thus, instead of (B3a), (B3b), we insert *q*_{x} = *A*_{H} [∂*b*′/∂*x* − *σ*(∂*H*/∂*x*)∂*b*′/∂*σ*] together with a similar formula for *q*_{y}. A calculation reveals that the error reverts to 0.08 cm s^{−1}, the SESK error. The question is, why not reformulate (B3a), (B3b)? The answer is external to this paper. The papers by Mellor and Blumberg (1985) and Mellor and Wang (1996) specifically model bottom boundary layers for which, of course, one needs to explicitly add vertical diffusion to (14). The problems treated in those papers were two-dimensional (*x, σ*) and the models ran without horizontal diffusion. We have returned to the code used in that paper and found that (B3a), (B3b)—reformulated so that diffusion is rotated to the horizontal and *A*_{H} = 2000 m^{2} s^{−1}—interacts unfavorably with the vertical diffusion so that the bottom boundary behavior differs importantly from the realistic boundary layer properties produced by Mellor and Wang with *A*_{H} = 0. The reason is that horizontal diffusion has a large component normal to a sloping bottom that competes deleteriously with the vertical diffusivity; physically, at the bottom, the normal component of momentum and scalar diffusion should be nil. Most of these findings were discussed in the paper by Mellor and Blumberg (1985), in which the point is made that the unhealthy interaction of horizontal diffusion and vertical diffusion near a sloping bottom is not uniquely a sigma coordinate problem but probably endemic to all models. Furthermore, we have seen that bottom flows like the Western Boundary Undercurrent (Ezer and Mellor 1997) and deep water formation (Zavatarelli and Mellor 1993) are well represented by the POM model even with crude resolution near the bottom since downslope Ekman transport is correctly modeled.

For the problem treated in this paper, POM is numerically tolerant of zero diffusion. We have also seen this to be the case in an Atlantic Ocean basin model. We do not know if zero diffusion will work in cases with, say, very crude resolution in which case one can set *A*_{H}/*A*_{M} to a small value; Fig. 7 suggests a value like 0.2. Another way of dealing with nonzero diffusion is to subtract a climatological density (or temperature and salinity) before computing fluxes as in (B3a), (B3b). Furthermore, the climatological density (or temperature and salinity) fields can be made to slowly adjust to the prognostic fields according to *dρ*_{clim}/*dt* = *α*(*ρ* − *ρ*_{clim}), where *α* is the inverse of a decay time. We have executed the problem of Fig. 6 with *A*_{H}/*A*_{M} = 1 and a decay time of 2 days and obtained very nearly the same results as those shown in Fig. 6 for *A*_{H}/*A*_{M} = 0.

## Footnotes

*Corresponding author address:* Prof. George L. Mellor, Program in Atmospheric and Oceanic Sciences, Princeton University, P.O. Box CN710, Sayre Hall, Princeton, NJ 08544-0710.

Email: glm@splash.princeton.edu