The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model

https://doi.org/10.1016/j.ocemod.2004.08.002Get rights and content

Abstract

The purpose of this study is to find a combination of optimal numerical algorithms for time-stepping and mode-splitting suitable for a high-resolution, free-surface, terrain-following coordinate oceanic model. Due to mathematical feedback between the baroclinic momentum and tracer equations and, similarly, between the barotropic momentum and continuity equations, it is advantageous to treat both modes so that, after a time step for the momentum equation, the computed velocities participate immediately in the computation of tracers and continuity, and vice versa, rather than advancing all equations for one time step simultaneously. This leads to a new family of time-stepping algorithms that combine forward–backward feedback with the best known synchronous algorithms, allowing an increased time step due to the enhanced internal stability without sacrificing its accuracy. Based on these algorithms we design a split-explicit hydrodynamic kernel for a realistic oceanic model, which addresses multiple numerical issues associated with mode splitting. This kernel utilizes consistent temporal averaging of the barotropic mode via a specially designed filter function to guarantee both exact conservation and constancy preservation properties for tracers and yields more accurate (up to second-order), resolved barotropic processes, while preventing aliasing of unresolved barotropic signals into the slow baroclinic motions. It has a more accurate mode-splitting due to redefined barotropic pressure-gradient terms to account for the local variations in density field, while maintaining the computational efficiency of a split model. It is naturally compatible with a variety of centered and upstream-biased high-order advection algorithms, and helps to mitigate computational cost of expensive physical parameterization of mixing processes and submodels.

Introduction

Realistic oceanic circulation models are usually based on Boussinesq, hydrostatic momentum and mass balances, material tracer conservation, seawater’s equation of state, and parameterized subgrid-scale transports. Their time integration is made with a decomposition of the 3D fields into barotropic (depth-averaged) and baroclinic (the residual) parts to facilitate the calculation of the pressure-gradient force (Bryan and Cox, 1969). The motivation to build a free-surface oceanic model is twofold. From a physical point of view, it is desirable to recapture processes lost or altered by the rigid-lid assumption. These include tidal motions, altered dispersion relations for the Rossby waves, etc. The other motivation comes from computational economics: as pointed out by Killworth et al. (1991), there is a natural physical ratio of phase speeds for the external and internal gravity-wave modes. Once the model time step is chosen from the CFL criterion based on the fastest baroclinic wave speed, the external mode has to be treated by either (i) a streamfunction method using rigid-lid approximation; or (ii) a two-dimensional (2D) pressure Poisson equation for pressure on the rigid-lid or due to free-surface elevation; or (iii) a special 2D barotropic submodel that uses a smaller time step chosen from a CFL criterion based on the barotropic speed. Approaches (i)–(ii) require solution of a 2D elliptic problem (Dukowicz, 1994, Dukowitz and Smith, 1994) at every time step that, with a conventional Successive Over-Relaxation (SOR) or similar method, requires a number of iterations on the order of the number of grid points in the longest direction of the computational domain. Since on this path the number of operations needed at every grid point at every time step tends to increase with the resolution, on finer grids the approach (iii) tends to be more efficient than the others with a threshold set by the ratio of phase speeds of the external and the fastest internal gravity waves compared to the number of grid points in the longest dimension,1|Vext||Vint|·max(Nx,Ny).

Despite the long-time existence of split-explicit versions for all three major classes of oceanic models––z-, sigma-, and density-coordinate––there are few published studies about the mathematical aspects: consistency, accuracy, and stability associated with mode splitting (e.g., Higdon and Bennett, 1996, Higdon and de Szoeke, 1997, Hallberg, 1997, Higdon, 2002 and an earlier theoretical work, Yanenko, 1971 and Skamarock and Klemp, 1992) more related to atmospheric modeling. Even within the limits of numerical stability based on the usual CFL condition, mode splitting may introduce additional sources of numerical instability not present in models with uniform time steps nor in rigid-lid models.

The purpose of the present study is to reconsider the computational kernel of an oceanic model, including the optimal choice of time-stepping algorithms for the barotropic and baroclinic momentum and tracer equations, and their mutual interaction. Here we advocate an integrated approach whose main focus is on the time-step limitation coming from the system as a whole that, as we will show, is typically more restrictive than the CFL limitations coming from each equation taken individually. We design a new finite-volume, finite-time-step discretization for the tracer equations to eliminate the conflict between integral conservation and constancy preservation properties associated with the variable free surface. We generalize the barotropic mode to take into account the nonuniform density. Collectively, these steps reduce the mode-splitting error and improve the stability, robustness, and efficiency of the model.

