Skip to content

Conversation

@matthewhoffman
Copy link

This pull request introduces a number of developments to get the subglacial hydrology model running reliably with channels active.

Most significant changes are:

  • A bug fix to eliminate uphill outflow at terrestrial margins
  • Disabling the calculation of iceThicknessHydro which was causing ice thickness artifacts near grounding lines
  • Adding a namelist option for a maximum channel area to prevent runaway channel growth

Less important changes are:

  • An option to disable terrestrial outflow entirely. This has little impact on AIS simulations, but I left it in the code as an option.
  • Some minor performance & timer adjustments
  • Changes to disable channel melt inside lakes that were subsequently removed. This introduced new problems so I removed the changes, but I left the commits in the history for posterity

Desirable for AIS or for the situation where we do not consider
thermodynamics.
When externally provided, these melt inputs will not change during the
SGH subcycling, so they don't need to be evaluated within the subcycling
loop.  This also moves one halo update out of the subcycle loop.
Useful for profiling where performance gains might occur
This commit fixes a bug that was introduced when downhill terrestrial
outflow was modified to ignore to the downhil bedtopo slope.  This
commit ensures that there is no outflow if the bedtopo external to the
grounded ice is a higher elevation than the bedtopo under the ice.  This
was previously occurring.
This seemed to be a failed experiment as it prevents channels from ever
getting very extensive.  I'm keeping the previous commits in the history
for posterity.
In AIS simulations, having iceThicknessHydro enabled was causing massive
lake formation and lack of channelization.  On further investigation, I
found there are thickened artifacts of increased thickness adjacent to
grounding lines.  This was in effect creating berms in the
hydropotential, making it harder for water to reach the ocean.  The
iceThicknessHydro calculations should be improved to eliminate this
issue, or, if that's not possible, removed from the code.  But for now I
am disabling it by default.

Note that I created an MPAS-Tools script to
preprocess ice thickness to eliminate depressions in the surface
elevation, and this greatly improves model stability around subglacial
lakes.  That preprocessing step should be used instead of
iceThicknessHydro.  It is unclear if that will be sufficient in the case
of evolving ice thickness.  If not, the surface-depression-filling
algorithm it employs should be ported to MALI and replace the existing
iceThicknessHydro algorithm.
@matthewhoffman matthewhoffman self-assigned this Nov 5, 2025
@matthewhoffman matthewhoffman added bug Something isn't working enhancement New feature or request clean up labels Nov 5, 2025
@matthewhoffman
Copy link
Author

matthewhoffman commented Nov 5, 2025

Test run with input files is here on perlmutter: /pscratch/sd/h/hoffman2/AIS_subglacial_hydro/dist_only/hager_high_res_params/output
Note the input file had the ice thickness preprocessed to fill surface depressions using script in MPAS-Dev/MPAS-Tools#520
Also note this run uses max channel area of 100km^2, but I subsequently set the default in the PR to 500 following https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2023JF007421. The run uses the hydro params from the high-res run in Hager et al. 2022 (noting a typo in the kq value listed on p.3582 - I pulled the files from the archive to confirm)

Mass balance & time series from plotting script in MPAS-Dev/MPAS-Tools#520 :
image

Comparison of water thickness within 1% of bump height (light shading) compared to radar specularity for Thwaites:
image (25)

Comparison for more of WAIS (obs from https://royalsocietypublishing.org/doi/full/10.1098/rsta.2014.0297), lakes in peach, water thickness within 1% of bump height in light shading:
image

Copy link

@alexolinhager alexolinhager left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for addressing these issues with the hydro model Matt! I made a few minor comments/suggestions, but otherwise should be ready to merge

channelMelt(iEdge) = 0.0_RKIND
channelPressureFreeze(iEdge) = 0.0_RKIND
endif
enddo

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We had previously tried something similar to this (E3SM-Project@c01e547), but had done so by disabling all channel evolution instead of just disabling channelMelt and channelPressureFreeze over lakes. The other attempt never got merged because we still couldn't get stable channels, but this seems like a better approach, especially when paired with the other commits in this PR

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reading two more commits down the line, I think the issue with this approach is familiar and seems to be
why we discontinued it previously


<nml_option name="config_SGH_use_iceThicknessHydro" type="logical" default_value=".true." units="unitless"
description="Option to use an altered ice thickness field called iceThicknessHydro that replaces local maxima/minima in upperSurface with a mean of the cells neighbors. This option has no significant effect on the behavior of the model but makes it more stable."
<nml_option name="config_SGH_use_iceThicknessHydro" type="logical" default_value=".false." units="unitless"

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we will eventually need to port your MPAS-Tools ice surface alteration code into MALI if we want to run with evolving dynamics. The primary reason for iceThicknessHydro was the development if single-cell surface depressions that formed when running coupled simulations. This was mostly an issue when running with the higher order velocity solver, which created speckled ice thicknesses along the domain edges as the run progressed.

It seems like the MPAS-Tools code should fix that issue, so I think we should try replacing the iceThicknessHydro code with what you wrote there. Would it be better to make that it's own PR, or tack it onto this one?

channelArea = channelChangeRate * deltatSGH + channelArea
! If sheet dissipation contributes to channel, there should be no need for a minimum channel size
channelArea = max(1.0e-8_RKIND, channelArea) ! make some tiny value when it goes negative
channelArea = min(config_SGH_chnl_max_area, channelArea) ! limit channel area to config-specified max value

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it be better to just limit channelMelt and channelPressureFreeze when channelArea exceeds config_SGH_chnl_max_area? That way the channel still wouldn't grow past the max value but we also wouldn't violate mass conservation if there was no actual unaccounted melting. Probably doesn't affect the outcome but could be cleaner to not have a cell where channelMelt far exceeds the channelArea

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working clean up enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants