Skip to content

Conversation

@vanroekel
Copy link
Collaborator

This is a draft design document for the upcoming submesoscale eddy parameterization. The compiled version is on nersc here. The document details two formulations, the original by Fox-Kemper et al (2011) and a new update by Bodner et al (2023). The latter formulation requires the boundary layer depth so cannot be implemented until KPP is implemented to Omega.

This PR also includes a change to conf.py to allow for equation numbering to restart in each individual design doc.

Checklist

  • Documentation:
  • Linting
  • Building
    • CMake build does not produce any new warnings from changes in this PR

@vanroekel
Copy link
Collaborator Author

@xylar a small question for you - I see diffs in the .precommit-config.yaml file that I didn't expect. Is this a rebase issue on my side?

The core of the FK08 scheme is the introduction of an eddy-induced velocity, $\vec{v}_e$, which advects buoyancy. The vertical component of the eddy-induced overturning streamfunction, $\Psi_e$, is parameterized as:

$$
\Psi_e = C_e \frac{H^2}{|f|} |\nabla_h b| \cdot \mu(z)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
\Psi_e = C_e \frac{H^2}{|f|} |\nabla_h b| \cdot \mu(z)
\Psi_e = C_e \frac{H^2}{|f|} |\nabla_h b_{ml}| \cdot \mu(z)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

changed



$$
\Psi(z) = C_e \frac{\Delta_s}{L_f} \frac{H^2 \nabla_h \overline{b_{ml}\times \hat{z}}}{\sqrt{f^2 + \tau^{-2}}} \mu(z)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
\Psi(z) = C_e \frac{\Delta_s}{L_f} \frac{H^2 \nabla_h \overline{b_{ml}\times \hat{z}}}{\sqrt{f^2 + \tau^{-2}}} \mu(z)
\Psi_e(z) = C_e \frac{\Delta_s}{L_f} \frac{H^2 \nabla_h \overline{b_{ml}\times \hat{z}}}{\sqrt{f^2 + \tau^{-2}}} \mu(z)