The topography-following vertical coordinate system implies that there is a transformationz=z(x,y,σ),where z is the Cartesian height and σ is the vertical distance from the surface measured as the fraction of the local water column thickness (i.e., −1  σ  0, σ = 0 corresponds to the free surface, z = ζ and σ = −1 corresponds to the oceanic bottom, z = h(x, y)). The resulting system of coordinates is nonorthogonal and leads to a set of chain rules for derivatives,xz=xσ-zxσ·z.In the case of the classical a-coordinate, (1.2) reduces toz=σ·h(x,y).

This may be combined with nonlinear stretching, S(σ)z(x,y,σ)=S(σ)·h(x,y)and further generalized into the S-coordinate of (Song and Haidvogel, 1994)—which in essence behaves like (1.4) in shallow regions and (1.5) in deep.

Past experience with σ-coordinate models and intercomparisons with z- and isopycnic-coordinate models (Beckmann, 1998, Willebrand et al., 2001) reveal that the solutions from (T-models exhibit stronger topographic sensitivity than the other two classes of models. This is attributed to the fact that the isosurfaces of the vertical coordinate intersect the isopycnals at some angle, even in the case of horizontally uniform stratification, which causes pressure-gradient error. One way to address this problem is to redesign the model algorithms making them less sensitive to such errors (Shchepetkin and McWilliams, 2003). It is also desirable to allow the possibility of a smooth transition from σ to z-coordinates, such that the top-most isosurfaces are nearly flat while the bottom-most are still aligned with topography. For example, one can chose a set of z-levels, {zk+12k=0,1,,N} where z12=-hmax is chosen to be the maximum depth and zN+12=0 is the unperturbed free surface.2 Then, starting from the bottom, for k = 0 setz12(x,y)=-h(x,y)and for each k = 1,  , N−1 setzk+12(x,y)=maxzk+12,zk-12(x,y)+Δzmin,where Δzmin is the chosen minimal vertical grid spacing (n.b., to avoid surfacing of coordinate isolines, Δzmin  hmin/N, where hmin is the minimal depth). In principle, Δzmin may be chosen as infinitely small, so the resultant system is equivalent to a z-coordinate with the necessity of handling the layers near the bottom as “mass-less” layers. Its disadvantage is the nonsmooth transition from z-level to topography-following regions. This nonsmoothness can be repaired by applying 2D diffusion to zk+12 with a variable diffusivity coefficient—zero for the bottom and increasing toward the surface. The resultant coordinate systems are shown in Fig. 1b–d for the cases of two different degrees of smoothing: case (b) is closer to the stretched σ-coordinate, while (c) retains more features of the original z-coordinate (d). Unlike the (σ-coordinate, in both cases the coordinate surfaces near the top are almost horizontal and have less resemblance to the topography.

Throughout this study we assume that our vertical system of coordinates is no longer separable in the sense that it cannot be generated by the simple relationships (1.5) where S(σ) is independent of horizontal coordinates, but involves a full three-dimensional (3D) transformation (1.2).

Consequently, the applicability of the methods developed here is not limited to just a σ- or S-coordinate class of models.

Discretization of vertical coordinate introduces a set of coordinate surfaces,{zk+12=zk+12(x,y),k=0,1,,N}.

If the ocean is at rest, the free-surface elevation is ζ = 0, hence zN+12=0, and the whole set corresponding to zero free-surface {zk+12(0)} is referred as an unperturbed coordinate system. In the case of a nonzero ζ, all zk+12 are displaced by a distance proportional to ζ and the distance from the bottom as the fraction of unperturbed local depthzk+12=zk+12(0)+ζ1+zk+12(0)h(recall that z12z12(0)-h and zN+12ζ). As a result the perturbed grid-box height, Δzkzk+12-zk-12, is related to the unperturbed height, Δzk(0)zk+12(0)-zk-12(0) according toΔzk=Δzk(0)1+ζh,where the multiplier (1 + ζ/h) is independent of the vertical coordinate. This choice is similar to Higdon and Bennett (1996) and Higdon and de Szoeke (1997), with the exception that they applied it for an isopycnic coordinate model. But it is different from Killworth et al. (1991) and Dukowitz and Smith (1994), where free-surface elevation affects only the top-most grid box, as well as from (Song and Haidvogel, 1994), where each grid box receives the same increment (hence Δzk=Δzk(0)+ζ/N) regardless of its unperturbed size Δzk(0). Later we show that (1.10) has several consequences, including the fact that vertical mass fluxes generated by a purely barotropic motion vanish identically at every interface, zk+12.

Combining the tracer equation in advective formqt+(u·)q=0with the nondivergence equation,(·u)=0,we derive the tracer equation in conservation form,qt+·(uq)=0.As a consequence of (1.11), if the tracer is specified as a spatially uniform field at the initial time, it remains so regardless of the velocity field. On the other hand, as a consequence of (1.13), the volume integral of the tracer concentration is conserved in the absence of incoming and outgoing fluxes across the domain boundary. The continuity equation (1.12) provides the compatibility condition between these two properties. Both properties are valuable and should be considered in constructing numerical oceanic models.

The discretization of (1.13) is usually done using a finite-volume approachΔVi,j,kn+1qi,j,kn+1=ΔVi,j,knqi,j,kn-Δtq˜i+12,j,kUi+12,j,k-q˜i-12,j,kUi-12,j,k+q˜i,j+12,kVi,j+12,k-q˜i,j-12,kVi,j-12,k+q˜i,j,k+12Wi,j,k+12-q˜i,j,k-12Wi,j,k-12,where qi,j,k is understood as a volume-averaged concentration over the grid-box ΔVi,j,k,qi,j,k=1ΔVi,j,knΔVi,j,knqdV.The q˜i+12,j,k (q with one index half-integer) are the interfacial values of tracer concentration. Uppercase3 Ui+12,j,k,,Vi,j+12,k, and Wi,j,k+12 are volumetric fluxes4 in the two horizontal and vertical directions. These are defined as velocity components multiplied by the contact area between two adjacent grid boxesUi+12,j,k=ui+12,j,kΔzi+12,j,kΔηi+12,jVi,j,+12,k=vi,j+12,kΔzi,j+12,kΔξi,j+12,where Δzi+12,j,k,Δηi,j+12, and Δzi,j+12,k,Δξi,j+12, are vertical and horizontal measures of the corresponding grid-box interfaces (Δε, Δη are assumed to be nonuniform because of curvilinear horizontal coordinates). The superscripts n + 1 and n denote new and old time steps. The time step for the flux variables in (1.14) is not specified yet (must be effectively at n + 1/2 to achieve the second-order temporal accuracy), but the flux form by itself guarantees exact conservation of the global volume integral of the advected quantity as long as there is no net flux across the domain boundary. Setting qi,j,k  1 in (1.14) yields the discretized continuity equation,ΔVi,j,kn+1=ΔVi,j,kn-Δt·Ui+12,j,k-Ui-12,j,k+Vi,j+12,k-Vi,j-12,k+Wi,j,k+12-Wi,j,k-12.Once it holds, the conservative form of the discrete tracer equation (1.14) also has the property of constancy preservation in addition to global content conservation.

In a hydrostatic model the discrete continuity equation (1.17) is needed to compute vertical velocity rather than grid-box volume ΔVi,j,kn+1. (The latter is entirely controlled by change of ζ via (1.10).) Hence,Wi,j,12=0attheseafloorandWi,j,k+12=-k=1kΔVi,j,kn+1-ΔVi,j,knΔt+Ui+12,j,k-Ui-12,j,k+Vi,j+12,k-Vi,j-12,kforallk=1,2,,N,which, in fact, defines the meaning of Wi,j,k+12 as a finite-volume flux across the moving grid-box interface zi,j,k+12. Vertical summation of (1.17) for different k leads to the equation for the free surface,ζi,jn+1=ζi,jn-ΔtΔAi,jU¯i+12,j-U¯i-12,j+V¯i,j+12-V¯i,j-12,where ΔAi,j is the horizontal area of the grid box i,j;U¯i+12,j=k=1NUi+12,j,k,V¯i,j+12=k=1NVi,j+12,jare vertically integrated (barotropic) volume fluxes; and we have used the identity(ζi,j+hi,j)·ΔAi,jk=1NΔVi,j,k,where hi,j is independent of time. Obviously, setting k = N in (1.19), consistently with (1.20), (1.21), (1.22) results inWi,j,N+12=0,as required by the kinematic boundary condition at the free surface.

Thus far we have assumed that the time step and time-stepping algorithm for the tracer (1.14) and for ζ (1.20) are the same. This would be the case if the barotropic and baroclinic components were advanced using the same small time step dictated by the stability criterion for the barotropic mode; or if the barotropic mode were treated implicitly with a special care to construct finite-volume fluxes Ui+12,j,k, Vi,j+12,k, and Wi,j,k+12 such that the (1.17) holds exactly and is compatible with (1.21), (1.23), (Dukowitz and Smith, 1994). In a split-explicit, free-surface model (cf., Blumberg and Mellor, 1987, Killworth et al., 1991), the equation for free-surface (1.20) and the vertically integrated (2D) momenta are advanced using a much smaller time step than the tracer equations. Each baroclinic time step starts with computation of the r.h.s of the 3D momentum equations. The r.h.s components are integrated vertically to provide forcing terms for the barotropic mode. During the barotropic time stepping, the free surface and the barotropic velocity components are averaged over the sequence of the barotropic steps and the averaged values are feed back into the 3D momenta. The averaging is needed to prevent temporal aliasing of the signals resolved by the barotropic, but not by the baroclinic step, and [in some models, cf., Nadiga et al., 1997, Hallberg, 1997, Higdon and de Szoeke, 1997, Higdon, 2002, where no dealiazing (averaging) of ζ-equation is actually performed] to provide vertically integrated fluxes consistent with finite-time-step “baroclinic time” free-surface equation (1.20). Then the 3D momenta are advanced to the baroclinic time step n + 1 (with violation of the external mode CFL criterion), and vertical integrals of the new fields are subtracted from the similar values from the barotropic submodel. The resultant differences is then uniformly distributed throughout the vertical column to make sure that the corrected 3D velocity components have the same vertical integrals as the barotropic ones. At the same time, free surface ζ at the new baroclinic step is assigned to its new state from the the barotropic submodel.

Perhaps the most delicate matter here is the replacement of free-surface ζ at n + 1 with its fast-time-averaged value: not doing so leaves room for aliasing error, while the replacement makes the “slow-time” discrete 2D continuity equation (1.20) hold only within the order of temporal accuracy, but no longer exactly (even thought it is exact at every fast time step). Consequently, it is no longer possible to reconstruct vertical velocity via (1.19) in such a way that the top kinematic boundary condition (1.23) is respected.5 As the result, a conservative update of the tracer fields (1.14) looses its constancy preservation property.

Section snippets

Accuracy and stability of time-stepping algorithms

Table 4 in Griffies et al. (2000), provides a comprehensive overview of time-stepping and mode-splitting algorithms for virtually all oceanic models currently in use. Despite the large diversity of models, the time-stepping algorithms are mainly limited to applications of classical methods––Leap-Frog (LF), Adams–Bashforth (AB), and Forward–Backward (FB, used almost exclusively for the barotropic mode). In this section we will show that, for oceanic modeling specifically where the time step is

Barotropic mode

In this section we address specific aspects of the barotropic mode as part of a coupled barotropic–baroclinic system.

A hybrid predictor–corrector for the baroclinic time step

Because of the mathematical similarity of their equations, time-stepping algorithms for the baroclinic mode are generally similar to the barotropic ones. The differences arise from the necessity for a conservative, constancy-preserving algorithm for tracers. We will show how this makes it necessary to update the velocities before the tracers. Similar to (3.48), the allowed time step is limited mainly by internal gravity waves, but the contrast with the Coriolis restriction is not so dramatic.

Time stepping the coupled baroclinic–barotropic system

We now summarize the time-stepping algorithm in ROMS, focusing on the discrete-time interactions between the modes (Fig. 18).

  • Stage 1: Compute the right-hand side for the 3D momentum equations at time step n (i.e., pressure-gradient, Coriolis, and advective terms only; no viscous terms are computed at this time). Apply this right-hand side to advance the 3D momenta using an LF step combined with a half-step, backward interpolation with AM3-like coefficients (the result is time-centered at n+12).

Conclusions

We have designed a robust computational kernel for a split-explicit, terrain-following-coordinate oceanic model, ROMS.

We use time-stepping algorithms with forward–backward feedback between the pairs of variables responsible for gravity wave propagation (surface or internal) that combine an extended range of stability with the temporal accuracy of the best known algorithms: in effect, generalizing the FB schemes to higher orders of accuracy. Among these schemes, the Euler step (2.13) can be

Acknowledgment

This study was supported by grant N00014-02-1-0236 from the Office of Naval Research.

References (41)

  • R. Asselin

    Frequency filter for time integations

    Month. Weath. Rev.

    (1972)
  • A. Beckmann

    The representation of bottom boundary layer processes in numerical ocean circulation models

  • R. Bleck et al.

    A wind-driven isopycnic coordinate model of the north and equatorial Atlantic Ocean: 1. Model development and supporting experiments

    J. Geophys. Res.

    (1990)
  • Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean circulation model. In: Heaps, N....
  • J.R. Brown et al.

    An economical time difference system for numerical weather prediction

    Month. Weath. Rev.

    (1978)
  • C. Canute et al.

    Spectral Methods in Fluid Mechanics

    (1988)
  • D.E. Dietrich et al.

    A high resolution numerical study of Gulf of Mexico fronts and Eddies

    Meteorol. Atmos. Phys.

    (1997)
  • J.K. Dukowitz et al.

    A reformulation and implementation of the Bryan–Cox–Semtner ocean model on the connection machine

    J. Atmos. Oceanic Technol.

    (1993)
  • J.K. Dukowitz et al.

    Implicit free-surface method for the Bryan–Cox–Semtner ocean model

    J. Geophys. Res.

    (1994)
  • D.R. Durran

    The third-order Adams–Bashforth method: an attractive alternative to leapfrog time differencing

    Month. Weath. Rev.

    (1991)
  • Cited by (0)

    View full text