High-order horizontal tracer advection produces incorrect results in time-stepping tests.

[hyungyukang]

I’ve been testing the time stepping PR (#392) across several test cases. During this process, I found that high-order horizontal tracer advection in Omega (HorzTracerFluxOrder: 3 and 4) appears to produce incorrect results compared to the second-order case when active tracers, such as temperature and salinity, are used to compute SpecVol through the EOS during time stepping.

The figure below shows the total kinetic energy from the QU240 spin-down test. In the legend, H denotes the horizontal tracer advection order, and V denotes the vertical advection order.

Since this test has no external forcing or eddy parameterization and includes only advection with strong viscosity (del2, del4), the total KE is expected to decrease and eventually saturate. However, unlike the other runs, including H2 in Omega and H2, H3, H4 in MPAS-Ocean, the H3 and H4 cases in Omega begin generating KE after about day 7, which is unexpected.

For reference, temperature field (k = 0) at day 30 is shown for Omega H4, Omega H2, MPAS-O H4, along with the difference field between Omega H4 and Omega H2.

Although only the horizontal tracer advection order is changed, errors in tracer advection can affect KE. Temperature and salinity are used to compute specific volume through the EOS, and the updated specific volume then feeds back onto the pressure gradient force in the momentum tendency through the following pathway:

Specific volume -> Pseudo layer thickness -> Pressure -> PGF in the momentum tendency

This issue does not show up in tracer transport tests, such as rotation_2d in Polaris, because those tests do not compute the PGF from temperature and salinity.

Please refer to the comment below for the updated results (#396 (comment)).
The figure below shows the total kinetic energy from the 4 km baroclinic channel test.

Even with a flat bottom, H3 and H4 in Omega still produce results that differ from all the other cases. This suggests that the high-order tracer advection issue in Omega is not specific to bottom topography.

Again, temperature field (k=0) at day 20:

[hyungyukang]

For anyone who wants to reproduce the issue using the QU240 spin-down test, I have prepared the test case on Frontier here:

/lustre/orion/cli115/world-shared/hgkang/data4others/Omega/spin_down_QU240

Please note:

1. Omega should be compiled with the time-stepping PR (Enable Omega time stepping with diagnostic computations #392).
2. A workaround for converting Geometric layer thickness to Pseudo thickness should be added to OceanInit.cpp. Please refer to this comment: Enable Omega time stepping with diagnostic computations #392 (comment)

[hyungyukang]

Update:
The previous MPAS-Ocean results for the baroclinic channel test were generated using incorrect configurations. The updated results are shown below. (QU240 was run correctly.)

For each horizontal advection order, Omega and MPAS-Ocean show very similar total KE trends through day 20, and they begin to diverge afterward. This behavior seems reasonable. The day-20 surface temperature fields for H4 from both models also now look very similar.

Based on these results, I can narrow down the high-order horizontal tracer advection issue to something related to bottom topography, more specifically vertical layer extension and masking.

[xylar]

Based on these results, I can narrow down the high-order horizontal tracer advection issue to something related to bottom topography, more specifically vertical layer extension and masking.

Just so I follow, you're asserting this because the baroclinic channel doesn't have topography? So things seem to be fine without it?

[hyungyukang]

@xylar , The high-order tracer advection looks fine with flat bottom topography like the baroclinic channel test. The issue still exists when the bottom topography is not flat like QU240. I will run the overflow test as well to verify this.

[hyungyukang]

I investigated this issue and think I found the cause.
The high-order masking array in the vertical direction was not constructed, which affected some grid points during high-order tracer advection. Near the boundaries, the scheme should stay at low order, but this was not happening with non-flat bottom topography.

To fix this issue, I created the vertical high-order masking array using BoundaryCell, following MPAS-Ocean, and updated the related computation in TendencyTerm.h.

The results after the fix are shown below (the fixed cases shown in red):

• QU240 spin-down test

After the fix, high-order tracer advection no longer produces unrealistic KE.

• Overflow test

Omega results with high-order tracer advection are now very close to the MPAS-Ocean results.

After verifying the fix with a few more experiments, I will open a PR to address this issue.