Skip to content

Hydro model channel improvements #156

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
+92 −31
Open

Hydro model channel improvements #156

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
95450f3
Add option to disallow terrestrial outflow
matthewhoffman Oct 16, 2025
d51a403
Move assignment of melt input outside of SGH subcycling loop
matthewhoffman Oct 17, 2025
f20a97b
Add more timers to hydrology solve
matthewhoffman Oct 21, 2025
c2d4944
Prevent uphill outflow at terrestrial margins
matthewhoffman Oct 27, 2025
ece5ce1
disable channel melt into/out of lakes
matthewhoffman Nov 2, 2025
6cd0d58
Allow channel melt out of lakes but not into or within
matthewhoffman Nov 4, 2025
3cd22f7
Remove code that disables channel melt in lakes
matthewhoffman Nov 5, 2025
b6581ba
Disable iceThicknessHydro option by default
matthewhoffman Nov 5, 2025
30b2d3d
Add max channel area option
matthewhoffman Nov 5, 2025
a89e8ef
Change from max channel area to area where melt shuts off
matthewhoffman Nov 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 10 additions & 2 deletions components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,8 @@
possible_values="positive real number"
/>

<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"
Copy link

@alexolinhager alexolinhager Nov 6, 2025

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?

Copy link
Author

@matthewhoffman matthewhoffman Nov 10, 2025

Choose a reason for hiding this comment

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

Thanks for reminding me about the noise along the boundaries of regional domains with the higher-order advection. I think change to iceThicknessHydro should occur in a separate PR. The depression-filling algorithm in my standalone script will be a decent amount of work to port from python to fortran, I think, because it currently relies on some python libraries for a lot of the work. Alternatively, we might be able to keep the existing iceThicknessHydro but disable modifications near the GL. And maybe there is a fix to the HO advection scheme to ensure it reduces to first-order near the boundaries of regional domains (cc: @trhille). With any of those options, I think we should tackle them separately from this PR.

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. NOTE: Current implementation causes unrealistic bumps in thickness near the grounding line in many places, which prevents realistic outflow. Enable this option with care."
possible_values=".true. or .false."
/>

Expand Down Expand Up @@ -112,6 +112,10 @@
description="width of sheet beneath/around channel that contributes to melt within the channel"
possible_values="positive real number"
/>
<nml_option name="config_SGH_chnl_area_shutoff" type="real" default_value="500.0" units="m^2"
description="channel area above which channel melting is disabled for stability reasons."
possible_values="positive real number"
/>
<nml_option name="config_SGH_include_pressure_melt" type="logical" default_value=".false." units="none"
description="whether to include the pressure melt term in the rate of channel opening"
possible_values=".true. or .false."
Expand All @@ -129,6 +133,10 @@
description="number of iterations to smooth waterPressure over when calculating waterPressureSlopeNormal. Used only to keep channelPressureFreeze stable and will not affect other aspects of the model that rely on waterPressure."
possible_values="positive integer or zero"
/>
<nml_option name="config_SGH_allow_terrestrial_outflow" type="logical" default_value=".true." units="none"
description="whether to allow discharge out at terrestrial margins"
possible_values=".true. or .false."
/>
</nml_record>

<!-- ======================================================================= -->
Expand Down
111 changes: 82 additions & 29 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -396,18 +396,6 @@ subroutine li_SGH_solve(domain, err)
block => block % next
end do


! =============
! =============
! =============
! subcycle hydro model until master dt is reached
! =============
! =============
! =============
timecycling: do while (timeLeft > 0.0_RKIND)
numSubCycles = numSubCycles + 1


! =============
! Calculate time-varying forcing, as needed
! =============
Expand Down Expand Up @@ -436,10 +424,6 @@ subroutine li_SGH_solve(domain, err)
err = ior(err, 1)
endif

! SHMIP forcing will override basalMeltInput.
call shmip_timevarying_forcing(block, err_tmp)
err = ior(err, err_tmp)

block => block % next
end do

Expand All @@ -453,6 +437,29 @@ subroutine li_SGH_solve(domain, err)
endif



! =============
! =============
! =============
! subcycle hydro model until master dt is reached
! =============
! =============
! =============
timecycling: do while (timeLeft > 0.0_RKIND)
numSubCycles = numSubCycles + 1


! apply shmip forcing if needed
block => domain % blocklist
do while (associated(block))
! SHMIP forcing will override basalMeltInput from above
call shmip_timevarying_forcing(block, err_tmp)
err = ior(err, err_tmp)

block => block % next
end do


! =============
! Update till water layer thickness
! =============
Expand Down Expand Up @@ -482,6 +489,7 @@ subroutine li_SGH_solve(domain, err)
! =============
! Calculate edge quantities and advective fluxes
! =============
call mpas_timer_start("calc edge quantities")
block => domain % blocklist
do while (associated(block))

Expand Down Expand Up @@ -518,12 +526,14 @@ subroutine li_SGH_solve(domain, err)

block => block % next
end do
call mpas_timer_stop("calc edge quantities")


