mom_hor_visc module reference

Calculates horizontal viscosity and viscous stresses.

Data Types

hor_visc_cs

Control structure for horizontal viscosity.

Functions/Subroutines

`horizontal_viscosity()`	Calculates the acceleration due to the horizontal viscosity.
`horizontal_viscosity_block()`	Calculates the acceleration due to the horizontal viscosity using an explicit k-block size.
`hor_visc_leith_grad()`	Calculates the magnitude of the vorticity and divergence gradients, and the Laplacian of vorticity, that are used within horizontal_viscosity by the Leith, modified Leith, Leith+E, and QG Leith viscosity schemes.
`hor_visc_backscatter_h()`	Adds the MEKE-based (EY24_EBT_BS) backscatter contribution to the diagonal term of the stress tensor at h-points.
`hor_visc_backscatter_q()`	Adds the MEKE-based (EY24_EBT_BS) backscatter contribution to the cross term of the stress tensor at q-points.
`hor_visc_leithy_ah()`	Computes the Leith+E antisymmetric-viscosity ratio m_leithy, and updates the biharmonic viscosity Ah (and Ah_h) to include the Leith+E contribution, optionally smoothing m_leithy and Ah in the process.
`hor_visc_gme_setup()`	Calculates the barotropic tension and shearing strain fields and the GME efficiency and isopycnal height diffusivity fields that are used within horizontal_viscosity when GME (CSuse_GME) is active.
`hor_visc_init()`	Allocates space for and calculates static variables used by horizontal_viscosity.
`leithy_taper_function()`	leithy_taper_function returns 1 if zc is shallower than leithy_depth; 0 if deeper than leithy_depth+leithy_width; and an interpolating cubic spline in between.
`hor_visc_vel_stencil()`	hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size
`align_aniso_tensor_to_grid()`	Calculates factors in the anisotropic orientation tensor to be align with the grid.
`smooth_gme()`	Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any horizontal two-grid-point noise.
`sum_5x5()`	Apply a 5x5 weighted sum.
`smooth_x9_h()`	Apply a 9-point smoothing filter twice to a field staggered at a thickness point to reduce horizontal two-grid-point noise.
`smooth_x9_uv()`	Apply a 9-point smoothing filter twice to a pair of velocity components to reduce horizontal two-grid-point noise.
`hor_visc_end()`	Deallocates any variables allocated in hor_visc_init.

Detailed Description

Horizontal viscosity in MOM

This module contains the subroutine horizontal_viscosity that calculates the effects of horizontal viscosity, including parameterizations of the value of the viscosity itself. Subroutine horizontal_viscosity calculates the acceleration due to some combination of a biharmonic viscosity and a Laplacian viscosity. Either or both may use a coefficient that depends on the shear and strain of the flow. All metric terms are retained. The Laplacian is calculated as the divergence of a stress tensor, using the form suggested by [77]. The biharmonic is calculated by twice applying the divergence of the stress tensor that is used to calculate the Laplacian, but without the dependence on thickness in the first pass. This form permits a variable viscosity, and indicates no acceleration for either resting fluid or solid body rotation.

The form of the viscous accelerations is discussed extensively in [33], and the implementation here follows that discussion closely. We use the notation of [78] with the exception that the isotropic viscosity is $\kappa_h$.

In general, the horizontal stress tensor can be written as

\[\begin{split}{\bf \sigma} = \begin{pmatrix} \frac{1}{2} \left( \sigma_D + \sigma_T \right) & \frac{1}{2} \sigma_S \\\\ \frac{1}{2} \sigma_S & \frac{1}{2} \left( \sigma_D - \sigma_T \right) \end{pmatrix}\end{split}\]

where $\sigma_D$, $\sigma_T$ and $\sigma_S$ are stresses associated with invariant factors in the strain-rate tensor. For a Newtonian fluid, the stress tensor is usually linearly related to the strain-rate tensor. The horizontal strain-rate tensor is

\[\begin{split}\dot{\bf e} = \begin{pmatrix} \frac{1}{2} \left( \dot{e}_D + \dot{e}_T \right) & \frac{1}{2} \dot{e}_S \\\\ \frac{1}{2} \dot{e}_S & \frac{1}{2} \left( \dot{e}_D - \dot{e}_T \right) \end{pmatrix}\end{split}\]

where $\dot{e}_D = \partial_x u + \partial_y v$ is the horizontal divergence, $\dot{e}_T = \partial_x u - \partial_y v$ is the horizontal tension, and $\dot{e}_S = \partial_y u + \partial_x v$ is the horizontal shear strain.

