Make landIceDraft diagnostic rather than a forcing variable

[cbegeman]

This PR changes the default way the sea-surface height gradient is computed near land ice in MPAS-Ocean for MPAS-SeaIce. In the absence of land ice, this PR has no effect [BFB]. A new ocean config option, config_land_ice_draft_mode, determines whether to use the new default method or the current method:

• [NCC] config_land_ice_draft_mode = 'pressure-dependent' is the new default whereby landIceDraft is computed from landIcePressure. The (spatially and temporally) constant ice density to use for this computation is provided by the new ocean config option config_land_ice_rho_ocean.
• [BFB] config_land_ice_draft_mode = 'data' preserves the current method whereby landIceDraft is provided as a forcing variable. It is applied to the following meshes
  • oQU240wLI
  • ECwISC30to60E1r2
  • SOwISC12to60E2r4
  • ECwISC30to60E2r1
  • IcoswISC30E3r5
  • FRISwISC08to60E3r1
  • FRISwISC04to60E3r1
  • FRISwISC02to60E3r1
  • FRISwISC01to60E3r1
  • RRSwISC6to18E3r5

config_land_ice_draft_mode = 'pressure-dependent' provides a more accurate determination of the SSH gradient near ice shelves when the SSH is adjusted during initialization. This is now done by default for new meshes #6375, which is why this option is the default here.

config_land_ice_draft_mode = 'data' provides a more accurate determination of the SSH gradient near ice shelves when the land ice pressure is adjusted during initialization. This config option is set for all meshes with ice shelf cavities employed in e3sm tests that were created prior to #6375.

The use of landIceDraft is:

• [NCC] Computing the pressure-adjusted SSH gradient for MPAS-Seaice. Formerly, landIceDraft was computed in compass using SHR_CONST_RHOSW. To increase accuracy, a representative ocean density near ice shelf fronts should be used. The separate config option config_land_ice_rho_ocean is used for this purpose.
• [Stealth, MPAS-Ocean standalone] For wetting-and-drying configurations in which there are ice shelf cavities, it is used to update the ice shelf area over which melting occurs. In coupled MALI configurations, this area will be determined in the coupler on the MALI mesh and this code will not be used.

[NML]
[NCC]

[cbegeman]

Copied from E3SM-Ocean-Discussion#125 (comment)

v3 LR piControl results

E3SM diags IICE against piControl: https://portal.nersc.gov/project/m2136/bin/iice/iice.cgi?url1=https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.cbegeman/analysis/v3.LR.piControl-landIceForcing/e3sm_diags/atm_monthly_180x360_aave/model_vs_obs_0051-0100/viewer&label1=landIceForcing&url2=https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.golaz/E3SMv3/v3.LR.piControl/e3sm_diags/atm_monthly_180x360_aave/model_vs_obs_0051-0100/viewer&label2=piControl&category=polar&plot=&diff=0&dconfig=&

This branch:
mpas_analysis: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.cbegeman/analysis/v3.LR.piControl-landIceForcing/mpas_analysis/ts_0001-0100_climo_0051-0100/

Baseline:
mpas_analysis: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.golaz/E3SMv3/v3.LR.piControl/mpas_analysis/ts_0001-0100_climo_0051-0100/

I did not find any significant differences.

[cbegeman]

@xylar Still planning to make that namelist change. Will post after running the e3sm_cryo_developer test suite.

[cbegeman]

Further discussion of this PR is located at E3SM-Ocean-Discussion#125

[cbegeman]

@jonbob I was able to run the e3sm_cryo_developer suite successfully with these namelist changes. When you have a chance, can you run your scripts and make sure I didn't miss anything?

[jonbob]

@cbegeman - done

[xylar]

@cbegeman, I'm rerunning the control analysis above so I can run a main vs. control (with the same version of Unified). I hope you don't mind. I'll post the results once I have them. I looked through the code and it looks great. I don't have any concerns there.

[xylar]

Just confirming that everything looks like ensemble-member level variability to me. Here is the main vs. control run:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis/FAnSSIE/v3.LR.piControl-landIceForcing_vs_ctrl/ts_0001-0100_climo_0051-0100/

[xylar]

Approving based on code inspection and comparison simulation showing non-climate-changing results.

[cbegeman]

Just confirming that everything looks like ensemble-member level variability to me. Here is the main vs. control run:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis/FAnSSIE/v3.LR.piControl-landIceForcing_vs_ctrl/ts_0001-0100_climo_0051-0100/

In the spirit of full disclosure, I wanted to point out that there are sea ice tendency differences due to transport as expected because we use a different ocean density to determine the SSH gradient near ice shelves.

https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis/FAnSSIE/v3.LR.piControl-landIceForcing_vs_ctrl/ts_0001-0100_climo_0051-0100/sea_ice/index.html#antarctic_extended_sh_tendency&gid=30&pid=1
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis/FAnSSIE/v3.LR.piControl-landIceForcing_vs_ctrl/ts_0001-0100_climo_0051-0100/sea_ice/index.html#antarctic_extended_sh_thick&gid=17&pid=1

[cbegeman]

A bit more context:

For land ice-ocean coupling, we need landIcePressure to be prescribed and not to be adjusted during the initialization SSH-adjustment step. This change was made here: #6375.

As a consequence, landIceDraft, used for adjusting sea ice SSH gradients near ice shelves, is no longer brought into hydrostatic equilibrium with the initial ocean density. (The reference v3 LR simulation above uses the old approach.) Instead we compute it based on a spatially and temporally constant reference ocean density config_land_ice_rho_ocean.

At init for the v3 LR mesh, the discrepancy between config_land_ice_rho_ocean and the local ocean density at ice shelf fronts introduces an SSH adjustment error O(10) cm. We believe this is acceptable, but in simulations where the ocean density diverges significantly from the initial state, users may consider updating config_land_ice_rho_ocean and developers may consider computing the local mean ocean density in situ to use for the pressure-adjusted SSH computation.

[cbegeman]

@darincomeau and @matthewhoffman Can you review when you have a chance? Thanks!

[matthewhoffman]

@cbegeman , I looked over the changes, and based on my limited understanding of the relevant code, the changes make sense based on the intent.

@xylar , @darincomeau , and I discussed the difference plots this morning and thought some further discussion about the potential impacts of this change would be good. In particular, the changes to sea ice near the Antarctic continent, while not huge, seem to systematic and probably deserve more consideration. Given the development cycle for v3, we'd also want to be careful about any changes that could inadvertently impact HR or SORRM development (as you are fully aware 😄 ). Let's plan to discuss at the FAnSSIE Ocean Team meeting next week, or we can schedule an alternate time earlier if you prefer.

[rljacob]

Does this impact coupled simulations that are not using land ice?

[xylar]

Thanks @cbegeman! I'm going to rerun the test that failed before and I'll approve if it passes.

[cbegeman]

@xylar Thanks for re-reviewing and testing!

[xylar]

Okay, with the latest changes, I my test passed and I think CI will as well. Thanks for putting up with the back-and-forth, @cbegeman!

[darincomeau]

@cbegeman do you expect any of these recent changes to affect my SORRM test (still in the queue). Do I need to pull these changes and start again?

[cbegeman]

@cbegeman do you expect any of these recent changes to affect my SORRM test (still in the queue). Do I need to pull these changes and start again?

It seems that only the gnu compiler had an issue with the associated check. If it passes a smoke test with your v3 SORRM config then I think you're fine. It would fail on the first time step. (If you haven't done a smoke test with this branch, then I would submit one to check.)

[jonbob]

@cbegeman -- is this still NCC? For all configurations?

[darincomeau]

10 year test with SORRM done (finally), MPAS-Analysis comparison vs. control: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.dcomeau/E3SMv3_dev/20250523.v3.SORRME3r3.CRYO1850-DISMF.alfred2-landicedraft.chrysalis/mpas_analysis/vs-ctrl/yrs1-10/

[xylar]

My take:

There are a few curves that hint at changes outside of ensemble-level variability, e.g. Southern Ocean temperature:

But there's nothing in the spatial patterns to suggest to me that this is related to a change in landIceDraft. Similarly, I didn't see anything in the sea-ice output that suggested major changes. There are a few plots with a halo of differences right at ice-shelf drafts, e.g.:

But not anything to suggest that these very local differences are causing discernible changes overall.

[cbegeman]

10 year test with SORRM done (finally), MPAS-Analysis comparison vs. control: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.dcomeau/E3SMv3_dev/20250523.v3.SORRME3r3.CRYO1850-DISMF.alfred2-landicedraft.chrysalis/mpas_analysis/vs-ctrl/yrs1-10/

This simulation shows the expected change: a decrease in sea ice transport away from ice shelf fronts (the red difference line right at ice shelf fronts)

Many of the larger differences in SO properties likely are due to the fact that the Weddell polynya takes longer to develop

Based on these results, this PR is implemented as intended.

[xylar]

Ah, yes, the Weddell polynya almost certainly explains the cooler SO temperatures.

[darincomeau]

Thanks for the comments @cbegeman @xylar ! Approved based on above testing.

[matthewhoffman]

I'm approving based on the code changes since the initial review and the testing and discussion above. @cbegeman , thanks for reworking this to retain the ability for backwards compatibility.

[njeffery]

Approved based on discussion and that it's bfb for HES simulations.

[darincomeau]

I'm now seeing this PR fail SMS_D_P1024.ne30pg2_r05_SOwISC12to30E3r3.CRYO1850.chrysalis_intel with

 286: forrtl: severe (408): fort: (7): Attempt to use pointer LANDICEDRAFTFORSSH when it is not associated with a target

My previous test was done with this branch up to 2d43229 (confirmed again SMS test passes here), so the suspect problem commit is 8019c15.

[cbegeman]

@darincomeau I found that @xylar's fail #7301 (review) was due to a duplicated pointer variable between the MLD AM routine and diagnostic variables.

This test now passes both:

• SMS_D_Ld1_P8.ne4pg2_oQU480.WCYCL2010NS.chrysalis_gnu
• SMS_D_P1024.ne30pg2_r05_SOwISC12to30E3r3.CRYO1850.chrysalis_intel

Now that I have removed that pointer, we can go back to checking for association of landIceDraftForSsh rather than checking the config option config_land_ice_flux_mode. The other change I made here was removing landIceDraftForSsh from the landIceFluxesPKG because it only needs to belong to landIcePressurePKG.

[darincomeau]

Approved based on @cbegeman testing after her fix - thanks!

[jonbob]

Passes:

• e3sm_cryo_developer on chrysalis
• SMS_D_Ld1.ne30pg2_r05_IcoswISC30E3r5.WCYCL1850.chrysalis_intel.allactive-wcprod

with expected NML DIFFs.

Merged to next

[jonbob]

merged to master and expected NML DIFFs blessed