Change shear mixing to LMD94 formulation and add ri smoothing#352
Change shear mixing to LMD94 formulation and add ri smoothing#352katsmith133 wants to merge 2 commits intoE3SM-Project:developfrom
Conversation
| BaseShearValue: 0.005 | ||
| RiCrit: 0.7 | ||
| Exponent: 3.0 | ||
| RiSmoothLoops: 3 |
There was a problem hiding this comment.
three is a lot for smoothing passes. I'd only do 2 here.
There was a problem hiding this comment.
Ok, I think 3 was the default in MPAS, so I was just going off of that, but happy to change it.
| 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. | ||
| number dependent shear mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear 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 `BruntVaisalaFreq` is needed, which |
There was a problem hiding this comment.
this is a little pedantic, but let's change shear mixing -> shear instability driven mixing throughout.
| BaseShearValue: 0.005 | ||
| RiCrit: 0.7 | ||
| Exponent: 3.0 | ||
| RiSmoothLoops: 3 |
| - `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) | ||
| - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) |
There was a problem hiding this comment.
| - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) | |
| - `BaseShearValue`: Base values of maximum viscosity and diffusivity for the LMD94 interior shear instability driven mixing formulation (Default: 0.005) |
| - `ShearAlpha`: Alpha parameter for P-P scheme (Default: 5.0) | ||
| - `ShearExponent`: Exponent parameter for P-P scheme (Default: 2.0) | ||
| - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) | ||
| - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) |
There was a problem hiding this comment.
| - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) | |
| - `RiCrit`: Critical Richardson number for the LMD94 formulation (Default: 0.7) |
| - `ShearExponent`: Exponent parameter for P-P scheme (Default: 2.0) | ||
| - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) | ||
| - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) | ||
| - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) |
There was a problem hiding this comment.
| - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) | |
| - `ShearExponent`: Exponent parameter number for the LMD94 formulation (Default: 3.0) |
| - `BaseShearValue`: Base values of shear for the LMD94 interior shear mixing formulation (Default: 0.005) | ||
| - `RiCrit`: Critical Richerson number for the LMB94 formulation (Default: 0.7) | ||
| - `ShearExponent`: Exponent parameter number for the LMB94 formulation (Default: 3.0) | ||
| - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) |
There was a problem hiding this comment.
| - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 3) | |
| - `RiSmoothLoops`: Number of 1-2-1 filter passes to apply to the gradient Richardson number smoothing (Default: 2) |
| 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. | ||
| dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear 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. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. |
There was a problem hiding this comment.
| 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. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear mixing calculation. | |
| 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. In addition to diffusivity and viscosity coefficients, the vertical mixing module calculates the gradient Richardson number and smooths that gradient Richardson number using a 1-2-1 filter before using it in the shear instability driven mixing calculation. |
| BaseShearValue: 0.005 # Base viscosity/diffusivity value | ||
| RiCrit: 0.7 # Critical Richardson number | ||
| Exponent: 3.0 # Richardson number exponent | ||
| RiSmoothLoops: 3 # Number of Richardson number smoothing loops |
There was a problem hiding this comment.
| RiSmoothLoops: 3 # Number of Richardson number smoothing loops | |
| RiSmoothLoops: 2 # Number of Richardson number smoothing loops |
|
|
||
| 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. | ||
| A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. |
There was a problem hiding this comment.
| A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior. | |
| A constant background mixing value that represents small-scale mixing processes not explicitly resolved or modeled. Typically, this is chosen to represent low values of vertical mixing happening in the ocean's stratified interior and is assumed roughly equivalent to the globally averaged interior mixing from all sources (e.g., from internal wave breaking). |
| $$ | ||
|
|
||
| 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). | ||
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $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. |
There was a problem hiding this comment.
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $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. | |
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear instability driven mixing. $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. |
There was a problem hiding this comment.
instead of 'surface layer' should we say 'at the surface'?
| 0.005; ///< Numerator of Pacanowski and Philander (1981) Eq (1). | ||
| Real ShearAlpha = 5.0; ///< Alpha value used in Pacanowski and Philander | ||
| ///< (1981) Eqs (1) and (2). | ||
| Real BaseShearValue = 0.005; ///< Base shear vertical viscosity and |
There was a problem hiding this comment.
| Real BaseShearValue = 0.005; ///< Base shear vertical viscosity and | |
| Real BaseShearValue = 0.005; ///< Base shear instability driven vertical viscosity and |
| ///< (1981) Eqs (1) and (2). | ||
| Real BaseShearValue = 0.005; ///< Base shear vertical viscosity and | ||
| ///< diffusivity (m^2 s^-1) of LMG94 | ||
| Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMG94 |
There was a problem hiding this comment.
| Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMG94 | |
| Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMD94 |
| Real ShearRiCrit = 0.7; ///< Critical Richardson number of LMG94 | ||
| Real ShearExponent = | ||
| 2.0; /// Exponent value used in Pacanowski and Philander (1981) Eqs (1). | ||
| 3.0; /// Exponent value used interior shear mixing calculation of LMG94 |
There was a problem hiding this comment.
| 3.0; /// Exponent value used interior shear mixing calculation of LMG94 | |
| 3.0; /// Exponent value used interior shear instability driven mixing calculation of LMG94 |
| Real ShearExponent = | ||
| 2.0; /// Exponent value used in Pacanowski and Philander (1981) Eqs (1). | ||
| 3.0; /// Exponent value used interior shear mixing calculation of LMG94 | ||
| I4 RiSmoothLoops = 3; ///< Number of smoothing loops for Richardson number |
There was a problem hiding this comment.
| I4 RiSmoothLoops = 3; ///< Number of smoothing loops for Richardson number | |
| I4 RiSmoothLoops = 2; ///< Number of smoothing loops for Richardson number |
| ShearExponent) * | ||
| BaseShearValue; | ||
| } else { | ||
| VertDiff(ICell, K) += 0.0_Real; |
There was a problem hiding this comment.
why do this? Do we really need this else chunk? We could just have a comment that says 'do nothing' or something like that
vanroekel
left a comment
vanroekel
left a comment
There was a problem hiding this comment.
Looks pretty good @katsmith133 just a few minor comments in the docs and one larger one in the main code.
|
Thanks for the review, @vanroekel! I will have a look at it today. Seems like I misspelled LMD94 A LOT, haha. |
|
@cbegeman I mostly included you on this review because the polaris omega_pr tests just hang and I want to figure out if this is a polaris problem or a this PR problem, especially after Youngsung mentioned that he saw these tests hanging on PM as well. Let me know your thoughts. |
|
@katsmith133 Thanks for the ping. I'd like to do polaris tests including the single column shear mixing one. Can you provide me with the appropriate mapping from Omega config options to MPAS-O config options (or indicate where there is no direct mapping? See https://github.com/E3SM-Project/polaris/pull/435/files#diff-d953af580e6050fd590b94c370ebc3401d6ed1b24e282d8d2d568536c35f8949R132 for an example. Is there anything else I should know about when comparing to MPAS-O and with the current Omega develop? |
|
|
||
| ### 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. |
There was a problem hiding this comment.
| Enhanced convective adjustment mixing that occurs in statically unstable regions of the water column to parameterize convection and homogenize properties. In Omega this 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. |
| 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. | ||
| dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit |
There was a problem hiding this comment.
For readability, I'd recommend hyphenating parts of this sentence (including previous line "Richardson-number-
| dependent shear instability driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit | |
| -dependent shear-instability-driven mixing value from the [Large et al (1994)](https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/94RG01872) or LMD94 interior shear parameterization. These are linearly additive and are describe a bit |
| @@ -51,20 +51,21 @@ | |||
|
|
|||
| ### 3. Shear Mixing | |||
There was a problem hiding this comment.
| ### 3. Shear Mixing | |
| ### 3. Shear-driven Mixing |
|
|
||
| $$ | ||
| \kappa = \frac{\nu}{(1+\alpha Ri)} + \kappa_b\,. | ||
| \nu_{shear}\,, \kappa_{shear} = = |
There was a problem hiding this comment.
| \nu_{shear}\,, \kappa_{shear} = = | |
| \nu_{shear} = = \kappa_{shear} = = |
| $$ | ||
|
|
||
| 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). | ||
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $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. |
There was a problem hiding this comment.
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 filter before being used to calculate the shear. $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. | |
| where $N^2$ is calculated by the EOS based on the ocean state and $\mathbf{U}$ is the magnitude of the horizontal velocity. $Ri$ is calculated by the vertical mixing module and then smoothed with a 1-2-1 (vertical) filter before being used to calculate the shear-instability-driven mixing. $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. |
|
@katsmith133 do you need more testing from me or just a code read-through? |
| 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` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreqSq`, as | ||
| adiabatically displaced to `K + KDisp` (where `K` counted positive downward, ie `K+1` is one layer below `K`). Note: `SpecVol` must be calculated before `BruntVaisalaFreq`, as |
There was a problem hiding this comment.
I am not sure why this change is there. Is that intentional or a version mix-up?
| \nu_{shear}\,, \kappa_{shear} = = | ||
| \begin{cases} | ||
| \nu_o \quad \text{ if } Ri_g < 0\\ | ||
| \nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ |
There was a problem hiding this comment.
| \nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ | |
| \nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 \leq Ri_g < Ri_{crit}\\ |
There was a problem hiding this comment.
@alicebarthel and @vanroekel This is actually a good point to discuss that I hadn't thought about. The LMD94 paper only has < and >, no <= or >=... I think I just naturally put the >= in for the code, but can't remember what my reasoning was. I am not too sure it matters too much which side the equals goes on when Ri_g ~ Ri_{crit}, since Ri_{crit} is a somewhat fudgey number anyway, but when Ri_g ~ 0 you are switching from a shear dominated regime to strong convection, so it might matter more which side 0 falls to? But maybe not, since this part of the code does not decide which is more dominant, it just decides how strong the shear is. Thoughts? I will take a look to see what cvmix does.
There was a problem hiding this comment.
Looking more at the form of the equation, I don't think it matters too much which side the equals is on since its asymptoting to those end values anyways.
There was a problem hiding this comment.
Yes that's right. It should be good to go either way. My recollection is the equals isn't necessary and the original could be fine too.
There was a problem hiding this comment.
My main point was to align the code and the documentation, and to make sure we don't create edge cases that we did not intend or do not test for. I have no preference as to where the equal goes.
| \begin{cases} | ||
| \nu_o \quad \text{ if } Ri_g < 0\\ | ||
| \nu_o \left[1 - \left( \frac{Ri_g}{Ri_{crit}} \right)^2 \right]^p \text{ if } 0 < Ri_g < Ri_{crit}\\ | ||
| 0.0 \quad \text{ if } Ri_{crit} < Ri_g |
There was a problem hiding this comment.
| 0.0 \quad \text{ if } Ri_{crit} < Ri_g | |
| 0.0 \quad \text{ if } Ri_{crit} \leq Ri_g |
| I4 K1 = K - 1; | ||
| I4 K2 = K; | ||
|
|
||
| if (K <= MinLayerCell(ICell)) { |
There was a problem hiding this comment.
A couple of clarifying questions here on the vertical indexing:
- can you confirm that maxLayerCell is below topography? I.e. the water column goes from minLayerCell (valid) to maxLayerCell-1?
- If K < MinLayerCell or K > MaxLayerCell, we are still at risk of using invalid/undefined velocities in the shear calculation. Looking at the loop, I don't think we are supposed to ever get there, but wouldn't it be safer to exit with an error (hopefully never exercised) in that case rather than risk using garbage? Is there a reason not to use something similar to the code above? e.g.
if (K < MinLayerCell(ICell) || K > MaxLayerCell(ICell))
continue;
alicebarthel
left a comment
alicebarthel
left a comment
There was a problem hiding this comment.
Thanks for your work on this @katsmith133!
I did a visual inspection and it looks good. I made minor suggestions, but let me know if it's helpful for me to run tests.
Question over testing in general: do we have any tests where we have more vertLayers than valid layers? I.e. do we ever test the interplay between nvertLayers and maxLayerCell, or minLayerCell? I suspect the former gets used when topography is used, but not the latter at this stage?
|
@cbegeman here is the mapping for Omega to MPAS-O for the VertMix variables: | Omega -> MPAS-O | These do not have direct mappings to anything in Omega, yet, but make sure they are set as following: These do not have direct mappings to anything in Omega, yet, and the default values should be fine to keep: Let me know if you need more than these variables! |
4 participants
This PR changes the shear vertical mixing formulation from Pacanowski and Philander (1981) to Large et al (1994), adds in vertical smoothing to the gradient Richardson number, and updates the tests and docs to reflect these changes.
Passes all CTests on PM-CPU and PM-GPU. Hangs on tests in polaris omega_pr test suite. Not sure how to address as I would not expect this PR to affect those tests.
Checklist
Testingwith the following:have been run on and indicate that are all passing.
has passed, using the Polaris
e3sm_submodules/Omegabaseline-pfor both the baseline (Polarise3sm_submodules/Omega) and the PR buildTesting
CTest unit tests
Polaris omega_pr regression suite
Tests added/modified/impacted
testGradRichNumandtestOneTwoOneFilter, and modified:testShearVertMixandtestTotalVertMix.Documentation
Documentation has be built locally here and charges are as expected: https://portal.nersc.gov/project/e3sm/katsmith/omega_docs/html/devGuide/VerticalMixingCoeff.html