Add vertical mixing coefficients module

components/omega/configs/Default.yml

@@ -59,6 +59,19 @@ Omega:
       DRhoDT: -0.2
       DRhoDS: 0.8
       RhoT0S0: 1000.0
+  VertMix:
+    Background:
+      Diffusivity: 1.0e-5
+      Viscosity: 1.0e-4
+    Convective:
+      Enable: true
+      Diffusivity: 1.0
+      TriggerBVF: 0.0
+    Shear:
+      Enable: true
+      NuZero: 0.005
+      Alpha: 5.0
+      Exponent: 2.0
   IOStreams:
     InitialVertCoord:
       UsePointerFile: false

components/omega/doc/devGuide/EOS.md

@@ -2,20 +2,27 @@
 
 # Equation of State (EOS)
 
-Omega includes an `Eos` class that provides functions that compute `SpecVol` and `SpecVolDisplaced`.
-Current EOS options are a linear EOS or an EOS computed using the TEOS-10 75 term expansion from
-[Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566). If
-`SpecVolDisplaced` is calculated with the linear EOS option, it will be equal to `SpecVol` as there
+Omega includes an `Eos` class that provides functions that compute `SpecVol`, `SpecVolDisplaced`,
+and `BruntVaisalaFreqSq`. Current EOS options are a linear EOS or an EOS computed using the TEOS-10
+75 term expansion from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566).
+If `SpecVolDisplaced` is calculated with the linear EOS option, it will be equal to `SpecVol` as there
 is no pressure/depth dependence for the linear EOS. `SpecVolDisplaced` computes specific volume
-adiabatically displaced to `K + KDisp`.
+adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreqSq`, as
+`SpecVol` is an input for the `BruntVaisalaFreqSq` calculation. If the linear EOS option is used, then the `BruntVaisalaFreqSq`
+is calculated using linear coefficients. If the TEOS-10 option is used, the `BruntVaisalaFreqSq` is calculated with non-linear
+coefficients according to the [TEOS-10 toolbox](https://www.teos-10.org/software.htm). Note: two assumption for ease of computation and efficiency have been made
+for the `BruntVaisalaFreqSq` TEOS-10 option that differ from how it is calculated in the TEOS-10 toolbox:
+(1) gravity is assumed to be constant and not a function of depth and latitude, and (2) the interface value of the specific volume is
+calculated as the average between two layer values, rather than being recalculated using the interface values of temperature,
+salinity, and pressure. Both of these assumptions incur less than a 1% error.
 
 ## Eos type
 
 An enumeration listing all implemented schemes is provided. It needs to be extended every time an
 EOS is added. It is used to identify which EOS method is to be used at run time.
 
 ```c++
-enum class EosType { Linear, Teos10Poly75t };
+enum class EosType { LinearEos, Teos10Eos };
 ```
 
 ## Initialization
@@ -49,9 +56,16 @@ and pressure arrays and displaced vertical index level, do
 Eos.computeSpecVolDisp(ConsrvTemp, AbsSalinity, Pressure, KDisp);
 ```
 
-where `KDisp` is the number of `k` levels you want to displace each specific volume level to.
+where `KDisp` is the number of `k` layers you want to displace each specific volume layer to.
 For example, to displace each level to one below, set `KDisp = 1`.
 
+To compute `BruntVaisalaFreqSq` for a particular set of temperature, salinity, pressure, and specific
+volume arrays, do
+
+```c++
+Eos.computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, Pressure, SpecVol);
+```
+
 ## Removal of Eos
 
 To clear the Eos instance do:

components/omega/doc/devGuide/VerticalMixingCoeff.md

