Update docs for horizontal viscosity and backscatter

docs/parameterizations_lateral.rst

@@ -27,15 +27,20 @@ See :cite:`jansen2015` and :cite:`marshall2010`.
 Backscatter
 -----------
 
-A parameterization of the upscale unresolved cascade utilizes MOM_MEKE
-and negative Laplacian viscosity in MOM_hor_visc.
+* A parameterization of the upscale unresolved cascade that utilizes MOM_MEKE
+  and negative Laplacian viscosity is available in MOM_hor_visc. For details see 
+  :cite:`jansen2015b,yankovsky2024`
+* The 2D biharmonic Leith+E backscatter scheme is also available via MOM_hor_visc. For details
+  see :cite:`grooms2023`.
+* The Stochastic GM+E backscatter scheme is also available via MOM_stochastics. For details
+  see :cite:`grooms2025`
 
 Mixed layer restratification by sub-mesoscale eddies
 ----------------------------------------------------
 
 Mixed layer restratification from :cite:`fox-kemper2008` and
 :cite:`fox-kemper2008-2` is implemented in MOM_mixed_layer_restrat,
-which now also contains the mixed layer restratication comes from :cite: Bodner2023.
+which now also contains the mixed layer restratication comes from :cite:`Bodner2023`.
 
     :ref:`namespacemom__mixed__layer__restrat_1section_mle`
 

docs/zotero.bib

@@ -1,5 +1,29 @@
 % This file is generated by zotero. Manual edits will be lost!
 