The trace of the stress tensor, $tr(\bf \sigma) = \sigma_D$, is usually absorbed into the pressure and only the deviatoric stress tensor considered. From here on, we drop $\sigma_D$. The trace of the strain tensor, $tr(\bf e) = \dot{e}_D$ is non-zero for horizontally divergent flow but only enters the stress tensor through $\sigma_D$ and so we will drop $\sigma_D$ from calculations of the strain tensor in the code. Therefore the horizontal stress tensor can be considered to be

\[\begin{split}{\bf \sigma} = \begin{pmatrix} \frac{1}{2} \sigma_T & \frac{1}{2} \sigma_S \\\\ \frac{1}{2} \sigma_S & - \frac{1}{2} \sigma_T \end{pmatrix} .\end{split}\]

The stresses above are linearly related to the strain through a viscosity coefficient, $\kappa_h$:

\[\begin{split}\begin{eqnarray*} \sigma_T & = & 2 \kappa_h \dot{e}_T \\\\ \sigma_S & = & 2 \kappa_h \dot{e}_S . \end{eqnarray*}\end{split}\]

The viscosity $\kappa_h$ may either be a constant or variable. For example, $\kappa_h$ may vary with the shear, as proposed by [77].

The accelerations resulting form the divergence of the stress tensor are

\[\begin{split}\begin{eqnarray*} \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) & = & \partial_x \left( \frac{1}{2} \sigma_T \right) + \partial_y \left( \frac{1}{2} \sigma_S \right) \\\\ & = & \partial_x \left( \kappa_h \dot{e}_T \right) + \partial_y \left( \kappa_h \dot{e}_S \right) \\\\ \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) & = & \partial_x \left( \frac{1}{2} \sigma_S \right) + \partial_y \left( - \frac{1}{2} \sigma_T \right) \\\\ & = & \partial_x \left( \kappa_h \dot{e}_S \right) + \partial_y \left( - \kappa_h \dot{e}_T \right) . \end{eqnarray*}\end{split}\]

The form of the Laplacian viscosity in general coordinates is:

\[\begin{split}\begin{eqnarray*} \hat{\bf x} \cdot \left( \nabla \cdot \sigma \right) & = & \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_T \right) + \partial_y \left( \kappa_h h \dot{e}_S \right) \right] \\\\ \hat{\bf y} \cdot \left( \nabla \cdot \sigma \right) & = & \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_S \right) - \partial_y \left( \kappa_h h \dot{e}_T \right) \right] . \end{eqnarray*}\end{split}\]

Laplacian viscosity coefficient

The horizontal viscosity coefficient, $\kappa_h$, can have multiple components. The isotropic components are: * A uniform background component, $\kappa_{bg}$.

A constant but spatially variable 2D map, $\kappa_{2d}(x,y)$.
A ‘’MICOM’’ viscosity, $U_\nu \Delta(x,y)$, which uses a constant velocity scale, $U_\nu$ and a measure of the grid-spacing $\Delta(x,y)^2 = \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}$.
A function of latitude, $\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n$.
A dynamic Smagorinsky viscosity, $\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}$.
A dynamic Leith viscosity, $\kappa_{Lth}(x,y,t) = C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}$.

A maximum stable viscosity, $\kappa_{max}(x,y)$ is calculated based on the grid-spacing and time-step and used to clip calculated viscosities.

The static components of $\kappa_h$ are first combined as follows:

\[\kappa_{static} = \min \left[ \max\left( \kappa_{bg}, U_\nu \Delta(x,y), \kappa_{2d}(x,y), \kappa_\phi(x,y) \right) , \kappa_{max}(x,y) \right]\]

and stored in the module control structure as variables Kh_bg_xx and Kh_bg_xy for the tension (h-points) and shear (q-points) components respectively.

The full viscosity includes the dynamic components as follows:

\[\kappa_h(x,y,t) = r(\Delta,L_d) \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right)\]

where $r(\Delta,L_d)$ is a resolution function.

The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each other.

Viscous boundary conditions

Free slip boundary conditions have been coded, although no slip boundary conditions can be used with the Laplacian viscosity based on the 2D land-sea mask. For a western boundary, for example, the boundary conditions with the biharmonic operator would be written as:

\[\partial_x v = 0 ; \partial_x^3 v = 0 ; u = 0 ; \partial_x^2 u = 0 ,\]

while for a Laplacian operator, they are simply

\[\partial_x v = 0 ; u = 0 .\]

These boundary conditions are largely dictated by the use of an Arakawa C-grid and by the varying layer thickness.

Anisotropic viscosity

[52] proposed enhancing viscosity in a particular direction and the approach was generalized in [78]. We use the second form of their two coefficient anisotropic viscosity (section 4.3). We also replace their $A^\prime$ and $D$ such that $2A^\prime = 2 \kappa_h + D$ and $\kappa_a = D$ so that $\kappa_h$ can be considered the isotropic viscosity and $\kappa_a=D$ can be consider the anisotropic viscosity. The direction of anisotropy is defined by a unit vector $\hat{\bf n}=(n_1,n_2)$.

