## Abstract

How to reduce the horizontal pressure gradient error is a key issue of using σ-coordinate ocean models, especially of using primitive equation models for coastal regions. The error is caused by the splitting of the horizontal pressure gradient term into two parts and the subsequent incomplete cancellation of the truncation errors of those parts. Due to the fact that the higher the order of the difference scheme, the less the truncation error and the more complicated the computation, a sixth-order difference scheme for the σ-coordinate ocean models is proposed in order to reduce error without increasing complexity of the computation. After the analytical error estimation, the Semi-spectral Primitive Equation Model is used to demonstrate the benefit of using this scheme. The stability and accuracy are compared with those of the second-order and fourth-order schemes in a series of calculations of unforced flow in the vicinity of an isolated seamount. The sixth-order scheme is shown to have error reductions by factors of 5 compared to the fourth-order difference scheme and by factors of 50 compared to the second-order difference scheme over a wide range of parameter space as well as a great parametric domain of numerical stability.

## 1. Introduction

In shallow-water prediction models the effects of bottom topography must be taken into account. This can be done by using a terrain-following σ-coordinate system, where the water column is divided into the same number of grid cells independent of depth. Let (x∗, y∗, z) denote Cartesian coordinates and (x, y, σ) sigma coordinates. In most sigma coordinate ocean models the relationship between the two coordinate systems is

where z and σ increase vertically upward such that z = 0, σ = 1 at the surface and z = −H and σ = −1 at the bottom; H = H(x, y) is the bottom topography.

We restrict our attention to two dimensions (x, σ) for simplicity. Horizontal gradients in the z- and σ-coordinates are related by

Most ocean models are hydrostatic balanced; that is,

where g is the acceleration due to gravity and ρ is the density of the fluid. Substitution of (3) into (2) and use of

yield

which is used for computing the horizontal pressure gradient.

A problem has long been recognized in computing the horizontal pressure gradient in the σ-coordinate system (e.g., Gary 1973; Haney 1991; Mellor et al. 1994; McCalpin 1994) the horizontal pressure gradient becomes a difference between two terms, which leads to a large truncation error at steep topography. We refer the reader to a comprehensive review by Haney (1991). Mellor et al. (1994) obtained an error estimation scheme and showed that the pressure gradient error is advectively eliminated after a long time integration. McCalpin (1994) employed a fourth-order difference scheme to eliminate the pressure gradient error by a factor of 10–20. Following McCalpin’s path, we propose a sixth-order difference scheme for a further reduction of pressure gradient error.

## 2. Second-, fourth-, and sixth-order schemes

In ocean models, the most common spatial difference scheme used to approximate the first derivative is a two point, centered, second-order staggered scheme such as the C-grid scheme (Hadivogel et al. 1991) and the B-grid scheme (Robinson 1993):

where Δ is the grid spacing. The truncation error of this scheme is O2).

For the fourth-order C-grid scheme,

the truncation error is O4). Comparing (6) with (5), the second-order difference, the error ratio between the fourth-order and second-order schemes is estimated by

Since the truncation error decreases with the increase of the order of the difference scheme, it might be beneficial to use an even higher-order difference scheme. Thus, we propose a sixth-order difference scheme

to compute the horizontal pressure gradient. For the two points nearest to the boundary, we use the lower-order schemes (taking one-side boundary as an example),

Comparing (8) with (6), the error ratio between the sixth-order and fourth-order schemes is estimated by

## 3. Seamount test case

### a. Model description

Suppose a seamount located inside a periodic f-plane (f0 = 10−4 s−1) channel with two solid, free-slip boundaries along constant y. Unforced flow over a seamount in the presence of resting, level isopycnals is an ideal test case for the assessment of pressure gradient errors in simulating stratified flow over topography. The flow is assumed to be reentrant (periodic) in the alongchannel coordinate (i.e., x axis). We use this seamount case of the Semi-spectral Primitive Equation Model (SPEM) version 3.9 to test the new difference scheme. The reader is referred to the original reference (Haidvogel et al. 1991) and the SPEM 3.9 User’s Manuel (Hedstrom 1995) for datailed information. In the horizontal directions the model uses the C-grid and the second-order finite difference discretization except for the horizontal pressure gradient for which the user has a choice of either second-order or fourth-order difference discretization (McCalpin 1994). In the vertical direction the model uses a boundary-fitted σ-coordinate system. The discretization is by spectral collocation using Chebyshev polynomials. Our model configuration is similar to that of Beckmann and Haidvogel (1993) and McCalpin (1994). The time step and grid size used here are

Δt = 675 s,  Δx = Δy = 5 km.

