Vertical S-coordinate

These equations are as given in the ROMS documentation

REMORA has a generalized vertical, terrain-following, coordinate system. Currently, the vertical transformation equation \(z=z \left( x,y, \sigma ,t \right)\), (Vtransform=2 from ROMS) is available which can support vertical stretching 1D-functions when several constraints are satisfied.

Transformation Equation

\[\begin{split}z\left( x,y,\sigma ,t\right) &=\zeta \left( x,y,t\right) +\left[ \zeta \left( x,y,t\right) +h\left( x,y\right) \right] S\left( x,y,\sigma \right) ,\\ S\left( x,y,\sigma \right) &=\frac{h_c\ \sigma +h \left( x,y\right) C\left( \sigma \right) }{h_c +h\left( x,y\right)}\end{split}\]

where \(S\left( x,y,\sigma \right)\) is a nonlinear vertical transformation functional, \(\zeta \left( x,y,t\right)\) is the time-varying free-surface, \(h\left( x,y\right)\) is the unperturbed water column thickness and \(z=-h\left( x,y\right)\) corresponds to the ocean bottom, \(\sigma\) is a fractional vertical stretching coordinate ranging from \(-1\leq \sigma \leq 0,\ C\left( \sigma \right)\) is a nondimensional, monotonic, vertical stretching function ranging from \(-1\leq C\left( \sigma \right) \leq 0\), and \(h_c\) is a positive thickness controlling the stretching.

We find it convenient to define:

\[H_z \equiv \frac{\partial z}{\partial \sigma }\]

where \(H_z = H_z \left( x,y,\sigma ,t\right)\) are the vertical grid thicknesses. In REMORA, \(H_z\) is computed discretely as \(\Delta z/\Delta \sigma\) since this leads to the vertical sum of \(H_z\) being exactly the total water column thickness \(D\).

Notice that,

\[\begin{split}S\left( x,y,\sigma \right) = \begin{cases} 0, & \text{if } \sigma = 0, & C\left( \sigma \right) = 0, & \text{at the free-surface;}\\ -1 , & \text{if } \sigma = -1, & C\left( \sigma \right) = -1, & \text{at the ocean bottom} \end{cases}\end{split}\]

This transformation offers several advantages:

Regardless of the design of \(C\left( \sigma \right)\), it behaves like equally-spaced sigma-coordinates in shallow regions, where \(h\left( x,y\right) \ll h_c\). This is advantageous because it avoids excessive resolution and associated CFL limitation in such areas.
Near-surface refinement behaves more or less like geopotential coordinates in deep regions (level thicknesses, \(H_z\), do not depend or weakly depend on bathymetry), while near-bottom like sigma coordinates (\(H_z\) is roughly proportional to depth). This reduces the extreme r-factors near the bottom and reduces pressure gradient errors.
The true sigma-coordinate system is recovered as \(h_c \rightarrow \infty\). This is useful when configuring applications with flat bathymetry and uniform level thickness. Practically, you can achieve this by setting tcline to 1.0d+16 in DataStruct.H. This will set \(h_c =1.0\times 10^{16}\). Although not necessary, we also recommend that you set \(\theta _S = 0\) and \(\theta _B =0\).

In an undisturbed ocean state, \(\zeta \equiv 0\), this transformation yields the following unperturbed depths, \(\hat{z}\),

\[\hat{z} \left( x,y,\sigma \right) \equiv h\left( x,y\right) S\left( x,y,\sigma \right) =h\left( x,y\right) \left[ \frac{h_c\ \sigma +h\left( x,y\right) C\left( \sigma \right)}{h_c +h\left( x,y\right) } \right]\]

and

\[d\hat{z} =d\sigma \ h\left( x,y\right) \left[ \frac{h_c}{h_c +h\left( x,y\right) } \right]\]

As a consequence, the uppermost grid box retains very little dependency from bathymetry in deep areas, where \(h_c \ll h\left( x,y\right)\). For example if \(h_c =250m\) and \(h\left( x,y\right)\) changes from \(2000\) to \(6000m\), the uppermost grid box changes only by a factor of 1.08 (less than 10%).

Vertical Stretching

The above generic vertical transformation design facilitates numerous vertical stretching functions \(C \left( \sigma \right)\) . This function is defined in terms of several parameters which are specified in the standard input file, inputs.

Stretching Function Properties

\(C\left( \sigma \right)\) is a dimensionless, nonlinear, monotonic function;
\(C\left( \sigma \right)\) is a continuous differentiable function, or a differentiable piecewise function with smooth transition;
\(C\left( \sigma \right)\) must be discretized in terms of the fractional stretched vertical coordinate \(\sigma\),

\[\begin{split}\sigma \left( k \right) = \begin{cases} \frac{k-N}{N}, & \text{at vertical }W\text{-points}, & k=0,\ldots ,N,\\ \frac{k-N-0.5}{N}, & \text{at vertical}\rho \text{-points}, & k=1,\ldots ,N \end{cases}\end{split}\]
\(C\left( \sigma \right)\) must be constrained by \(-1 \leq C\left( \sigma \right) \leq 0\), that is,

\[\begin{split}C\left( \sigma \right) = \begin{cases} 0, & \text{if } \sigma = 0, & C\left( \sigma \right) = 0, & \text{at the free-surface};\\ -1, & \text{if } \sigma = -1, & C\left( \sigma \right) = -1, & \text{at the ocean bottom}. \end{cases}\end{split}\]

Available Stretching Functions

(A. Shchepetkin (2010) Vstretching=4 from ROMS). \(C\left( \sigma \right)\) is defined as a continuous, double stretching function:

Surface refinement function:

\[C\left( \sigma \right) = \frac{1-\cosh \left( \theta _S \sigma \right) }{\cosh \left( \theta _S \right) -1}, \quad \text{for } \theta _S > 0, \quad C\left( \sigma \right) = - \sigma ^2, \quad \text{for }\theta _S \leq 0\]

Bottom refinement function:

\[C\left( \sigma \right) = \frac{\exp \left( \theta _B C\left( \sigma \right) \right) -1}{1-\exp \left( -\theta _B \right) }, \quad \text{for }\theta _B >0\]

Notice that the bottom function is the second stretching of an already stretched transform. The resulting stretching function is continuous with respect to \(\theta _S\) and \(\theta _B\) as their values approach zero. The range of meaningful values for \(\theta _S\) and \(\theta _B\) are:

\[0\leq \theta _S \leq 10 \quad \text{and} \quad 0\leq \theta _B \leq 4\]

However, users need to pay attention to extreme r-factor (rx1) values near the bottom.