Fixup ISOMIP+ land ice pressure #723

cbegeman · 2023-10-25T20:56:37Z

This PR changes the calculation of land ice pressure for ISOMIP+ thin-film test cases so that it is a function of the land ice thickness rather than land ice draft. The land ice pressure is not perfectly identical between Ocean0 and thin_film_Ocean0 test cases, mostly because in the former case we smooth the draft and in the latter we smooth ice thickness. This PR also fixes a bug wherein the thin film region was culled from the domain.

Checklist

User's Guide has been updated
API documentation in the Developer's Guide (api.rst) has any new or modified class, method and/or functions listed
Documentation has been built locally and changes look as expected
Document (in a comment titled Testing in this PR) any testing that was used to verify the changes

cbegeman · 2023-10-25T20:57:36Z

Testing

Ocean0 and thin_film_Ocean0 test cases have been run on chrys with intel, openmpi. thin_film_wetting_Ocean0 has been run through initial_state step to verify time-varying forcing.

cbegeman · 2023-10-25T20:58:54Z

@xylar Can you let me know what you think? I ended up doing something a little different from you to make the Ocean0 and thin_film_Ocean0 land ice pressure fields more similar. Open to other ideas though.

compass/ocean/iceshelf.py

compass/ocean/tests/isomip_plus/initial_state.py

xylar · 2023-10-26T07:46:42Z

compass/ocean/tests/isomip_plus/initial_state.py

-        ref_density = constants['SHR_CONST_RHOSW']
-        landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft(
-            ssh=ds.ssh, modify_mask=ds.ssh < 0., ref_density=ref_density)
+        land_ice_draft = ds.ssh


The reason I don't do this for the thin film case in Polaris is because the ocean won't want to have the SSH at the flotation height that BISICLES wants to give it. This is at least partly because BISICLES used a different reference density for the ocean (1028 instead of 1026) and partly because it has been my experience that computing the ice draft (SSH) from the land-ice pressure and the reference ocean density tends to produce the most reliable starting point for the SSH adjustment step.

I'm not sure if this will help with the state validation issues during SSH adjustment for the thin-film case but I think it might.

I discovered that I was missing a check in my Polaris computation of land_ice_draft:

land_ice_draft = np.maximum(ds_topo.bedrockTopography, land_ice_draft)

We obviously don't what to allow the ice draft to go deeper than the bedrock in the thin film region.

If that's the case, then wouldn't we want to assign the ssh to the recomputed land ice draft for both thin film and non-thin film cases?

We could certainly do that. It would dovetail better with the coupling strategy.

The reason I had been working with SSH/draft as the known quantity and land-ice pressure as the one to adjust was that this is more how we operate with static cavities with realistic topography (e.g. from BedMachine). Without wetting and drying, you can end up with too thin a water column if you are using land-ice pressure to determine SSH. But we could either adjust the grounding line if this happens. We can't adjust the minimum column thickness unless we also change the land-ice pressure because otherwise the pressure will just slam the column back down.

I modify the bottom depth to maintain min water column thickness here:

compass/compass/ocean/tests/isomip_plus/initial_state.py

Lines 156 to 167 in b977252

# Deepen the bottom depth to maintain the minimum water-column

# thickness

min_column_thickness = section.getfloat('min_column_thickness')

min_layer_thickness = section.getfloat('min_layer_thickness')

min_levels = section.getint('minimum_levels')

min_column_thickness = max(min_column_thickness,

min_levels * min_layer_thickness)

min_depth = -ds.ssh + min_column_thickness

ds['bottomDepth'] = np.maximum(ds.bottomDepth, min_depth)

print(f'Adjusted bottomDepth for '

f'{np.sum(ds.bottomDepth.values<min_depth.values)} cells '

f'to achieve minimum column thickness of {min_column_thickness}')

But I only do this for the initial state so it wouldn't guarantee for a time-varying pressure that the water column thickness would stay above the threshold without wetting and drying.

Right. In my experience, the only to have things behave robustly without wetting-and-drying has been to treat the water-column thickness as the thing you know and the land-ice pressure as the thing you can adjust (in the poorly named "ssh adjustment" step). Even if you thicken the water column in the initial condition, if the land-ice pressure is too high, it will just slam right back down and try to evacuate the water column.

cbegeman · 2023-10-26T16:21:26Z

@xylar Does the latest commit look like what you have in mind? It results in these water column thicknesses:
Ocean0:

thin_film_Ocean0:

xylar · 2023-10-26T17:05:17Z

Does the latest commit look like what you have in mind?

Yes exactly, that looks great to me! Are you happy with it thought? I definitely want to make sure we're on the same page.

cbegeman · 2023-10-26T19:59:57Z

