Add submesoscale eddy parametrization by mwarusz · Pull Request #438 · E3SM-Project/Omega

mwarusz · 2026-06-20T05:40:22Z

This PR adds the submesoscale eddy parametrization in the form presented in Fox-Kemper et al. 2011. The parametrization class provides three main methods that

determine mixed layer depth
compute buoyancy gradient
compute submesoscale induced velocity

respectively. A new auxiliary variable named NormalTransportVelocity has been added. This velocity variable is used to advect pseudo-thickness and tracers, and can optionally include the submesocale induced velocity.

Checklist

mwarusz · 2026-06-20T05:58:13Z

The parametrization unit tests are passing. However, after turning it on in the baroclinic channel polaris test with the default parameters, I see NaNs after about 3 hours of simulation time, so there is still some debugging to be done. Moreover, user and developer documentation needs to be added.

Copilot

Pull request overview

This PR introduces a new submesoscale eddy parameterization module (Fox-Kemper
et al. 2011) into Omega’s ocean model, wires it into initialization/finalization
and tendency/auxiliary-state workflows, and adds a dedicated unit test plus
broad test-harness updates to support new analytic-field helper signatures.

Changes:

Add SubmesoEddies module (config-driven) and integrate it into ocean module
init/finalize and auxiliary-state computation.
Introduce TransportAuxVars::NormalTransportVelocity and use it in thickness
and tracer tendency computations as the transport velocity (optionally
including submesoscale-induced velocity).
Update unit test utilities and many tests to accommodate updated analytic
field setters; add SubmesoEddiesTest.

Reviewed changes

Copilot reviewed 21 out of 21 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
components/omega/test/timeStepping/TimeStepperTest.cpp	Initialize/destroy `SubmesoEddies` in the time stepper test harness.
components/omega/test/ocn/TendencyTermsTest.cpp	Update analytic-field lambdas to match updated `setScalar`/`setVectorEdge` signatures.
components/omega/test/ocn/TendenciesTest.cpp	Initialize/destroy `SubmesoEddies`; update analytic-field lambdas and `setScalar`/`setVectorEdge` call signatures.
components/omega/test/ocn/SubmesoEddiesTest.cpp	New unit test covering mixed-layer depth, buoyancy gradient, and eddy velocity computations.
components/omega/test/ocn/OceanTestCommon.h	Extend `setScalar`/`setVectorEdge` helpers to pass element indices and optional Z coordinates.
components/omega/test/ocn/HorzOperatorsTest.cpp	Update analytic-field lambdas to match updated helper signatures.
components/omega/test/ocn/AuxiliaryVarsTest.cpp	Update analytic-field lambdas to match updated helper signatures.
components/omega/test/ocn/AuxiliaryStateTest.cpp	Initialize/destroy `SubmesoEddies`; update analytic-field lambdas and helper call signatures.
components/omega/test/CMakeLists.txt	Register new `SUBMESOEDDIES_TEST`.
components/omega/src/timeStepping/ForwardBackwardStepper.cpp	Pass tracer array into thickness tendency computation.
components/omega/src/ocn/Tendencies.h	Extend thickness tendency API to accept a tracer array (for aux-state computations).
components/omega/src/ocn/Tendencies.cpp	Route thickness/tracer tendencies through new aux-state compute routines and use `NormalTransportVelocity`.
components/omega/src/ocn/SubmesoEddies.h	New `SubmesoEddies` public interface.
components/omega/src/ocn/SubmesoEddies.cpp	New implementation and IO field registration for `SubmesoEddies`.
components/omega/src/ocn/OceanInit.cpp	Initialize `SubmesoEddies` during ocean module init.
components/omega/src/ocn/OceanFinal.cpp	Destroy `SubmesoEddies` during ocean finalize.
components/omega/src/ocn/auxiliaryVars/TransportAuxVars.h	New aux-var container for `NormalTransportVelocity`.
components/omega/src/ocn/auxiliaryVars/TransportAuxVars.cpp	Register/unregister IO field for `NormalTransportVelocity`.
components/omega/src/ocn/AuxiliaryState.h	Add `TransportAuxVars` member and new aux compute entry points.
components/omega/src/ocn/AuxiliaryState.cpp	Compute Brunt–Väisälä frequency, compute transport velocity (optionally + submeso), and use it for vertical pseudo-velocity.
components/omega/configs/Default.yml	Add `Omega.Submeso` configuration group and defaults.

+       KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+          const int KMin   = MinLayerEdgeBot(IEdge);
+          const int KMax   = MaxLayerEdgeTop(IEdge);
+          const int KRange = vertRangeChunked(KMin, KMax);
+
+          parallelForInner(
+              Team, KRange, INNER_LAMBDA(int KChunk) {
+                 LocPseudoThicknessAux.computeVarsOnEdge(
+                     IEdge, KChunk, PseudoThick, NormalVelEdge);
+              });


+/// Destroy instance of SubmesoEddies
+void SubmesoEddies::destroyInstance() {
+   delete Instance;
+   Instance = nullptr;
+}


mwarusz · 2026-07-15T23:22:59Z

Testing

CTest unit tests

Machine: pm-cpu / pm-gpu / frontier
Compiler: gnu / gnugpu / craygnu-mphipcc
Build type: Release
Result: All tests passed (on all machines)

Polaris `omega_pr` suite

Baseline workdir: /pscratch/sd/m/mwarusz/omega-pr-testing/submeso/baseline-gnu
Baseline build: /pscratch/sd/m/mwarusz/omega-pr-testing/submeso/baseline-gnu/build
PR build: /pscratch/sd/m/mwarusz/omega-pr-testing/submeso/build-gnu
PR workdir: /pscratch/sd/m/mwarusz/omega-pr-testing/submeso/pr-gnu
Machine: pm-cpu
Compiler: gnu
Build type: Release
Log: /pscratch/sd/m/mwarusz/omega-pr-testing/submeso/pr-gnu/polaris_omega_pr.o55955301
Result: All tests passed

mwarusz · 2026-07-15T23:29:40Z

The baroclinic channel test is now running stably with the submesoscale eddy parametrization turned on. This PR is ready for review. However, I haven't checked that the parametrization has the intended effect, and I'm not sure that I'm the best person to do that. @vanroekel could you help with this kind of testing ?

vanroekel · 2026-07-23T00:01:41Z

@mwarusz I ran a couple tests with baroclinic channel and it's too hard to tell what the effects must be. We need to add a mixed layer. Let me look and see how straight forward that might be.

vanroekel · 2026-07-24T22:23:49Z

update - I have a mixed layer variant of baroclinic channel I have a test run in the queue on frontier. Hopefully will have results by Monday

vanroekel · 2026-08-11T21:38:03Z

+                 PseudoThickML, GradBuoyML, BVFreqML);
+
+             GradBuoyML /= PseudoThickML;
+             BVFreqML /= PseudoThickML;


@mwarusz for small ML this could be be problematic, it's not clear above in the MLD code if there is a floor. maybe it's not an issue though

I agree that we should probably add a tiny value here. MPAS-O has it.

vanroekel · 2026-08-11T21:38:50Z

