S-coordinate models (SPEM/SCRUM)

Design philosophy

The desire to avoid spurious effects associated with discontinuous (i.e., step-wise) representation of bathymetry and sidewall geometry suggests the use of a numerical grid which smoothly fits the irregular shape of the model domain. This has led in both meteorology and oceanography to classes of models with topography-following vertical and/or orthogonal curvilinear horizontal coordinates. In these systems, coordinate transformations map the ocean (or atmosphere) to a rectangular computational domain. This ensures efficient use of computer resources, since all grid points lie within the fluid.

The terrain-following (or sigma) coordinate was first introduced in atmospheric modeling (Phillips, 1957) and has since become a standard alternative in ocean modeling. One ocean model which uses a terrain-following vertical coordinate is SPEM (the s-coordinate primitive equation model) originally developed by Haidvogel et al. (1991) for high-resolution process studies. In its earliest version, a mixed finite difference (in the horizontal) and Galerkin (in z) solution procedure was employed, along with a single state variable (density). The most recent version, which incorporates two thermodynamic variables and S and a nonlinear equation of state, a special masking technique for islands and promontories, and a staggered finite difference treatment in the vertical coordinate, is suitable to a wider range of regional and basin-wide studies. In addition to the rigid lid version, a separate free sea surface version (SCRUM, the s-coordinate Rutgers University model) is available (Song and Haidvogel, 1994). Both programs are written in a modular fashion, and run very efficiently on vector computers; parallel versions for different computer architectures exist as well. Detailed user's guides for SPEM and SCRUM are maintained by Hedstrom (1994, 1997).

System of equations

The vertical topography-fitting transformations employed are either

for the rigid lid version (SPEM) or, in the case of free sea surface (SCRUM),

In the case of a linear relationship between s and z, the transformation is equivalent to the traditional sigma coordinate and the hydrostatic relation, therefore, becomes

and the pressure gradient terms in are

Note that the pressure gradient splits into the sum of two terms, the gradient along the coordinate lines and a hydrostatic correction.

Both SPEM and SCRUM are written in horizontal curvilinear coordinates. The transformed coordinates are defined by the relations

where m(x,y) and n(x,y) are the scale factors which relate the differential distances to the physical arc lengths . (See, for example, Arakawa and Lamb, 1977.) Note that this general formulation of curvilinear coordinates includes Cartesian coordinates (by setting m = n = constant) as well as spherical coordinates (usually employed in large-scale ocean modeling) with

where is the geographical latitude and

Under these horizontal and vertical transformation, the hydrostatic primitive equations become


where are the (X, Y, s) components of the velocity vector , h = h(X,Y) is the fluid depth, and P(X, Y, s, t) is the dynamic pressure . All other variables are used in their standard notation. As before, D and F are dissipative and forcing terms, if any. The "vertical velocity" in this coordinate system,

includes both "upwelling" and "upsloping" components of the vertical motion. In the transformed coordinate system, the kinematic boundary conditions at the surface (s = 0) and the bottom(s = -1) becomes

a significant simplification of the lower boundary condition for vertical velocity. Lastly, the full UNESCO equation of state is implemented in a vectorized form, with modified coefficients for use with potential rather than in-situ temperatures.

Depth-integrated flow

The elliptic solver for the streamfunction is based on a multi-grid algorithm (MUD2 by Adams, 1991 and Adams et al. 1992), which is more efficient and accurate than traditional SOR schemes. A potential drawback is the restriction to certain combinations of horizontal numbers of grid points. For masked areas of the model domain, additional modifications are required. In SPEM, a capacitance matrix method is utilized (Wilkin et al., 1995). The explicit free surface code (SCRUM) is algorithmically simpler, but requires, for computational economy, special treatment of the fast barotropic mode and its coupling to the slower baroclinic modes. This is accomplished via split-explicit time-stepping, and subsequent temporal filtering of the barotropic mode. The latter step involves time averaging of the vertically integrated velocities and the sea surface height, before they replace those values found by the three-dimensional equations on a longer time-step.

Spatial discretization, grids and topography

Standard sigma coordinates subdivide the vertical into an equal number of points, according to s = z/h(X,Y). Figure 2 shows four examples of a terrain-following discretization at the continental margin. In all cases, the vertical resolution in shallow areas is extremely fine. This high vertical resolution over the continental margins and other areas of shallow bathymetry is a potential virtue of terrain-following coordinate models. Note however that this may cause a time-step limitation due to the very fine vertical grid spacing, unless implicit methods are used. Higher resolution near the surface can be achieved by non-uniform placement of the levels in sigma space (e.g., Barnier et al., 1998; Ezer and Mellor, 1997, see figure below).

Figure 2: Vertical discritization for terrain-following (s) coordinate model with 20 levels. (a) standard (equidistant) sigma coordinate, (b) discretization with higher resolution near the surface (quadratic in s-space).

As an extension to standard terrain-following coordinate transformation, a nonlinear stretching of the vertical coordinate can be applied that depends on local water depth (Song and Haidvogel, 1994). This option can be used to generate a more uniform vertical resolution near the surface and consequently a better representation of the mixed layer and thermocline. The strong (and sometimes undesirable) coupling between s-coordinate surfaces and local water depth is also reduced. The transformation used in SPEM/SCRUM is


with , and . For large , the coordinate lines crowd near the surface (Figure 2c); additionally, if approaches 1, resolution at the bottom boundary is enhanced (Figure 2d). This latter figure shows a discretization similar to the collocation points of Chebyshev polynomials (used in the original version of SPEM), and has been found to be particularly well suited in process studies (Haidvogel et al., 1993; Beckmann and Haidvogel, 1997).

