Horizontal Mixing

Horizontal Mixing

Deviatory Stress Tensor (Horizontal Viscosity)

The horizontal components of the divergence of the stress tensor ( Wajsowicz, 1993) in nondimesional, orthogonal curvilinear coordinates ($\xi$, $\eta$, $s$) with dimensional, spatially-varying metric factors (${\frac {1}{m}}$, ${\frac {1}{n}}$, $H_{z}$) and velocity components ($u$, $v$, $\omega H_{z}$) are given by:

$${\begin{aligned}F^{u}\equiv {\widehat {\xi }}\cdot \left(\nabla \cdot {\vec {\sigma }}\right)={\frac {mn}{H_{z}}}{\Biggl [}&{\frac {\partial }{\partial \xi }}{\Biggl (}{\frac {H_{z}{\sigma }_{\xi \xi }}{n}}{\Biggr )}+{\frac {\partial }{\partial \eta }}{\Biggl (}{\frac {H_{z}{\sigma }_{\xi \eta }}{m}}{\Biggr )}+{\frac {\partial }{\partial s}}{\Biggl (}{\frac {{\sigma }_{\xi s}}{mn}}{\Biggr )}+\\&H_{z}{\sigma }_{\xi \eta }{\frac {\partial }{\partial \eta }}\left({\frac {1}{m}}\right)-H_{z}{\sigma }_{\eta \eta }{\frac {\partial }{\partial \xi }}\left({\frac {1}{n}}\right)-{\frac {1}{n}}{\sigma }_{ss}{\frac {\partial H_{z}}{\partial \xi }}{\Biggr ]}\end{aligned}}$$

(1)

$${\begin{aligned}F^{v}\equiv {\widehat {\eta }}\cdot \left(\nabla \cdot {\vec {\sigma }}\right)={\frac {mn}{H_{z}}}{\Biggl [}&{\frac {\partial }{\partial \xi }}{\Biggl (}{\frac {H_{z}{\sigma }_{\eta \xi }}{n}}{\Biggr )}+{\frac {\partial }{\partial \eta }}{\Biggl (}{\frac {H_{z}{\sigma }_{\eta \eta }}{m}}{\Biggr )}+{\frac {\partial }{\partial s}}{\Biggl (}{\frac {{\sigma }_{\eta s}}{mn}}{\Biggr )}+\\&H_{z}{\sigma }_{\eta \xi }{\frac {\partial }{\partial \xi }}\left({\frac {1}{n}}\right)-H_{z}{\sigma }_{\xi \xi }{\frac {\partial }{\partial \eta }}\left({\frac {1}{m}}\right)-{\frac {1}{m}}{\sigma }_{ss}{\frac {\partial H_{z}}{\partial \eta }}{\Biggr ]}\end{aligned}}$$

(2)

where

$${\begin{aligned}{\sigma }_{\xi \xi }&=\left(A_{M}+\nu \right)e_{\xi \xi }+\left(\nu -A_{M}\right)e_{\eta \eta }\,,\\{\sigma }_{\eta \eta }&=\left(\nu -A_{M}\right)e_{\xi \xi }+\left(A_{M}+\nu \right)e_{\eta \eta }\,,\\{\sigma }_{ss}&=2\,\nu \,e_{ss}\,,\\{\sigma }_{\xi \eta }&={\sigma }_{\eta \xi }=2\,A_{M}\,e_{\xi \eta }\,,\\{\sigma }_{\xi s}&=2\,K_{M}\,e_{\xi s}\,,\\{\sigma }_{\eta s}&=2\,K_{M}\,e_{\eta s}\,,\end{aligned}}$$

(3)

and the strain field is:

$${\begin{aligned}e_{\xi \xi }&=m{\frac {\partial u}{\partial \xi }}+mnv{\frac {\partial }{\partial \eta }}\left({\frac {1}{m}}\right)\,,\\e_{\eta \eta }&=n\;{\frac {\partial v}{\partial \eta }}+mnu{\frac {\partial }{\partial \xi }}\left({\frac {1}{n}}\right)\,,\\e_{ss}&={\frac {1}{H_{z}}}{\frac {\partial \left(\omega H_{z}\right)}{\partial s}}+{\frac {m}{H_{z}}}u{\frac {\partial H_{z}}{\partial \xi }}+{\frac {n}{H_{z}}}v{\frac {\partial H_{z}}{\partial \eta }}\,,\\2\,e_{\xi \eta }&={\frac {m}{n}}{\frac {\partial \left(nv\right)}{\partial \xi }}+{\frac {n}{m}}{\frac {\partial \left(mu\right)}{\partial \eta }}\,,\\2\,e_{\xi s}&={\frac {1}{mH_{z}}}{\frac {\partial \left(mu\right)}{\partial s}}+mH_{z}{\frac {\partial \omega }{\partial \xi }}\,,\\2\,e_{\eta s}&={\frac {1}{nH_{z}}}\;{\frac {\partial \left(nv\right)}{\partial s}}\;+n\;H_{z}{\frac {\partial \omega }{\partial \eta }}\,.\end{aligned}}$$

