From 95450f3cf3e2218a80cb6567024579b7b6bd11a9 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 16 Oct 2025 15:43:27 -0700 Subject: [PATCH 1/9] Add option to disallow terrestrial outflow Desirable for AIS or for the situation where we do not consider thermodynamics. --- .../src/Registry_subglacial_hydro.xml | 4 +++ .../mode_forward/mpas_li_subglacial_hydro.F | 30 +++++++++++-------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml index 024f23a9bb49..2976c518595f 100644 --- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml +++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml @@ -129,6 +129,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" /> + diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index bf9616b4339b..a4b6e0263de2 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -879,6 +879,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 @@ -909,6 +910,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) @@ -972,19 +974,21 @@ 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 + 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 + else + 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 From d51a403927ae947bb7ce70741d574cb7c46c7082 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Thu, 16 Oct 2025 20:15:12 -0700 Subject: [PATCH 2/9] Move assignment of melt input outside of SGH subcycling loop 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. --- .../mode_forward/mpas_li_subglacial_hydro.F | 39 +++++++++++-------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index a4b6e0263de2..dcb51a05c92d 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -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 ! ============= @@ -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 @@ -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 ! ============= From f20a97b1cc7b07aff5cd8583634da3bb80223b01 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Tue, 21 Oct 2025 12:13:55 -0700 Subject: [PATCH 3/9] Add more timers to hydrology solve Useful for profiling where performance gains might occur --- .../src/mode_forward/mpas_li_subglacial_hydro.F | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index dcb51a05c92d..7c0a0306c517 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -489,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)) @@ -525,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)) @@ -546,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 @@ -562,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) @@ -598,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)) @@ -618,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)) @@ -631,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) @@ -667,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 @@ -686,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)) @@ -716,6 +731,7 @@ subroutine li_SGH_solve(domain, err) block => block % next end do + call mpas_timer_stop("calc pressure") deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH From c2d49443fdcf7f95b819b84b4c025c73d63b319f Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 27 Oct 2025 12:52:56 -0700 Subject: [PATCH 4/9] Prevent uphill outflow at terrestrial margins 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. --- .../mode_forward/mpas_li_subglacial_hydro.F | 24 +++++++++++++++---- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index 7c0a0306c517..91c1f8fcc579 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -1001,13 +1001,27 @@ subroutine calc_edge_quantities(block, err) 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) + 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 - hydropotentialBaseSlopeNormal(iEdge) = waterPressure(cell2) / dcEdge(iEdge) - hydropotentialSlopeNormal(iEdge) = (rho_water * gravity * waterThickness(cell2) + waterPressure(cell2)) / dcEdge(iEdge) + 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 + else ! if disallow terrestrial outflow hydropotentialBaseSlopeNormal(iEdge) = 0.0_RKIND hydropotentialSlopeNormal(iEdge) = 0.0_RKIND waterPressureSlopeNormal(iEdge) = 0.0_RKIND From ece5ce180b9965b9e54997ee2ee5e048fad4a20a Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Sat, 1 Nov 2025 20:22:23 -0700 Subject: [PATCH 5/9] disable channel melt into/out of lakes --- .../src/mode_forward/mpas_li_subglacial_hydro.F | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index 91c1f8fcc579..c32abac14d7b 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -1885,6 +1885,7 @@ 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), dimension(:), pointer :: channelArea real (kind=RKIND), dimension(:), pointer :: channelMelt @@ -1903,6 +1904,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 @@ -1926,6 +1928,7 @@ 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_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve) call mpas_pool_get_array(hydroPool, 'channelArea', channelArea) @@ -1949,6 +1952,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 @@ -1994,6 +1998,16 @@ subroutine update_channel(block, err) (channelDischarge + waterFlux * config_SGH_incipient_channel_width) & * waterPressureSlopeNormal / latent_heat_ice + ! disable channel melt into/out of lakes + do iEdge = 1, nEdgesSolve + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + if ((waterThickness(cell1) > bedRoughMax) .or. (waterThickness(cell2) > bedRoughMax)) then + channelMelt(iEdge) = 0.0_RKIND + channelPressureFreeze(iEdge) = 0.0_RKIND + endif + enddo + if (config_SGH_include_pressure_melt) then channelOpeningRate = (channelMelt - channelPressureFreeze) / rhoi else From 6cd0d58d4155546f151fc1bbd2866bfb8044a4d5 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Mon, 3 Nov 2025 19:33:26 -0800 Subject: [PATCH 6/9] Allow channel melt out of lakes but not into or within --- .../src/mode_forward/mpas_li_subglacial_hydro.F | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index c32abac14d7b..923d68295236 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -1998,11 +1998,21 @@ subroutine update_channel(block, err) (channelDischarge + waterFlux * config_SGH_incipient_channel_width) & * waterPressureSlopeNormal / latent_heat_ice - ! disable channel melt into/out of lakes + ! disable channel melt into and in lakes do iEdge = 1, nEdgesSolve cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) - if ((waterThickness(cell1) > bedRoughMax) .or. (waterThickness(cell2) > bedRoughMax)) then + if ((waterThickness(cell1) > bedRoughMax) .and. (waterThickness(cell2) > bedRoughMax)) then + ! no melt in a lake + channelMelt(iEdge) = 0.0_RKIND + channelPressureFreeze(iEdge) = 0.0_RKIND + elseif ((channelDischarge(iEdge) > 0.0_RKIND) .and. (waterThickness(cell2) <= bedRoughMax) .and. (waterThickness(cell1) > bedRoughMax)) then + ! no melt for channel feeding into lake (channel flowing from cell2 to cell1) + ! channelDischarge defined as positive for flow from cell2 to cell1 + channelMelt(iEdge) = 0.0_RKIND + channelPressureFreeze(iEdge) = 0.0_RKIND + elseif ((channelDischarge(iEdge) < 0.0_RKIND) .and. (waterThickness(cell2) > bedRoughMax) .and. (waterThickness(cell1) <= bedRoughMax)) then + ! no melt for channel feeding into lake (channel flowing from cell1 to cell2) channelMelt(iEdge) = 0.0_RKIND channelPressureFreeze(iEdge) = 0.0_RKIND endif From 3cd22f7dc09f03dba930570656f76c60bf8e9d65 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 5 Nov 2025 09:41:12 -0800 Subject: [PATCH 7/9] Remove code that disables channel melt in lakes 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. --- .../mode_forward/mpas_li_subglacial_hydro.F | 20 ------------------- 1 file changed, 20 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index 923d68295236..300942794586 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -1998,26 +1998,6 @@ subroutine update_channel(block, err) (channelDischarge + waterFlux * config_SGH_incipient_channel_width) & * waterPressureSlopeNormal / latent_heat_ice - ! disable channel melt into and in lakes - do iEdge = 1, nEdgesSolve - cell1 = cellsOnEdge(1, iEdge) - cell2 = cellsOnEdge(2, iEdge) - if ((waterThickness(cell1) > bedRoughMax) .and. (waterThickness(cell2) > bedRoughMax)) then - ! no melt in a lake - channelMelt(iEdge) = 0.0_RKIND - channelPressureFreeze(iEdge) = 0.0_RKIND - elseif ((channelDischarge(iEdge) > 0.0_RKIND) .and. (waterThickness(cell2) <= bedRoughMax) .and. (waterThickness(cell1) > bedRoughMax)) then - ! no melt for channel feeding into lake (channel flowing from cell2 to cell1) - ! channelDischarge defined as positive for flow from cell2 to cell1 - channelMelt(iEdge) = 0.0_RKIND - channelPressureFreeze(iEdge) = 0.0_RKIND - elseif ((channelDischarge(iEdge) < 0.0_RKIND) .and. (waterThickness(cell2) > bedRoughMax) .and. (waterThickness(cell1) <= bedRoughMax)) then - ! no melt for channel feeding into lake (channel flowing from cell1 to cell2) - channelMelt(iEdge) = 0.0_RKIND - channelPressureFreeze(iEdge) = 0.0_RKIND - endif - enddo - if (config_SGH_include_pressure_melt) then channelOpeningRate = (channelMelt - channelPressureFreeze) / rhoi else From b6581baf1e57cd6570dac6e3993ff4958944baf3 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 5 Nov 2025 09:44:31 -0800 Subject: [PATCH 8/9] Disable iceThicknessHydro option by default 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. --- .../mpas-albany-landice/src/Registry_subglacial_hydro.xml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml index 2976c518595f..43c644b04608 100644 --- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml +++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml @@ -78,8 +78,8 @@ possible_values="positive real number" /> - From 30b2d3d6de1ad6676c535f8d82adb0a7b79211ef Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 5 Nov 2025 09:55:02 -0800 Subject: [PATCH 9/9] Add max channel area option --- .../mpas-albany-landice/src/Registry_subglacial_hydro.xml | 4 ++++ .../src/mode_forward/mpas_li_subglacial_hydro.F | 3 +++ 2 files changed, 7 insertions(+) diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml index 43c644b04608..bd365ff1ef75 100644 --- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml +++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml @@ -112,6 +112,10 @@ description="width of sheet beneath/around channel that contributes to melt within the channel" possible_values="positive real number" /> +