### b. Topography

The domain is a periodic channel, 320 km long and 320 km wide. The channel walls are solid (no normal flow) with free-slip viscous boundary conditions. The channel has a far-field depth hmax and in the center includes an isolated Gaussian-shape seamount with a width L and an amplitude hs (Fig. 1),

where (x0, y0) are the longitude and latitude of the seamount center. The far-field depth (hmax) is fixed as 5000 m. But the seamount amplitude (hs) changes from 500 to 4500 m, and the lateral scale of the seamount (L) varies from 6 to 40 km for the study.

### c. Density field

Suppose that the background fluid is at rest and with an exponentially stratified mean density

where z is the vertical coordinate, and Hρ = 1000 m. Here a constant density, 1000 kg m−3, has been subtracted for error reduction. When a density disturbance ρi = ρ̂ exp(z/Hρ) was introduced initially,

ρi = ρ̄(z) + ρi,
(12)

the fluid was slightly less stably stratified than the reference field for each computation. Here ρ̂ is called the density anomaly. The larger the ρ̂, the less stable the fluid will be. In this study, ρ̂ varies from 0.1 to 1 kg m−3.

Following Beckmann and Haidvogel (1993) and McCalpin (1994), we subtract the mean density field ρ̄(z) before integrating the density field to obtain pressure from the hydrostatic equation.

### d. Lateral viscosity

Ideally, the new difference scheme should be tested with no lateral diffusion of density. This is due to the fact that the density diffusion along σ surfaces generates horizontal gradients wherever the σ surfaces are not flat, and then produces horizontal pressure gradients that drive currents in much the same way as the pressure gradient errors (McCalpin 1994). Unfortunately, the absence of the horizontal diffusion keeps the small-scale pressure disturbances generated by topographic-scale density advection and the induced small-scale velocity fields, which in turn cause the computational instability problem. Thus, some lateral viscosity on σ surfaces is required in the momentum equations to maintain stability. A constant coefficient (AH) biharmonic formulation is used here for the lateral viscosity, which varies from 1010 to 108 m4 s−1 in this study.

## 4. Standard case

### a. Model parameters

In this section, we assess the sixth-order scheme (8) by the cumulative effects of the pressure gradient errors. A variety of configuration were used for the test. These will be described in next section. At first, we set up a standard test case to compare errors among three different (second, fourth, and sixth) schemes:

The SPEM3.9 standard was integrated for 20 days for the standard test case using ordinary and compact fourth- and sixth-order difference schemes for computing horizontal pressure gradient. Figure 2 displays errors in the streamfunction and velocity fields after performing 5 days of integration using the sixth-order scheme. The mass transport streamfunction has a large-scale eight-lobe pattern centered on the seamount. This symmetric structure can be found in all the fields. After 5 days of integration, the model generates spurious currents of O(0.001 cm s−1).

Fig. 2.

Instantaneous error pattern after 5 days of integration for the unforced (exponential stratification) experiment with sixth-order difference scheme: (a) u, (b) υ, (c) w, and (d) mass transport streamfunction. Here u, υ, w are evaluated at the slice (facing upchannel) through the center of the seamount.

Fig. 2.

Instantaneous error pattern after 5 days of integration for the unforced (exponential stratification) experiment with sixth-order difference scheme: (a) u, (b) υ, (c) w, and (d) mass transport streamfunction. Here u, υ, w are evaluated at the slice (facing upchannel) through the center of the seamount.

### b. Temporal variations of peak error velocity

Owing to a very large number of calculations performed, we discuss the results exclusively in terms of the maximum absolute value of the spurious velocity (called peak error velocity) generated by the pressure gradient errors. Figure 3 shows the time evolution of the peak error velocity for the first 20 days of integration with the second-, fourth-, and sixth-order schemes. The peak error velocity fluctuates rapidly during the first few days integration. After the 5 days of integration, the peak error velocity shows the decaying inertial oscillation superimposed on mean values: 0.4 cm s−1 (the second order), 0.042 cm s−1 (the fourth order), and 0.01 cm s−1 (the sixth order). The steady approach of the peak error velocities to these values for the three schemes indicates the stability of the computation. Furthermore, the time-mean peak error velocity for the second-order case is roughly 40–50 times that of the sixth-order case.

### c. Performance of the second-, fourth-, and sixth-order schemes

