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

 
formula

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 δxH 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 δyH. 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 |(σδxH)/Hδσ)| > 1; we do not believe the concept is particularly meaningful. The error in (1) is actually proportional to (δσ)2σ2(δxH/H)2 so that |(σδxH)/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.

Fig. 1.

The seamount geometry. The grid is stretched so that theresoluton is highest at the center (after Beckman and Haidvogel 1993).

Fig. 1.

The seamount geometry. The grid is stretched so that theresoluton is highest at the center (after Beckman and Haidvogel 1993).

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

 
Hx, yH0[1.0 − 0.90e−(x2+y2)/L2]⁠
(2)

A reference stratification was imposed according to

 
ρr−3zρz
(3)

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

 
ρρrρρ−3z
(4)

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

 
formula

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

 
formula

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 rmax = 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.

Fig. 2.

From Beckman and Haidvogel (1993; their Fig. 6a). Ten-day time series of the maximum velocity (m s−1) from the Burger number survey for exponential stratification.

Fig. 2.

From Beckman and Haidvogel (1993; their Fig. 6a). Ten-day time series of the maximum velocity (m s−1) from the Burger number survey for exponential stratification.

Table 1.

Maximum velocity (cm s−1) after 10 days of integration of the unforced exponential stratification. The dash denotes the case where the computation could not be continued to day 10. Decimated table from Beckman and Haidvogel (1993).

Maximum velocity (cm s−1) after 10 days of integration of the unforced exponential stratification. The dash denotes the case where the computation could not be continued to day 10. Decimated table from Beckman and Haidvogel (1993).
Maximum velocity (cm s−1) after 10 days of integration of the unforced exponential stratification. The dash denotes the case where the computation could not be continued to day 10. Decimated table from Beckman and Haidvogel (1993).

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

 
formula