(4)

Here, $A_{M}(\xi ,\eta )$ and $K_{M}(\xi ,\eta ,s)$ are the spatially varying horizontal and vertical viscosity coefficients, respectively, and $\nu$ 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, $H_{z}={\frac {\partial z}{\partial s}}$. Also in these models, the vertical velocity is computed as ${\frac {\omega H_{z}}{mn}}$ and has units of ${\hbox{m}}^{3}/{\hbox{s}}$.

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

$${\begin{aligned}H_{z}F^{u}&={n^{2}}m{\partial \over \partial \xi }\left({\frac {H_{z}F^{u\xi }}{n}}\right)+{m^{2}}n{\partial \over \partial \eta }\left({\frac {H_{z}F^{u\eta }}{m}}\right)\\H_{z}F^{v}&={n^{2}}m{\partial \over \partial \xi }\left({\frac {H_{z}F^{v\xi }}{n}}\right)+{m^{2}}n{\partial \over \partial \eta }\left({\frac {H_{z}F^{v\eta }}{m}}\right)\end{aligned}}$$

(5)

where

$${\begin{aligned}F^{u\xi }&={\frac {1}{n}}\;A_{M}\left[{\frac {m}{n}}{\frac {\partial \left(nu\right)}{\partial \xi }}\;-{\frac {n}{m}}{\frac {\partial \left(mv\right)}{\partial \eta }}\right]\,,\\F^{u\eta }&={\frac {1}{m}}A_{M}\left[{\frac {n}{m}}{\frac {\partial \left(mu\right)}{\partial \eta }}+{\frac {m}{n}}{\frac {\partial \left(nv\right)}{\partial \xi }}\;\right]\,,\\F^{v\xi }&={\frac {1}{n}}\;A_{M}\left[{\frac {m}{n}}{\frac {\partial \left(nv\right)}{\partial \xi }}\;+{\frac {n}{m}}{\frac {\partial \left(mu\right)}{\partial \eta }}\right]\,,\\F^{v\eta }&={\frac {1}{m}}A_{M}\left[{\frac {n}{m}}{\frac {\partial \left(mv\right)}{\partial \eta }}-{\frac {m}{n}}{\frac {\partial \left(nu\right)}{\partial \xi }}\;\right]\,.\end{aligned}}$$

(6)

Notice the flux form of (5) and the symmetry between the $F^{u\xi }$ and $F^{v\eta }$ terms which are defined at density points on a C-grid. Similarly, the $F^{u\eta }$ and $F^{v\xi }$ 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 $H_{z}$ 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 $s$-surfaces. The transverse stress tensor rotated along geopotentials (constant depth) is, then, given by

$${\begin{aligned}H_{z}R^{u}&={n^{2}}m{\frac {\partial }{\partial \xi }}{\Biggl (}{\frac {H_{z}R^{u\xi }}{n}}{\Biggr )}+{m^{2}}n{\frac {\partial }{\partial \eta }}{\Biggl (}{\frac {H_{z}R^{u\eta }}{m}}{\Biggr )}+{\frac {\partial }{\partial s}}{\Biggl (}R^{us}{\Biggr )}\\H_{z}R^{v}&={n^{2}}m{\frac {\partial }{\partial \xi }}{\Biggl (}{\frac {H_{z}R^{v\xi }}{n}}{\Biggr )}+{m^{2}}n{\frac {\partial }{\partial \eta }}{\Biggl (}{\frac {H_{z}R^{v\eta }}{m}}{\Biggr )}+{\frac {\partial }{\partial s}}{\Biggl (}R^{vs}{\Biggr )}\end{aligned}}$$

(7)

where

$${\begin{aligned}R^{u\xi }=&{\frac {1}{n}}\;A_{M}\left[{\frac {1}{n}}\;\left(m{\frac {\partial \left(nu\right)}{\partial \xi }}-m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nu\right)}{\partial s}}\right)-{\frac {1}{m}}\left(n{\frac {\partial \left(mv\right)}{\partial \eta }}-n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mv\right)}{\partial s}}\right)\right]\,,\\R^{u\eta }=&{\frac {1}{m}}A_{M}\left[{\frac {1}{m}}\left(n{\frac {\partial \left(mu\right)}{\partial \eta }}-n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mu\right)}{\partial s}}\right)+{\frac {1}{n}}\;\left(m{\frac {\partial \left(nv\right)}{\partial \xi }}-m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nv\right)}{\partial s}}\right)\right]\,,\\R^{us}=&m{\frac {\partial z}{\partial \xi }}A_{M}\left[{\frac {1}{n}}\;\left(m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nu\right)}{\partial s}}-m{\frac {\partial \left(nu\right)}{\partial \xi }}\right)-{\frac {1}{m}}\left(n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mv\right)}{\partial s}}-n{\frac {\partial \left(mv\right)}{\partial \eta }}\right)\right]+\\&n\;{\frac {\partial z}{\partial \eta }}A_{M}\left[{\frac {1}{m}}\left(n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mu\right)}{\partial s}}-n{\frac {\partial \left(mu\right)}{\partial \eta }}\right)+{\frac {1}{n}}\;\left(m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nv\right)}{\partial s}}-m{\frac {\partial \left(nv\right)}{\partial \xi }}\right)\right]\,,\end{aligned}}$$