Feasibility of using the sixth-order scheme is twofold: 1) drastic error reduction and 2) no drastic CPU time increase. We use four quantities to compare the errors: the maximum peak error velocity during the 10-day integration (V1), the averaged peak error velocity during the 5–10 day integration (V2), the maximum peak error pressure gradient during the 10-day integration (PG1), and averaged peak error pressure gradient during the 5–10 day integration (PG2). Here V1 (or PG1) indicates the maximum error during the geostrophic adjustment stage, and V2 (or PG2) shows the mean error after the adjustment. Usually, V1 > V2, and PG1 > PG2. The performance of each scheme is listed in Table 1. Both V1 and V2 reduce their magnitude by factors around 10 from the second-order scheme to the fourth-order scheme (V1 from 7.4 × 10−1 to 8.17 × 10−2 cm s−1, V2 from 4.61 × 10−1 to 4.43 × 10−2 cm s−1), and by factors around 5 from the fourth-order scheme to the sixth-order scheme (V1 from 8.17 × 10−2 to 1.61 × 10−2 cm s−1, V2 from 4.43 × 10−2 to 9.75 × 10−3 cm s−1). The pressure gradient error reduces even more when the high-order schemes are used: PG1 reduces from 9.7 × 10−1 to 8.3 × 10−2 N m−3 from the second-order scheme to the fourth-order scheme (a factor of 11.7) and from 8.3 × 10−2 to 1.27 × 10−2 N m−3 from the fourth-order scheme to the sixth-order scheme (a factor of 6.5). PG2 reduces from 7.1 × 10−3 to 4.96 × 10−4 N m−3 from the second-order scheme to the fourth-order scheme (a factor of 14.3) and from 4.96 × 10−4 to 8.75 × 10−5 N m−3 from the fourth-order scheme to the sixth-order scheme (a factor of 5.7). The CPU time increases 10% from the second-order to fourth-order and sixth-order scheme. Notice that there is no CPU increase from the fourth-order to the sixth-order scheme, which shows the great potential of using the sixth-order scheme.

Table 1.

Ratios of CPU time and various errors between the second-order and the fourth (sixth)-order schemes. ## 5. Sensitivity studies

We performed various sensitivity study to show the dependence of model results on the four parameters listed in (13). For each sensitivity study, we vary only one parameter and then integrate the model for 10 days using the second-, fourth-, and sixth-order schemes, respectively.

### a. Seamount width

Seamount widths (L) of 40, 20, 15, 12, and 10 km were used. A summary of V1 and V2 after 10 days of integration for each case is presented in Fig. 4. The seamount with L larger than 10 km can be resolved by the grid spacing used here (5 km). The figure clearly shows a great improvement of using high-order schemes. This improvement enhances as L increases. At L = 10 km, V1 reduces by a factor of 1.45 from the second-order scheme to the fourth-order scheme and a factor of 1.18 from the fourth-order scheme to the sixth-order scheme. At L = 40 km, V1 reduces by a factor of 9.06 from the second-order scheme to the fourth-order scheme and a factor of 5.08 from the fourth-order scheme to the sixth-order scheme. Similar reduction of V2 was found (Fig. 4) from the second-order scheme to the fourth-order scheme, and from the fourth-order scheme to the sixth-order scheme.

### b. Seamount amplitude

Seamount amplitudes (hs) of 4500, 4000, 3500, 3000, 2000, 1000, and 500 m were used. A summary of V1 and V2 after 10 days of integration for each case is presented in Fig. 5. The figure clearly shows a great improvement using high-order schemes. This improvement enhances as hs increases. At hs = 500 m (i.e., a small seamount with a slope of 0.0125), V1 reduces by a factor of 2.55 from the second-order scheme to the fourth-order scheme and a factor of 1.10 from the fourth-order scheme to the sixth-order scheme. At hs = 4500 m (i.e., a large seamount with a slope of 0.1125), V1 reduces by a factor of 9.06 from the second-order scheme to the fourth-order scheme and of 5.08 from the fourth-order scheme to the sixth-order scheme. Similar reduction of V2 was found (Fig. 5) from the second-order scheme to the fourth-order scheme, and from the fourth-order scheme to the sixth-order scheme.

### c. Density anomaly

Density anomalies (ρ̂) of 1, 0.8, 0.6, 0.4, 0.2, and 0.1 kg m−3 were used. A summary of V1 and V2 after 10 days of integration for each case is presented in Fig. 6. The figure clearly shows a great improvement using high-order schemes. This improvement is not sensitive to the change of ρ̂. At ρ̂ = 0.1 kg m−3 (i.e., a small density anomaly), V1 reduces by a factor of 9.05 from the second-order scheme to the fourth-order scheme and a factor of 5.08 from the fourth-order scheme to the sixth-order scheme. At ρ̂ = 1 kg m−3 (i.e., a large density anomaly), V1 reduces by a factor of 9.10 from the second-order scheme to the fourth-order scheme and of 5.07 from the fourth-order scheme to the sixth-order scheme. Similar reduction of V2 was found (Fig. 6) from the second-order scheme to the fourth-order scheme, and from the fourth-order scheme to the sixth-order scheme.

