Review of MOM section in 25km config doc

[dougiesquire]

Contributes to #904

---

Current rendered version here.

[dougiesquire]

I've tried to clarify the text, including removing a few things that I wasn't sure were correct. Happy to add them back in if we are confident about them.

I haven't really touched the sections on isopycnal mixing because I'm not familiar enough with MEKE and GM in MOM6 without spending a fair bit of time. I'm hoping @AndyHoggANU and/or @aekiss may be able to check over it. This is what the sections currently say:

Isopycnal mixing

At 25km resolution, the model begins to resolve some mesoscale eddies, but parameterisation is still needed for the unresolved part. The configuration uses a hybrid parameterisation for mesoscale eddies, combining neutral diffusion [@redi1982oceanic] and a dynamic Gent-McWilliams (GM) scheme [@gent1990isopycnal] based on an eddy kinetic energy budget.

Isopycnal thickness diffusion (Gent McWilliams diffusion)

GM is turned on via THICKNESSDIFFUSE = True. Instead of using a fixed GM thickness diffusivity (KHTH = 0.0), the Mesoscale Eddy Kinetic Energy (MEKE) scheme (USE_MEKE = True) is turned on. MEKE activates a prognostic equation for eddy kinetic energy (EKE) and a spatially varying GM streamfunction. The MEKE parameterisation is based on the work of [@jansen2015parameterization], where an EKE budget is solved. The model converts that EKE into an eddy diffusivity (GM diffusivity) via mixing-length theory. In practice, this means the thickness diffusion coefficient is not a fixed number but evolves according to local conditions. Our configuration does not feed external EKE data (EKE_SOURCE = "prog"), so the model instability growth provides the source of EKE. MEKE_BGSRC = 1.0E-13 prevents EKE from decaying to zero in very quiet regions. It serves as a floor to aid numerical stability and is analogous to a background diffusivity but in energy form. MEKE_GMCOEFF = 1.0 means the scheme converts eddy potential energy to eddy kinetic energy with 100% efficiency for the GM effect. MEKE_KHTR_FAC = 0.5 and MEKE_KHTH_FAC = 0.5 map some of the eddy energy to tracer diffusivity and lateral thickness diffusivity, respectively. So the configuration actually uses MEKE to the job of GM: flatterning isopycnals to remove available potential energy, but in a physically informed way using a local EKE prognostic variable. We use KHTH_USE_FGNV_STREAMFUNCTION = True which solves a 1D boundary value problem so the GM streamfunction is automatically smooth in the vertical and vanishes at the surface and bottom [@ferrari2010boundary]. FGNV_FILTER_SCALE = 0.1 is used to damp the eddy field noise.

By using MEKE, the model is effectively resolution-aware, as resolution increases and resolves more eddies, the diagnostic EKE and hence GM coefficient naturally reduces. At the same time , in coarser areas or higher latitudes where eddies are still under-resolved, MEKE ramps up the eddy mixing. This avoids the need for ad-hoc spatial maps of GM coefficients. By using FGNV[@ferrari2010boundary], it ensures a robust energetically consistent vertical structure.

We set RES_SCALE_MEKE_VISC = False, meaning the viscosity is not explicitly scaled by MEKE.

Isopycnal tracer mixing (Redi)

Neutral tracer diffusion is turned on with USE_NEUTRAL_DIFFUSION = True, which means that tracers are mixed primarily along surfaces of constant density, which greatly reduces spurious diapycnal mixing in stratified oceans. The coefficient for along-isopycnal tracer diffusion is set to KHTR = 50.0. This number is adopted from GFDL OM4_05 configuration. In addition, we also use USE_STORED_SLOPES = True and keep NDIFF_CONTINUOUS = True.

[chrisb13]

Suggested change

	Code-formatted text in the following sections gives the parameter values set in the `MOM_input` configuration file.
	Code-formatted text (`example`) in the following sections gives the parameter values set in the `MOM_input` configuration file.

[chrisb13]

Suggested change

There are 75 vertical layers (`NK = 75`). MOM6 supports an arbitrary Lagrangian Euler (ALE) vertical coordinate scheme [@griffies2020primer] that allows for completely general vertical coordinates. A number of coordinates are supported "out-of-the-box", including isopycnal, geopotential, terrain-following, HYCOM1 etc. This configuration uses ALE with a stretched geopotential (z-star) vertical coordinate (`USE_REGRIDDING = True`, `REGRIDDING_COORDINATE_MODE = "ZSTAR"`). The layer spacing specified via an input file (`ALE_COORDINATE_CONFIG = "FILE:ocean_vgrid.nc,interfaces=zeta"`). The deepest ocean depth is set to `MAXIMUM_DEPTH = 6000.0`.

