Spectral Element Ocean Model (SEOM)
The numerical algorithms utilized in the models described in the previous sections have relied on first- and second-order finite difference approximations and structured horizontal grids to discretize the governing equations. Their virtue has been the simplicity of their numerical foundation and their robustness. Their disadvantages are a limited geometrical flexibility, low convergence rates as resolution is increased, and good but limited scalability on parallel computers. The spectral element model uses alternate numerical methods that seek to address these design issues.
SEOM is based on the spectral element method (Patera, 1984). The spectral element method can be most concisely described as an h-p type Galerkin finite element method which approximates the solution within each element with a high-order polynomial interpolant. It offers several desirable properties for ocean simulations: geometrical flexibility with a spatial discretization based on unstructured grids, high-order convergence rates, and dense computations at the elemental level leading to extremely good scalability characteristics on parallel computers (Fischer, 1989; Curchitser et al., 1996).
Spectral elements have most of the advantages of spherical harmonics including exponential converge under p-refinement, and the avoidance of polar singularities. In addition, spectral elements allow for local refinement of the mesh, high-order algebraic convergence under h-refinement, and complete parallelism. Trade-offs for these useful properties are a greater computational cost per gridpoint/time-step (as expected for a higher-order Galerkin scheme) and a lower degree of implicit smoothing (thus necessitating more clever subgridscale viscous operators). The use of widely variable elemental sizes also places a premium on spatially dependent smoothers.
The governing equations (in vector notation) are:
Note that, although formulated in geopotential coordinates, the solution procedure is essentially bathymetry-following since the three-dimensional elements are isoparametrically mapped to the bathymetry and continental sidewalls as part of the solution procedure.
The velocity boundary conditions at the top and bottom surfaces are the kinematic boundary conditions of no-normal flow:
and the dynamic boundary conditions specifying the stresses at the sea surface (
)
where 
is the prescribed wind stress, 
is the outward unit normal, 
is the resting depth, and 
is the linear drag coefficient.
Surface gravity waves are isolated in a set of two-dimensional equations which are advanced in time separately from the three-dimensional equations. The two-dimensional equations are similar to the shallow water equations and can be obtained by vertical integration of the three-dimensional momentum and continuity equations:
where 
is the depth-average velocity:
and C is the coupling term between the two and three-dimensional momentum equations:
The first and second terms on the right-hand-side of the above equation contain the contributions to the vertically averaged flow of the baroclinic velocity and of the baroclinic pressure gradient, respectively. The surface and bottom boundary conditions on the velocity have been used in the derivation of the depth-averaged equations.
Spatial discretization, grids and topography
Since the sea surface is moving, the domain V occupied by the fluid is time-dependent. In order to simplify the calculations, the domain is mapped to a steady computational space 
where
The height 
is thus mapped into the interval 
.
Under this transformation, the horizontal gradient, material derivative and divergence operators become, e.g.,
where 
, 
and 
denote the x and y derivatives, and the horizontal gradient operator along constant Z-lines. 
is the vertical velocity in the mapped space:
The Jacobian of the mapping between the unsteady physical domain and the steady computational domain is
which is set equal to 1 since 
in basin-scale applications. This approximation results in a tremendous savings in computational overhead since the three-dimensional mass matrix discussed below becomes steady and need not be updated at every time-step.
The spatial discretization in SEOM relies on isoparametric and conforming hexahedral elements (cubes with curved surfaces) with 
continuity. The three-dimensional variables in each element are approximated with a high-order Lagrangian interpolant whose collocation points are the Gauss-Lobatto roots of the Legendre polynomials. (See, e.g., Boyd (1989).) Hence, for example, 
in element e is interpolated as
where 
and 
are the horizontal interpolation functions, 
is the vertical interpolation function, 
is the local coordinate system in element e, N is the number of collocation points in each horizontal direction, M is the number of collocation points in the vertical, and 
is the value of 
at the 
Gauss-Lobatto-Legendre roots 
. Notice that the orders of the interpolation in the vertical and horizontal directions are independent, so they can be adjusted individually to achieve the best h-p balance.
The two-dimensional variables are interpolated according to
where 
refers to the surface displacement interpolation function whose collocation points are given by the (N-2) Gauss-Lobatto-Legendre roots. The 
collocation points are staggered with respect to those of 
in order to eliminate spurious pressure oscillations when the barotropic flow is almost divergence-free (Iskandarani et al., 1995).
Given the complicated bathymetry of ocean basins, and the lack of appropriate three-dimensional grid generation software, the elemental partitioning in the vertical is simplified by restricting it to a structured discretization; i.e., the number of vertical elements does not vary horizontally. A three-dimensional spectral element grid can thus be produced by stacking vertically a number of two-dimensional grids. The figure below gives an example of two possible vertical grids in a region of strong topographic contrast. Note that the positioning of the elemental boundaries in the vertical is quite general. For example, a single element can be assigned to follow the sea surface, providing uniform vertical resolution in the surface mixed layer irrespective of local water depth and the fact that the underlying elemental structure is terrain-following.
The local 
coordinate lines are restricted to be vertical in order to decouple the vertical and horizontal calculations. This decoupling reduces the operation count tremendously if the vertical diffusion operator has to be integrated implicitly. These two restrictions are not severe since it is still possible to fit the three-dimensional elements to the bathymetry of the basin, to distribute the vertical resolution according to a priori dynamical considerations, and to retain the unstructured nature of the grid in the horizontal direction.
Below is a sample spectral element grid encompassing the global ocean, including the Arctic Ocean and emphasizing the North Eastern Pacific Ocean.
The interpolations of the two previous equations are substituted into the variational formulation of the HPE to yield a system of ordinary differential equation, after settings
and evaluating the resulting integrals numerically. The variational formulation of the governing equation
where 
is the three-dimensional test function associated with the velocity and tracer fields, 
refers to the free surface boundary, and 
refers to the sea-bed boundary. The boundary integrals weakly enforce the Neumann and Robin (mixed) boundary conditions imposed on the velocity field. Essential boundary conditions are enforced by zeroing the test function 
on those portions of the boundary where the boundary conditions are of Dirichlet type, and setting 
equal to the imposed boundary value, e.g., zero on no-slip walls.
Likewise, the variational formulation of
where 
and 
are surface inputs of heat and salt, respectively. The continuity equation is differentiated in the vertical before integration in order to ease the imposition of surface and bottom boundary conditions. The variational formulation becomes
The variational formulation of the barotropic equations is a little more elaborate since the 
field uses a staggered grid with respect to the velocity. It is
where 
and 
refer to the two-dimensional test functions associated with 
and 
, respectively. The divergence term in the continuity equation has been integrated by parts to allow the imposition of inflow as weak boundary conditions, and to produce a symmetric positive-definite matrix equation for 
in case an implicit integration scheme is adopted. The q in the top equation refers to the inflow per unit width. Notice that the integrals in the above equations are performed over a horizontal surface.
The prognostic variables are the horizontal velocity 
, the temperature T, the salinity S, the surface displacement 
and the barotropic velocity 
. The vertical velocity, density, and baroclinic pressure are calculated diagnostically from the continuity equation, the equation of state, and the hydrostatic equation, respectively.
The temporal discretization of
is based on a semi-implicit integration scheme. The nonlinear advection term, the Coriolis term, the horizontal viscous term, and the baroclinic pressure gradient term are integrated with a third-order Adams-Bashforth scheme. The vertical viscous term is integrated with an implicit second-order Crank-Nicholson scheme for enhanced numerical stability, since the vertical viscosity 
may need to be increased during calculations to simulate rapid vertical mixing associated with convective overturning events.
The tracer equations are also integrated in time with a semi-implicit scheme: third-order Adams-Bashforth for the advection and horizontal diffusion terms, and second-order Crank-Nicholson for the vertical diffusion term. In matrix notation, the tracer equations become:
where c is the weight factor of the implicit scheme (c = 1 / 2 for the Crank-Nicholson scheme).
The SEOM code supports two methods of integration for the barotropic equations: a synchronous integration scheme for circumstances (e.g., in shallow fluids) in which the permissible barotropic and baroclinic time-steps are nearly equal, and a split-explicit scheme for the usual circumstance in which the time-stepping restriction associated with surface gravity waves is dominant. The split-explicit scheme proceeds in two stages. First, the barotropic equations are integrated explicitly using a small time-step, 
, that respects the CFL stability limit of the gravity waves. The barotropic equations take the form:
If the ratio between the long and short time-step is denoted by L, then the above integration is carried out from s=0 to s=2L. The surface displacement thus calculated is time-averaged to filter out the fast time scales:
The three-dimensional equations are then integrated, using the filtered surface displacement, on a long time-step whose size can be chosen based on the time scales of the baroclinic dynamics. The main advantage of this procedure is that it avoids the solution of an implicit system of equations, and is simple to implement. Its main disadvantage is that it complicates the evaluation of the coupling term, C, between the barotropic and three-dimensional equations. In the current implementation of the model, this term is fixed in time at the previous time level and is updated only after the integration of the three-dimensional equations. Moreover, the time-filtering exacts a computational penalty: in order to center the time-averaging procedure around the next long time level, 
, the barotropic equations must be integrated until time level 
before the averaging procedure is applied.
The parallelization of SEOM relies on a domain decomposition approach and a SPMD programming paradigm. An optimal decomposition would minimize the amount of communication required while maximizing the density of local computations. A breadth-first-search algorithm has been devised with these constraints in mind, and has been presented in Curchitser et al. (1996). Given the short extent of the computational domain in the vertical and the tight coupling of the equations in that direction, it is best to allocate all vertically aligned elements to the same processor. The three-dimensional sub-domains are hence simple extensions of the two-dimensional sub-domains and can be viewed as columns of elements. Since the 3D grid is structured in the vertical, load-balancing is automatically guaranteed if the two-dimensional sub-domains have an equal number of elements.
Procedures for the two-way coupling (nesting) of multiple nonconforming grids (i.e., those having differing h-p truncations) have recently been implemented (Levin et al., 1999). Based on the mortar element method, the resulting procedures make SEOM the only ocean model at this time to offer a unified approach to non-uniform horizontal resolution and two-way grid nesting.
Lastly, a variety of special approaches to subgridscale mixing and horizontal smoothing have been developed for use with SEOM. Examples include spectral filtering procedures to emulate bi-harmonic and higher-order lateral mixing and flow-dependent mixing coefficients based (e.g.) on grid Reynolds numbers (Levin et al., 1997).
The model is written in f90 and uses mks units throughout.
Another finite element ocean circulation model, based on a low-order triangular decomposition is QUODDY (Ip and Lynch, 1994).
Copyright ©2002-2007 Ocean-Modeling.org. All rights reserved. Webmaster