@@ -0,0 +1,82 @@
+(omega-dev-vertmix) =
+
+# Vertical Mixing Coefficients
+
+Omega includes a `VertMix` class that provides functions that compute `VertDiff` and `VertVisc`, the
+vertical diffusivity and viscosity, where both are defined at the center of the cell and top of the layer.
+Currently the values of `VertDiff` and `VertVisc` are calculated using the linear combination of three options: (1) a
+constant background mixing value, (2) a convective instability mixing value, and (3) a Richardson
+number dependent shear mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These options are linearly additive. In the future, additional additive options will be implemented, such as the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872). For both the convective and shear mixing values `BruntVaisalaFreqSq` is needed, which
+is calculated by the `EOS` class.
+
+## Initialization and Usage
+
+The primary class `VertMix` is implemented using the Singleton pattern to ensure a single instance manages all vertical mixing calculations.
+
+```c++
+// Initialize VertMix
+VertMix* VMix = VertMix::getInstance();
+
+// Compute mixing coefficients
+VMix->computeVertMix(NormalVelocity, TangentialVelocity, BruntVaisalaFreqSq);
+```
+
+## Configuration
+
+The initialization process reads parameters from the yaml configuration file with the following structure and
+default values:
+
+```yaml
+VertMix:
+  Background:
+    Viscosity: 1e-4
+    Diffusivity: 1e-5
+  Convective:
+    Enable: true
+    Diffusivity: 1.0
+    TriggerBVF: 0.0
+  Shear:
+    Enable: true
+    NuZero: 0.005
+    Alpha: 5.0
+    Exponent: 2.0
+```
+
+## Class Structure
+
+### Core Data Members
+
+- `VertDiff`: 2D array storing vertical diffusivity coefficients (m²/s)
+- `VertVisc`: 2D array storing vertical viscosity coefficients (m²/s)
+
+### Mixing Parameters
+
+1. Background Mixing:
+   - `BackDiff`: Background vertical diffusivity (m²/s; Default: 1e-5)
+   - `BackVisc`: Background vertical viscosity (m²/s: Default: 1e-4)
+
+2. Convective Mixing:
+   - `EnableConvMix`: Flag to enable/disable convective mixing (Default: True)
+   - `ConvDiff`: Convective mixing coefficient (m²/s; Default: 1.0)
+   - `ConvTriggerBVF`: Trigger threshold for convective mixing (Default: 0.0)
+
+3. Shear Mixing:
+   - `EnableShearMix`: Flag to enable/disable shear mixing (Default: True)
+   - `ShearNuZero`: Base coefficient for Pacanowski-Philander scheme (Default: 0.005)
+   - `ShearAlpha`: Alpha parameter for P-P scheme (Default: 5.0)
+   - `ShearExponent`: Exponent parameter for P-P scheme (Default: 2.0)
+
+## Core Functionality (Vertical Mixing Coefficient Calculation)
+
+The main computation is handled by:
+
+```cpp
+void computeVertMix(const Array2DReal &NormalVelocity,
+                   const Array2DReal &TangentialVelocity,
+                   const Array2DReal &BruntVaisalaFreqSq);
+```
+
+This method combines the effects of:
+- Background mixing (constant coefficients)
+- Convective mixing (triggered by static instability)
+- Shear instability driven mixing (Pacanowski-Philander scheme; to be changed to [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) shear mixing scheme in a later development)

components/omega/doc/index.md

@@ -50,6 +50,7 @@ userGuide/Tracers
 userGuide/TridiagonalSolvers
 userGuide/VertCoord
 userGuide/Timing
+userGuide/VerticalMixingCoeff
 ```
 
 ```{toctree}
@@ -92,6 +93,7 @@ devGuide/Tracers
 devGuide/TridiagonalSolvers
 devGuide/VertCoord
 devGuide/Timing
+devGuide/VerticalMixingCoeff
 ```
 
 ```{toctree}

components/omega/doc/userGuide/EOS.md