There are 75 vertical layers (`NK = 75`). MOM6 supports an arbitrary Lagrangian Euler (ALE) vertical coordinate scheme [@griffies2020primer] that allows for completely general vertical coordinates. A number of coordinates are supported "out-of-the-box", including isopycnal, geopotential, terrain-following, HYCOM1 etc. This configuration uses ALE (`USE_REGRIDDING = True`) with a stretched geopotential (z-star) vertical coordinate (`REGRIDDING_COORDINATE_MODE = "ZSTAR"`). The layer spacing is specified via an input file (`ALE_COORDINATE_CONFIG = "FILE:ocean_vgrid.nc,interfaces=zeta"`). The maximum ocean depth is capped at 6000m (`MAXIMUM_DEPTH = 6000.0`).

[chrisb13]

Is it possible to be more explicit here?

The layer spacing is specified via an input file (ALE_COORDINATE_CONFIG = "FILE:ocean_vgrid.nc,interfaces=zeta")

Kind of thing that always gets discussed in a paper. I imagine @ezhilsabareesh8 will have some comments. I think we used Kial's tools? Maybe even just a link to the grid information?

[ezhilsabareesh8]

The vertical spacing is based on the following function

$\Delta z(z) = \Delta z_{\max}\tanh\left(-\frac{2\pi}{S_h} H_{\max}\right) + \varepsilon$

Where, $\varepsilon = 10^{-3}m$, based on

"Stewart, K.D., Hogg, A.M., Griffies, S.M., Heerdegen, A.P., Ward, M.L., Spence, P. and England, M.H., 2017. Vertical resolution of baroclinic modes in global ocean models. Ocean Modelling, 113, pp.50-65. http://dx.doi.org/10.1016/j.ocemod.2017.03.012"

[chrisb13]

Thanks @ezhilsabareesh8. Can I suggest you add that detail?

I note that currently we have:

Vertical grid
In the ocean model, we use a 75 level vertical grid unchanged from many OM2 configurations, following Stewart et al.(2017)(Stewart et al., 2017)2.

I think putting it here^ makes more sense? Then just adding a link from here. If you do want to add it here, please do it on this PR:
#912
(Just add your commit to the branch.)

Is it possible to even list what the layers are (as in what they start at)?

@dougiesquire has just been looking at how this varies with z* -- a comment in that sense might be nice too but not necessary.

[ezhilsabareesh8]

Thanks @chrisb13 , I have updated the vertical grids section in Grids.md in this PR 912

[chrisb13]

Thanks @ezhilsabareesh8 can you please add a reference to that here too? "Further details of the vertical grid can be found..." or similar.

[ezhilsabareesh8]

Thanks @chrisb13 , I added in this commit cd4944a

[chrisb13]

Suggested change

The configuration uses the "ROQUET_RHO" equation of state for seawater thermodynamic properties (`EQN_OF_STATE = "ROQUET_RHO"`). This scheme is based on TEOS-10, but uses a 75-term polynomial to compute in-situ density as a function of conservative temperature and absolute salinity, closely approximating the full TEOS-10 results [@roquet2015accurate]. The `_RHO` variant specifically fits density rather than specific volume and is well-suited for layered models. Prognostic temperature and salinity are conservative temperature and absolute salinity (`USE_CONTEMP_ABSSAL = True`), consistent with the equation of state.

The configuration uses the "ROQUET_RHO" equation of state for seawater thermodynamic properties (`EQN_OF_STATE = "ROQUET_RHO"`). This scheme is based on TEOS-10 using a 75-term polynomial to compute in-situ density as a function of conservative temperature and absolute salinity, closely approximating the full TEOS-10 results [@roquet2015accurate]. The `_RHO` variant specifically fits density rather than specific volume and is well-suited for layered models. Prognostic temperature and salinity are conservative temperature and absolute salinity (`USE_CONTEMP_ABSSAL = True`), consistent with the equation of state.

[chrisb13]

Currently we have:
https://access-om3-configs.access-hive.org.au/pr-preview-912/inputs/Initial-conditions-generation/

The workflow produces Conservative Temperature (CT) and Absolute Salinity (SA), the prognostic variables required by the TEOS-10 equation of state (EOS) for use in ACCESS-OM3.

I think it would be worth mentioning or linking to that here.

.. Here's a suggestion:
65eefba

(Part of this PR)

[chrisb13]

