Several fixes to global ocean meshes with MALI topography

[xylar]

This merge:

• Fixes MALI ice draft calculation for grounded ice. The ice draft is determined by flotation in floating regions but is equal to the bed elevation in grounded regions. We also include sea_level (currently zero) in the calculation of the ice draft.
• Fixes renormalization of MALI topography. It is first masked to below sea level, then remapped, then renormalized by the ocean fraction (area below sea level).
• Uses the correct target scrip file for topo remapping. We need to use the appropriate unsmoothed or smoothed scrip file, whereas we were previously creating a new MPAS-Ocean scrip file that was always without smoothing.

Checklist

• Document (in a comment titled Testing in this PR) any testing that was used to verify the changes

[xylar]

Testing

I ran the ocean/global_ocean/SOwISC12to30/MALI_topo_AIS_4to20km/mesh test case after:

• setting ocean_includes_grounded=True in two places in compass/ocean/tests/global_ocean/mesh/__init__.py
• setting expand_factor = 2 in mesh.cfg

I now see, as expected:

• the draft is exactly equal to the bed both before and after remapping of the MALI topography in fully grounded regions
• the draft is never below the bed (the water column thickness is non-negative)
• the MALI topography is masked by oceanFrac = bed < 0, as expected
• the remapped MALI topography is renormalized by the remapped oceanFrac as expected: The topography does not drop toward zero at locations where oceanFrac < 1
• the remapped MALI topography in the remap_topo_smoothed step is, in fact, visibly smoother than the remapped topography in the unsmoothed step.
• the combined topography doesn't show any visible jumps at the boundary where landIceFraction has fractional values (used to interpolate between the two topography deatasets)
• the culled mesh removes masked areas of zero bathymetry and ice draft (locations where bathymetry was above sea level before masking)
• on the culled mesh, the 0.27 < oceanFracObserved < 1; areas where oceanFracObserved < 0.5 are few and likely are explained by topological requirements (e.g. adding cells to the ocean to prevent isolated vertices)

For paraview verification, see xdmf files in:

/lcrc/group/e3sm/ac.xylar/compass_1.9/chrysalis/test_20260113/fix-mali-draft-with-smoothing3/ocean/global_ocean/SOwISC12to30/MALI_topo_AIS_4to20km/mesh/*/xdmf_*/

[xylar]

@cbegeman, could you take a look at this when you have time? It likely makes sense to rebase #931 onto this branch.

[cbegeman]

@xylar Can you let me know whether this is what you expect and whether you intend that we would correct the negative water column thicknesses (red in figure below) before the digging phase?

All figures below are plotting the difference between land ice draft and bed elevation

In remap_topography_unsmoothed:

mali_topography_mod, landIceDraftObserved minus bed_elevation:

mali_topography_remapped, landIceDraftObserved minus bed_elevation:

topography_remapped, landIceDraftObserved minus bed_elevation:

[xylar]

The values before remapping in the masked region aren't meaningful. They get culled so negative values are just irrelevant noise. So I think this is expected. I didn't see any negative results after culling in my testing, so I think this is as expected. Thanks for checking.

[cbegeman]

The values before remapping in the masked region aren't meaningful. They get culled so negative values are just irrelevant noise. So I think this is expected. I didn't see any negative results after culling in my testing, so I think this is as expected. Thanks for checking.

I'm a bit confused by your comment. I agree that the red values in the topography_remapped.nc's (landIceDraftObserved minus bed_elevation) are noise, but my interpretation is that this was generated during the remapping to topography_ncremap.nc. And at least some of these negative values in the ice sheet interior wouldn't be culled, we would just dig the bed to accommodate the lower landIceDraft. I don't think it's likely to be too much lower than bed_elevation so I think it's practically not an issue. I was just wondering whether we might want to limit the landIceDraft to bed_elevation here

compass/compass/ocean/mesh/remap_topography.py

Line 316 in b253579

def _modify_remapped_bathymetry(self):

in a similar way as you've done in compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py

[cbegeman]

Here's what I mean: this is topography_culled.nc

As above, legend title is not correct, this is actually landIceDraftObserved minus bed_elevation

[xylar]

@cbegeman, sorry, I clearly misunderstood. It was late and I was trying to get in a response (on my phone) so you wouldn't have to wait until next week but that wasn't very helpful of me.

I plotted both of the source datasets that get blended together:

• bedmachine_topography_remapped.nc -- this is what the parent class RemapTopography produces. It gets renamed here so we don't clobber it in RemapMaliTopography.
• mali_topography_remapped.nc

In both cases, I see exactly zero difference between landIceDraftObserved and bed_elevation in fully grounded regions as expected. But, like you, I see nonzero, noisy values (both positive and negative) in topography_remapped.nc. This suggests to me that the problem arises in combining the bed_elevation datasets here:
https://github.com/xylar/compass/blob/8931918748dbf6523db2962bed7601db35139b63/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py#L301-L304
We don't combine the landIceDraftObserved data set in the same way, we just take the MALI version:
https://github.com/xylar/compass/blob/8931918748dbf6523db2962bed7601db35139b63/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py#L284-L286
The problem could arise from either of these. But I agree that the easiest solution is to make sure the draft can't be deeper than the bed. That would happen here:
https://github.com/xylar/compass/blob/8931918748dbf6523db2962bed7601db35139b63/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py#L273
rather than in RemapTopography. I'll add that.

[xylar]

This seems to have worked. There is some noise in the water-column thickness where it should be grounded and zero but the differences are now always positive.

[xylar]

The results looked as expected -- landiceFloatingFracObserved was never higher than oceanFracObserved and:

landIceFloatingFracObserved + landIceGroundedFracObserved - landIceFracObserved

was close to zero, though it was noisy ~+/-1e-8 around grounding lines for some reason. I don't really understand why, but it doesn't seem worth spending too much time on, since landIceGroundedFracObserved is never used.

[cbegeman]

@xylar I think this is ready to go in. Thanks for iterating with me on this. I've looked at the resulting output for the mesh task for SOwISC12to30 and IcoswISC with grounded ice areas included in the mesh. Have you looked at the output when grounded ice areas are excluded?

[xylar]

Have you looked at the output when grounded ice areas are excluded?

No I haven't looked that that. It would be good to run that quickly before I merge. Thanks for the review!

[xylar]

I ran the SOwISC12to30 mesh again, excluding the grounded region. The culled-mesh topography fields look good to me in ParaView.

[cbegeman]

@xylar Thanks!