Numerical Solution Technique
Horizontal Discretization
In the horizontal, the REMORA governing equations are discretized over a boundary-fitted, orthogonal curvilinear coordinates \(\left(\xi,\eta\right)\) grid. The general formulation of the curvilinear coordinates system allows Cartesian, polar, and spherical coordinates applications. The transformation of any of these coordinates to REMORA \(\left(\xi,\eta\right)\) grid is specified in the metric terms (pm, pn).
The model state variables are staggered using an Arakawa C-grid. As illustrated below, the free-surface (zeta), density (rho), and active/passive tracers (t) are located at the center of the cell whereas the horizontal velocity (u and v) are located at the west/east and south/north edges of the cell, respectively. That is, the density is evaluated between points where the currents are evaluated.
Staggered Horizontal Grid
In REMORA all the state arrays are dimensioned the same size to facilitate parallelization. However, the computational ranges for all the state variables are:
Grid Cell
Variable |
Interior Range |
Full Range |
|---|---|---|
\(\rho\text{-type}\) |
1: |
0: |
\(\psi\text{-type}\) |
2: |
1: |
\(\text{u-type}\) |
2: |
1: |
\(\text{v-type}\) |
1: |
0: |
Vertical Discretization
The REMORA governing equations are discretized over variable topography using a stretched, terrain-following, vertical coordinate. As a result, each grid cell may have different level thickness (Hz) and volume. The model state variables are vertically staggered so that horizontal momentum (u, v), (rho), and active/passive tracers (t) are located at the center of the grid cell. The vertical velocity (omega, w) and vertical mixing variables (Akt, Akv, etc) are located at the bottom and top faces of the cell. See diagram below.
Vieste-Dubrovnik Transect
Staggered Vertical Grid
In this diagram, indices are 1-indexed (as in ROMS), while the indices in REMORA are 0-indexed.
The total thickness of the water column is \(\zeta\left(i,j\right)+h\left(i,j\right)\). The bathymetry (h) is usually time invariant whereas the free-surface (zeta) evolves in time. At input and output, the bathymetry is always a positive quantity. However, the depths z_r(i,j,k) and z_w(i,j,k) are negative for all locations below the mean sea level.
Grid Variables
Variable |
Variable in Code |
Definition |
Location |
Origin |
|---|---|---|---|---|
m |
x-dir coordinate |
corners |
||
n |
y-dir coordinate |
corners |
||
\(\xi\) |
transformed x-dir orthogonal curvilinear coordinate |
corners |
||
\(\eta\) |
transformed y-dir orthogonal curvilinear coordinate |
corners |
||
\(\zeta\) |
free-surface |
center |
||
\(\rho\) |
|
density |
center |
equation of state |
\(t\) |
active/passive tracers |
center |
||
\(u\) |
|
x-dir horizontal velocity |
west/east faces |
|
\(v\) |
|
y-dir horizontal velocity |
north/south faces |
|
\(\overline{u}\) |
|
x-dir vertically integrated momentum |
west/east faces |
|
\(\overline{v}\) |
|
y-dir vertically integrated momentum |
west/east faces |
|
\(\psi\) |
corners |
|||
\(H_z\) |
|
level thickness |
center |
|
\(\omega\) |
vertical velocity |
bottom/top faces |
||
\(w\) |
vertical velocity |
bottom/top faces |
||
Akt |
|
vertical mixing |
bottom/top faces |
|
Akv |
vertical mixing |
bottom/top faces |
||
\(h\) |
|
bathymetry (always positive) |
||
\(z_{r\left(i,j,k\right)}\) |
depth (negative below sea level) |
center |
||
\(z_{w\left(i,j,k\right)}\) |
depth (negative below sea level) |
bottom/top faces |
||
\(T\) |
|
temperature |
||
\(S\) |
|
salinity |
||
\(D_{crit}\) |
critical depth |
|||
\(D\) |
total water depth |
|||
\(C\) |
concentration |
|||
\(\delta_{\xi}\) |
centered finite-difference approximation of \(\Delta\xi\) |
|||
\(\delta_{\eta}\) |
centered finite-difference approximation of \(\Delta\eta\) |
|||
\(\delta_{\sigma}\) |
centered finite-difference approximation of \(\Delta\sigma\) |
|||
\(\Delta_{\xi}\) |
transformed x-dir differential distance |
|||
\(\Delta_{\eta}\) |
transformed y-dir differential distance |
|||
\(\Delta_{\sigma}\) |
|
transformed z-dir differential distance |
||
\(\Delta_{z}\) |
vertical distance from one \(\rho\) to another |
|||
\(\overline{\left(\qquad\right)}^{\xi}\) |
average taken over \(\Delta\xi\) |
|||
\(\overline{\left(\qquad\right)}^{\eta}\) |
average taken over \(\Delta\eta\) |
|||
\(\overline{\left(\qquad\right)}^{\sigma}\) |
average taken over \(\Delta\sigma\) |
|||
\(\Delta V\) |
grid box volume |
Conservation Properties
From Shchepetkin and McWilliams (2005), we have a tracer concentration equation in advective form:
and also a tracer concentration equation in conservation form:
The continuity equation:
can be used to get from one tracer equation to the other. As a consequence of eq. (1), if the tracer is spatially uniform, it will remain so regardless of the velocity field (constancy preservation). On the other hand, as a consequence of (2), the volume integral of the tracer concentration is conserved in the absence of internal sources and fluxes through the boundary. Both properties are valuable and should be retained when constructing numerical ocean models.
The semi-discrete form of the tracer equation is:
Here \(\delta_{\xi},\delta_{\eta}\) and \(\delta_{\sigma}\) denote simple centered finite-difference approximations to \(\partial/\partial\xi,\partial/\partial\eta\) and \(\partial/\partial\sigma\) with the differences taken over the distances \(\Delta\xi,\Delta\eta\) and \(\Delta\sigma\), respectively. \(\Delta z\) is the vertical distance from one \(\rho\) point to another. \(\overline{\left(\qquad\right)}^{\xi}, \overline{\left(\qquad\right)}^{\eta}\) and \(\overline{\left(\qquad\right)}^{\sigma}\) represent averages taken over the distances \(\Delta\xi,\Delta\eta\) and \(\Delta\sigma\).
The finite volume version of the same equation is no different, except that a quantity \(C\) is defined as the volume-averaged concentration over the grid box \(\Delta V\):
The quantity \(\left(\dfrac{u\overline{H_z}^{\xi}\overline{C}^{\xi}}{\overline{n}^{\xi}}\right)\) is the flux through an interface between adjacent grid boxes.
This method of averaging was chosen because it internally conserves first moments in the model domain, although it is still possible to exchange mass and energy through the open boundaries. The method is similar to that used in Arakawa and Lamb; though their scheme also conserves enstrophy. Instead, we will focus on (nearly) retaining constancy preservation while coupling the barotropic (depth-integrated) equations and the baroclinic equations.
The timestep in eq. (4) is assumed to be from time \(n\) to time \(n+1\), while the other terms being evaluated at time \(n+1/2\) for second-order accuracy. Setting \(C\) to \(1\) everywhere reduces eq. (4) to:
If this equation holds true for the step from time \(n\) to time \(n+1\), then constancy preservation will hold.
In a hydrostatic model such as REMORA, the discrete continuity equation is needed to compute vertical velocity rather than grid-box volume \(\dfrac{H_z}{mn}\) (the latter is controlled by changes in \(\zeta\) in the barotropic mode computations). Here, \(\dfrac{H_z\Omega}{mn}\) is the finite-volume flux across the moving grid-box interface, vertically on the \(w\) grid.
The vertical integral of the continuity equation (5), using the vertical boundary conditions on \(\Omega\), is:
where \(\zeta\) is the surface elevation, \(D=h+\zeta\) is the total depth, and \(\overline{u},\overline{v}\) are the depth-integrated horizontal velocities. This equation and the corresponding 2-D momentum equations are timestepped on a shorter timestep than eq.(4) and the other 3-D equations. Due to the details in the mode coupling, it is only possible to maintain constancy preservation to the accuracy of the barotropic timesteps.
Depth-Integrated Equations
The depth average of a quantity \(A\) is given by:
where the overbar indicates a vertically averaged quantity and
is the total depth of the water column. The vertical integral of the momentum equations are:
and
where \(\phi_2\) includes the \(\frac{\partial z}{\partial\xi}\) term, \(\overline{\mathcal{D}}_{h_u}\) is the horizontal viscosity, and the vertical viscosity only contributes through the upper and lower boundary conditions. We also need the vertical integral of the continuity equation, shown above as eq. (6). The presence of a free surface introduces waves which propagate at a speed of \(\sqrt{gh}\). These waves usually impose a more severe time-step limit than any of the internal processes. We have therefore chosen to solve the full equations by means of a split time step. In other words, the depth integrated equations (9), (10), and (6) are integrated using a short time step and the values of \(\overline{u}\) and \(\overline{v}\) are used to replace those found by integrating the full equations on a longer time step. A diagram of the barotropic time stepping is shown here:
Some of the terms in equations (9) and (10) are updated on the short time step while others are not. The contributions from the slow terms are computed once per long time step and stored. If we call these terms \(R_{u_{\text{slow}}}\) and \(R_{v_{\text{slow}}}\), equations (9) and (10) become:
and
When time stepping the model, we compute the right-hand-sides for the full 3-D momentum equations as well as the right-hand-sides for equations (11) and (12). The vertical integral of the 3-D right-hand-sides are obtained and then the 2-D right-hand-sides are subtracted. The resulting fields are the slow forcings \(R_{u_{\text{slow}}}\) and \(R_{v_{\text{slow}}}\). This was found to be the easiest way to retain the baroclinic contributions of the non-linear terms such as \(\overline{uu}-\overline{u}\,\overline{u}\). The model is time stepped from time \(n\) to time \(n+1\) by using short time steps on equations (11), (12) and (6). Equation (6) is time stepped first, so that an estimate of the new \(D\) is available for the time rate of change terms in equations (11) and (12). A third-order predictor-corrector time stepping is used. In practice, we actually time step all the way to time \(\left(n+\textbf{dtfast}\times M^*\right)\) and while maintaining weighted averages of the values of \(\overline{u},\overline{v}\) and \(\zeta\). The averages are used to replace the values at time \(n+1\) in both the baroclinic and barotropic modes, and for recomputing the vertial grid spacing \(H_z\). The following figure shows one option for how these weights might look:
The primary weights, \(a_m\), are used to compute \(\langle\zeta\rangle^{n+1}\equiv\overunderset{M^*}{m=1}{\sum}a_m\zeta^m\). There is a related set of secondary weights \(b_m\), used as \(\langle\!\langle\overline{u}\rangle\!\rangle^{n+\frac{1}{2}}\equiv\overunderset{M^*}{m=1}{\sum}b_m\overline{u}^m\). In order to maintain constancy preservation, this relation must hold:
Shchepetkin and McWilliams (2005) introduce a range of possible weights, but the ones used here have a shape function:
where \(p,q\) are parameters and \(A_0,\tau_0\), and \(r\) are chosen to satisfy normalization, consistency, and second-order accuracy conditions,
using Newton iterations. \(\tau^*\) is the upper limit of \(\tau\) with \(A\left(\tau\right)\geq0\). In practice we initially set
compute \(A\left(\tau\right)\) using eq.(14), normalize using:
and adjust \(r\) iteratively to satisfy the \(n=2\) condition of (15). We are using values of \(p=2,q=4\), and \(r=0.284\). This form allows some negative weights for small \(m\), allowing \(M^*\) to be less than \(1.5M\).
Pressure Gradient Terms in Mode Coupling
Equation (11) contains the term \(R_{u_{\text{slow}}}\), computed as the difference between the 3-D right-hand-side and the 2-D right-hand-side. The pressure gradient therefore has the form:
where the term in square brackets is the mode coupling term and is held fixed over all the barotropic steps and
is the vertically integrated pressure gradient. The latter is a function of the bathymetry, free surface gradient, and the free surface itself, as well as the vertical distribution of density.
The disadvantage of this approach is that after the barotropic time stepping is complete and the new free surface is substituted into the full baroclinic pressure gradient, its vertical integral will no longer be equal to the sum of the new surface slope term and the original coupling term based on the old free surface. This is one form of mode-splitting error which can lead to trouble because the vertically integrated pressure gradient is not in balance with the barotropic mass flux.
Instead, let us define the following:
Changing the vertical coordinate to \(\sigma\) yields:
which implies that \(\overline{\rho}\) and \(\rho^*\) are actually independent of \(\zeta\) as long as the density profile \(\rho=\rho\left(\sigma\right)\) does not change. The vertically integrated pressure gradient becomes:
In the case of uniform density \(\rho_0\), we obtain \(\rho^*\equiv\overline{\rho}\equiv\rho_0\), but we otherwise have two new terms. The accuracy of these terms depends on an accurate vertical integration of the density, as described in Shchepetkin and McWilliams (2005).
Horizontal and Vertical Advection
The advection of a tracer \(C\) has an equation of the form
where we have introduced the advective fluxes:
Fourth-order Centered
The barotropic advection is centered fourth-order. Create gradient terms:
The fluxes now become:
Third-order Upwind
There is a class of third-order upwind advection schemes, both one-dimensional (Leanord, 1979) and two-dimensional (Rasch, 1994 and Shchepetkin and McWilliams, 1998). This scheme is known as UTOPIA (Uniformly Third-Order Polynomial Interpolation Algorithm). Applying flux limiters to UTOPIA is explored in Thuburn (1995), although it is not implemented in REMORA. The two-dimensional formulation in Rasch contains terms of order \(u^2C\) and \(u^3C\), including cross terms (\(uvC\)). The terms which are nonlinear in velocity have been dropped in REMORA, leaving one extra upwind term in the computation of the advective fluxes:
The second derivative terms are centered on a \(\rho\) point in the grid, but are needed at a \(u\) or \(v\) point in the flux. The upstream value is used:
The value of \(\gamma\) in the model is \(\frac{1}{8}\) while that in Rasch(1994) is \(\frac{1}{6}\).
Because the third-order upwind scheme is designed to be two-dimensional, it is not used in the vertical (though one might argue that we are simply performing one-dimensional operations here). Instead we use a centered fourth-order scheme in the vertical when the third-order upwind option is turned on:
One advantage of UTOPIA over MPDATA is that it can be used on variables having both negative and positive values. Therefore, it can be used on velocity as well as scalars (is there a reference for this?). For the \(u\)-velocity, we have:
while for the \(v\)-velocity we have:
In all these terms, the second derivatives are evaluated at an upstream location.
Vertical Velocity
Having obtained a complete specification of the \(u,v,T\), and \(S\) fields at the next time level by the methods outlined above, the vertical velocity and density fields can be calculated. the vertical velocity is obtained by rewriting equation (5) as:
and combining it with equation (6) to obtain:
Solving for \(H_z\Omega/mn\) and using the semi-discrete notation we obtain:
The integral is actually computed as a sum from the bottom upwards and also as a sum from the top downwards. The value used is a linear combination of the two, weighted so that the surface down value is used near the surface while the other is used near the bottom. [Is this still done?]
Equation of State
The density is obtained from temperature and salinity via an equation of state. REMORA provides a choice of a nonlinear equation of state \(\rho=\rho\left(T,S,z\right)\) or a linear equation of state \(\rho=\rho\left(T\right)\). The nonlinear equation of state has been modified and now corresponds to the UNESCO equation of state as derived by Jackett and McDougall (1995). It computes in situ density as a function of potential temperature, salinity and pressure.
Warning: although we have used it quite extensively in the past, McDougall (personal communication) claims that the single-variable \(\left(\rho=\rho\left(T\right)\right)\) equation of state is not dynamically appropriate as is. He has worked out the extra source and sink terms required, arising from vertical motions and the compressibility of water. They are quite complicated and we have not implemented them to see if they alter the flow.