+@article{grooms2025,
+author = {Grooms, I. and Agarwal, N. and Marques, G. and Pegion, P. J. and Yassin, H.},
+title = {The Stochastic GM + E Closure: A Framework for Coupling Stochastic Backscatter With the Gent and McWilliams Parameterization},
+journal = {Journal of Advances in Modeling Earth Systems},
+volume = {17},
+number = {5},
+pages = {e2024MS004560},
+doi = {https://doi.org/10.1029/2024MS004560},
+url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2024MS004560},
+year = {2025}
+}
+
+@article{yankovsky2024,
+author = {Yankovsky, Elizabeth and Bachman, Scott and Smith, K. Shafer and Zanna, Laure},
+title = {Vertical Structure and Energetic Constraints for a Backscatter Parameterization of Ocean Mesoscale Eddies},
+journal = {Journal of Advances in Modeling Earth Systems},
+volume = {16},
+number = {7},
+pages = {e2023MS004093},
+doi = {https://doi.org/10.1029/2023MS004093},
+url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2023MS004093},
+year = {2024}
+}
+
 @article{agarwal2023,
   title={Impact of Stochastic Ocean Density Corrections on Air-Sea Flux Variability},
   author={Agarwal, Niraj and Small, R Justin and Bryan, Frank O and Grooms, Ian and Pegion, Philip J},
@@ -317,6 +341,18 @@ @article{farneti2015
 	pages = {84--120}
 }
 
+@article{grooms2023,
+    title = {Backscatter in energetically-constrained Leith parameterizations},
+    journal = {Ocean Modelling},
+    volume = {186},
+    pages = {102265},
+    year = {2023},
+    issn = {1463-5003},
+    doi = {https://doi.org/10.1016/j.ocemod.2023.102265},
+    url = {https://www.sciencedirect.com/science/article/pii/S1463500323001051},
+    author = {Ian Grooms}
+}
+
 @article{griffies2000,
 	title = {Biharmonic {Friction} with a {Smagorinsky}-{Like} {Viscosity} for {Use} in {Large}-{Scale} {Eddy}-{Permitting} {Ocean} {Models}},
 	volume = {128},
@@ -362,6 +398,19 @@ @article{jansen2015
 	pages = {28--41}
 }
 
+@article{jansen2015b,
+    title={Energy budget-based backscatter in an eddy permitting primitive equation model},
+    author={Jansen, Malte F and Held, Isaac M and Adcroft, Alistair and Hallberg, Robert},
+    journal={Ocean Modelling},
+    volume={94},
+    pages={15--26},
+    year={2015},
+    issn = {1463-5003},
+    doi = {https://doi.org/10.1016/j.ocemod.2015.07.015},
+    url = {https://www.sciencedirect.com/science/article/pii/S1463500315001341},
+    month = oct
+}
+
 @article{stern2016,
 	title = {The effects of {Antarctic} iceberg calving-size distribution in a global climate model},
 	volume = {121},

src/core/MOM.F90

@@ -357,8 +357,8 @@ module MOM
                                 !! This is only valid when interp_p_surf is true.
   real, dimension(:,:), pointer :: &
     p_surf_prev  => NULL(), &   !< surface pressure [R L2 T-2 ~> Pa] at end  previous call to step_MOM
-    p_surf_begin => NULL(), &   !< surface pressure [R L2 T-2 ~> Pa] at start of step_MOM_dyn_...
-    p_surf_end   => NULL()      !< surface pressure [R L2 T-2 ~> Pa] at end   of step_MOM_dyn_...
+    p_surf_begin => NULL(), &   !< surface pressure [R L2 T-2 ~> Pa] at start of step_MOM_dyn_???
+    p_surf_end   => NULL()      !< surface pressure [R L2 T-2 ~> Pa] at end   of step_MOM_dyn_???
 
   ! Variables needed to reach between start and finish phases of initialization
   logical :: write_IC           !< If true, then the initial conditions will be written to file
@@ -525,7 +525,7 @@ module MOM
 contains
 
 !> This subroutine orchestrates the time stepping of MOM.  The adiabatic
-!! dynamics are stepped by calls to one of the step_MOM_dyn_...routines.
+!! dynamics are stepped by calls to one of the step_MOM_dyn_??? routines.
 !! The action of lateral processes on tracers occur in calls to
 !! advect_tracer and tracer_hordiff.  Vertical mixing and possibly remapping
 !! occur inside of diabatic.

src/core/MOM_dynamics_split_RK2.F90

@@ -1973,7 +1973,7 @@ end subroutine end_dyn_split_RK2
 !!  initialize_dyn_split_RK2 initializes the cpu clocks that are
 !!  used in this module.  For largely historical reasons, this module
 !!  does not have its own control structure, but shares the same
-!!  control structure with MOM.F90 and the other MOM_dynamics_...
+!!  control structure with MOM.F90 and the other MOM_dynamics_???
 !!  modules.
 
 end module MOM_dynamics_split_RK2

src/core/MOM_dynamics_split_RK2b.F90

@@ -1782,7 +1782,7 @@ end subroutine end_dyn_split_RK2b
 !!  initialize_dyn_split_RK2b initializes the cpu clocks that are
 !!  used in this module.  For largely historical reasons, this module
 !!  does not have its own control structure, but shares the same
-!!  control structure with MOM.F90 and the other MOM_dynamics_...
+!!  control structure with MOM.F90 and the other MOM_dynamics_???
 !!  modules.
 
 end module MOM_dynamics_split_RK2b

src/core/MOM_dynamics_unsplit.F90

@@ -33,7 +33,7 @@ module MOM_dynamics_unsplit
 !*  initialize_dyn_unsplit  initializes the cpu clocks that are        *                                      *
 !*  used in this module.  For largely historical reasons, this module  *
 !*  does not have its own control structure, but shares the same       *
-!*  control structure with MOM.F90 and the other MOM_dynamics_...      *
+!*  control structure with MOM.F90 and the other MOM_dynamics_???      *
 !*  modules.                                                           *
 !*                                                                     *
 !*  Macros written all in capital letters are defined in MOM_memory.h. *

src/core/MOM_dynamics_unsplit_RK2.F90

@@ -31,7 +31,7 @@ module MOM_dynamics_unsplit_RK2
 !*  initialize_dyn_unsplit_RK2 initializes the cpu clocks that are     *
 !*  used in this module.  For largely historical reasons, this module  *
 !*  does not have its own control structure, but shares the same       *
-!*  control structure with MOM.F90 and the other MOM_dynamics_...      *
+!*  control structure with MOM.F90 and the other MOM_dynamics_???      *
 !*  modules.                                                           *
 !*                                                                     *
 !*  Macros written all in capital letters are defined in MOM_memory.h. *

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -3622,15 +3622,15 @@ end subroutine hor_visc_end
 !! 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 \cite Smagorinsky1993. The biharmonic
+!! a stress tensor, using the form suggested by :cite:`Smagorinsky1993`. 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 \cite griffies2000,
+!! The form of the viscous accelerations is discussed extensively in :cite:`griffies2000`,
 !! and the implementation here follows that discussion closely.
-!! We use the notation of \cite Smith2003 with the exception that the
+!! We use the notation of :cite:`Smith2003` with the exception that the
 !! isotropic viscosity is \f$\kappa_h\f$.
 !!
 !! In general, the horizontal stress tensor can be written as
@@ -3680,7 +3680,7 @@ end subroutine hor_visc_end
 !! \f}
 !!
 !! The viscosity \f$\kappa_h\f$ may either be a constant or variable. For example,
-!! \f$\kappa_h\f$ may vary with the shear, as proposed by \cite Smagorinsky1993.
+!! \f$\kappa_h\f$ may vary with the shear, as proposed by :cite:`Smagorinsky1993`.
 !!
 !! The accelerations resulting form the divergence of the stress tensor are
 !! \f{eqnarray*}{
@@ -3723,42 +3723,85 @@ end subroutine hor_visc_end
 !! The horizontal viscosity coefficient, \f$\kappa_h\f$, can have multiple components.
 !! The isotropic components are:
 !!   - A uniform background component, \f$\kappa_{bg}\f$.
-!!   - A constant but spatially variable 2D map, \f$\kappa_{2d}(x,y)\f$.
+!!   - A steady but spatially variable 2D map, \f$\kappa_{2d}(x,y)\f$.
 !!   - A ''MICOM'' viscosity, \f$U_\nu \Delta(x,y)\f$, which uses a constant
 !! velocity scale, \f$U_\nu\f$ and a measure of the grid-spacing \f$\Delta(x,y)^2 =
 !! \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}\f$.
 !!   - A function of
 !! latitude, \f$\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n\f$.
-!!   - A dynamic Smagorinsky viscosity, \f$\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}\f$.
-!!   - A dynamic Leith viscosity, \f$\kappa_{Lth}(x,y,t) =
+!!   - A Smagorinsky viscosity, \f$\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}\f$.
+!!   - A Leith viscosity, \f$\kappa_{Lth}(x,y,t) =
 !!                                    C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}\f$.
-!!
-!! A maximum stable viscosity, \f$\kappa_{max}(x,y)\f$ is calculated based on the
-!! grid-spacing and time-step and used to clip calculated viscosities.
+!!   - A MEKE-based viscosity \f$\kappa_{MEKE}\f$; see the :ref:`meke_viscosity_section`.
 !!
 !! The static components of \f$\kappa_h\f$ are first combined as follows:
 !! \f[
-!! \kappa_{static} = \min \left[ \max\left(
+!! \kappa_{static} = \max\left(
 !! \kappa_{bg},
 !! U_\nu \Delta(x,y),
 !! \kappa_{2d}(x,y),
 !! \kappa_\phi(x,y)
 !! \right)
-!! , \kappa_{max}(x,y) \right]
 !! \f]
 !! and stored in the module control structure as variables <code>Kh_bg_xx</code> and
 !! <code>Kh_bg_xy</code> for the tension (h-points) and shear (q-points) components
 !! respectively.
 !!
-!! The full viscosity includes the dynamic components as follows:
+!! The full viscosity includes the flow-aware components as follows:
 !! \f[
 !! \kappa_h(x,y,t) = r(\Delta,L_d)
-!! \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right)
+!! \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth}, \kappa_{MEKE} \right)
+!! \f]
+!! where \f$r(\Delta,L_d)\f$ is a resolution function. The flow-aware components can optionally
+!! be added to the static components instead. A minimum can be also enforced after resolution
+!! scaling. The parameters of the resolution function for viscosity are controlled by
+!! <code>VISC_RES_SCALE_COEF</code> and <code>VISC_RES_FN_POWER</code>.
+!!
+!! A maximum stable viscosity, \f$\kappa_{max}(x,y)\f$ is calculated based on the
+!! grid-spacing and time-step and used to clip calculated viscosities. The maximum is applied
+!! after all other components have been combined.
+!!
+!! \subsection section_biharmonic_viscosity_coefficient Biharmonic viscosity coefficient
+!!
+!! The horizontal biharmonic viscosity coefficient, \f$A_h\f$, can have multiple components.
+!! The isotropic components are:
+!!   - A uniform background component, \f$A_{bg}\f$.
+!!   - A ''MICOM'' viscosity, \f$U_\nu \Delta(x,y)^3\f$, which uses a constant
+!! velocity scale, \f$U_\nu\f$ and a measure of the grid-spacing \f$\Delta(x,y)\f$ as above.
+!!   - A component based on the grid scale and a time scale \f$Delta(x,y)^4 / T\f$.
+!!   - A component that keeps the grid-Reynolds number fixed \f$\sqrt{KE} Delta(x,y)^3 / Re\f$.
+!!   - A Smagorinsky-like biharmonic viscosity,
+!!     \f$A_{Sm}(x,y,t) = C_{Sm} \Delta^4 \sqrt{\dot{e}_T^2 + \dot{e}_S^2} / 8\f$.
+!!   - A Leith viscosity, \f$A_{Lth}(x,y,t) =
+!!                                    C_{Lth} \Delta^r63 \sqrt{|\nabla^2 \zeta|^2}\f$.
+!!   - A MEKE-based viscosity \f$A_{MEKE}\f$; see the :ref:`meke_viscosity_section`.
+!!
+!! The static components of \f$A_h\f$ are first combined as follows:
+!! \f[
+!! A_{static} = \max\left(
+!! A_{bg},
+!! U_\nu \Delta(x,y)^3,
+!! \Delta(x,y)^4 / T
+!! \right)
+!! \f]
+!! and stored in the module control structure as variables <code>Kh_bg_xx</code> and
+!! <code>Kh_bg_xy</code> for the tension (h-points) and shear (q-points) components
+!! respectively.
+!!
+!! The full viscosity includes the flow-aware components as follows:
+!! \f[
+!! \kappa_h(x,y,t) =
+!! \max \left( A_{static}, A_{Sm}, A_{Lth} \right) + A_{MEKE}.
 !! \f]
-!! where \f$r(\Delta,L_d)\f$ is a resolution function.
 !!
-!! The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each
-!! other.
+!! The grid-Reynolds number scheme is exclusive of all other schemes including the
+!! static background.
+!!
+!! A maximum stable viscosity, \f$A_{max}(x,y)\f$ is calculated based on the
+!! grid-spacing and time-step and used to clip calculated viscosities. The maximum is applied
+!! after all other components have been combined.
+!!
+!! The 2D biharmonic Leith+E scheme from :cite:`grooms2023` is also available.
 !!
 !! \subsection section_viscous_boundary_conditions Viscous boundary conditions
 !!
@@ -3778,8 +3821,8 @@ end subroutine hor_visc_end
 !!
 !! \subsection section_anisotropic_viscosity Anisotropic viscosity
 !!
-!! \cite Large2001 proposed enhancing viscosity in a particular direction and the
-!! approach was generalized in \cite Smith2003. We use the second form of their
+!! :cite:`Large2001` proposed enhancing viscosity in a particular direction and the
+!! approach was generalized in :cite:`Smith2003`. We use the second form of their
 !! two coefficient anisotropic viscosity (section 4.3). We also replace their
 !! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and
 !! \f$\kappa_a = D\f$ so that \f$\kappa_h\f$ can be considered the isotropic
@@ -3899,6 +3942,9 @@ end subroutine hor_visc_end
 !! Monthly Weather Review, 128(8), 2935-2946.
 !! https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2
 !!
+!! Grooms, I., 2023. Backscatter in energetically-constrained Leith parameterizations.
+!! Ocean Modelling, 186, p.102265. https://doi.org/10.1016/j.ocemod.2023.102265
+!!
 !! 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.

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -2196,8 +2196,12 @@ end subroutine VarMix_end
 !! r(\Delta,L_d) = \frac{1}{1+(\alpha R)^p}
 !! \f]
 !!