Suggested change

	Frazil formation in the ocean is turned on (`FRAZIL = TRUE`). This scheme works upwards through each water column, transferring heat downwards from the layer above as needed to prevent the in-situ temperature falling below the local freezing point in each layer in turn. If the top layer is below freezing, heat is extracted from the sea ice model, which grows frazil ice in response. More details are [here](https://mom6.readthedocs.io/en/main/api/generated/pages/Frazil_Ice.html).
	Frazil ice typically consists of individual particles that appear like a coarse slush. Frazil ice can contribute to the creation of first-year sea ice, evolving to form a soupy ‘grease ice’ layer on the ocean surface, which may then develop into pancake ice. Frazil formation in the ocean is turned on (`FRAZIL = TRUE`, `ENABLE_THERMODYNAMICS = True`), Frazil ice forms in the model when the in situ temperature drops below the local freezing point, taking into account the in situ salinity and pressure. This scheme works upwards through each water column, transferring heat downwards from the layer above as needed to prevent the in-situ temperature falling below the local freezing point in each layer in turn. If the top layer is below freezing, heat is extracted from the sea ice model, which grows frazil ice in response. In these configurations, we also have:
	1. `MASK_SRESTORE_UNDER_ICE = False`, i.e. disables SSS restoring under sea-ice based on a frazil criteria;
	1. `PRESSURE_DEPENDENT_FRAZIL = True` a pressure dependent freezing temperature is used when making frazil;
	1. `RECLAIM_FRAZIL = True` use any frazil heat deficit to cool any overlying layers down to the freezing point.
	More details are [here](https://mom6.readthedocs.io/en/main/api/generated/pages/Frazil_Ice.html.

A bit of background taken from here. @ezhilsabareesh8 seen that study? Alexander Babanin was involved.

[github-actions]

PR Preview
🚀 Preview of PR head commit fdd757d deployed to https://access-om3-configs.access-hive.org.au/pr-previews/911
2026-01-28 09:11 AEDT
Preview generated through the Deploy to GitHub Pages workflow run 21416215645.

[chrisb13]

Suggested change

There are 75 vertical layers (`NK = 75`). MOM6 supports an arbitrary Lagrangian Euler (ALE) vertical coordinate scheme [@griffies2020primer] that allows for completely general vertical coordinates. A number of coordinates are supported "out-of-the-box", including isopycnal, geopotential, terrain-following, HYCOM1 etc. This configuration uses ALE with a stretched geopotential (z-star) vertical coordinate (`USE_REGRIDDING = True`, `REGRIDDING_COORDINATE_MODE = "ZSTAR"`). The layer spacing specified via an input file (`ALE_COORDINATE_CONFIG = "FILE:ocean_vgrid.nc,interfaces=zeta"`). The deepest ocean depth is set to `MAXIMUM_DEPTH = 6000.0`. Further details of the vertical grid can be found [here](https://access-om3-configs.access-hive.org.au/inputs/Grids/#vertical-grid).

There are 75 vertical layers (`NK = 75`). MOM6 supports an arbitrary Lagrangian Euler (ALE) vertical coordinate scheme [@griffies2020primer] that allows for completely general vertical coordinates. A number of coordinates are supported "out-of-the-box", including isopycnal, geopotential, terrain-following, HYCOM1 etc. This configuration uses ALE with a stretched geopotential (z-star) vertical coordinate (`USE_REGRIDDING = True`, `REGRIDDING_COORDINATE_MODE = "ZSTAR"`). The layer spacing specified via an input file (`ALE_COORDINATE_CONFIG = "FILE:ocean_vgrid.nc,interfaces=zeta"`). The deepest ocean depth is set to `MAXIMUM_DEPTH = 6000.0`. Further details of the vertical grid can be found [here](../../inputs/Grids/#vertical-grid).

[chrisb13]

Thanks @ezhilsabareesh8. FYI, please use relative links in future.
e.g. here --> here

[chrisb13]

Suggested change

Sea surface salinity is restored toward a monthly climatological dataset calculated from the World Ocean Atlas 2023 available at [NOAA NCEI](https://www.ncei.noaa.gov/products/world-ocean-atlas) (`RESTORE_SALINITY = True`, `SALT_RESTORE_FILE = "salt_sfc_restore.nc"`). A piston velocity of 0.11 $m/day$ is applied to control the strength of the salinity relaxation (`FLUXCONST = 0.11`). The restoring is implemented as a virtual salt flux (`SRESTORE_AS_SFLUX = True`). This approach conserves salt overall, balanced globally by subtracting the mean flux to avoid altering global salinity (`ADJUST_NET_SRESTORE_TO_ZERO = True`). No effective limit is applied to the salinity restoring flux (`MAX_DELTA_SRESTORE = 999`). More discussion can be found in [issues/350](https://github.com/ACCESS-NRI/access-om3-configs/issues/350), [issues/325](https://github.com/ACCESS-NRI/access-om3-configs/issues/325), [issues/257](https://github.com/ACCESS-NRI/access-om3-configs/issues/257).

Sea surface salinity is restored toward a monthly climatological dataset calculated from the World Ocean Atlas 2023 available at [NOAA NCEI](https://www.ncei.noaa.gov/products/world-ocean-atlas) (`RESTORE_SALINITY = True`, `SALT_RESTORE_FILE = "salt_sfc_restore.nc"`). A piston velocity of 0.11 $m/day$ is applied to control the strength of the salinity relaxation (`FLUXCONST = 0.11`). The restoring is implemented as a virtual salt flux (`SRESTORE_AS_SFLUX = True`). This approach conserves salt overall, balanced globally by subtracting the mean flux to avoid altering global salinity (`ADJUST_NET_SRESTORE_TO_ZERO = True`). No effective limit is applied to the salinity restoring flux (`MAX_DELTA_SRESTORE = 999`). More discussion can be found in [issues/350](https://github.com/ACCESS-NRI/access-om3-configs/issues/350), [issues/325](https://github.com/ACCESS-NRI/access-om3-configs/issues/325), [issues/257](https://github.com/ACCESS-NRI/access-om3-configs/issues/257).

On the More discussion part, given we use virtual salt flux, I think it's less relevant but we could add:
#920

Possibly this too:
https://forum.access-hive.org.au/t/salt-restoring-flux-from-file-in-mom5-and-mom6/4610/17

[chrisb13]

Would it be possible to give a link(s) or some information as to which variables are output on which coordinates?

[chrisb13]

Suggested change

	An ideal age tracer is configured (`USE_IDEAL_AGE_TRACER = True`). This tracer ages at a rate of 1/year once it is isolated from the surface and is useful for understanding water mass ventilation and residence times.
	An ideal age tracer is initialised (`USE_IDEAL_AGE_TRACER = True`). This tracer, once isolated from the surface, ages at a rate of 1 per year. It is useful for understanding water mass ventilation and residence times.

[chrisb13]

Suggested change

The configuration handles the vertical mixing in the ocean surface boundary layer with the ePBL scheme rather the the traditional KPP. The ePBL scheme is an energy-based 1D turbulence closure approach that integrates a boundary layer energy budget to determine mixing coefficients. It was developed by [@reichl2018simplified] to improve upon KPP for climate simulations by including the effect of turbulent kinetic energy input and wind-driven mixing in a more physically constrained way. Relevant discussion can be found in [issues/465](https://github.com/ACCESS-NRI/access-om3-configs/issues/465), [issues/426](https://github.com/ACCESS-NRI/access-om3-configs/issues/426), [issues/373](https://github.com/ACCESS-NRI/access-om3-configs/issues/373).

The energetic Planetary Boundary Layer (ePBL) mixing scheme is an energy-based 1D turbulence closure approach that integrates a boundary layer energy budget to determine mixing coefficients. It was developed by [@reichl2018simplified] to improve upon KPP for climate simulations by including the effect of turbulent kinetic energy input and wind-driven mixing in a more physically constrained way. Relevant discussion can be found in [issues/465](https://github.com/ACCESS-NRI/access-om3-configs/issues/465), [issues/426](https://github.com/ACCESS-NRI/access-om3-configs/issues/426), [issues/373](https://github.com/ACCESS-NRI/access-om3-configs/issues/373).

[chrisb13]

I guess we could say how the effects of wave coupling are explored here:
#950

[chrisb13]

The inclusion of wave effects is expected to reduce warm SST biases by enhancing mixing under strong winds

is that what we've found in mcw @ezhilsabareesh8 @NoahDay ?

[ezhilsabareesh8]

" the Langmuir effect is parameterised via a predetermined enhancement factor using MOM6 default values. The inclusion of wave effects is expected to reduce warm SST biases by enhancing mixing under strong winds, as found in studies of Langmuir turbulence (e.g., USE_LA_LI2016 = True from [@li2016langmuir])"

This section applies only to KPP. The Langmuir mixing formulation of Li et al. (2016) is explicitly based on the KPP scheme, and recent test runs confirm a clear sensitivity of mixed-layer depths to the inclusion of wave effects. In KPP, the Langmuir effect is parameterised through the VR12-MA enhancement factor using WW3 inputs, consistent with MOM6 default settings. In contrast, wave effects in ePBL are represented through surface Stokes drift velocities from WW3 (when coupled), with non-breaking wave–induced mixing incorporated via additive or rescaling approaches rather than an explicit Langmuir enhancement factor.

When WW3 is not coupled, still ePBL does not use a Langmuir enhancement factor; instead, Langmuir mixing is represented through surface Stokes drift velocities derived from a Pierson–Moskowitz wind–wave spectrum refer here, where significant wave height is parameterised from the 10 m wind speed refer here, which is ultimately utilised in the calculation of m* by modifying the surface Langmuir number, as implemented through EPBL_LANGMUIR_SCHEME = "ADDITIVE".

[ezhilsabareesh8]

When WW3 is fully coupled to ePBL vertical mixing, we find that the resulting mixed-layer depths are unrealistic; therefore, we disable wave-induced mixing in ePBL (i.e., no wave to ocean coupling) to obtain reasonable mixed-layer depths.

[chrisb13]

@atteggiani and @anton-seaice the pr-preview on this seems to have broken:
https://access-nri.github.io/access-om3-configs/pr-preview-911

Perhaps that's related to the infra changes you've both been looking at?

(fyi @dougiesquire I notice there's some MC_25KM.md conflicts too)

[anton-seaice]

The link changed during this PR, its lost in the comments above.

Try

https://access-om3-configs.access-hive.org.au/pr-previews/911/

Not sure if youll also need to merge with latest main

[chrisb13]

Suggested change

Shear-driven mixing is parameterised use the Jackson-Hallberg-Legg shear mixing scheme [@jackson2008parameterization] using the MOM6 default critical Richardson number (Ri) and shear mixing rate (`USE_JACKSON_PARAM = True`, `RINO_CRIT = 0.25`, `SHEARMIX_RATE = 0.089`). This scheme targets mixing in stratified shear zones, effectively adding interior vertical diffusivity when the local Ri is less than the critical value. The shear is computed at cell vertices to better capture shear between adjacent grid cells (`VERTEX_SHEAR = True`). The Jackson-Hallberg-Legg parameterisation is energetically constrained: it iteratively finds a diffusivity such that the energy extracted from the mean flow equals the energy used in mixing plus that lost to dissipation (`MAX_RINO_IT = 25`, inherited from the [GFDL OM5 configuration](https://github.com/NOAA-GFDL/MOM6-examples/blob/3c1de3512e2200bfc10d9e5150715c9df76dbd30/ice_ocean_SIS2/Baltic_OM5_025/MOM_parameter_doc.all#L2088)).

Shear-driven mixing is parameterised using the Jackson-Hallberg-Legg shear mixing scheme [@jackson2008parameterization] with the MOM6 default critical Richardson number (Ri; `RINO_CRIT = 0.25`) and shear mixing rate (`USE_JACKSON_PARAM = True` and `SHEARMIX_RATE = 0.089`). This scheme targets mixing in stratified shear zones, effectively adding interior vertical diffusivity when the local Ri is less than the critical value. The shear is computed at cell vertices to better capture shear between adjacent grid cells (`VERTEX_SHEAR = True`). The Jackson-Hallberg-Legg parameterisation is energetically constrained: it iteratively finds a diffusivity such that the energy extracted from the mean flow equals the energy used in mixing plus that lost to dissipation (`MAX_RINO_IT = 25`, inherited from the [GFDL OM5 configuration](https://github.com/NOAA-GFDL/MOM6-examples/blob/3c1de3512e2200bfc10d9e5150715c9df76dbd30/ice_ocean_SIS2/Baltic_OM5_025/MOM_parameter_doc.all#L2088)). See the [MOM6 documentation for further details](https://mom6.readthedocs.io/en/main/api/generated/pages/Internal_Vert_Mixing.html#shear-driven-mixing-in-jackson).

[chrisb13]

@atteggiani, the old pr-preview isn't working anymore but the new link is. Can I delete this post?

[atteggiani]

@atteggiani, the old pr-preview isn't working anymore but the new link is. Can I delete this post?

Yes thank you

[chrisb13]

Thanks, done.

[chrisb13]

Suggested change

Ocean interior background mixing is also configured based on the [GFDL OM5 configuration](https://github.com/NOAA-GFDL/MOM6-examples/blob/3c1de3512e2200bfc10d9e5150715c9df76dbd30/ice_ocean_SIS2/Baltic_OM5_025/MOM_parameter_doc.all#L2004). A weak constant background diapycnal diffusivity is set for diapycnal mixing (`KD = 1.5E-05`), with a floor for numerical stability (`KD_MIN = 2.0e-6`). Vertical mixing is enhanced in the salt-fingering regime (`DOUBLE_DIFFUSION = True`) and Henyey-type internal wave scaling is configured with default MOM6 parameters (`HENYEY_IGW_BACKGROUND = True`, `HENYEY_N0_2OMEGA = 20.0`, `HENYEY_MAX_LAT = 73.0`). The total diffusivity increment from TKE-based schemes is capped to prevent unbounded growth of shear-based or convective mixing (`KD_MAX = 0.1`). This is a large upper bound that would only take effect in extremely unstable cases.

Ocean interior background mixing is also configured based on the [GFDL OM5 configuration](https://github.com/NOAA-GFDL/MOM6-examples/blob/3c1de3512e2200bfc10d9e5150715c9df76dbd30/ice_ocean_SIS2/Baltic_OM5_025/MOM_parameter_doc.all#L2004). A weak constant background diapycnal diffusivity is set for diapycnal mixing (`KD = 1.5E-05`; m^2 s-1), with a floor for numerical stability (`KD_MIN = 2.0e-6`; m^2 s-1). Vertical mixing is enhanced in the salt-fingering regime (`DOUBLE_DIFFUSION = True`) and Henyey-type internal wave scaling is configured with default MOM6 parameters (`HENYEY_IGW_BACKGROUND = True`, `HENYEY_N0_2OMEGA = 20.0`, `HENYEY_MAX_LAT = 73.0`). The total diffusivity increment from TKE-based schemes is capped to prevent unbounded growth of shear-based or convective mixing (`KD_MAX = 0.1`; m^2 s-1). This is a large upper bound that would only take effect in extremely unstable cases.

[chrisb13]

Feel free to use LaTeX if you remember the syntax.

[chrisb13]

Suggested change

	- no constant background viscosity (`AH = 0.0`)
	- no constant background viscosity (`AH = 0.0` default value; `m^2 s-1`)

[chrisb13]

Suggested change

	- no constant background viscosity (`KH = 0.0`)
	- no constant background viscosity (`KH = 0.0` default value; `m^4 s-1`)

[chrisb13]

Suggested change

	- no file-base background viscosity (`USE_KH_BG_2D = False`)
	- a 2-d spatial map of background viscosity is not applied (`USE_KH_BG_2D = False`)

[chrisb13]

Suggested change

	The Lapacian viscosity includes:
	The [Lapacian viscosity](https://mom6.readthedocs.io/en/main/api/generated/modules/mom_hor_visc.html#laplacian-viscosity-coefficient) includes:

[chrisb13]

This should be on a new line but it isn't on the rendered version and currently look like a run-on sentence.

Maybe this will help?

Suggested change

	The Laplacian and biharmonic coefficients are both limited locally to guarantee stability (`BOUND_KH = True`, `BETTER_BOUND_KH = True`, `BOUND_AH = True`, `BETTER_BOUND_AH = True`).
	The Laplacian and biharmonic coefficients are both limited locally to guarantee stability (`BOUND_KH = True`, `BETTER_BOUND_KH = True`, `BOUND_AH = True`, `BETTER_BOUND_AH = True`).

[chrisb13]

Suggested change

Baroclinic instability converts available potential energy (APE) stored in sloping isopycnals into eddy kinetic energy (EKE). In an eddy-permitting 25km configuration, this conversion is only partly resolved hence the model does not fully capture the eddy processes that naturally flatten isopycnals and release APE. As a result, isopycnal slopes remain steeper than they should be unless the unresolved eddy effects are parameterised. More details can be found in [@mitgcm_gmredi].

Baroclinic instability converts available potential energy (APE) stored in sloping isopycnals into eddy kinetic energy (EKE). In an eddy-permitting 25km configuration, this conversion is only partly resolved [1] hence the model does not fully capture the eddy processes that naturally flatten isopycnals and release APE. As a result, isopycnal slopes remain steeper than they should be unless the unresolved eddy effects are parameterised. More details can be found in [@mitgcm_gmredi].

Don't commit the above! Rather suggesting a reference insert at [1]
[1] the standard ref: https://www.sciencedirect.com/science/article/abs/pii/S1463500313001601

Also on this

More details can be found in [@mitgcm_gmredi]

Details of what exactly?

[minghangli-uni]

Sorry for the confusion. I meant details of the physical basis and formulation of GM.

Suggested change

Baroclinic instability converts available potential energy (APE) stored in sloping isopycnals into eddy kinetic energy (EKE). In an eddy-permitting 25km configuration, this conversion is only partly resolved hence the model does not fully capture the eddy processes that naturally flatten isopycnals and release APE. As a result, isopycnal slopes remain steeper than they should be unless the unresolved eddy effects are parameterised. More details can be found in [@mitgcm_gmredi].

Baroclinic instability converts available potential energy (APE) stored in sloping isopycnals into eddy kinetic energy (EKE). In an eddy-permitting 25km configuration, this conversion is only partly resolved hence the model does not fully capture the eddy processes that naturally flatten isopycnals and release APE. As a result, isopycnal slopes remain steeper than they should be unless the unresolved eddy effects are parameterised. More details on the physical basis and mathematical formulation of the GM parameterisation can be found in [@mitgcm_gmredi].

[chrisb13]

I found these two sentences a little confusing. I think you are referring to "the real ocean" and not the parameterisation? In the second sentence the use of the word diabatic feels counter to what you are saying about a parameterisation that is focused on "along neutral surfaces".

The GM parameterization was originally derived to apply to the adiabatic interior ocean

source

[minghangli-uni]

Good catch! The use of diabatic is indeed not correct and misleading!

Suggested change

	Mesoscale eddies also induce irreversible mixing of tracers but primarily along neutral density surfaces rather than vertically. This diabatic component needs to be represented through an isopycnal diffusion parameterisation that diffuses tracers along neutral surfaces while avoiding spurious diapycnal mixing.
	Mesoscale eddies also induce irreversible mixing of tracers but primarily along neutral density surfaces rather than vertically. In ocean models, this effect is represented using an isopycnal diffusion parameterisation, which diffuses tracers along neutral surfaces while minimising spurious diapycnal mixing.

[chrisb13]

To represent both the flattening of isopycnals and the along-isopycnal tracer mixing

On my previous comment, I think it is the "adiabatic flattening of isopycnals"? (Please correct me if I'm wrong!)

[minghangli-uni]

You are correct! But strictly quasi-adiabatic! :)

[chrisb13]

Suggested change

In MOM5, these coefficients were prescribed constants or latitude-dependent maps, but these choices are ad-hoc and not dynamically constrained. Hence MOM6 also offers the Mesoscale Eddy Kinetic Energy (`MEKE`) scheme [@mom6_mom_meke] which provides a flow-dependent, scale-aware `GM` diffusivity. `MEKE` prognostically computes an eddy kinetic energy field $E(x,y,z,t)$ from which it derives an eddy velocity scale,

In MOM5, these coefficients were prescribed constants or latitude-dependent maps, but these choices are ad-hoc and not dynamically constrained. Hence MOM6 also offers the Mesoscale Eddy Kinetic Energy (`MEKE`) scheme [@mom6_mom_meke] which provides a flow-dependent, scale-aware `GM` diffusivity. `MEKE` prognostically computes an eddy kinetic energy field $E(x,y,z,t)$ ($m^2 s^-2$) from which it derives an eddy velocity scale:

[chrisb13]

I think we should use the MOM6 notation here (as well), as in κ_o see here

Not sure if khtr has a symbol?

[minghangli-uni]

I'd prefer to use the code-level parameter name here (KHTH), so we dont need to worry about differences in the MOM6 notation.

[chrisb13]

Thanks for all your work on this @minghangli-uni! I found the new version much easier to read and I've learnt a lot from reading it. 😃

I've left some additional comments on top of my older comments above (it seems some have appeared as one line reviews and others got collated).

I briefly chatted to @AndyHoggANU late last week and his GM sensitivity work is done and he'll share it soon. @AndyHoggANU feel free to have a look at my comments in the Isopycnal mixing section. Once we've chatted, I'd still like to chase down a few of the code paths and add that information in.

So we can get this out the door (as discussed at last week's OSIT)... @minghangli-uni, once you've resolved the suggestions. Can I suggest we move the Isopycnal mixing section (and sub-sections) to a new PR. We can then proceed with a release as @anton-seaice suggested.

One minor thing, I noticed there's a tendency to flip flop between referencing the MOM docs/MIT with a hyperlink versus sometimes with a footnote. I think we should be consistent. Although, if we are going to use footnoted ones we should probably use it correctly (e.g. date accessed etc).

[chrisb13]

The docs use U_e I think it helps if we are consistent.

[minghangli-uni]

Suggested change

	U_{eddy}=\sqrt{2E}
	U_{e}=\sqrt{2E}

[chrisb13]

Do you have a reference for this?

See earlier comment about notation for khth.

[minghangli-uni]

https://mom6.readthedocs.io/en/main/api/generated/modules/mom_hor_visc.html#namespacemom-hor-visc-1section-horizontal-viscosity

[chrisb13]

MOM6 adds a diffusive flux of layer thickness to the thickness (mass continuity) equation

When you say "mass continuity", in a Boussinesq Approximation model, I think of volume?

e.g.

This is also called the incompressibility equation. In the Boussinesq approximation, volume conservation replaced mass conservation.

source

We are running mom6 in Boussinesq I imagine? (MOM6 supports a non-Boussinesq formulation)

[minghangli-uni]

mass continuity is kind of like a fluid mechanics terminology but yes in Boussinesq, that is also volume continuity.

[chrisb13]

This configuration uses MEKE (USE_MEKE = True) to provide the GM diffusivity.

I think @AndyHoggANU suggested it would be more like an upper bound? Or does it actually specify the amount of diffusion to be done by GM?

[minghangli-uni]

My understanding was MEKE computes the GM diffusivity and not just simply an upper bound.

[chrisb13]

Suggested change

This configuration uses `MEKE` (`USE_MEKE = True`) to provide the `GM` diffusivity. Because we do not supply an external EKE field (`EKE_SOURCE = "prog"`), EKE is generated internally through instability growth. `MEKE_BGSRC = 1.0E-13` prevents `EKE` from decaying to zero in very quiet regions. It serves as a floor to aid numerical stability and is analogous to a background diffusivity but in energy form. `MEKE_GMCOEFF = 1.0` means the scheme converts eddy potential energy to eddy kinetic energy with 100% efficiency for the `GM` effect. `MEKE_KHTR_FAC = 0.5` and `MEKE_KHTH_FAC = 0.5` map some of the eddy energy to tracer diffusivity and thickness diffusivity, respectively. In practice, `MEKE` supplies the dynamically varying `GM` coefficient that flattens isopycnals and removes APE in a physically informed way. We use `KHTH_USE_FGNV_STREAMFUNCTION = True`, which solves a 1-D boundary-value problem to ensure the `GM` streamfunction is smooth in the vertical and vanishes at the surface and bottom [@ferrari2010boundary]. `FGNV_FILTER_SCALE = 0.1` applies light spatial filtering to reduce noise in the diagnosed streamfunction. We set `RES_SCALE_MEKE_VISC = False`, meaning viscosity is not explicitly scaled by `MEKE`.

This configuration uses `MEKE` (`USE_MEKE = True`) to provide the `GM` diffusivity. Because we do not supply an external EKE field (`EKE_SOURCE` is not equal to "file"`), EKE is generated internally through instability growth, i.e., solving the EKE equation (`EKE_SOURCE = "prog"`). `MEKE_BGSRC = 1.0E-13` prevents `EKE` from decaying to zero in very quiet regions. It serves as a floor to aid numerical stability and is analogous to a background diffusivity but in energy form. `MEKE_GMCOEFF = 1.0` means the scheme converts eddy potential energy to eddy kinetic energy with 100% efficiency for the `GM` effect. `MEKE_KHTR_FAC = 0.5` and `MEKE_KHTH_FAC = 0.5` map some of the eddy energy to tracer diffusivity and thickness diffusivity, respectively. In practice, `MEKE` supplies the dynamically varying `GM` coefficient that flattens isopycnals and removes APE in a physically informed way. We use `KHTH_USE_FGNV_STREAMFUNCTION = True`, which solves a 1-D boundary-value problem to ensure the `GM` streamfunction is smooth in the vertical and vanishes at the surface and bottom [@ferrari2010boundary]. `FGNV_FILTER_SCALE = 0.1` applies light spatial filtering to reduce noise in the diagnosed streamfunction. We set `RES_SCALE_MEKE_VISC = False`, meaning viscosity is not explicitly scaled by `MEKE`.

[chrisb13]

removing the need for ad-hoc spatial maps of GM coefficients

perhaps re-phrase without the use of the word "ad-hoc". Perhaps less dynamically aware?

[minghangli-uni]

How about this?

Suggested change

	With `MEKE`, MOM6 becomes explicitly resolution-aware: as horizontal resolution increases and more eddy processes are partially resolved, the diagnosed `GM` coefficient naturally decreases; conversely, `GM` strengthens where eddies remain under-resolved (e.g. high latitudes), removing the need for ad-hoc spatial maps of `GM` coefficients.
	With `MEKE`, MOM6 becomes explicitly resolution-aware: as horizontal resolution increases and more eddy processes are partially resolved, the diagnosed `GM` coefficient naturally decreases; conversely, `GM` strengthens where eddies remain under-resolved (e.g. high latitudes), reducing the need for prescribed spatially varying `GM` coefficients.

[chrisb13]

The along-isopycnal diffusivity is set to KHTR = 50.0

Some of @AndyHoggANU's sensitivity work suggested that this might not be occurring so we should double check.

[minghangli-uni]

Comments on this @AndyHoggANU ?

[dougiesquire]

@chrisb13, you added a link to the MOM6-examples/ice_ocean_SIS2/OM4_025 configuration to a comment about our configuration.

[chrisb13]

Oh sorry, my mistake, I think I meant to point to ours.

[dougiesquire]

Can I suggest we move the Isopycnal mixing section (and sub-sections) to a new PR. We can then proceed with a release as @anton-seaice suggested.*

@chrisb13 can you please clarify what you mean here. Do you mean:

• remove the section altogether for the release and open a new PR to add it with the suggestions from this PR?
• leave the section as it currently is (prior to this PR) and open a new PR with the suggested changes from this PR?

@minghangli-uni, all but three of @chrisb13's review comments are to changes I made in this PR (all but the Isopycnal mixing section). I suggest I respond to these. We can decide what to do with the comments on the Isopycnal mixing section pending @chrisb13's response to my question above.

[chrisb13]

can you please clarify what you mean here

Based on the summary, I think you might have missed my previous line? Emphasised below:

@minghangli-uni, once you've resolved the suggestions. Can I suggest we move the Isopycnal mixing section (and sub-sections) to a new PR.

As in, I'm suggesting @minghangli-uni can implement or not the suggestions and then it gets moved to a new PR. Make sense?

[dougiesquire]

As in, I'm suggesting @minghangli-uni can implement or not the suggestions and then it gets moved to a new PR. Make sense?

What gets moved to a new PR?

[chrisb13]

Can I suggest we move the Isopycnal mixing section (and sub-sections) to a new PR.

This was just a means of getting a release out as discussed last week. Happy to chat if still unclear?

[dougiesquire]

@chrisb13 and I just had a chat. The plan is:

• @dougiesquire respond to review comments on his changes
• @minghangli-uni respond to review comments on his changes
• We all re-review
• We copy the Isopycnal mixing section into a new branch, delete it from this branch and merge
• We open a new PR to add the Isopycnal mixing section

[minghangli-uni]

@minghangli-uni respond to review comments on his changes

I’ve responded to the comments relevant to my changes.

We copy the Isopycnal mixing section into a new branch, delete it from this branch and merge
We open a new PR to add the Isopycnal mixing section

Happy to do this once the above items are addressed.