! =============
! Update Channel fields
! =============
if (config_SGH_chnl_active) then
call mpas_timer_start("update channel")
block => domain % blocklist
do while (associated(block))

Expand All @@ -539,13 +549,16 @@ subroutine li_SGH_solve(domain, err)
call mpas_dmpar_field_halo_exch(domain, 'channelMelt')
call mpas_dmpar_field_halo_exch(domain, 'channelPressureFreeze')
call mpas_timer_stop("halo updates")
call mpas_timer_stop("update channel")
endif


! =============
! Calculate adaptive time step
! =============
call mpas_timer_start("calc SGH dt")
call check_timestep(domain, timeLeft, numSubCycles, err_tmp)
call mpas_timer_stop("calc SGH dt")
err = ior(err, err_tmp)
if (err > 0) then
exit timecycling
Expand All @@ -555,6 +568,7 @@ subroutine li_SGH_solve(domain, err)
! =============
! Compute flux divergence
! =============
call mpas_timer_start("compute flux div")
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
Expand Down Expand Up @@ -591,12 +605,14 @@ subroutine li_SGH_solve(domain, err)
call mpas_timer_start("halo updates")
call mpas_dmpar_field_halo_exch(domain, 'divergence')
call mpas_timer_stop("halo updates")
call mpas_timer_stop("compute flux div")


! =============
! Update channel area now that we have time step
! =============
if (config_SGH_chnl_active) then
call mpas_timer_start("evolve channel")
block => domain % blocklist
do while (associated(block))

Expand All @@ -611,11 +627,13 @@ subroutine li_SGH_solve(domain, err)
call mpas_dmpar_field_halo_exch(domain, 'channelMeltInputCell')
call mpas_dmpar_field_halo_exch(domain, 'channelArea')
call mpas_timer_stop("halo updates")
call mpas_timer_stop("evolve channel")
endif

! =============
! Calculate total grounding line discharges
! =============
call mpas_timer_start("calc SGH GL flux")
block => domain % blocklist
do while (associated(block))

Expand All @@ -624,10 +642,12 @@ subroutine li_SGH_solve(domain, err)

block => block % next
end do
call mpas_timer_stop("calc SGH GL flux")

! =============
! Update water layer thickness
! =============
call mpas_timer_start("update water layer thickness")
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
Expand Down Expand Up @@ -660,6 +680,7 @@ subroutine li_SGH_solve(domain, err)

block => block % next
end do
call mpas_timer_stop("update water layer thickness")

! =============
! Calculate pressure field
Expand All @@ -679,6 +700,7 @@ subroutine li_SGH_solve(domain, err)
! ordering choice is only needed to support channels. Without it, runs with
! channels would use an out of date waterThickness in the hydropotential and
! do not restart correctly due to the order of operations mismatch.
call mpas_timer_start("calc pressure")
block => domain % blocklist
do while (associated(block))

Expand Down Expand Up @@ -709,6 +731,7 @@ subroutine li_SGH_solve(domain, err)

block => block % next
end do
call mpas_timer_stop("calc pressure")

deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH

Expand Down Expand Up @@ -879,6 +902,7 @@ subroutine calc_edge_quantities(block, err)
real (kind=RKIND), pointer :: bedRoughMax
real (kind=RKIND) :: conduc_coeff_wtd
character (len=StrKIND), pointer :: config_SGH_tangent_slope_calculation
logical, pointer :: config_SGH_allow_terrestrial_outflow
real (kind=RKIND), pointer :: config_sea_level
real (kind=RKIND), pointer :: rhoo
integer, pointer :: nEdges
Expand Down Expand Up @@ -909,6 +933,7 @@ subroutine calc_edge_quantities(block, err)
call mpas_pool_get_config(liConfigs, 'config_SGH_conduc_coeff_drowned', conduc_coeff_drowned)
call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness_max', bedRoughMax)
call mpas_pool_get_config(liConfigs, 'config_SGH_tangent_slope_calculation', config_SGH_tangent_slope_calculation)
call mpas_pool_get_config(liConfigs, 'config_SGH_allow_terrestrial_outflow', config_SGH_allow_terrestrial_outflow)
call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)

Expand Down Expand Up @@ -972,19 +997,35 @@ subroutine calc_edge_quantities(block, err)
! This one-sided implementation also creates outflowing conditions at terrestrial boundary
do iEdge = 1, nEdges
if (hydroTerrestrialMarginMask(iEdge) == 1) then
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
if (li_mask_is_grounded_ice(cellMask(cell1))) then ! cell2 is the icefree cell - replace phi there with cell1 Phig

hydropotentialBaseSlopeNormal(iEdge) = - waterPressure(cell1) / dcEdge(iEdge)
hydropotentialSlopeNormal(iEdge) = - (rho_water * gravity * waterThickness(cell1) + waterPressure(cell1)) / dcEdge(iEdge)

else ! cell1 is the icefree cell - replace phi there with cell2 Phig