+                 Team, Range{MinLyrEdgeBot, IndexMLEdge},
+                 INNER_LAMBDA(int K, Real &AccumThick, Real &AccumGradBuoy,
+                              Real &AccumBVFreq) {
+                    const Real PseudoThickKm1 =


If your function starts at MinLyrEdgeBot can this lead to values lower than min layer edge?

No, this variable is set to

((K - 1) >= MinLyrEdgeBot ? MeanPseudoThickEdge(IEdge, K - 1) : 0)

which is 0 for K = MinLyrEdgeBot. In the end an average of PseudoThickKm1 and PseudoThickK is used. This means that for K = MinLyrEdgeBot you get

PseudoThickAvg = 0.5 * (0 + MeanPseudoThickEdge(IEdge, MinLyrEdgeBot));

This is a kinda hacky treatment of the first half layer. I will add a comment to explain this.

vanroekel · 2026-08-11T21:43:11Z

@mwarusz in testing I'm getting NaN for EddyVelocity quickly. CoPilot added some guards on the ML average and zero initializing and it fixed the NaN issue. I can point you at what the changes made were if you'd like

mwarusz · 2026-08-11T22:05:33Z

@vanroekel

@mwarusz in testing I'm getting NaN for EddyVelocity quickly. CoPilot added some guards on the ML average and zero initializing and it fixed the NaN issue. I can point you at what the changes made were if you'd like

Please do.

vanroekel · 2026-08-11T22:15:43Z

here it is - https://github.com/vanroekel/E3SM/tree/temp-submeso-branch

There is some extra stuff that is adding more diagnostics as well.

mwarusz · 2026-08-11T22:39:20Z

Hmm, Copilot added this check

if (PseudoThickML > Tiny && MLDepthEdge > Tiny)

but afterwards it is also explicitly checking for NaNs and skipping calculations if it finds any. It would good to see if the first check is enough.

sbrus89 · 2026-09-08T16:06:49Z

+
+   // Compute mixed layer index and  depth based on the density difference
+   // criterion
+   void computeDenMixLayerDepth(const Array2DReal &SpecVol);


@mwarusz, you were going to move this to an aux variable, correct?

Yes. I haven't worked on this PR in a while. I was planning to get back to it after the split-explicit is merged, since the transport velocity refactoring that I did there affects this PR in non-trivial ways.

Sounds good. Thanks!

mwarusz · 2026-09-21T21:46:13Z

To test this branch after rebasing on develop and fixing a couple of bugs I ran ocean/spherical/realistic_global/EC30to60E2r2/analysis_members_test for about one year using Omega and MPAS-O with the following changes to both models:

enable submesoscale eddy parametrization
use split-explicit time stepping (RK2 with 45 min time step for Omega, AB2 with 30 min time step for MPAS-O)
use monotonic third-order horizontal advection

Here are some plots comparing temperature, mixed layer depth, and eddy velocity kinetic energy at the final time. @vanroekel Do you think these look reasonable ?

Omega

MPAS-O

vanroekel · 2026-09-29T04:31:11Z

+- `Tau`: MLI timescale parameter (s).
+- `Ce`: nondimensional efficiency coefficient.
+- `LfMin`: minimum frontal width limiter (m).
+- `DsMax`: maximum grid-length limiter used in the closure (m).


somewhere, either here or in the SubmesoEddy class we should include reasonable ranges of these parameters.

vanroekel · 2026-09-30T16:12:29Z

+      int NDims = 2;
+      std::vector<std::string> DimNames(NDims);
+      DimNames[0] = "NEdges";
+      DimNames[1] = "NVertLayers";


if this is interface - should this be NVertLayersP1?

These lines are for EddyVelocity, which is not located at interfaces.

vanroekel · 2026-09-30T16:18:22Z

+
+   parallelForOuter(
+       LaunchConfig({Mesh->NEdgesAll},
+                    TeamScratch<Real>(3 * NVertLayers + 3 * NVertLayersP1)),


Could you explain the 3 * here? I don't understand where they come from

This is the number of scratch arrays used in this loop. Three scratch arrays for layer data and three for interface data. I added a comment and changed the code slightly to hopefully make it more clear.

vanroekel · 2026-09-30T16:20:37Z

+                 const int K      = MaxLyrEdgeTop + 1;
+                 const int JCell0 = CellsOnEdge(IEdge, 0);
+                 const int JCell1 = CellsOnEdge(IEdge, 1);
+                 BVFSqEdge(K)     = 0.5_Real * (BruntVaisalaFreqSq(JCell1, K) +


how is this part different than just above on L183?

It's not. In addition to BVF, the loop above interpolates fields that are located at layer mid-points (SpecVol), so it only loops over layers. This code handles the last interface interpolation for BVF. I rephrased the comment above to make it more clear.

vanroekel · 2026-09-30T16:30:15Z

@mwarusz this looks pretty good to me. Most of my comments are clarification questions. Your results look pretty good between the two. I hope the MLD difference is without submeso - with. Submeso should shallow MLD everywhere.

katsmith133 · 2026-09-30T17:44:47Z

@mwarusz this looks pretty good to me. Most of my comments are clarification questions. Your results look pretty good between the two. I hope the MLD difference is without submeso - with. Submeso should shallow MLD everywhere.

Where is this MLD difference plot? I don't see it. I see the MPAS and Omega MLD plots, but they seem like they are just straight up MLD, not differences... unless I am mistaken. I think a difference plot would be useful to see. Maybe just a difference between an average of the last 3 months of the run.

vanroekel · 2026-09-30T17:48:47Z

Oh yes, sorry @katsmith133 I thought the middle plot was a diff, but the positive didn't make sense. Now it does since it is MLD.

I agree, and would add a difference plot of MLE on and off would be super helpful in each model in addition to a MPAS-Omega diff.

mwarusz · 2026-09-30T17:55:58Z

Thanks for your suggestions @katsmith133 and @vanroekel . I will create such plots.

Co-authored-by: Luke Van Roekel <lvanroekel@lanl.gov>

katsmith133

Nice work on this @mwarusz! I just have a few questions and I think the docs still need to be updated now that the MLD calculation has been moved.

katsmith133 · 2026-09-30T18:01:46Z

   EosInstance->computeDepthMeanSpecificVolume(PseudoThickCell);

+   // compute Brunt-Vaisala freqency squared
+   EosInstance->computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, PressureMid,


I guess this is a part of Omega I don't completely understand. I see here that BVF is being computed again, but with pressure at the mid layer, does that then overwrite the previous BVF that was calculated at the interface in a different part of the code? And which one gets outputted?

I think at every time step there is only this calculation of BVF and one in implicit vertical mixing, that happens after the state update. Which one gets outputted depends on the config options, and is a known issue that we should fix (#532). As for the use of PressureMid, I assumed this is what the BVF code expects, based on the following lines

Omega/components/omega/src/ocn/Eos.h

Lines 861 to 862 in 73469eb

Real PInt =

0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1));

I see now that the call in implicit vertical mixing uses PressureInterface

Omega/components/omega/src/ocn/VertMix.cpp

Lines 717 to 718 in 73469eb

EqState->computeBruntVaisalaFreqSq(

ConservTemp, AbsSalinity, VCoord->PressureInterface, EqState->SpecVol);

Is this a bug ?

katsmith133 · 2026-09-30T19:02:42Z

+
+// Generic linear interpolation routine. This should be in OmegaMath.h or
+// something like that once it exists.
+KOKKOS_INLINE_FUNCTION Real linearInterp(Real x, Real y1, Real x1, Real y2,


@alicebarthel's frazil PR (#466) also includes a function like this that should be put into something like OmegaMath.h. Should we make a separate PR for this or just include the new file in this PR or #466?

katsmith133 · 2026-09-30T21:08:28Z

+   }
+}
+
+void SubmesoEddies::computeTimeScale() {


Does this need to be done as a separate NEdgesAll loop? Seems like this could be calculated when its needed and save the extra compute time? Maybe I am missing something

Or is calculating it once in its own loop less computation than having to calculate it every time in another loop with other things? I'm guessing that is the case since this is only called once in the init.

Yes, exactly.

katsmith133 · 2026-09-30T22:33:54Z

+                 Team, NVertLayersP1,
+                 INNER_LAMBDA(int K) { StreamFunction(K) = 0; });
+
+             const Real MLDepthEdge = Kokkos::min(DenMixLayerDepth(JCell0),


In MPAS we are taking the average of these two cells, but here its the minimum... what is the reason for the change?

MPAS-O takes the minimum for the mixed-layer index here

Omega/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F

Line 127 in de99dc4

indMLDedge = min(indMLD(cell1),indMLD(cell2))

but then, like you said, uses the average depth

Omega/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F

Line 146 in de99dc4

zMLD = 0.5_RKIND*(dThreshMLD(cell1)+dThreshMLD(cell2))

Me and @vanroekel looked at the MPAS-O code a while back and we thought it was better to be consistent here.

Co-authored-by: Luke Van Roekel <lvanroekel@lanl.gov>

mwarusz requested a review from Copilot June 20, 2026 05:59

Copilot started reviewing on behalf of mwarusz June 20, 2026 05:59 View session

Copilot AI reviewed Jun 20, 2026

View reviewed changes

mwarusz force-pushed the omega/submesoscale-eddies branch 2 times, most recently from 1d3f283 to 9825eda Compare July 14, 2026 21:14

mwarusz requested review from katsmith133 and vanroekel July 15, 2026 23:29

vanroekel reviewed Aug 11, 2026

View reviewed changes

mwarusz force-pushed the omega/submesoscale-eddies branch 2 times, most recently from 3660296 to a7bd36d Compare August 24, 2026 22:33

sbrus89 reviewed Sep 8, 2026

View reviewed changes

mwarusz and others added 7 commits September 18, 2026 17:09

Refactor setScalar and setVector

8b95f48

Add submesocale eddies

ef4d790

Add Omega user and developer docs for submesoscale eddies

c1f27e4

Add computation of Brunt-Vaisala freq in vert aux

6967128

Add submesoscale contribution to transport velocity

70eeedc

Clamp mixed layer depth

a0ab018

Make mixed layer depth an auxiliary variable

0bf23ad

mwarusz force-pushed the omega/submesoscale-eddies branch from a7bd36d to 0bf23ad Compare September 18, 2026 23:10

vanroekel reviewed Sep 26, 2026

View reviewed changes

Comment thread components/omega/doc/userGuide/SubmesoEddies.md Outdated

vanroekel reviewed Sep 26, 2026

View reviewed changes

Comment thread components/omega/doc/userGuide/SubmesoEddies.md Outdated

vanroekel reviewed Sep 26, 2026

View reviewed changes

Comment thread components/omega/doc/userGuide/SubmesoEddies.md Outdated

vanroekel reviewed Sep 28, 2026

View reviewed changes

Comment thread components/omega/src/ocn/auxiliaryVars/MixedLayerAuxVars.h Outdated

vanroekel reviewed Sep 29, 2026

View reviewed changes

vanroekel reviewed Sep 30, 2026

View reviewed changes

Update tau description in submesoscale user docs

af17239

Co-authored-by: Luke Van Roekel <lvanroekel@lanl.gov>

katsmith133 reviewed Sep 30, 2026

View reviewed changes

mwarusz and others added 4 commits September 30, 2026 16:51

Update docs to reflect move of DenMixLayerDepth to aux

7c48bb5

Return early in mixed layer depth computation for shallow columns

8ba02d1

Expand some comments in submesoscale code

a75be84

Update transport velocity description in submesoscale user docs

0f3a088

Co-authored-by: Luke Van Roekel <lvanroekel@lanl.gov>

	Real PInt =
	0.5_Real * (Pressure(ICell, K) + Pressure(ICell, K - 1));

	EqState->computeBruntVaisalaFreqSq(
	ConservTemp, AbsSalinity, VCoord->PressureInterface, EqState->SpecVol);

Conversation

mwarusz commented Jun 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mwarusz commented Jun 20, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mwarusz commented Jul 15, 2026

Testing

CTest unit tests

Polaris omega_pr suite

Uh oh!

mwarusz commented Jul 15, 2026

Uh oh!

vanroekel commented Jul 23, 2026

Uh oh!

vanroekel commented Jul 24, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mwarusz Aug 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

vanroekel commented Aug 11, 2026

Uh oh!

mwarusz commented Aug 11, 2026

Uh oh!

vanroekel commented Aug 11, 2026

Uh oh!

mwarusz commented Aug 11, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mwarusz commented Sep 21, 2026

Omega

MPAS-O

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

vanroekel commented Sep 30, 2026

Uh oh!

katsmith133 commented Sep 30, 2026

Uh oh!

vanroekel commented Sep 30, 2026

Uh oh!

mwarusz commented Sep 30, 2026

mwarusz commented Jun 20, 2026 •

edited

Loading

mwarusz commented Jun 20, 2026 •

edited

Loading

Polaris `omega_pr` suite

mwarusz Aug 11, 2026 •

edited

Loading