Deviatory Stress Tensor (Horizontal Viscosity)
The horizontal components of the divergence of the stress tensor ( Wajsowicz, 1993) in nondimesional, orthogonal curvilinear coordinates (, , ) with dimensional, spatially-varying metric factors (, , ) and velocity components (, , ) are given by:
and the strain field is:
Here, and are the spatially varying horizontal and vertical viscosity coefficients, respectively, and is another (very small, often neglected) horizontal viscosity coefficient. Notice that because of the generalized terrain-following vertical coordinates of ROMS, we need to transform the horizontal partial derivatives from constant z-surfaces to constant s-surfaces. And the vertical metric or level thickness is the Jacobian of the transformation, . Also in these models, the vertical velocity is computed as and has units of .
Transverse Stress Tensor
Assuming transverse isotropy, as in Sadourny and Maynard (1997) and Griffies and Hallberg (2000), the deviatoric stress tensor can be split into vertical and horizontal sub-tensors. The horizontal (or transverse) sub-tensor is symmetric, it has a null trace, and it possesses axial symmetry in the local vertical direction. Then, transverse stress tensor can be derived from (1) and (2) yielding
Notice the flux form of (5) and the symmetry between the and terms which are defined at density points on a C-grid. Similarly, the and terms are symmetric and defined at vorticity points. These staggering positions are optimal for the discretization of the tensor; it has no computational modes and satisfy first-moment conservation.
The biharmonic friction operator can be computed by applying twice the tensor operator (5), but with the squared root of the biharmonic viscosity coefficient (Griffies and Hallberg, 2000). For simplicity and momentum balance, the thickness appears only when computing the second harmonic operator as in Griffies and Hallberg (2000).
Rotated Transverse Stress Tensor
In some applications with tall and steep topography, it will be advantageous to reduce substantially the contribution of the stress tensor (5) to the vertical mixing when operating along constant -surfaces. The transverse stress tensor rotated along geopotentials (constant depth) is, then, given by
Notice that transverse stress tensor remains invariant under coordinate transformation. The rotated tensor (7) retains the same properties as the unrotated tensor (5). The additional terms that arise from the slopes of -surfaces along geopotentials are discretized using a modified version of the triad approach of Griffies et al. (1998).
The Laplacian of a scalar in curvilinear coordinates is:
In ROMS, this term is multiplied by and becomes
where is any tracer. This form guarantees that the term does not contribute to the volume-integrated equations.
The biharmonic operator is ; the corresponding term is computed using a temporary variable :
where is once again any tracer and is the square root of the input value so that it can be applied twice.
Rotated mixing tensors
Both the Laplacian and biharmonic terms above operate on surfaces of constant and can contribute substantially to the vertical mixing. However, the oceans are thought to mix along constant density surfaces so this is not entirely satisfactory. Therefore, the option of using rotated mixing tensors for the Laplacian and biharmonic operators has been added. Options exist to diffuse on constant surfaces (MIX_GEO_TS) and constant potential density surfaces (MIX_ISO_TS).
The horizontal Laplacian diffusion operator is computed by finding the three components of the flux of the quantity . The and components are locally horizontal, rather than along the surface. The diffusive fluxes are:
and there is some trickery such that the computations depend on the sign of and of . No flux boundary conditions are easily imposed by setting
Finally, the flux divergence is calculated and is added to the right-hand-side term for the field being computed:
The biharmonic rotated mixing tensors are computed much as the non-rotated biharmonic mixing. We define a temporary variable based on equation (16):
We then build up fluxes of as in equations (14). We then apply equation (16) to these fluxes to obtain the biharmonic mixing tensors. Again, the value of is the square root of that read in so that it can be applied twice.
Guidelines for Coefficient Values
The horizontal viscosity and diffusion coefficients are scalars which are read in from ocean.in. Several factors to consider when choosing these values are:
- spindown time The spindown time on wavenumber is for the Laplacian operator and for the biharmonic operator. The smallest wavenumber corresponds to the length and is , leading to
This time should be short enough to damp out the numerical noise which is being generated but long enough on the larger scales to retain the features you are interested in. This time should also be resolved by the model timestep.
- boundary layer thickness The western boundary layer has a thickness proportional to:
for the Laplacian and biharmonic viscosity, respectively. We have found that the model typically requires the boundary layer to be resolved with at least one grid cell. This leads to coarse grids requiring large values of .
We have chosen anything from zero to the value of the horizontal viscosity for the horizontal diffusion coefficient. One common choice is an order of magnitude smaller than the viscosity.