Leveraging turbulent thermal wind balance [McWilliams et al, 2015](https://journals.ametsoc.org/view/journals/phoc/45/8/jpo-d-14-0211.1.xml) and leveraging the ePBL mixing closure [Reichl and Hallberg, 2018](https://www.sciencedirect.com/science/article/pii/S1463500318301069), [](#submeso-definition) can be rewritten as

$$
\Psi(z) = C_r \frac{\Delta_s \left|f\right| h H^2 \nabla b_{ml} \times \hat{z} }{(m_* u_u^3 + n_* w_*^3)^{2/3}}\mu(z)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
\Psi(z) = C_r \frac{\Delta_s \left|f\right| h H^2 \nabla b_{ml} \times \hat{z} }{(m_* u_u^3 + n_* w_*^3)^{2/3}}\mu(z)
\Psi(z) = C_r \frac{\Delta_s \left|f\right| h H^2 \nabla b_{ml} \times \hat{z} }{(m_* u_*^3 + n_* w_*^3)^{2/3}}\mu(z)

\left|\rho(p) - \rho(p_0)\right| \geq \Delta \rho_t
$$ (mld-definition)

In [](#mld-definition), $p_0$ is a reference pressure, which is set to $10^5$ Pa ($\approx 10$m), and $\Delta \rho_t$ is a density threshold criterion, which is typically set to 0.03 kg m$^{-3}$. When the density threshold occurs between layers, the mixed layer depth is found via a linear fit between the bounding layers.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
In [](#mld-definition), $p_0$ is a reference pressure, which is set to $10^5$ Pa ($\approx 10$m), and $\Delta \rho_t$ is a density threshold criterion, which is typically set to 0.03 kg m$^{-3}$. When the density threshold occurs between layers, the mixed layer depth is found via a linear fit between the bounding layers.
In [](#mld-definition), $p_0$ is a reference pressure, which is set to $10^5$ Pa ($\approx 10$m), $\rho(p_0)$ is the local density at that pressure level, and $\Delta \rho_t$ is a density threshold criterion, which is typically set to 0.03 kg m$^{-3}$. When the density threshold occurs between layers, the mixed layer depth is found via a linear fit between the bounding layers.


In [](#mld-definition), $p_0$ is a reference pressure, which is set to $10^5$ Pa ($\approx 10$m), and $\Delta \rho_t$ is a density threshold criterion, which is typically set to 0.03 kg m$^{-3}$. When the density threshold occurs between layers, the mixed layer depth is found via a linear fit between the bounding layers.

In future versions of Omega, ice shelf cavities slightly complicate the determination of the mixed layer depth. The pressure must be referenced to the land ice pressure. This document will be updated once ice shelf cavities are implemented in Omega.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be a safe guess to just build land ice pressure in additively: $p_0 = p_{ice} + p_t$ where $p_t = 10^5$ Pa

```
## 5. Testing

### Test: Convergence to an exact solution

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds like this would be a CTest?


#### 4.2.4 Submesoscale eddy velocity

The horizontal buoyancy gradient averaged over the mixed layer is then used in [](#submeso-definition) to compute the submesoscale eddy stream function. The vertical divergence of the stream function is then used to create a `NormalVelocityMLE`. This term will be added to the `NormalVelocity` for use in the advection tendencies.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be helpful to specify whether buoyancy gradient and the streamfunction are originally computed on cells, edges or vertices?

Copy link

@katsmith133 katsmith133 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks great, @vanroekel! Thanks for this work. Most of my comments are minor and can be taken or left depending on how you feel they fit.

since $g$ is constant, we can write

$$
b_{ml} = g - \frac{1}{H}\int_{-H}^0 g \frac{\rho_0}{\rho} d \tilde{z} = g - \frac{1}{H}\int_{-H}^0 g \rho_0 \alpha

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
b_{ml} = g - \frac{1}{H}\int_{-H}^0 g \frac{\rho_0}{\rho} d \tilde{z} = g - \frac{1}{H}\int_{-H}^0 g \rho_0 \alpha
b_{ml} = g - \frac{1}{H}\int_{-H}^0 g \frac{\rho_0}{\rho} d \tilde{z} = g - \frac{1}{H}\int_{-H}^0 g \rho_0 \alpha d \tilde{z}

\nabla b = \nabla_h b + \frac{\partial b}{\partial \tilde{z}} \nabla_h \tilde{z}
$$

where $\nabla_h$ is a gradient along a fixed model layer in Omega. The vertical derivative of buoyancy will utilize the Brunt Vaisala Frequency variable from the vertical mixing calculation.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
where $\nabla_h$ is a gradient along a fixed model layer in Omega. The vertical derivative of buoyancy will utilize the Brunt Vaisala Frequency variable from the vertical mixing calculation.
where $\nabla_h$ is a gradient along a fixed model layer in Omega. The vertical derivative of buoyancy will utilize the Brunt Vaisala Frequency variable calculated in the Equation of State.


where $\nabla_h$ is a gradient along a fixed model layer in Omega. The vertical derivative of buoyancy will utilize the Brunt Vaisala Frequency variable from the vertical mixing calculation.

### 2.3 Requirement: The parameterization must be configurable in the traditiona FK11 form and the B23 update

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
### 2.3 Requirement: The parameterization must be configurable in the traditiona FK11 form and the B23 update
### 2.3 Requirement: The parameterization must be configurable in the traditional FK11 form and the B23 update

\tau = \sqrt{f^2 + \tau_{min}^{-2}}
$$

Here, $f$ is the Coriolis parameter and $\tau_{min}$ is a chosen minimum frontal slumping timescale to prevent excessive values of $\tau$ near the equator.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we have a typical/suggested value for $\tau_{min}$ or is this something we will need to tune? Maybe also note that this differs from the FK11 formula.


### 3.1 The original FK11 formulation

The basic definition of the submesoscale parameterization is given in [](#submeso-definition). In this equation, $L_f$, $C_e$, and $\tau$ are parameters. Some models, such as MOM6 merge $\frac{\Delta_s}{L_f}$ into a single parameter. This is not done in Omega to allow for variations in the grid size ($\Delta_s$) common to many target unstructured meshes. As in MPAS-Ocean, the characteristic frontal width ($L_f$) is given by

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The basic definition of the submesoscale parameterization is given in [](#submeso-definition). In this equation, $L_f$, $C_e$, and $\tau$ are parameters. Some models, such as MOM6 merge $\frac{\Delta_s}{L_f}$ into a single parameter. This is not done in Omega to allow for variations in the grid size ($\Delta_s$) common to many target unstructured meshes. As in MPAS-Ocean, the characteristic frontal width ($L_f$) is given by
The basic definition of the submesoscale parameterization is given in [](#submeso-definition). In this equation, $C_e$, $L_f$, $\Delta_s$, and $\tau$ are parameters. Some models, such as MOM6 merge $\frac{\Delta_s}{L_f}$ into a single parameter. This is not done in Omega to allow for variations in the grid size ($\Delta_s$) common to many target unstructured meshes. $C_e$ is defined above and is typically hard-coded to a value like 0.06. As in MPAS-Ocean, the characteristic frontal width ($L_f$) is taken as the maximum of three length scales


## 4. Design

The most complex part of this design is the spatial dependence implicit in the averaging of the buoyancy. The necessary `parallelReduce` calculation requires reduction over a spatially variable inner range. This will be accomplished with the hierarchical parallelism wrappers described in CITE HP design doc. In the initial implementation, only the FK11 form will be included. The B23 formulation requires the boundary layer depth, which can not be included until KPP [Large et al, 1994](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/94RG01872) is implemented for vertical mixing.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The most complex part of this design is the spatial dependence implicit in the averaging of the buoyancy. The necessary `parallelReduce` calculation requires reduction over a spatially variable inner range. This will be accomplished with the hierarchical parallelism wrappers described in CITE HP design doc. In the initial implementation, only the FK11 form will be included. The B23 formulation requires the boundary layer depth, which can not be included until KPP [Large et al, 1994](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/94RG01872) is implemented for vertical mixing.
The most complex part of this design is the spatial dependence implicit in the averaging of the buoyancy. The necessary `parallelReduce` calculation requires reduction over a spatially variable inner range. This will be accomplished with the hierarchical parallelism wrappers described in CITE HP design doc. In the initial implementation, only the FK11 form will be included. The B23 formulation requires the boundary layer depth, which can not be included until KPP ([Large et al, 1994](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/94RG01872) is implemented for vertical mixing).

- `Ce` submesoscale eddy efficiency (0.06-0.08)
- `RiCrit` critical Richardson number for the B23 form (~0.25)
- `DsMax` maximum grid length to include in the submesoscale calculation (~100 km)
- `DrhoCrit` critical threshold for the density MLD calculation (0.03 kg m$^-2$)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Two additional parameters:

  • Something like EnableSubmesos (we need to be able to turn them off completely as well, no?)
  • $\rho_0$, but maybe note that this can be grabbed from the vertical coord section its defined in in the config Yaml.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep, thanks. I've made the changes


#### 4.2.2 Density MLD calculation

The public `computeDenMLD` mthod will utilize the specific volume frmo the TEOS-10 calculation to compute a mixed layer depthy based on the threshold criterion. Given Omega has easy access to `SpecVol` and not density, we rearrange the standard density threshold calculation

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The public `computeDenMLD` mthod will utilize the specific volume frmo the TEOS-10 calculation to compute a mixed layer depthy based on the threshold criterion. Given Omega has easy access to `SpecVol` and not density, we rearrange the standard density threshold calculation
The public `computeDenMLD` method will utilize the specific volume from the TEOS-10 calculation to compute a mixed layer depth based on the threshold criterion. Given Omega has easy access to `SpecVol` and not density, we rearrange the standard density threshold calculation

b = g \frac{\rho - \rho_0}{\rho}
$$

where $\rho_0$ is the reference density used in the pseudo height coordinate definition. This can be easily rearranged to

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
where $\rho_0$ is the reference density used in the pseudo height coordinate definition. This can be easily rearranged to
where $\rho_0$ is the reference density used in the pseudo height coordinate definition. Again, since Omega calculates `SpecVol` and not density itself, this can be easily rearranged to

b = g \left(1 - \alpha \rho_0 \right)
$$

This form is preferred since Omega calculates `SpecVol` and not density itself. Since the buoyancy is only used for a horizontal gradient the reference pressure is not a consideration.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
This form is preferred since Omega calculates `SpecVol` and not density itself. Since the buoyancy is only used for a horizontal gradient the reference pressure is not a consideration.
Since the buoyancy is only used for a horizontal gradient the reference pressure is not a consideration.

Copy link

@katsmith133 katsmith133 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks great, @vanroekel! Thanks for this work. Most of my comments are minor and can be taken or left depending on how you feel they fit.

@vanroekel
Copy link
Collaborator Author

Thanks @katsmith133 and @cbegeman for the great comments. I believe I've addressed them all. I'd appreciate a rereview when you have a chance. @alicebarthel I'd also appreciate a review from you as well.

Copy link

@cbegeman cbegeman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your changes look good to me. Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants