How to reduce the computational error is a key issue in numerical modeling and simulation. The higher the order of the difference scheme, the less the truncation error and the more complicated the computation. For compromise, a simple, three-point accuracy progressive (AP) finite-difference scheme for numerical calculation is proposed. The major features of the AP scheme are three-point, high-order accuracy, and accuracy progressive. The lower-order scheme acts as a “source” term in the higher-order scheme. This treatment keeps three-point schemes with high accuracy. The analytical error estimation shows the sixth-order accuracy that the AP scheme can reach. The Fourier analysis of errors indicates the accuracy improvement from lower-order to higher-order AP schemes. The Princeton Ocean Model (POM) implemented for the Japan/East Sea (JES) is used to evaluate the AP scheme. Consider a horizontally homogeneous and stably stratified JES with realistic topography. Without any forcing, initially motionless ocean will keep motionless forever; that is to say, there is a known solution (V = 0). Any nonzero model velocity can be treated as an error. The stability and accuracy are compared with those of the second-order scheme in a series of calculations of unforced flow in the JES. The three-point sixth-order AP scheme is shown to have error reductions by factors of 10–20 compared to the second-order difference scheme. Due to their three-point grid structure, the AP schemes can be easily applied to current ocean and atmospheric models.
Most ocean models are hydrostatically balanced. Such a balance sets up a lowest limit, (Δx, Δy)L, for the grid spacings. If the model resolution is finer than that limit, nonhydrostatic ocean models should be used. If the hydrostatic balance is kept unchanged, the grid spacing cannot be artificially small and the reduction of discretization errors within the grid spacing limitation becomes important. A logical approach is to use high-order difference schemes.
High-order schemes can be applied to the σ-coordinate ocean models for error reduction. Let (x∗, y∗, z) denote Cartesian coordinates and (x, y, σ) σ coordinates. The conventional relationships between the two coordinates are
where H is the water depth, and z and σ increase vertically upward such that z = σ = 0 at the surface and σ = −1 and z = −H at the bottom. The horizontal pressure gradient becomes a difference between two terms using the σ coordinates,
which leads to a so-called hydrostatic inconsistency with a large truncation error at a steep topography (Gary 1973; Haney 1991). Mellor et al. (1994) argue that if the pressure gradient force is computed as a vertical integral of discretized Jacobian of density and z coordinate, the problem of hydrostatic inconsistency can be avoided and the spurious residual pressure gradient is drastically reduced. Furthermore, Song (1998) shows that a modified Jacobian formulation on the discredited pressure gradient term with a specially designed nonlinear weighting procedure drastically reduces the residue errors.
Thus, we have two approaches to reduce errors of the σ-coordinate ocean models: 1) bringing certain symmetries of the continuous forms into the discrete level to ensure cancellations of these terms (e.g., Mellor et al. 1994; Song 1998; Song and Wright 1998), and 2) increasing numerical accuracy (e.g., McCalpin 1994; Chu and Fan 1997, 1998, 1999, 2000). Note that the use of high-order schemes does not conflict with the use of symmetric form. We may combine the two approaches (symmetric forms with high-order schemes) to further diminish the computational errors.
Usually, the higher the accuracy of the scheme, the more the grid points are needed. For example, a C-grid fourth-order ordinary scheme needs four grid points (McCalpin 1994) and a C-grid sixth-order ordinary scheme needs six grid points (Chu and Fan 1997) for computing the first derivative. However, the current ocean models generally use three-point second-order schemes. Conversion of the existing ocean models from second- to high-order (fourth or sixth order) schemes becomes easier if the high-order scheme preserves the three-point grid.
Concurrently, Hirsh (1975) and Kreiss (1975) have proposed Hermitian compact techniques using fewer grid points (three instead of five) at each grid point to solve partial differential equations (PDEs). Later on, as pointed out by Adam (1977), the truncation errors are usually four to six times smaller than the same order noncompact schemes. Since then, much work has been done in developing compact schemes for various applications, such as an implicit compact fourth-order algorithm (Navon and Riphagen 1979), a fourth-order compact difference scheme for nonuniform grids (Goedheer and Potters 1985), and fourth-order and sixth-order compact difference schemes for the stagged grid (Chang and Shirer 1985). Recently, Chu and Fan (1998, 1999, 2000) proposed a three-point sixth-order combined compact difference (CCD) scheme for uniform, nonuniform, and staggered grids. This scheme follows earlier work on use of second derivatives in compact differences (such as Rubin and Khosla 1977),
which is referred as the Hermitian formula. Here, f is a dependent variable. The CCD schemes use implicit solvers to calculate various differences (twin-tridiagonal solver) and to solve PDEs (triple-tridiagonal solver). Such a complicated procedure may limit the application of the CCD schemes to the ocean models.
In this study, we propose a simple, three-point accuracy progressive (AP) finite-difference scheme for numerical calculation. The major features of the AP scheme are three point, high-order accuracy, and accuracy progressive. The lower-order scheme acts as a “source” term in the higher-order scheme. This treatment keeps in three-point grid with high accuracy and is explicit in time integration.
2. Progress of accuracy
Let the dependent variable f(x) be defined on the interval, 0 ≤ x ≤ L. Use a uniform grid, 0 = x1 < x2 < x3 < · · · < xN < xN+1 = L with a spacing h = xi+1 − xi = L/N. Let the dependent variable f(x) at any grid point xi and two neighboring points xi−1 and xi+1 be given by fi, fi−1, and fi+1 and let its derivatives at the two neighboring points xi−1 and xi+1 be given by f′i−1, f″i−1, … , f(k)i−1, and f′i+1, f″i+1, … , f(k)i+1.
A new approach is to relate higher-order schemes to lower-order schemes. For example, in the standard Pade (fourth order) scheme,
the second-order differences of f′ and f″ at the grid point xi,
are the source terms of the standard Pade (fourth order) scheme in the right-hand sides of (1) and (2).
Following this path, we propose a series of three-point AP finite-difference schemes with f′, f″, … , f(k) computed by a lower-order scheme as source terms in higher-order schemes. They have the following features: three-point, accuracy progressive, and only one tridiagonal solver needed. The second-order AP scheme is the same as the central difference scheme, and the fourth-order AP scheme is the same as the standard Pade scheme.
3. AP schemes
a. Standard grid
Let (Δf/Δx, Δ2f/Δx2), (δf/δx, δ2f/δx2), and (Df/Dx, D2f/Dx2) be the second-order, fourth-order, and sixth-order differences for (f′, f″), respectively. The accuracy progress from the second-order scheme into the fourth-order AP scheme is given by
Taylor series expansion shows that if α1 = 1/4, β1 = 3/2, α2 = 1/10, and β2 = 6/5, the schemes (5) and (6) have the fourth-order accuracy.
The accuracy progress from the fourth-order AP scheme into the sixth-order AP scheme can be defined by
When α1 = 7/16, β1 = 15/8, γ1 = 1/16, α2 = −1/8, β2 = 3, and γ2 = −9/4, the schemes (7) and (8) have the sixth-order accuracy (see appendix).
b. Staggered C grid
The second-order AP scheme on a staggered C grid is given by
The fourth-order scheme is given by
Taylor series expansion shows that if α1 = 1/22, β1 = 12/11, α2 = 1/10, and β2 = 6/5, the schemes (11) and (12) have the fourth-order accuracy.
The sixth-order AP scheme on a staggered C grid can be defined by
When α1 = −7/254, β1 = 120/127, γ1 = −17/254, α2 = −5/94, β2 = −102/47, and γ2 = 144/47, the schemes (13) and (14) have the sixth-order accuracy.
4. Fourier analysis of errors
Fourier analysis of errors is commonly used to evaluate various difference schemes, described extensively in Orzag (1971), Kreiss and Oliger (1972), Lele (1992), and Chu and Fan (1998). As pointed out by Lele (1992), Fourier analysis provides an effective way to quantify the resolution characteristics of differencing approximations.
For the purpose of Fourier analysis, the dependent variable f(x) is assumed to be periodic over the domain [0, L] of the independent variable, that is, f1 = fN+1 and h = L/N. The dependent variable may be decomposed into Fourier series,
where i = −1. It is convenient to introduce a scaled wavenumber w = 2πkh/L = 2πk/N, and a scaled coordinate s = x/h. The Fourier modes in terms of these are simply exp(iws). The exact first-order and second-order derivatives of (15) generate a function with exact Fourier coefficients:
However, the Fourier coefficients of the derivatives obtained from the differencing scheme might not be the same as the exact Fourier coefficients, that is,
where w′ = w′(w) and w″ = w″(w) are the modified wavenumber (both real numbers) for the first-order and second-order differencing. The smaller the difference between the exact and modified wavenumbers, the better the difference scheme.
According to Lele (1992), if the modified wavenumbers of the current generalized difference schemes as
the modified wavenumbers will be
For the second-order AP scheme (central difference scheme), we have
for the fourth-order AP scheme (standard Pade' scheme), we have
and for the sixth-order AP schemes (7) and (8), we have
Among various difference schemes, the modified wavenumbers of the first-order differencing w′ (Fig. 1a) and of the second-order differencing w″ (Fig. 1b) of the AP scheme are closest to the exact wavenumber w. Besides, we computed the errors of the modified wavenumbers for the first-order derivative (Fig. 1c)
and for the second-order derivative (Fig. 1d)
respectively. We clearly see the accuracy increase of the AP schemes from the second-order to the fourth-order, and from the fourth-order to the sixth-order. Furthermore, the performance of the sixth-order AP scheme is very close to the CCD scheme.
In multidimensional problems the phase error of first-order differencing scheme appear in the form of anisotropy:
Figure 1e shows polar plots of phase speed anisotropy of various schemes for first derivative approximations. The phase speeds for wavenumbers (magnitude) w/π = 1/50, 5/50, … , 45/50, 50/50 are plotted. Here, we also see that the sixth-order AP scheme shows improvement.
5. Boundary treatment
Due to the progressive accuracy, the AP schemes need special treatment near the boundaries x1 and xN+1. For the standard AP schemes, the derivatives can be computed at (x2, x3, … , xN) using the second-order scheme, at (x3, x4, … , xN−1) using the fourth-order scheme, and at (x4, x5, … , xN−2) using the sixth-order scheme (Fig. 2a). Thus, we use the second-order scheme for the grid points x2 and xN (next to the boundaries), the fourth-order scheme for the grid points x3 and xN−1, and the sixth-order scheme for the grid points (x4, x5, … , xN−2). For the sixth-order AP scheme on the staggered C grid, we have the similar treatment (Fig. 2b).
6. Evaluation using the Princeton Ocean Model
a. Known solution
Accuracy of any numerical scheme should be evaluated by either an analytical or known solution. For a realistic ocean model, the analytical solution is hard to find. Consider a horizontally homogeneous and stably stratified ocean with realistic topography. Without forcing, initially motionless ocean will keep motionless forever; that is to say, we have a known solution (V = 0). Any nonzero model velocity can be treated as an error.
b. JES geography
The Japan Sea, known as the East Sea in Korea, has a steep bottom topography that makes it a unique semi-enclosed ocean basin overlaid by a pronounced monsoon surface wind. The JES covers an area of 106 km2. It has a maximum depth in excess of 3700 m, and is isolated from open oceans except for small (narrow and shallow) straits, which connect the JES with the North Pacific through the Tsushima/Korean and Tsugaru Straits and with the Okhotsk Sea through the Soya and Tatar Straits. The smoothed JES bathymetry is shown in Fig. 3.
c. Model description
The POM is a time-dependent, σ-coordinate, primitive equation circulation model on a three-dimensional grid that includes realistic topography and a free surface (Blumburg and Mellor 1987). The model contains 94 × 100 × 15 horizontally fixed grid points. The horizontal spacing of 10′ latitude and longitude (approximately 11.54–15.18 km in the zonal direction and 18.53 km in the latitudinal direction) and 15 vertical σ-coordinate levels. The model domain is from 35.0° to 51.5°N, and from 127.0° to 142.5°E. The bottom topography is obtained from the smoothed Naval Oceanographic Office Digital Bathymetry Data Base 5 min by 5 min resolution. The horizontal diffusivities are modeled using the Smagorinsky (1963) form with the coefficient chosen to be 0.2 for this application. No atmospheric forcing is applied to the model.
Closed lateral boundaries, that is, the modeled ocean bordered by land, were defined using a free-slip condition for velocity and a zero-gradient condition for temperature and salinity. No advective or diffusive heat, salt, or velocity fluxes occur through these boundaries. At open boundaries, we use the radiative boundary condition with zero volume transport.
d. Initial conditions
The model was integrated with all three components of velocity (u, υ, w) initially set to zero, and with the same potential density used by Mellor et al. (1994) for investigating the horizontal pressure gradient error,
which can be treated as an approximation to the area-averaged, vertical density distribution. Here z is the vertical coordinate, Hρ = 1000 m, and ρ0 = 103 kg m−3.
e. Mode splitting
For computational efficiency, the mode splitting technique (Blumburg and Mellor 1987) is applied with a barotropic time step of 25 s, based on the Courant–Friederichs–Levy (CFL; 1928) computational stability condition and the external wave speed; and a baroclinic time step of 900 s, based on the CFL condition and the internal wave speed.
f. Symmetric cancellation
The POM contains the symmetric cancellation scheme (Mellor et al. 1994). Without the cancellation scheme, we compute the scalar pressure field first by vertical integration of density and then its horizontal gradient via straightforward discretization. With the cancellation scheme, we compute the pressure gradient force as the vertical integral of discretized Jacobian of density and z coordinate. With this treatment, the problem of hydrostatic inconsistency can be avoided and the spurious residual gradient is drastically reduced (Mellor et al. 1994).
7. Horizontal pressure gradient force error
When the density field (23) is used, the horizontal pressure gradient should be zero, theoretically; however, it is usually nonzero in the σ-coordinate system due to the computational error. We may use the horizontal pressure gradient force on the three-dimensional grid cell
at day 20 to evaluate the performance of the three schemes in computing the horizontal pressure gradient vector for the JES POM: the second-order scheme without the symmetric cancellation (Fig. 4a), the second-order scheme with the symmetric cancellation (Fig. 4b), and the sixth-order AP scheme with the symmetric cancellation (Fig. 4c) at four different σ levels (1, 5, 9, and 14). Here, Δx, Δy, and Δσ are the grid spacings in the x, y, and σ axes. The σ levels 1 and 14 are the surface and bottom, respectively. We use the cubic-spline interpolation to compute ∂p/∂σ. Using the same second-order difference scheme, we find an order of magnitude reduction in the horizontal pressure gradient error with the symmetric cancellation (Fig. 4b) compared to without the symmetric cancellation (Fig. 4a). Using the same symmetric cancellation (Mellor et al. 1994), we find another order of magnitude reduction in the horizontal pressure gradient error with the sixth-order AP scheme (Fig. 4c) compared to the second-order scheme (Fig. 4b).
We compute horizontally averaged and maximum values of |A|, over a σ level, to compare the performance of these schemes (Fig. 5). The symmetric cancellation scheme (Mellor et al. 1994) has the capability to drastically increase the accuracy using the same second-order scheme. For example, the horizontally averaged (maximum) values of |A| generally reduce by factors of 5–32 (8–23) with the symmetric cancellation scheme compared to without the symmetric cancellation scheme. Moreover, using the symmetric cancellation we see a further reduction in the horizontally averaged (maximum) values of |A| by factors of 6–9 (4–9) with the sixth-order AP scheme compared to the second-order scheme.
8. Error estimation using the JES POM
We integrated the diagnostic JES POM with Mellor et al.'s (1994) symmetric cancellation scheme from zero velocity and horizontally integrated density field (23) for 20 days with the second-order scheme and the sixth-order AP scheme, respectively. Without computational error, the model ocean should keep motionless. With the computational error, the model generates false motion. The absolute value of the velocity is taken as the error velocity. The smaller the error velocity, the more accurate the numerical scheme would be.
a. Volume transport streamfunction
Volume transport streamfunction (Ψ) on days 5 and 20 are presented in Fig. 6 for the second-order scheme and in Fig. 7 for the sixth-order AP scheme. All the figures show a multi-eddy structure in coincidence with the steep bottom topography (Fig. 3). The difference between two Ψ isolines indicates the volume transport at that location. The smaller the difference, the less the volume transport. The contour interval of Ψ isolines is five times smaller in Fig. 7 than in Fig. 6. This indicates the error reduction by a factor of 5 using the AP sixth-order scheme compared to using the second-order scheme.
b. Temporal variations of global and peak error velocities
Two quantities are commonly used for identifying the error evolution: peak and global error velocities. The peak error velocity is the maximum absolute value the spurious velocity generated by the pressure gradient errors. Figure 8a shows the time evolution of the peak error velocity for the first 20 days of integration with the second-order and sixth-order AP schemes. The peak error velocities increase to a maximum value and then decay as the inertial oscillation superimposed into some mean values: 4.80 × 10−3 m s−1 (the second-order scheme), and 0.62 × 10−3 m s−1 (the sixth-order AP scheme). The maximum peak error velocity is 5.38 × 10−3 m s−1 (0.76 × 10−3 m s−1) at day 12 (day 7.5) for the second-order (sixth-order AP) scheme. The steady approach of the peak error velocities to these values for the three schemes indicates the stability of the computation. Furthermore, the peak error velocity for the second-order case is roughly eight times that of the sixth-order AP case.
At each time step, we compute the spatially averaged error velocity over the whole JES (called the global error velocity) to represent the overall accuracy of the schemes. Figure 8b shows the time evolution of the global error velocity for the first 20 days of integration with the second-order and the sixth-order AP schemes. It increases with time at a much slower pace using the sixth-order AP scheme than the second-order scheme. On day 20 the global error velocity is 2.02 × 10−4 m s−1 for the second-order scheme and 0.34 × 10−4 m s−1 for the sixth-order AP scheme. The global error velocity for the second-order scheme is roughly six times that of the sixth-order AP scheme.
c. CPU time
We integrate the POM for 20 days on the SGI Origin-200 machine to identify the cost of pressure gradient schemes. The central processing unit (CPU) times are 21.067 h for the second-order scheme and 21.167 h for the sixth-order AP scheme, respectively. There is no noticeable CPU time increase from the second-order scheme to the sixth-order AP scheme.
9. Open boundary influence
a. Coastal geometry
One of the difficult problems in shallow water modeling is the uncertainty of the open boundary condition (OBC). At open boundaries where the numerical grid ends, the fluid motion should be unrestricted. Ideal open boundaries are transparent to motions. Two approaches, local type and inverse type, are available for determining OBC (Chu et al. 1997) with low-order conditions. Before converting any ocean model from a second-order scheme to a sixth-order scheme, it is important to verify if a high-order interior plus low-order boundary conditions (such as open boundary conditions) would degrade the interior solution. Consider a horizontally homogeneous and stably stratified coastal ocean with a longitudinal and straight coastline and three open boundaries in the south (20°N), the north (33.5°N), and the east (x = xE = 320 km). Choose coordinates so that the y axis coincides with the coast, positive x pointing offshore. The coastal water is divided into three parts (Fig. 9a): shelf (0 ≤ x ≤ x0), slope (x0 < x ≤ x1 = 200 km), and deep water (x1 < x ≤ xE).
Here, the water depth (h) is only a function of offshore distance x. At the coast (x = 0), the water has a minimum depth hmin (Fig. 9). An analytical bottom topography is proposed such that shelf and slope are arcs of two circles. The shelf has a smaller radius (r), and the slope has a larger radius (R). The two arcs are connected in such a way that the tangent of the bottom topography, dh/dx, is continuous at the shelf break (x = x0). This requirement is met using the same expanding angle (θ) for both arcs (Fig. 9b). Thus, this bottom topography has 5 degrees of freedom: hmin, hmax, r, R, and θ.
Off-shore distance x and bottom topography h are nondimensionalized by
and in turn r and R are nondimensionalized by
The bottom topography (Fig. 9) can be represented analytically by
As in the previous section, initially motionless coastal ocean with an exponentially stratified density (23) and without forcing will keep motionless forever, that is, V = 0. Any nonzero model velocity can be treated as an error. We use the POM implemented for such a coastal region to test the open boundary influence for the high-order AP scheme. At open boundaries, we use the radiative boundary condition with zero volume transport.
The parameters are chosen in this study by
This implies that the width of the shelf break (x0) is 33.33 km.
In the numerical integration, the space is discretized by
and barotropic and baroclinic time steps are chosen as 6 and 180 s in order to satisfy the CFL condition.
b. Volume transport streamfunction
Volume transport streamfunction (Ψ) on day 10 is presented in Fig. 10a for the second-order scheme and in Fig. 10b for the sixth-order AP scheme. Both figures show a multi-eddy structure in coincidence with the continental slope (Fig. 9). The difference between two Ψ isolines indicates the volume transport at that location. The smaller the difference, the less the volume transport. The contour interval of Ψ isolines is around 200 times smaller in Fig. 10b than in Fig. 10a. This indicates the error reduction by a factor of 200 using the AP sixth-order scheme compared to using the second-order scheme.
c. Horizontal pressure gradient force error
We compute horizontally averaged and maximum values of |A| at day 10 over a σ level, to compare the performance of the sixth-order AP and second-order schemes (Fig. 11). The horizontally averaged (solid) and maximum (dashed) |A| values are quite small at the upper ocean for both the second-order scheme (Fig. 11a) and the sixth-order AP scheme (Fig. 11b). Using the second-order scheme, the errors (|A| values) increase monotonically versus σ, with largest errors next to the ocean bottom (σ = 14): 2943 N for the averaged |A| and 18 582 N for the maximum |A|. Using the sixth-order AP scheme (Fig. 11b), the errors (|A| values) are much smaller than the second-order scheme. The largest errors occur at one level higher than the second-order scheme (σ = 13): 44 N for the averaged |A| and 977 N for the maximum |A|.
Error ratios between second-order and sixth-order schemes show drastic reductions by factors of 30–514 in the horizontally averaged |A| values (solid) and by factors of 43–1160 in the horizontally maximum |A| values (dashed) (Fig. 11c). The error reduction becomes very evident in the middle level (σ = 6) of the ocean.
d. CPU time
We integrate the POM for this open boundary test case for 10 days on the SGI Origin-200 machine to identify the cost of pressure gradient schemes. The CPU times are 501.07 min for the second-order scheme and 504.19 min for the sixth-order AP scheme, respectively. There is no noticeable CPU time change from the second-order scheme to the sixth-order AP scheme.
1) This study shows a process of building a series of accuracy progress difference schemes with the lower-order scheme as the “source” term in the higher-order scheme. Such schemes have the following features: three-point, accuracy progressive, and only tridiagonal solver needed. The second-order AP scheme is the same as the central difference scheme, the fourth-order AP scheme is the same as the standard Pade scheme, and the sixth-order AP scheme has the sixth-order accuracy.
2) Fourier analysis shows the progressive error decrease in modified wavenumber from the second-order AP scheme to the fourth-order scheme, and then from the fourth-order scheme to the sixth-order scheme.
3) The Princeton Ocean Model implemented for the Japan/East Sea is used to demonstrate the benefit of using the sixth-order AP scheme. Two calculations of unforced flow are performed with the second-order scheme and the sixth-order AP scheme, respectively. The results show that the sixth-order AP scheme has error reductions by factors of 5–8 compared to the second-order difference scheme.
4) Using the sixth-order AP scheme for computing the pressure gradient does not require much more CPU time. Taking POM for the Japan/East Sea as an example, the CPU time for the sixth-order AP scheme is almost the same as for the second-order scheme.
5) Both standard and staggered AP schemes were developed. The sixth-order AP scheme on a staggered C grid can be easily implemented into current ocean models.
6) The high-order AP scheme has the capability to handle the open boundary condition. There is no degradation of the interior solution using high-order AP schemes and low-order boundary conditions.
7) Future studies include applying AP schemes to nonuniform grid systems as well as designing even higher-order schemes such as an eighth-order AP scheme.
This work was funded by the Office of Naval Research, the Naval Oceanographic Office, and the Naval Postgraduate School.
Error Estimation of AP Schemes on a Standard Grid
The truncation error is usually obtained from Taylor series expansion. The first and second derivatives at the grid point xi are represented by
which shows the second-order accuracy of Δf/Δx, Δ2f/Δx2. Also, f′ and f″ at the grid point xi can be represented by (standard Pade scheme)
which indicates the fourth-order accuracy; that is, the difference between the first-order derivative and its fourth-order difference counterpart at two grid points xi+1 and xi−1 should be
The difference between the second-order derivative and its fourth-order difference at two grid points xi+1 and xi−1 should be
Chu and Fan (1998) provide the following error estimation:
Elimination of (f″i+1 − f″i−1) from (A5) and (A6) leads to
Corresponding author address: Dr. Peter Chu, Dept. of Oceanography, Naval Postgraduate School, Monterey, CA 93943. Email: firstname.lastname@example.org