-!! The resolution function can be applied independently to thickness diffusion \(module mom_thickness_diffuse\),
-!! tracer diffusion \(mom_tracer_hordiff\) lateral viscosity \(mom_hor_visc\).
+!! The resolution function can be applied to thickness diffusion \(module mom_thickness_diffuse\),
+!! tracer diffusion \(mom_tracer_hordiff\), and lateral viscosity \(mom_hor_visc\).
+!! <code>VISC_RES_SCALE_COEF</code> and <code>VISC_RES_FN_POWER</code> control the \f$\alpha\f$ and \f$p\f$
+!! parameters of the resolution function for viscosity and tracer diffusion, while
+!! <code>KH_RES_SCALE_COEF</code> and <code>KH_RES_FN_POWER</code> control the parameters of the
+!! resolution function for thickness diffusion.
 !!
 !! Robert Hallberg, 2013: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.
 !! Ocean Modelling, 71, pp 92-103.  http://dx.doi.org/10.1016/j.ocemod.2013.08.007
@@ -2208,10 +2212,10 @@ end subroutine VarMix_end
 !! | -                     | <code>RESOLN_SCALED_KH</code> |
 !! | -                     | <code>RESOLN_SCALED_KHTH</code> |
 !! | -                     | <code>RESOLN_SCALED_KHTR</code> |