(8)

$${\begin{aligned}R^{v\xi }=&{\frac {1}{n}}\;A_{M}\left[{\frac {1}{n}}\;\left(m{\frac {\partial \left(nv\right)}{\partial \xi }}-m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nv\right)}{\partial s}}\right)+{\frac {1}{m}}\left(n{\frac {\partial \left(mu\right)}{\partial \eta }}-n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mu\right)}{\partial s}}\right)\right]\,,\\R^{v\eta }=&{\frac {1}{m}}A_{M}\left[{\frac {1}{m}}\left(n{\frac {\partial \left(mv\right)}{\partial \eta }}-n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mv\right)}{\partial s}}\right)-{\frac {1}{n}}\;\left(m{\frac {\partial \left(nu\right)}{\partial \xi }}-m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nu\right)}{\partial s}}\right)\right]\,,\\R^{vs}=&m{\frac {\partial z}{\partial \xi }}A_{M}\left[{\frac {1}{n}}\;\left(m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nv\right)}{\partial s}}-m{\frac {\partial \left(nv\right)}{\partial \xi }}\right)+{\frac {1}{m}}\left(n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mu\right)}{\partial s}}-n{\frac {\partial \left(mu\right)}{\partial \eta }}\right)\right]+\\&n\;{\frac {\partial z}{\partial \eta }}A_{M}\left[{\frac {1}{m}}\left(n{\frac {\partial z}{\partial \eta }}{\frac {1}{H_{z}}}{\frac {\partial \left(mv\right)}{\partial s}}-n{\frac {\partial \left(mv\right)}{\partial \eta }}\right)-{\frac {1}{n}}\;\left(m{\frac {\partial z}{\partial \xi }}{\frac {1}{H_{z}}}{\frac {\partial \left(nu\right)}{\partial s}}-m{\frac {\partial \left(nu\right)}{\partial \xi }}\right)\right]\,.\end{aligned}}$$

(9)

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 $s$-surfaces along geopotentials are discretized using a modified version of the triad approach of Griffies et al. (1998).

Horizontal Diffusion

Laplacian

The Laplacian of a scalar $C$ in curvilinear coordinates is:

$$\nabla ^{2}C=\nabla \cdot \nabla C=mn\left[{\partial \over \partial \xi }\!\!\left({m \over n}{\partial C \over \partial \xi }\right)+{\partial \over \partial \eta }\!\!\left({n \over m}{\partial C \over \partial \eta }\right)\right]$$

(10)

In ROMS, this term is multiplied by ${\nu _{2}H_{z} \over mn}$ and becomes

$$\left[{\partial \over \partial \xi }\!\!\left({\nu _{2}H_{z}m \over n}{\partial C \over \partial \xi }\right)+{\partial \over \partial \eta }\!\!\left({\nu _{2}H_{z}n \over m}{\partial C \over \partial \eta }\right)\right]$$

(11)

where $C$ is any tracer. This form guarantees that the term does not contribute to the volume-integrated equations.

Biharmonic

The biharmonic operator is $\nabla ^{4}=\nabla ^{2}\nabla ^{2}$; the corresponding term is computed using a temporary variable $Y$:

$$Y={mn \over H_{z}}\left[{\partial \over \partial \xi }\!\!\left({\nu _{4}H_{z}m \over n}{\partial C \over \partial \xi }\right)+{\partial \over \partial \eta }\!\!\left({\nu _{4}H_{z}n \over m}{\partial C \over \partial \eta }\right)\right]$$

(12)

and is

$$-\left[{\partial \over \partial \xi }\!\!\left({\nu _{4}H_{z}m \over n}{\partial Y \over \partial \xi }\right)+{\partial \over \partial \eta }\!\!\left({\nu _{4}H_{z}n \over m}{\partial Y \over \partial \eta }\right)\right]$$

(13)

where $C$ is once again any tracer and $\nu _{4}$ 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 $s$ 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 $z$ 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 $C$. The $\xi$ and $\eta$ components are locally horizontal, rather than along the $s$ surface. The diffusive fluxes are:

$${\begin{aligned}F^{\xi }&=\nu _{2}\left[m{\frac {\partial C}{\partial \xi }}-\underbrace {\left(m{\frac {\partial z}{\partial \xi }}\underbrace {+S_{x}} _{{\text{MIX}}{\_}{\text{ISO}}}\right){\frac {\partial C}{\partial z}}} _{{\text{MIX}}{\_}{\text{GEO}}}\right]\\F^{\eta }&=\nu _{2}\left[n{\partial C \over \partial \eta }-\underbrace {\left[n{\partial z \over \partial \eta }\underbrace {+S_{y}} _{{\text{MIX}}{\_}{\text{ISO}}}\right){\partial C \over \partial z}} _{{\text{MIX}}{\_}{\text{GEO}}}\right]\\F^{s}&=-\underbrace {{1 \over H_{z}}\left(m{\partial z \over \partial \xi }\underbrace {+S_{x}} _{{\text{MIX}}{\_}{\text{ISO}}}\right)F^{\xi }} _{{\text{MIX}}{\_}{\text{GEO}}}-\underbrace {{1 \over H_{z}}\left(n{\partial z \over \partial \eta }\underbrace {+S_{y}} _{{\text{MIX}}{\_}{\text{ISO}}}\right)F^{\eta }} _{{\text{MIX}}{\_}{\text{GEO}}}\end{aligned}}$$

(14)

where

$${\begin{aligned}S_{x}&={{\partial \rho \over \partial x} \over {\partial \rho \over \partial z}}={\left[m{\partial \rho \over \partial \xi }-{m \over H_{z}}{\partial z \over \partial \xi }{\partial \rho \over \partial s}\right] \over {1 \over H_{z}}{\partial \rho \over \partial s}}\\S_{y}&={{\partial \rho \over \partial y} \over {\partial \rho \over \partial z}}={\left[n{\partial \rho \over \partial \eta }-{n \over H_{z}}{\partial z \over \partial \eta }{\partial \rho \over \partial s}\right] \over {1 \over H_{z}}{\partial \rho \over \partial s}}\end{aligned}}$$

(15)

and there is some trickery such that the computations depend on the sign of $\partial z \over \partial \xi$ and of $\partial z \over \partial \eta$. No flux boundary conditions are easily imposed by setting

$${\begin{aligned}F^{\xi }=0\,\,&{\text{at}}\,\xi \,{\text{walls}}\\F^{\eta }=0\,\,&{\text{at}}\,\eta \,{\text{walls}}\\F^{s}=0\,\,&{\text{at}}\,s=-1,0\end{aligned}}$$

Finally, the flux divergence is calculated and is added to the right-hand-side term for the field being computed:

$${\partial \over \partial \xi }\left({H_{z}F^{\xi } \over n}\right)+{\partial \over \partial \eta }\left({H_{z}F^{\eta } \over m}\right)+{\partial \over \partial s}\left({H_{z}F^{s} \over mn}\right)$$

(16)

The biharmonic rotated mixing tensors are computed much as the non-rotated biharmonic mixing. We define a temporary variable $Y$ based on equation (16):

$$Y={mn \over H_{z}}\left[{\partial \over \partial \xi }\left({\nu _{4}H_{z}F^{\xi } \over n}\right)+{\partial \over \partial \eta }\left({\nu _{4}H_{z}F^{\eta } \over m}\right)+{\partial \over \partial s}\left({\nu _{4}H_{z}F^{s} \over mn}\right)\right]\,.$$

(17)

We then build up fluxes of $Y$ as in equations (14). We then apply equation (16) to these $Y$ fluxes to obtain the biharmonic mixing tensors. Again, the value of $\nu _{4}$ 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 roms.in. Several factors to consider when choosing these values are:

• spindown time The spindown time on wavenumber $k$ is ${1 \over k^{2}\nu _{2}}$ for the Laplacian operator and ${1 \over k^{4}\nu _{4}}$ for the biharmonic operator. The smallest wavenumber corresponds to the length $2\Delta x$ and is $k={\pi \over \Delta x}$, leading to

$$\Delta t<t_{damp}={\frac {\Delta x^{2}}{\pi ^{2}\nu _{2}}}\,\,\,\,{\text{or}}\,\,\,\,{\frac {\Delta x^{4}}{\pi ^{4}\nu _{4}}}$$

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:

$$\Delta x<L_{BL}=\left({\frac {\nu _{2}}{\beta }}\right)^{\frac {1}{3}}\,\,\,\,{\text{and}}\,\,\,\,\left({\frac {\nu _{4}}{\beta }}\right)^{\frac {1}{5}}$$

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 $\nu$.

Horizontal Diffusion

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.