@@ -4,7 +4,7 @@
 
 The equation of state (EOS) for the ocean describes the relationship between specific volume of seawater (in $\textrm{m}^3/\textrm{kg}$; the reciprocal of density) and temperature (in $^{\circ}\textrm{C}$), salinity (in $\textrm{g/kg}$), and pressure (in $\textrm{dbar}$). Through the hydrostatic balance (which relates density/specific volume gradients to pressure gradients), the equation of state provides a connection between active tracers (temperature and salinity) and the fluid dynamics.
 
-Two choices of EOS are provided by Omega: a linear EOS and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. However, this option is only recommended for simpler idealized test cases as its accuracy is not sufficient for real ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf) , but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega.
+Two choices of EOS are provided by Omega: a linear EOS and a TEOS-10 EOS. The linear EOS simplifies the relationship by excluding the influence of pressure and using constant expansion/contraction coefficients, making the specific volume a simple linear function of temperature and salinity. However, this option should only be used for simpler idealized test cases as it is not appropriate for realistic ocean simulations. The TEOS-10 EOS is a 75-term polynomial expression from [Roquet et al. 2015](https://www.sciencedirect.com/science/article/pii/S1463500315000566) that approximates the [Thermodynamic Equation of Seawater 2010](https://www.teos-10.org/pubs/TEOS-10_Manual.pdf), but in a less complex and more computationally efficient manner, and is the preferred EOS for real ocean simulations in Omega.
 
 The user-configurable options are: `EosType` (choose either `Linear` or `Teos-10`), as well as the parameters needed for the linear EOS.
 
@@ -19,4 +19,12 @@ Eos:
 
 where `DRhoDT` is the thermal expansion coefficient ($\textrm{kg}/(\textrm{m}^3 \cdot ^{\circ}\textrm{C})$), `DRhoDS` is the saline contraction coefficient ($\textrm{kg}/\textrm{m}^3$), and `RhoT0S0` is the reference density at (T,S)=(0,0) (in $\textrm{kg}/\textrm{m}^3$).
 
-In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` is also calculated by the EOS. This calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels, capturing the effects of pressure/depth changes. This is primarily used to calculate quantities for determining the water column stability (i.e. the stratification) and the vertical mixing coefficients (viscosity and diffusivity). Note: when using the linear EOS, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth.
+In addition to `SpecVol`, the displaced specific volume `SpecVolDisplaced` and `BruntVaisalaFreqSq` are also calculated by the EOS.
+
+## Displaced Specific Volume
+
+The `Eos` class calculates the density of a parcel of fluid that is adiabatically displaced by a relative `k` levels (`k` counted positive downward), capturing the effects of pressure/depth changes. This is primarily used to calculate quantities for determining the water column stability (i.e. the stratification) and the vertical mixing coefficients (viscosity and diffusivity). Note: when using the `Linear` EOS option, `SpecVolDisplaced` will be the same as `SpecVol` since the specific volume calculation is independent of pressure/depth.
+
+## Squared Brunt Vaisala Frequency
+
+The `Eos` class also calculates the squared Brunt Vaisala Frequency, which is used by the `VertMix` class to determine water column stability for both the convective adjustment and shear-instability driven mixing. When using the `Linear` EOS option, the `BruntVaisalaFreqSq` is calculated using linear coefficients. When using the `Teos-10` EOS option, the `BruntVaisalaFreqSq` is calculated with non-linear coefficients according to the TOES-10. For both options, `SpecVol` must be calculated prior to calculating `BruntVaisalaFreqSq`, as it is an input to the `BruntVaisalaFreqSq` calculation.

components/omega/doc/userGuide/VerticalMixingCoeff.md

@@ -0,0 +1,70 @@
+(omega-user-vertmix)=
+
+# Vertical Mixing Coefficients
+
+The vertical mixing module in Omega handles the parameterization of unresolved vertical mixing
+processes in the ocean. It calculates vertical diffusivity and viscosity coefficients that
+determine how properties (like momentum, heat, salt, and biogeochemical tracers) mix vertically
+in the ocean model. Currently, Omega offers three different mixing processes within the water column: (1) a constant
+background mixing value, (2) a convective instability mixing value, and (3) a Richardson number
+dependent shear instability driven mixing value from the [Pacanowski and Philander (1981)](https://journals.ametsoc.org/view/journals/phoc/11/11/1520-0485_1981_011_1443_povmin_2_0_co_2.xml) parameterization. These are linearly additive and are describe a bit
+more in detail below. Other mixing processes and parameterizations, such as the the K Profile Parameterization [(KPP; Large et al., 1994)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872) will be added in the future.
+
+The user-configurable options are the following parameters in the yaml configuration file:
+
+```yaml
+VertMix:
+  Background:
+    Viscosity: 1.0e-4    # Background vertical viscosity (m²/s)
+    Diffusivity: 1.0e-5  # Background vertical diffusivity (m²/s)
+  Convective:
+    Enable: true         # Enables the convective-induced mixing option
+    Diffusivity: 1.0     # Convective mixing coefficient (m²/s)
+    TriggerBVF: 0.0      # Squared Brunt-Vaisala frequency threshold
+  Shear:
+    Enable: true         # Enables the shear-instability driven mixing option
+    NuZero: 1.0e-2       # Base viscosity coefficient
+    Alpha: 5.0           # Stability parameter
+    Exponent: 2          # Richardson number exponent
+```
+
+## Vertical Mixing Processes/Types
+
+### 1. Background Mixing
+
+A constant background mixing value that represents small-scale mixing processes not explicitly resolved by the model. Typically, this is chosen to represent low values of vertical mixing
+happening in the ocean's stratified interior.
+
+### 2. Convective Mixing
+
+Enhanced convective adjustment mixing that occurs in statically unstable regions of the water column to parameterize convection and homogenize properties. In Omega this is mixing is defaulted to occur when the squared Brunt Vaisala Frequency is less than 0.0 (unstable), and is off when equal to (neutral) or greater than (stable) 0.0.
+
+$$
+\kappa =
+\begin{cases}
++\kappa_b + \kappa_{conv} \quad \text{ if } N^2 < N^2_{crit}\\
++\kappa_b \quad \text{ if } N^2 \geq N^2_{crit}
+\end{cases}
+$$
+
+This is different than some current implementations (i.e. in MPAS-Ocean and the CVMix package), where convective adjustment occurs both with unstable and neutral conditions ($N^2 \leq N^2_{crit}$). $\kappa_{conv}$ and $N^2_{crit}$ are constant parameters set in the `VertMix` section of the yaml file (`Diffusivity` and `TriggerBVF` under the `Convective` header).
+
+### 3. Shear Mixing
+
+Mixing induced by vertical velocity shear, implemented using the Pacanowski-Philander scheme, through the gradient Richardson number (ratio of buoyancy to shear).
+
+$$
+\nu = \frac{\nu_o}{(1+\alpha Ri)^n} + \nu_b\,,
+$$
+
+$$
+\kappa = \frac{\nu}{(1+\alpha Ri)} + \kappa_b\,.
+$$
+
+where $Ri$ is defined as:
+
+$$
+Ri = \frac{N^2}{\left|\frac{\partial \mathbf{U}}{\partial z}\right|^2}\,,
+$$
+
+where $\nu_o$, $\alpha$, $n$, $\nu_b$, $\kappa_b$ are constant parameters set in the `VertMix` section of the yaml file (`NuZero`, `Alpha`, `Exponent` under the `Shear` header and `Viscosity`, `Diffusivity` under the `Background` header). $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $N^2$, $\partial \mathbf{U}}{\partial z}\right|^2$ and $Ri$ of `K` are all defined at the cell center, top interface of layer `K`.  $N^2$, $\nu_{shear}$ and $\kappa_{shear}$ are set to zero for the surface layer. In a later development, the shear mixing option will be changed to the interior shear mixing formulation in [Large et al., 1994](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94rg01872).