The contributions to the stress tensor are

\[\begin{split}\begin{pmatrix} \sigma_T \\\\ \sigma_S \end{pmatrix} = \left[ \begin{pmatrix} 2 \kappa_h + \kappa_a & 0 \\\\ 0 & 2 \kappa_h \end{pmatrix} + 2 \kappa_a n_1 n_2 \begin{pmatrix} - 2 n_1 n_2 & n_1^2 - n_2^2 \\\\ n_1^2 - n_2^2 & 2 n_1 n_2 \end{pmatrix} \right] \begin{pmatrix} \dot{e}_T \\\\ \dot{e}_S \end{pmatrix}\end{split}\]

Dissipation of kinetic energy requires $\kappa_h \geq 0$ and $2 \kappa_h + \kappa_a \geq 0$. Note that when anisotropy is aligned with the x-direction, $n_1 = \pm 1$, then $n_2 = 0$ and the cross terms vanish. The accelerations in this aligned limit with constant coefficients become

\[\begin{split}\begin{eqnarray*} \hat{\bf x} \cdot \nabla \cdot {\bf \sigma} & = & \partial_x \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right) + \partial_y \left( \kappa_h \dot{e}_S \right) \\\\ & = & \left( \kappa_h + \kappa_a \right) \partial_{xx} u + \kappa_h \partial_{yy} u - \frac{1}{2} \kappa_a \partial_x \left( \partial_x u + \partial_y v \right) \\\\ \hat{\bf y} \cdot \nabla \cdot {\bf \sigma} & = & \partial_x \left( \kappa_h \dot{e}_S \right) - \partial_y \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right) \\\\ & = & \kappa_h \partial_{xx} v + \left( \kappa_h + \kappa_a \right) \partial_{yy} v - \frac{1}{2} \kappa_a \partial_y \left( \partial_x u + \partial_y v \right) \end{eqnarray*}\end{split}\]

which has contributions akin to a negative divergence damping (a divergence enhancement?) but which is weaker than the enhanced tension terms by half.

Discretization

The horizontal tension, $\dot{e}_T$, is stored in variable sh_xx and discretized as

\[\dot{e}_T = \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} u \right) - \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} v \right) .\]

The horizontal divergent strain, $\dot{e}_D$, is stored in variable div_xx and discretized as

\[\dot{e}_D = \frac{1}{h A} \left( \delta_i \left( \overline{h}^i \Delta y \, u \right) + \delta_j \left( \overline{h}^j\Delta x \, v \right) \right) .\]

Note that for expediency this is the exact discretization used in the continuity equation.

The horizontal shear strain, $\dot{e}_S$, is stored in variable sh_xy and discretized as

\[\dot{e}_S = v_x + u_y\]

where

\[\begin{split}\begin{align*} v_x &= \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} v \right) \\\\ u_y &= \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} u \right) \end{align*}\end{split}\]

which are calculated separately so that no-slip or free-slip boundary conditions can be applied to $v_x$ and $u_y$ where appropriate.

The tendency for the x-component of the divergence of stress is stored in variable diffu and discretized as

\[\hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) = \frac{1}{A \overline{h}^i} \left( \frac{1}{\Delta y} \delta_i \left( h \Delta y^2 \kappa_h \dot{e}_T \right) + \frac{1}{\Delta x} \delta_j \left( \tilde{h}^{ij} \Delta x^2 \kappa_h \dot{e}_S \right) \right) .\]

The tendency for the y-component of the divergence of stress is stored in variable diffv and discretized as

\[\hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) = \frac{1}{A \overline{h}^j} \left( \frac{1}{\Delta y} \delta_i \left( \tilde{h}^{ij} \Delta y^2 A_M \dot{e}_S \right) - \frac{1}{\Delta x} \delta_j \left( h \Delta x^2 A_M \dot{e}_T \right) \right) .\]

References

Griffies, S.M., and Hallberg, R.W., 2000: Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Monthly Weather Review, 128(8), 2935-2946. https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2

Large, W.G., Danabasoglu, G., McWilliams, J.C., Gent, P.R. and Bryan, F.O., 2001: Equatorial circulation of a global ocean climate model with anisotropic horizontal viscosity. Journal of Physical Oceanography, 31(2), pp.518-536. https://doi.org/10.1175/1520-0485(2001)031%3C0518:ECOAGO%3E2.0.CO;2

Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear viscosities. Large eddy simulation of complex engineering and geophysical flows, 1, 69-106.

Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for ocean models. Ocean Modelling, 5(2), 129-156. https://doi.org/10.1016/S1463-5003(02)00016-1

Type Documentation

type mom_hor_visc/hor_visc_cs