Figure 2: Vertical discretization for terrain-following (s) coordinate model with 20 levels. (c stretched coordinate () for higher resolution near the surface, (d) stretching () for boundary layer resolution.

The terrain-following coordinate system, when discretized, is subject to a special systematic error, associated with the possibility of significant errors in the horizontal pressure gradients (Haney, 1991; Deleersnijder and Beckers, 1992; Beckmann and Haidvogel, 1993; Mellor et al., 1994). These errors arise due to the splitting of the pressure gradient term into an "along-coordinate surface" component and a "hydrostatic correction". The latter term is there to remove that portion of the pressure variations along the sloping s-surfaces which are due to vertical hydrostatic pressure changes (which do not vary horizontally and hence do not accelerate the fluid). Unfortunately, both terms are large, and cancellation of the hydrostatic resting pressure is not exact, due to imbalanced numerical truncation errors in the two terms.

The pressure gradient errors that result depend on the steepness of the topography, both the horizontal and vertical resolution, and the strength of the stratification. Various methods have been proposed to reduce them, e.g., by applying the pressure gradient force only to the dynamically active deviation from a domain-wide reference profile (Gary, 1973), and by using higher-order finite difference approximations in the horizontal (McCalpin, 1994; Chu and Fan, 1997). While these work successfully in idealized process studies with relatively smooth topography, the improvement is less significant in coarser-resolution applications.

Idealized forms of bathymetry for process studies can be chosen such that gradients are well resolved. In most cases, however, realistic topography requires some smoothing before it can be used in terrain-following models. For SPEM/SCRUM, a useful parameter is found to be (Beckmann and Haidvogel, 1993)

Empirical studies have shown that reliable results are obtained if r does not significantly exceed 0.2.

Several alternate forms of the s-coordinate pressure gradient terms can be formulated, each with slightly different numerical properties. For example, the form

conserves the net bottom torque; it can be shown to give the correct JEBAR. The properties of the discretized version have been explored and a residual error remains in most cases. It seems advisable to test the magnitude of the resulting error field in resting stratification test runs (see Beckmann and Haidvogel, 1993; Barnier et al., 1998), where a horizontally uniform stratification is prescribed and any developing flow is known to be erroneous.

While the current versions of SPEM/SCRUM use second-order finite differences in all three spatial directions, earlier versions were based on a spectral approach in the vertical (using Chebyshev polynomials) and had a non-staggered grid with collocation levels at the surface and the bottom. The advantages were highly accurate vertical operations and a more direct application of the vertical boundary conditions. However, numerical stability was found to require shorter time-steps (as expected for a higher-order Galerkin-based scheme) and more care in the selection of viscous parameters.

In the horizontal (X,Y), a centered second-order finite-difference approximation is adopted (an Arakawa "C" grid; Arakawa and Lamb, 1977). Lateral boundary conditions are usually free-slip (the more "natural" form for the "C" grid, as the derivative of the tangential velocity has to be specified) but an approximation to no-slip boundary conditions is also available. It should be noted that a similar restriction to equation 4.10 applies to the curvilinear horizontal grid (see, e.g., Figure 4, which should vary gradually rather than abruptly.) Measures of this variation are

(in both horizontal directions) which should be of moderate size to exclude spurious effects on the variable grid (phase and amplitude errors during partial reflection of waves, discontinuous transport of tracers, etc.).

Semi-discrete equations

We define , and as centered finite-difference approximations to , and with the differences taken of the distances , and , respectively. Similarly, , and represent averages taken over the distances , and . represents a second-order vertical integral computed as a sum from level s to the surface at s = 0. We also define

Then, the semi-discrete form of the dynamic equations is

The vertical velocity is computed as

Finally, in the case of a rigid lid (SPEM), the barotropic vorticity is

and the barotropic velocities are computed from

For the free sea surface () code (SCRUM), we have

Note that on a "C" grid, there is no averaging of the variables for the barotropic mode.

Temporal Discretization

The time-stepping scheme follows standard procedures for primitive equation models. The baroclinic (depth-varying) modes of the velocity distribution are obtained by direct time-stepping of the momentum equations, having removed their depth-averaged component. The temperature and salinity equations are similarly advanced in time using a leapfrog time-stepping scheme with an occasional trapezoidal correction. In the presence of a rigid lid, the barotropic (depth-averaged) components of the horizontal velocity are determined by first taking the curl of the depth-integrated momentum equations. The resulting vorticity equation for the transport streamfunction may then be time-stepped, subject to closed or cyclic boundary conditions in one or both directions. The free surface version SCRUM integrates the surface elevation explicitly with a much smaller time step (split-explicit time-stepping). Implicit vertical diffusion is available for cases where the water depth is very shallow or where convection is parameterized by locally increased vertical diffusion/viscosity.

Additional features

Barnier et al. (1998) describe a set of open boundaries for SPEM, based on an Orlanski radiation condition with sponge layer damping. SCRUM offers a wide variety of similar options. A horizontal Shapiro filter (Shapiro, 1970) that eliminates grid scale noise is available, and found helpful for the epineutral mixing option, where it is used to smooth the rotation angle.

The task of horizontal grid generation can be a complicated, often iterative procedure; a separate package, called gridpak, is provided for convenient set-up and modification of an orthogonal horizontal grid in an irregular domain (Wilkin and Hedstrom, 1998). Figure 4 shows an example of a curvilinear grid (including masked areas) for the Denmark Strait region between Greenland and Iceland. The horizontal curvilinear grid requires a flexible way to produce forcing and initial data from regularly gridded datasets. For this purpose, an objective analysis package is available for SCRUM as a preprocessing support program for model set-up.

Figure 3: Horizontalcurvilinear orthogonal grid (97 x 65 points) for a regional model of the Denmark Strait region. Minimum and maximum grid spacings are approximately 5 and 40 km, respectively. Masking is applied to some parts of the very irregular Icelandic penninsula.

Diagnostics for process studies are mostly hand-tailored for each application. For example, averaging along one coordinate or in one direction may be required, second moments are needed, dynamic balances have to be evaluated (see, e.g., Haidvogel et al., 1991), terms need to be transformed into cylindrical coordinates (see, e.g., Haidvogel et al., 1993; Beckmann and Haidvogel, 1997). To date, many analysis tools exist in the SPEM/SCRUM user community.

Unlike other models of this complexity, the SPEM/SCRUM family has a standardized set of plotting included in the code. The routines are based on NCAR graphics (Clare and Kennison, 1989) and provide a variety of instantaneous model fields for monitoring of the model run. An interface is available to write NetCDF binary output which then can be read in by other state-of-the-art visualization programs.

The model is written in f77 and uses mks units throughout. Individual selections can be specified via cpp options. There are makefiles for various computer platforms. Parallel model versions exist for SGI and Cray T3E. They can be obtained from


under the SPEM and SCRUM keywords.

Concluding Remarks

Other models based on terrain-following coordinates include:

These models have been mainly used for coastal ocean studies (see also Haidvogel and Beckmann, 1998). They are all formulated on a "C" grid but differ in aspects of the numerical algorithms, available parameterizations and coding.

Copyright ©2002-2007 Ocean-Modeling.org. All rights reserved. Webmaster