where br(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 −∂br/∂z = N2 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 (Fx, Fy) 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

 
formula

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

 
formula

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

 
formula

where ωwσuH/∂xσvH/∂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η/∂xH 0σ [Hb′/∂xσ′(∂H/∂x)(∂b′/∂σ′)] ′ 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

 
formula

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 rmax = 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, AM = 2000 m2 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 AM 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 rmax. 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.

Fig. 3.

The maximum velocity as a function of time for the case rmax = 0.211. These are fully prognostic results. If run diagnostically (b′ held constant), then Vmax would be identical to the case for S = 0; this would also be true for other values of rmax. The singular behavior of the S = 0 case is explained in appendix A.

Fig. 3.

The maximum velocity as a function of time for the case rmax = 0.211. These are fully prognostic results. If run diagnostically (b′ held constant), then Vmax would be identical to the case for S = 0; this would also be true for other values of rmax. The singular behavior of the S = 0 case is explained in appendix A.

Table 2.

Maximum velocities (cm s−1) after 100 days of integration of the POM model. The viscosity is AM = 2000 m2 s−1 and the diffusivity is nil. The asterisk indicates cases in which the calculation becomes unstable. For the stable cases, the values at 10 days and 100 days are nearly the same.

Maximum velocities (cm s−1) after 100 days of integration of the POM model. The viscosity is AM = 2000 m2 s−1 and the diffusivity is nil. The asterisk indicates cases in which the calculation becomes unstable. For the stable cases, the values at 10 days and 100 days are nearly the same.
Maximum velocities (cm s−1) after 100 days of integration of the POM model. The viscosity is AM = 2000 m2 s−1 and the diffusivity is nil. The asterisk indicates cases in which the calculation becomes unstable. For the stable cases, the values at 10 days and 100 days are nearly the same.

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 rmax 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. For S > 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.

So long as NoS > 0, the velocity error, according to appendix A, should decay completely. We label this decaying error, that associated with Eq. (1), as an SEFK.

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 ∂2p*/∂x*∂y* − ∂2p*/∂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

 
formula

where, from (13a), (13b)

 
formula

(The first i, j location spans ui,j in Fig. 4; the second spans vi,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,

 
formula

The vertical integral of (18), 0−1 Γp dσ, for S = 4, rmax = 0.211, and t = 0 is plotted in Fig. 5. This error is labeled an SESK.

Fig. 4.

The stencil for a C grid. (a) The solid box denotes the principal grid cell. (b) The dashed box denotes the offset cell used to finite difference vorticity. The variables H and b are collocated with p.

Fig. 4.

The stencil for a C grid. (a) The solid box denotes the principal grid cell. (b) The dashed box denotes the offset cell used to finite difference vorticity. The variables H and b are collocated with p.

Fig. 5.

SESK, the vertical integral of Eq. (18). The central region of the computational domain is shown. Here S = 4, rmax = 0.211, and t = 0. Solid lines are positive; dashed lines are negative. The straight lines are zero contours. Here CI = 5 × 10−13 m2 s−2.

Fig. 5.

SESK, the vertical integral of Eq. (18). The central region of the computational domain is shown. Here S = 4, rmax = 0.211, and t = 0. Solid lines are positive; dashed lines are negative. The straight lines are zero contours. Here CI = 5 × 10−13 m2 s−2.

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 rmax 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 m2 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 AM.

Table 3.

Maximum velocities (cm s−1) after 10 days of integration of the POM model. The viscosity is AM = 200 m2 s−1 and the diffusivity is nil.

Maximum velocities (cm s−1) after 10 days of integration of the POM model. The viscosity is AM = 200 m2 s−1 and the diffusivity is nil.
Maximum velocities (cm s−1) after 10 days of integration of the POM model. The viscosity is AM = 200 m2 s−1 and the diffusivity is nil.

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 Hi,jHi,j−1 = (δyH)o + δx(δyH)/2, etc.; then

 
formula

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

 
formula

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 δxyH = (∂2H/∂xy) δx δy, etc., and from (20), Γp = F(r) cosθ sinθ (cos2θ − sin2θ). 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 δxyH 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δx4 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 DH + η and we define the operator

 
formula

then the equations of motion are

 
formula

Note that, consistent with the previous linear system, ρ′ = ρρr is used in the baroclinic integrals of (23a), (23b). Here, Fx and Fy 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 AH = 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 m2 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 AH = 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.

Fig. 6.

The full nonlinear model with a throughflow of 0.20 m s−1 applied uniformly at x = 0. The diffusivity is nil. Shown are the differences between the densities at day 2 and the initial densities at a depth of 2000 m. Contour interval = 20 × 10−4 kg m−3. The arrows represent particle trajectories for 2-day flights based on instantaneous Eulerian velocities. The model ran successfully for 30 days with no sign of numerical instability.

Fig. 6.

The full nonlinear model with a throughflow of 0.20 m s−1 applied uniformly at x = 0. The diffusivity is nil. Shown are the differences between the densities at day 2 and the initial densities at a depth of 2000 m. Contour interval = 20 × 10−4 kg m−3. The arrows represent particle trajectories for 2-day flights based on instantaneous Eulerian velocities. The model ran successfully for 30 days with no sign of numerical instability.

Fig. 7.

The maximum velocity error as a function of the horizontal viscosity AM (m2 s−1) and the ratio AH/AM. Here S = 4 and t = 10 days. Computations were made at the circled points.

Fig. 7.

The maximum velocity error as a function of the horizontal viscosity AM (m2 s−1) and the ratio AH/AM. Here S = 4 and t = 10 days. Computations were made at the circled points.

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

REFERENCES
Baringer, M. O., and J. F. Price, 1997a: Mixing and spreading of the Mediterranean outflow. J. Phys. Oceanogr., 27, 1654–1677.
——, and ——, 1997b: Momentum and energy balance of the Mediterranean outflow. J. Phys. Oceanogr., 27, 1678–1692.
Beckman, A., and D. Haidvogel, 1993: Numerical simulation of flow around a tall isolated seamount. J. Phys. Oceanogr., 23, 1736–1753.
Blumberg, A. F., and G. L. Mellor, 1987: A description of a three-dimensional coastal ocean circulation model. Three-Dimensional Coastal Ocean Models, N. Heaps, Ed., Coastal Estuarine Science, Vol. 4, Amer. Geophys. Union, 1–16.
Ezer, T., and G. L. Mellor, 1997: Simulations of the Atlantic Ocean with a free surface sigma coordinate model. J. Geophys. Res., 102, 15 647–15 657.
Gary, J. M., 1973: Estimate of truncation error in transformed coordinate, primitive equation atmospheric models. J. Atmos. Sci., 30, 223–233.
Haidvogel, D. B., J. L. Wilkins, and R. E. Young, 1991: A semi-spectral primitive equation ocean circulation model using vertical and orthogonal curvilinear horizontal coordinates. J. Comput. Phys., 94, 151–185.
Haney R. L., 1991: On the pressure gradient force over steep topography in sigma coordinate models. J. Phys. Oceanogr., 21, 610–619.
Hedstom, K. S., 1990: User’s manual for a semi-spectral primitive equation regional ocean-circulation model, version 3.0. Institute for Naval Oceanography Tech. Note FY90-2, 90 pp. [Available from IMCS, 71 Dudley Rd., New Brunswick, NJ 08901.]
.
Jungclaus, J., and J. O. Backhaus, 1994: Application of a transient reduced gravity plume model to the Denmark Strait overflow. J. Geophys. Res., 99, 12 375–12 396.
McCalpin, J. D., 1994: A comparison of second-order and fourth-order pressure gradient algorithms in a σ-coordinate ocean model. Int. J. Numer. Methods Fluids, 18, 361–383.
Mellor, G. L., 1996: User’s guide for a three-dimensional, primitive equation, numerical ocean model. Princeton University Rep. Princeton University, Princeton, NJ, 40 pp. [Available from Princeton University, Princeton, NJ 08540.]
.
——, and A. F. Blumberg, 1985: Modeling vertical and horizontal diffusivities with the sigma coordinate system. Mon. Wea. Rev., 113, 1379–1383.
——, and X. H. Wang, 1996: Pressure compensation and the bottom boundary layer. J. Phys. Oceanogr., 26, 2214–2222.
——, T. Ezer, and L.-Y. Oey, 1994: The pressure gradient conundrum of sigma coordinate ocean models. J. Atmos. Oceanic Technol., 11, 1126–1134.
Rousseau, D., and H. L. Pham, 1971: Premiers resultat d’un modele de prevision numerique à courte escheance sur l’Europe. Meteorologie, 20, 1–12.
Zavatarelli, M., and G. L. Mellor, 1995: A numerical study of the Mediterranean Sea circulation. J. Phys. Oceanogr., 25, 1384–1414.

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

 
formula

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) = bo(z), ∂bo/∂x = ∂bo/∂y = ∂bo/∂t = 0 so that bo(z) can be subtracted from b′ and added back after the time-dependent solution has been obtained. Thus the initial conditions are

 
formula

For S = N2 = 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, up, vp, wp, bp = 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

 
formula

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

 
upvpwpbpEnei(lx+my)
(A6)

where En(z) will subsequently be related to ɛ and Wn(z).

Essentially the problem is to find bh so that b′ = bp + bh = 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

 
formula

where υ ≡ (ℓ2 + m2)υ. Substitution of the above into Eq. (A1) yields

 
formula

where the dependencies σ = σn(ℓ, m), are understood. Subject to the boundary conditions, Wn(0) = Wn(−H) = 0, Wn(z), and λ2n may be determined as an eigenvalue problem (for N2 = constant = N2o, Wn ∝ sinλnz, and λn = /H for n = 1, 2, 3, . . .).

Now the dispersion relation (A9) may be written

 
ω3iAω2MA2ωiAM
(A10)

where

 
formula

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

 
formula

The imaginary parts of σ = fω, when inserted into (A5a)–(A5d), yields a zero or negative exponent for eiσt.

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

 
formula

where we have set a2 = −a1 = i/2. Since the homogeneous solutions may be multiplied by any constant, to satisfy (A4d), we have

 
formula

and En(z) = N2(z)Wn(z). If we invert

 
formula

to determine An(ℓ, m), then b′ and w can be determined as well as u and v from (A7a), (A7b). In particular,

 
formula

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 anticipation of appendix B, we add diffusion to (A3) so that

 
formula

In steady state it may be seen directly, for N2 ≠ 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

 
formula

Upn, Vpn, and Bpn may be similarly related to Wpn or En. 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 + m2)(υ + α/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

 
formula

and

 
formula

where

 
formula

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, rmax = 0.21, and, now, AH = AM = 2000 m2 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 AH and the inverse Prandtl number AH/AM. A lesson is that one might profitably set AH < AM.

How about rotating the diffusion so that it is horizontal? Thus, instead of (B3a), (B3b), we insert qx = AH [∂b′/∂xσ(∂H/∂x)∂b′/∂σ] together with a similar formula for qy. 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 AH = 2000 m2 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 AH = 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 AH/AM 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 clim/dt = α(ρρclim), where α is the inverse of a decay time. We have executed the problem of Fig. 6 with AH/AM = 1 and a decay time of 2 days and obtained very nearly the same results as those shown in Fig. 6 for AH/AM = 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.