### d. Horizontal viscosity

Horizontal viscosities (AH) of 1010, 109, and 108 m4 s−1 were used. A summary of V1 and V2 after 10 days of integration for each case is presented in Tables 2 and 3, which clearly shows a great improvement using high-order schemes. This improvement is not very sensitive to the change of AH. At AH = 1010 m4 s−1, V1 reduces by a factor of 9.06 from the second-order scheme to the fourth-order scheme and a factor of 5.08 from the fourth-order scheme to the sixth-order scheme. At AH = 108 m4 s−1, V1 reduces by a factor of 9.19 from the second-order scheme to the fourth-order scheme and of 5.10 from the fourth-order scheme to the sixth-order scheme. Similar reduction of V2 was found (Table 3) from the second-order scheme to the fourth-order scheme, and from the fourth-order scheme to the sixth-order scheme.

Table 2.

Sensitivity of V1 (cm s−1) on AH for the second-, fourth-, and sixth-order schemes. Table 3.

Sensitivity of V2 (cm s−1) on AH for the second-, fourth-, and sixth-order schemes. ## 6. Conclusions

1) The σ-coordinate, pressure gradient error depends on the choice of difference schemes. By choosing an optimal scheme, we may reduce the error a great deal without increasing the horizontal resolution. Analytical analysis shows that the truncation error of the fourth-order scheme may be 1–2 orders of magnitude smaller than the second-order scheme, and the truncation error of the sixth-order scheme may be 1–2 orders of magnitude smaller than the fourth-order scheme.

2) The SPEM Version 3.9 is used to demonstrate the benefit of using the sixth-order scheme. A series of calculations of unforced flow in the vicinity of an isolated seamount are performed. The results show that the sixth-order scheme has error reductions by factors of 5 compared to the fourth-order difference scheme, and by factors of 50 compared to the second-order difference scheme over a wide range of parameter space as well as a great parametric domain of numerical stability.

3) Using the sixth-order scheme does not require much more CPU time. Taking SPEM3.9 as an example, the CPU time for the sixth-order scheme is almost the same as for the fourth-order scheme, and 10% more than for the second-order scheme.

4) Since the fourth-order different scheme has error reductions by factors of 10 compared to the second-order difference scheme, there is no real advantage to going to a higher-order scheme if the bottom topography is not too complicated. The need for more accuracy will go up with increasingly complex bottom topography on small scales, so one might expect that future demand for accuracy will increase as models strive for more realism.

## Acknowledgments

The authors wish to thank Dale Haidvogel and Kate Hedstrom of the Rutgers University for most kindly proving us with a copy of the SPEM code. This work was funded by the Office of Naval Research NOMP Program, the Naval Oceanographic Office, and the Naval Postgraduate School.

## REFERENCES

REFERENCES
Beckmann, A., and D. B. Haidvogel, 1993: Numerical simulation of flow around a tall isolated seamount. Part I: Problem formulation and model accuracy. J. Phys. Oceanogr., 23, 1736–1753.
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. Wikin, and R. Young, 1991: A semi-spectral primitive equation model using vertical sigma and orthogonal curvilinear coordinates. J. Comput. Phys., 94, 151–185.
Haney, R. L., 1991: On the pressure gradient force over steep topography in sigma coordinate ocean models. J. Phys. Oceanogr., 21, 610–619.
Hedstrom, K., 1995: User’s Manual for a Semi-Spectral Primitive Equation Ocean Circulation Model Version 3.9. Rutgers University, New Jersey, 131 pp. [Available from Institute of Marine and Coastal Sciences, Rutgers University, New Brunswick, NJ 08903.]
.
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., T. Ezer, and L.-Y. Oey, 1994: The pressure gradient conundrum of sigma coordinate ocean models. J. Atmos. Oceanic Technol., 11, 1126–1134.
Robinson, A. R., 1993: Physical processes, field estimation and interdisciplinary ocean modeling. Harvard Open Ocean Model Reports, Harvard University, 71 pp. [Available from Harvard University, 29 Oxford St., Cambridge, MA 02138.]
.

## Footnotes

Corresponding author address: Dr. Peter C. Chu, Department of Oceanography, Naval Postgraduate School, Code OC/CU, Monterey, CA 93943-5000.