-!! | \f$ \alpha \f$        | <code>KH_RES_SCALE_COEF</code> (for thickness and tracer diffusivity) |
-!! | \f$ p \f$             | <code>KH_RES_FN_POWER</code> (for thickness and tracer diffusivity) |
-!! | \f$ \alpha \f$        | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity) |
-!! | \f$ p \f$             | <code>VISC_RES_FN_POWER</code> (for lateral viscosity) |
+!! | \f$ \alpha \f$        | <code>KH_RES_SCALE_COEF</code> (for thickness diffusivity) |
+!! | \f$ p \f$             | <code>KH_RES_FN_POWER</code> (for thickness diffusivity) |
+!! | \f$ \alpha \f$        | <code>VISC_RES_SCALE_COEF</code> (for lateral viscosity and tracer diffusivity) |
+!! | \f$ p \f$             | <code>VISC_RES_FN_POWER</code> (for lateral viscosity and tracer diffusivity) |
 !! | -                     | <code>GILL_EQUATORIAL_LD</code> |
 !!
 !!

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -2560,6 +2560,11 @@ end subroutine thickness_diffuse_end
 !! (enabled with <code>USE_VARIABLE_MIXING=True</code> and the term \f$<SN>\f$ is the vertical average slope
 !! times the buoyancy frequency prescribed by \cite visbeck1996.
 !!
+!! Whether the resolution function is applied is controlled by <code>RESOLN_SCALED_KHTH</code>.
+!! When <code>INTERPOLATE_RES_FN = True</code>, the parameters of the resolution function for GM are
+!! controlled by <code>VISC_RES_SCALE_COEF</code> and <code>VISC_RES_FN_POWER</code>. Otherwise they are
+!! controlled by <code>KH_RES_SCALE_COEF</code> and <code>KH_RES_FN_POWER</code>.
+!!
 !! The result of the above expression is subsequently bounded by minimum and maximum values, including an upper
 !! diffusivity consistent with numerical stability (\f$ \kappa_{cfl} \f$ is calculated internally).
 !! \f[

src/tracer/MOM_tracer_hor_diff.F90

@@ -1849,5 +1849,11 @@ end subroutine tracer_hor_diff_end
 !!  diffusion if Khtr is defined and positive.  The tracer diffusion
 !!  can use a suitable number of iterations to guarantee stability
 !!  with an arbitrarily large time step.
+!!
+!! Whether the resolution function is applied is controlled by <code>RESOLN_SCALED_KHTR</code>.
+!! The parameters of the resolution function for tracer horizontal diffusion are controlled by
+!! <code>VISC_RES_SCALE_COEF</code> and <code>VISC_RES_FN_POWER</code>. They are not controlled
+!! by the parameters of the GM resolution function, <code>KH_RES_SCALE_COEF</code> and
+!! <code>KH_RES_FN_POWER</code>.
 
 end module MOM_tracer_hor_diff