hydropotentialBaseSlopeNormal(iEdge) = waterPressure(cell2) / dcEdge(iEdge)
hydropotentialSlopeNormal(iEdge) = (rho_water * gravity * waterThickness(cell2) + waterPressure(cell2)) / dcEdge(iEdge)

endif ! which cell is icefree
if (config_SGH_allow_terrestrial_outflow) then
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
if (li_mask_is_grounded_ice(cellMask(cell1))) then ! cell2 is the icefree cell - replace phi there with cell1 Phig
if (bedTopography(cell2) < bedTopography(cell1)) then
hydropotentialBaseSlopeNormal(iEdge) = - waterPressure(cell1) / dcEdge(iEdge)
hydropotentialSlopeNormal(iEdge) = - (rho_water * gravity * waterThickness(cell1) + waterPressure(cell1)) / dcEdge(iEdge)
else
! If bare ground elevation is higher than elevation beneath the ice, no flux
hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND
hydropotentialSlopeNormal(iEdge) = 0.0_RKIND
waterPressureSlopeNormal(iEdge) = 0.0_RKIND
endif
else ! cell1 is the icefree cell - replace phi there with cell2 Phig
if (bedTopography(cell1) < bedTopography(cell2)) then
hydropotentialBaseSlopeNormal(iEdge) = waterPressure(cell2) / dcEdge(iEdge)
hydropotentialSlopeNormal(iEdge) = (rho_water * gravity * waterThickness(cell2) + waterPressure(cell2)) / dcEdge(iEdge)
else
! If bare ground elevation is higher than elevation beneath the ice, no flux
hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND
hydropotentialSlopeNormal(iEdge) = 0.0_RKIND
waterPressureSlopeNormal(iEdge) = 0.0_RKIND
endif
endif ! which cell is icefree
else ! if disallow terrestrial outflow
hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND
hydropotentialSlopeNormal(iEdge) = 0.0_RKIND
waterPressureSlopeNormal(iEdge) = 0.0_RKIND
endif ! if allow terrestrial outflow
endif ! if edge of grounded ice
end do

Expand Down Expand Up @@ -1844,7 +1885,9 @@ subroutine update_channel(block, err)
real (kind=RKIND), pointer :: creep_coeff
real (kind=RKIND), pointer :: rhoi
real (kind=RKIND), pointer :: config_SGH_incipient_channel_width
real (kind=RKIND), pointer :: bedRoughMax
logical, pointer :: config_SGH_include_pressure_melt
real (kind=RKIND), pointer :: config_SGH_chnl_area_shutoff
real (kind=RKIND), dimension(:), pointer :: channelArea
real (kind=RKIND), dimension(:), pointer :: channelMelt
real (kind=RKIND), dimension(:), pointer :: channelPressureFreeze
Expand All @@ -1862,6 +1905,7 @@ subroutine update_channel(block, err)
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: channelDiffusivity
real (kind=RKIND), dimension(:), pointer :: waterThicknessEdgeUpwind
real (kind=RKIND), dimension(:), pointer :: waterThickness
integer, dimension(:), pointer :: waterFluxMask
integer, dimension(:), pointer :: hydroMarineMarginMask
integer, dimension(:), pointer :: edgeMask
Expand All @@ -1885,6 +1929,8 @@ subroutine update_channel(block, err)
call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_creep_coefficient', creep_coeff)
call mpas_pool_get_config(liConfigs, 'config_SGH_incipient_channel_width', config_SGH_incipient_channel_width)
call mpas_pool_get_config(liConfigs, 'config_SGH_include_pressure_melt', config_SGH_include_pressure_melt)
call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness_max', bedRoughMax)
call mpas_pool_get_config(liConfigs, 'config_SGH_chnl_area_shutoff', config_SGH_chnl_area_shutoff)
call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
call mpas_pool_get_array(hydroPool, 'channelArea', channelArea)
Expand All @@ -1908,6 +1954,7 @@ subroutine update_channel(block, err)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(hydroPool, 'channelDiffusivity', channelDiffusivity)
call mpas_pool_get_array(hydroPool, 'waterThicknessEdgeUpwind', waterThicknessEdgeUpwind)
call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
! Calculate terms needed for opening (melt) rate

Expand Down Expand Up @@ -1949,6 +1996,12 @@ subroutine update_channel(block, err)
channelMelt = (abs(channelDischarge * hydropotentialSlopeNormal) & ! channel dissipation
+ abs(waterFlux * hydropotentialSlopeNormal * config_SGH_incipient_channel_width) & !some sheet dissipation
) / latent_heat_ice
! disable channel melting above area limit
do iEdge = 1, nEdgesSolve
if (channelArea(iEdge) > config_SGH_chnl_area_shutoff) then
channelMelt = 0.0_RKIND
endif
enddo
channelPressureFreeze = -1.0_RKIND * iceMeltingPointPressureDependence * cp_freshwater * rho_water * &
(channelDischarge + waterFlux * config_SGH_incipient_channel_width) &
* waterPressureSlopeNormal / latent_heat_ice
Expand Down