@xylar I'm fine with that. I can add the same thing to the time-varying forcing. I assume we keep the floating fraction the same, to mimic a case where the GL is different in MALI from MPAS-O. Do you agree?

xylar · 2023-10-26T20:07:52Z

Right, I agree. We want to treat the floating fraction as something that is prescribed or supplied by MALI, so it doesn't get updated based on whether MPAS-Ocean thinks the cell is wet or dry.

cbegeman · 2023-10-31T17:19:40Z

@xylar Following up on your question earlier: Wetting and drying is on for thin-film cases, forward and ssh-adjustment steps, so that is not the issue with ssh_adjustment failures.

cbegeman · 2023-11-08T15:49:30Z

@xylar Just wanted to make sure this was still on your radar to review. No particular rush.

xylar · 2023-11-08T16:23:37Z

Yes, I think this was where we were discussing going back to one strategy where we trust landIcePressure and dig when necessary to prevent drying (for non-W&D cases). Did you want to have this reviewed as it is or for me to make that change here?

cbegeman · 2023-11-08T17:23:42Z

Yes, I think this was where we were discussing going back to one strategy where we trust landIcePressure and dig when necessary to prevent drying (for non-W&D cases). Did you want to have this reviewed as it is or for me to make that change here?

@xylar Right, that's worth chatting about. I was assuming we merge this PR and tackle the digging strategy in a separate PR.
If you want to tackle it in this PR, that's fine by me. It would be good to have you review at least the thin film approach here because that's what I'll use in my W&D testing in the very near future.

xylar · 2023-11-08T17:51:13Z

Okay, let's just focus on the thin film here and I'll make a follow-up PR to change the strategy for the non-thin-film cases. I could do a rush review tonight but it would be better if it's on my list for tomorrow.

xylar

I ran all of the 5 km test cases. I'm seeing a failure in initial_state in:

ocean/isomip_plus/planar/5km/z-star/time_varying_Ocean0

because of the missing landIceDraft variable noted below.

I'm seeing state validation failures in the performance steps of:

ocean/isomip_plus/planar/5km/z-star/thin_film_time_varying_Ocean0
ocean/isomip_plus/planar/5km/sigma/thin_film_time_varying_Ocean0
ocean/isomip_plus/planar/5km/sigma/thin_film_wetting_Ocean0
ocean/isomip_plus/planar/5km/single_layer/thin_film_drying_Ocean0
ocean/isomip_plus/planar/5km/single_layer/thin_film_time_varying_Ocean0
ocean/isomip_plus/planar/5km/single_layer/thin_film_wetting_Ocean0

This is causing the test to hang, which is especially annoying but unrelated to this PR, of course.

I'm not sure how many of these need investigation as part of this PR and which we just let be for now and investigate going forward. That's partly because there are now a lot of tests and I'm not sure I have a good handle on when they're all supposed to be working.

compass/ocean/tests/isomip_plus/initial_state.py

cbegeman · 2024-02-01T00:12:52Z

@xylar FYI: I'm finding that the run hangs at the performance step of time_varying_Ocean0. I'll ping you when I figure out why that is.

xylar · 2024-02-01T12:23:27Z

@cbegeman, I think the hanging itself isn't related to this PR. It seems to happen sometimes when MPAS-Ocean fails with state validation errors but I don't understand why this leads to hanging as opposed to MPAS-Ocean returning an error.

But it's worth trying to fix the state validation issues so the hanging doesn't happen.

cbegeman · 2024-07-05T16:56:26Z

Oops! comment posted on the wrong PR

cbegeman · 2024-10-24T16:12:57Z

@xylar I haven't yet tested all the cases, but this seems to be working now

cbegeman self-assigned this Oct 25, 2023

cbegeman added bug Something isn't working ocean labels Oct 25, 2023

xylar reviewed Oct 26, 2023

View reviewed changes

cbegeman marked this pull request as ready for review October 26, 2023 23:39

xylar reviewed Nov 9, 2023

View reviewed changes

compass/ocean/tests/isomip_plus/initial_state.py Outdated Show resolved Hide resolved

cbegeman force-pushed the ocn-fix-grounded-land-ice-pressure branch from 3902ae8 to 0a21b83 Compare November 27, 2023 19:06

xylar reviewed Nov 27, 2023

View reviewed changes

compass/ocean/tests/isomip_plus/initial_state.py Outdated Show resolved Hide resolved

xylar reviewed Nov 27, 2023

View reviewed changes

compass/ocean/tests/isomip_plus/initial_state.py Outdated Show resolved Hide resolved

cbegeman force-pushed the ocn-fix-grounded-land-ice-pressure branch from 0a21b83 to 73e20b9 Compare January 30, 2024 23:45

cbegeman force-pushed the ocn-fix-grounded-land-ice-pressure branch from 73e20b9 to aea2ec7 Compare July 5, 2024 15:30

cbegeman closed this Jul 5, 2024

cbegeman reopened this Jul 5, 2024

cbegeman added 3 commits October 23, 2024 14:02

Reorg land ice pressure function

ccd1b15

Add new land ice pressure function

9be0b8b

Use new land ice pressure function in isomip_plus

70c4d83

cbegeman added 13 commits October 23, 2024 14:02

fixup thin film culling

c3e2fc4

Remove ice density calculation

7ea8c50

Extend treatment to time-varying forcing

f4d89d5

Update docs

1fac8a5

Add landIceDraft to initial_state

2f6b1d0

Apply suggestions from code review

b547b48

Fixup time varying IC

a8de443

Fix time varying forcing for thin film cases

28d3901

fixup viz

c12e272

Make compute_land_ice_draft_from_pressure more flexible

90a3d09

Plot pressure adjustment in grounded ice region

a28a73a

Make isomip plotting more flexible

824648b

Only visualize init for thin film cases

e285b25

cbegeman force-pushed the ocn-fix-grounded-land-ice-pressure branch from aea2ec7 to e285b25 Compare October 24, 2024 16:12

cbegeman added 6 commits October 24, 2024 11:20

Add comments

e4e00fc

Add new thin film namelist options

bdcc274

Add more output vars for wetdry

5833135

Plot wetdry factor as scatter

4d550a2

Fix viz when land ice fluxes off

42a4ed1

Add isomip+ rk4 suite

82f3df5

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixup ISOMIP+ land ice pressure #723

Fixup ISOMIP+ land ice pressure #723

cbegeman commented Oct 25, 2023 •

edited

Loading

cbegeman commented Oct 25, 2023 •

edited

Loading

cbegeman commented Oct 25, 2023

xylar Oct 26, 2023

xylar Oct 26, 2023

xylar Oct 26, 2023

cbegeman Oct 26, 2023

xylar Oct 26, 2023

cbegeman Oct 26, 2023

xylar Oct 26, 2023

cbegeman commented Oct 26, 2023

xylar commented Oct 26, 2023

cbegeman commented Oct 26, 2023

xylar commented Oct 26, 2023

cbegeman commented Oct 31, 2023

cbegeman commented Nov 8, 2023

xylar commented Nov 8, 2023

cbegeman commented Nov 8, 2023

xylar commented Nov 8, 2023

xylar left a comment

cbegeman commented Feb 1, 2024

xylar commented Feb 1, 2024

cbegeman commented Jul 5, 2024

cbegeman commented Oct 24, 2024

	# Deepen the bottom depth to maintain the minimum water-column
	# thickness
	min_column_thickness = section.getfloat('min_column_thickness')
	min_layer_thickness = section.getfloat('min_layer_thickness')
	min_levels = section.getint('minimum_levels')
	min_column_thickness = max(min_column_thickness,
	min_levels * min_layer_thickness)
	min_depth = -ds.ssh + min_column_thickness
	ds['bottomDepth'] = np.maximum(ds.bottomDepth, min_depth)
	print(f'Adjusted bottomDepth for '
	f'{np.sum(ds.bottomDepth.values<min_depth.values)} cells '
	f'to achieve minimum column thickness of {min_column_thickness}')

Fixup ISOMIP+ land ice pressure #723

Are you sure you want to change the base?

Fixup ISOMIP+ land ice pressure #723

Conversation

cbegeman commented Oct 25, 2023 • edited Loading

cbegeman commented Oct 25, 2023 • edited Loading

Testing

cbegeman commented Oct 25, 2023

xylar Oct 26, 2023

Choose a reason for hiding this comment

xylar Oct 26, 2023

Choose a reason for hiding this comment

xylar Oct 26, 2023

Choose a reason for hiding this comment

cbegeman Oct 26, 2023

Choose a reason for hiding this comment

xylar Oct 26, 2023

Choose a reason for hiding this comment

cbegeman Oct 26, 2023

Choose a reason for hiding this comment

xylar Oct 26, 2023

Choose a reason for hiding this comment

cbegeman commented Oct 26, 2023

xylar commented Oct 26, 2023

cbegeman commented Oct 26, 2023

xylar commented Oct 26, 2023

cbegeman commented Oct 31, 2023

cbegeman commented Nov 8, 2023

xylar commented Nov 8, 2023

cbegeman commented Nov 8, 2023

xylar commented Nov 8, 2023

xylar left a comment

Choose a reason for hiding this comment

cbegeman commented Feb 1, 2024

xylar commented Feb 1, 2024

cbegeman commented Jul 5, 2024

cbegeman commented Oct 24, 2024

cbegeman commented Oct 25, 2023 •

edited

Loading

cbegeman commented Oct 25, 2023 •

edited

Loading