diff --git a/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml b/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml
old mode 100644
new mode 100755
index 622d33933cc0..50ad79158354
--- a/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml
+++ b/components/mpas-albany-landice/src/analysis_members/Registry_global_stats.xml
@@ -100,6 +100,41 @@
                 <var name="basalSpeedMax" type="real" dimensions="Time" units="m yr^{-1}"
                         description="maximum basal speed in the domain"
                 />
+                <!-- stats related to subglacial hydrology -->
+                <var name="totalSubglacialWaterVolume" type="real" dimensions="Time" units="m^3"
+                        description="total subglacial water volume"
+                />
+                <var name="totalSubglacialLakeVolume" type="real" dimensions="Time" units="m^3"
+                        description="total volume of subglacial lakes, defined as water volume exceeding bed bump height"
+                />
+                <var name="totalSubglacialLakeArea" type="real" dimensions="Time" units="m^2"
+                        description="total area of subglacial lakes"
+                />
+                <var name="totalBasalMeltInput" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total basal meltwater contributing to the subglacial hydrologic system"
+                />
+                <var name="totalExternalWaterInput" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total external meltwater contributing to the subglacial hydrologic system"
+                />
+                <var name="totalChannelMelt" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total melt rate in the subglacial hydrologic system "
+                />
+                <var name="totalDistWaterFluxMarineMargin" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total distributed subglacial water flux across marine boundaries (grounding lines or grounded marine margins)"
+                />
+                <var name="totalDistWaterFluxTerrestrialMargin" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total distributed subglacial water flux across terrestrial margins"
+                />
+                <var name="totalChnlWaterFluxMarineMargin" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total channelized subglacial water flux across marine boundaries (grounding lines or grounded marine margins)"
+                />
+                <var name="totalChnlWaterFluxTerrestrialMargin" type="real" dimensions="Time" units="kg s^{-1}"
+                        description="total channelized subglacial water flux across terrestrial margins"
+                />
+                <var name="avgFlotationFraction" type="real" dimensions="Time" units="none"
+                        description="area-weighted average of the flotation fraction under grounded ice"
+                />
+
  	</var_struct>
 	<streams>
 		<stream name="globalStatsOutput" type="output"
diff --git a/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F b/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F
old mode 100644
new mode 100755
index d7a8123bf371..2e1d05cbc390
--- a/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F
+++ b/components/mpas-albany-landice/src/analysis_members/mpas_li_global_stats.F
@@ -163,11 +163,13 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       type (mpas_pool_type), pointer :: globalStatsAM
       type (mpas_pool_type), pointer :: geometryPool
       type (mpas_pool_type), pointer :: velocityPool
+      type (mpas_pool_type), pointer :: hydroPool
 
       ! arrays, vars needed from other pools for calculations here
       real (kind=RKIND), pointer ::  deltat
       real (kind=RKIND), dimension(:), pointer ::  areaCell
       real (kind=RKIND), dimension(:), pointer ::  dvEdge
+      real (kind=RKIND), dimension(:), pointer ::  dcEdge
       real (kind=RKIND), dimension(:), pointer ::  thickness
       real (kind=RKIND), dimension(:), pointer ::  bedTopography
       real (kind=RKIND), dimension(:), pointer ::  sfcMassBalApplied
@@ -181,7 +183,18 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       real (kind=RKIND), dimension(:), pointer ::  basalSpeed
       real (kind=RKIND), dimension(:), pointer ::  fluxAcrossGroundingLine
       real (kind=RKIND), dimension(:), pointer ::  groundedToFloatingThickness
+
+      real (kind=RKIND), dimension(:), pointer ::  waterThickness
+      real (kind=RKIND), dimension(:), pointer ::  basalMeltInput
+      real (kind=RKIND), dimension(:), pointer ::  externalWaterInput
+      real (kind=RKIND), dimension(:), pointer ::  channelMelt
+      real (kind=RKIND), dimension(:), pointer ::  waterFlux
+      real (kind=RKIND), dimension(:), pointer ::  channelDischarge
+      real (kind=RKIND), dimension(:), pointer ::  waterPressure
+
       integer, dimension(:), pointer :: cellMask
+      integer, dimension(:), pointer :: hydroMarineMarginMask
+      integer, dimension(:), pointer :: hydroTerrestrialMarginMask
       integer, pointer :: nCellsSolve
       integer, pointer :: nEdgesSolve
 
@@ -189,9 +202,11 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       real (kind=RKIND), pointer :: config_sea_level
       real (kind=RKIND), pointer :: rhoi ! config_ice_density
       real (kind=RKIND), pointer :: rhow ! config_ocean_density
+      real (kind=RKIND), pointer :: bedBumpMax ! config_SGH_bed_roughness_max
+      logical, pointer :: config_SGH
 
       ! Local counters
-      integer :: k, iCell, iEdge
+      integer :: k, iCell, iEdge, nCellsGrounded
 
       ! scalars to be calculated here from global reductions
       real (kind=RKIND), pointer ::  totalIceArea, totalIceVolume
@@ -211,6 +226,17 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       real (kind=RKIND), pointer ::  groundingLineMigrationFlux
       real (kind=RKIND), pointer ::  surfaceSpeedMax
       real (kind=RKIND), pointer ::  basalSpeedMax
+      real (kind=RKIND), pointer ::  totalSubglacialWaterVolume
+      real (kind=RKIND), pointer ::  totalBasalMeltInput
+      real (kind=RKIND), pointer ::  totalExternalWaterInput
+      real (kind=RKIND), pointer ::  totalChannelMelt
+      real (kind=RKIND), pointer ::  totalSubglacialLakeVolume
+      real (kind=RKIND), pointer ::  totalSubglacialLakeArea
+      real (kind=RKIND), pointer ::  totalDistWaterFluxMarineMargin
+      real (kind=RKIND), pointer ::  totalDistWaterFluxTerrestrialMargin
+      real (kind=RKIND), pointer ::  totalChnlWaterFluxMarineMargin
+      real (kind=RKIND), pointer ::  totalChnlWaterFluxTerrestrialMargin
+      real (kind=RKIND), pointer ::  avgFlotationFraction
 
       ! scalar reductions over all blocks on this processor
       real (kind=RKIND) ::  blockSumIceArea, blockSumIceVolume
@@ -227,11 +253,23 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       real (kind=RKIND) ::  blockMaxBasalSpeed
       real (kind=RKIND) ::  blockGLflux
       real (kind=RKIND) ::  blockGLMigrationFlux
+      real (kind=RKIND) ::  blockSumSubglacialWaterVolume
+      real (kind=RKIND) ::  blockSumBasalMeltInput
+      real (kind=RKIND) ::  blockSumExternalWaterInput
+      real (kind=RKIND) ::  blockSumChannelMelt
+      real (kind=RKIND) ::  blockSumLakeVolume
+      real (kind=RKIND) ::  blockSumLakeArea
+      real (kind=RKIND) ::  blockSumGLMeltFlux
+      real (kind=RKIND) ::  blockSumTerrestrialMeltFlux
+      real (kind=RKIND) ::  blockSumChannelGLMeltFlux
+      real (kind=RKIND) ::  blockSumChannelTerrestrialMeltFlux
+      real (kind=RKIND) ::  blockSumFlotationFraction
 
       ! Local variables for calculations
+      real (kind=RKIND) ::  totalFlotationFraction
 
       ! variables for processing stats
-      integer, parameter :: kMaxVariables = 24 ! Increase if number of stats increase
+      integer, parameter :: kMaxVariables = 35 ! Increase if number of stats increase
       integer :: nVars
       real (kind=RKIND), dimension(kMaxVariables) :: reductions, sums, mins, maxes
 
@@ -257,6 +295,17 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       blockSumFaceMeltingFlux = 0.0_RKIND
       blockGLflux = 0.0_RKIND
       blockGLMigrationFlux = 0.0_RKIND
+      blockSumSubglacialWaterVolume = 0.0_RKIND
+      blockSumBasalMeltInput = 0.0_RKIND
+      blockSumExternalWaterInput = 0.0_RKIND
+      blockSumChannelMelt = 0.0_RKIND
+      blockSumLakeVolume = 0.0_RKIND
+      blockSumLakeArea = 0.0_RKIND
+      blockSumGLMeltFlux = 0.0_RKIND
+      blockSumTerrestrialMeltFlux = 0.0_RKIND
+      blockSumChannelGLMeltFlux = 0.0_RKIND
+      blockSumChannelTerrestrialMeltFlux = 0.0_RKIND
+      blockSumFlotationFraction = 0.0_RKIND
 
       ! initialize max to 0, min to large number
       blockThickMin = 100000.0_RKIND
@@ -268,6 +317,8 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
       call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi)
       call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhow)
+      call mpas_pool_get_config(liConfigs, 'config_SGH_bed_roughness_max', bedBumpMax)
+      call mpas_pool_get_config(liConfigs, 'config_SGH', config_SGH)
 
       ! loop over blocks
       block => domain % blocklist
@@ -278,6 +329,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_subpool(block % structs, 'globalStatsAM', globalStatsAMPool)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
          call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
+         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
 
          ! get values and arrays from standard pools
          call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
@@ -285,6 +337,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(meshPool, 'deltat', deltat)
          call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
          call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
+         call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
@@ -299,6 +352,17 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(velocityPool, 'surfaceSpeed', surfaceSpeed)
          call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed)
          call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine)
+         if (config_SGH) then
+            call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
+            call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
+            call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
+            call mpas_pool_get_array(hydroPool, 'channelMelt', channelMelt)
+            call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
+            call mpas_pool_get_array(hydroPool, 'hydroTerrestrialMarginMask', hydroTerrestrialMarginMask)
+            call mpas_pool_get_array(hydroPool, 'waterFlux', waterFlux)
+            call mpas_pool_get_array(hydroPool, 'channelDischarge', channelDischarge)
+            call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
+         endif
 
          ! loop over cells
          do iCell = 1,nCellsSolve
@@ -314,11 +378,13 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
 
             blockSumGroundedIceArea = blockSumGroundedIceArea + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) &
                  * areaCell(iCell)
+
             blockSumGroundedIceVolume = blockSumGroundedIceVolume + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) &
                 * areaCell(iCell) * thickness(iCell)
 
             blockSumFloatingIceArea = blockSumFloatingIceArea + real(li_mask_is_floating_ice_int(cellMask(iCell)),RKIND) &
                 * areaCell(iCell)
+
             blockSumFloatingIceVolume = blockSumFloatingIceVolume + real(li_mask_is_floating_ice_int(cellMask(iCell)),RKIND) &
                 * areaCell(iCell) * thickness(iCell)
 
@@ -335,6 +401,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
             blockSumGroundedSfcMassBal = blockSumGroundedSfcMassBal + areaCell(iCell) * groundedSfcMassBalApplied(iCell) * scyr
             blockSumFloatingSfcMassBal = blockSumFloatingSfcMassBal + &
                (sfcMassBalApplied(iCell) - groundedSfcMassBalApplied(iCell)) * areaCell(iCell) * scyr
+
             ! BMB (kg yr-1)
             blockSumBasalMassBal = blockSumBasalMassBal + areaCell(iCell) * basalMassBalApplied(iCell) * scyr
             blockSumGroundedBasalMassBal = blockSumGroundedBasalMassBal + areaCell(iCell) * groundedBasalMassBalApplied(iCell) * scyr
@@ -362,15 +429,63 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
             blockGLMigrationFlux = blockGLMigrationFlux + groundedToFloatingThickness(iCell) * areaCell(iCell) &
                * rhoi / (deltat / scyr)  ! convert from m to kg/yr
 
+            !! Subglacial Hydrology Calculations
+            if (config_SGH) then
+
+               ! Subglacial Water Volume
+               blockSumSubglacialWaterVolume = blockSumSubglacialWaterVolume + waterThickness(iCell) * areaCell(iCell)
+
+               ! Basal melt input
+               blockSumBasalMeltInput = blockSumBasalMeltInput + real(li_mask_is_grounded_ice_int(cellMask(iCell)),RKIND) * &
+                                        basalMeltInput(iCell) * areaCell(iCell)
+
+               ! External water input
+               blockSumExternalWaterInput = blockSumExternalWaterInput + externalWaterInput(iCell) * areaCell(iCell)
+
+               ! Lake Volume
+               if (waterThickness(iCell) > bedBumpMax) then
+                  blockSumLakeVolume = blockSumLakeVolume + (waterThickness(iCell) - bedBumpMax) * areaCell(iCell)
+               endif
+
+               ! Lake Area
+               if (waterThickness(iCell) > bedBumpMax) then
+                   blockSumLakeArea = blockSumLakeArea + areaCell(iCell)
+               endif
+
+               ! Area-weighted flotation fraction for grounded ice
+               if (li_mask_is_grounded_ice(cellMask(iCell))) then
+                   blockSumFlotationFraction = blockSumFlotationFraction + ( waterPressure(iCell) / rhoi / gravity / thickness(iCell) ) * areaCell(iCell)
+               endif
+            endif
+
+
          end do ! end loop over cells
 
-         ! Loop over edges
-         do iEdge = 1, nEdgesSolve
-            ! Flux across GL, units = kg/yr
-            blockGLflux = blockGLflux + fluxAcrossGroundingLine(iEdge) * dvEdge(iEdge) &
-                      * scyr * rhoi ! convert from m^2/s to kg/yr
-         end do ! end loop over edges
+         if (config_SGH) then
+            ! Loop over edges
+            do iEdge = 1, nEdgesSolve
+
+               ! Flux across GL, units = kg/yr
+               blockGLflux = blockGLflux + fluxAcrossGroundingLine(iEdge) * dvEdge(iEdge) &
+                         * scyr * rhoi ! convert from m^2/s to kg/yr
 
+               ! Channel Melt
+               blockSumChannelMelt = blockSumChannelMelt + abs(channelMelt(iEdge) * dcEdge(iEdge))
+
+               ! Meltwater Flux across the grounding line
+               blockSumGLMeltFlux = blockSumGLMeltFlux + abs(hydroMarineMarginMask(iEdge) * waterFlux(iEdge) * dvEdge(iEdge) * rho_water)
+
+               ! Meltwater Flux across terrestrial margins
+               blockSumTerrestrialMeltFlux = blockSumTerrestrialMeltFlux + abs(hydroTerrestrialMarginMask(iEdge) * waterFlux(iEdge) * dvEdge(iEdge) * rho_water)
+
+               ! Meltwater Discharge in channels across grounding line
+               blockSumChannelGLMeltFlux = blockSumChannelGLMeltFlux + abs(hydroMarineMarginMask(iEdge) * channelDischarge(iEdge) * rho_water)
+
+               ! Meltwater discharge in channels across terrestrial margin
+               blockSumChannelTerrestrialMeltFlux = blockSumChannelTerrestrialMeltFlux + abs( hydroTerrestrialMarginMask(iEdge) * channelDischarge(iEdge) * rho_water)
+
+            end do ! end loop over edges
+         endif
          block => block % next
       end do    ! end loop over blocks
 
@@ -409,16 +524,30 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
       sums(7) = blockSumSfcMassBal
       sums(8) = blockSumGroundedSfcMassBal
       sums(9) = blockSumFloatingSfcMassBal
-
       sums(10) = blockSumBasalMassBal
       sums(11) = blockSumGroundedBasalMassBal
       sums(12) = blockSumFloatingBasalMassBal
       sums(13) = blockSumCalvingFlux
       sums(14) = blockSumFaceMeltingFlux
       sums(15) = blockSumVAF
-      sums(16) = blockGLflux
-      sums(17) = blockGLMigrationflux
-      nVars = 17
+      if (config_SGH) then
+         sums(16) = blockGLflux
+         sums(17) = blockGLMigrationflux
+         sums(18) = blockSumSubglacialWaterVolume
+         sums(19) = blockSumBasalMeltInput
+         sums(20) = blockSumExternalWaterInput
+         sums(21) = blockSumChannelMelt
+         sums(22) = blockSumLakeVolume
+         sums(23) = blockSumLakeArea
+         sums(24) = blockSumGLMeltFlux
+         sums(25) = blockSumTerrestrialMeltFlux
+         sums(26) = blockSumChannelGLMeltFlux
+         sums(27) = blockSumChannelTerrestrialMeltFlux
+         sums(28) = blockSumFlotationFraction
+         nVars = 28
+      else
+        nVars = 15
+      endif
 
       call mpas_dmpar_sum_real_array(dminfo, nVars, sums(1:nVars), reductions(1:nVars))
 
@@ -445,8 +574,21 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(globalStatsAMPool, 'avgSubshelfMelt', avgSubshelfMelt)
          call mpas_pool_get_array(globalStatsAMPool, 'totalCalvingFlux', totalCalvingFlux)
          call mpas_pool_get_array(globalStatsAMPool, 'totalFaceMeltingFlux', totalFaceMeltingFlux)
-         call mpas_pool_get_array(globalStatsAMPool, 'groundingLineFlux', groundingLineFlux)
-         call mpas_pool_get_array(globalStatsAMPool, 'groundingLineMigrationFlux', groundingLineMigrationFlux)
+         if (config_SGH) then
+            call mpas_pool_get_array(globalStatsAMPool, 'groundingLineFlux', groundingLineFlux)
+            call mpas_pool_get_array(globalStatsAMPool, 'groundingLineMigrationFlux', groundingLineMigrationFlux)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalSubglacialWaterVolume', totalSubglacialWaterVolume)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalBasalMeltInput', totalBasalMeltInput)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalExternalWaterInput', totalExternalWaterInput)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalChannelMelt', totalChannelMelt)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalSubglacialLakeVolume', totalSubglacialLakeVolume)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalSubglacialLakeArea', totalSubglacialLakeArea)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalDistWaterFluxMarineMargin',totalDistWaterFluxMarineMargin)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalDistWaterFluxTerrestrialMargin', totalDistWaterFluxTerrestrialMargin)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalChnlWaterFluxMarineMargin',totalChnlWaterFluxMarineMargin)
+            call mpas_pool_get_array(globalStatsAMPool, 'totalChnlWaterFluxTerrestrialMargin', totalChnlWaterFluxTerrestrialMargin)
+            call mpas_pool_get_array(globalStatsAMPool, 'avgFlotationFraction', avgFlotationFraction)
+         endif
 
          totalIceArea = reductions(1)
          totalIceVolume = reductions(2)
@@ -463,8 +605,21 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          totalCalvingFlux = reductions(13)
          totalFaceMeltingFlux = reductions(14)
          volumeAboveFloatation = reductions(15)
-         groundingLineFlux = reductions(16)
-         groundingLineMigrationFlux = reductions(17)
+         if (config_SGH) then
+            groundingLineFlux = reductions(16)
+            groundingLineMigrationFlux = reductions(17)
+            totalSubglacialWaterVolume = reductions(18)
+            totalBasalMeltInput = reductions(19)
+            totalExternalWaterInput = reductions(20)
+            totalChannelMelt = reductions(21)
+            totalSubglacialLakeVolume = reductions(22)
+            totalSubglacialLakeArea = reductions(23)
+            totalDistWaterFluxMarineMargin = reductions(24)
+            totalDistWaterFluxTerrestrialMargin = reductions(25)
+            totalChnlWaterFluxMarineMargin = reductions(26)
+            totalChnlWaterFluxTerrestrialMargin = reductions(27)
+            totalFlotationFraction = reductions(28)
+         endif
 
          if (totalIceArea > 0.0_RKIND) then
             iceThicknessMean = totalIceVolume / totalIceArea
@@ -473,16 +628,27 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
             iceThicknessMean = 0.0_RKIND
             avgNetAccumulation = 0.0_RKIND
          endif
+
          if (groundedIceArea > 0.0_RKIND) then
             avgGroundedBasalMelt = -1.0_RKIND * totalGroundedBasalMassBal / groundedIceArea / rhoi
          else
             avgGroundedBasalMelt = 0.0_RKIND
          endif
+
          if (floatingIceArea > 0.0_RKIND) then
             avgSubshelfMelt = -1.0_RKIND * totalFloatingBasalMassBal / floatingIceArea / rhoi
          else
             avgSubshelfMelt = 0.0_RKIND
          endif
+         if (config_SGH) then
+            if (groundedIceArea > 0.0_RKIND) then
+                avgFlotationFraction = totalFlotationFraction / groundedIceArea
+            else
+                avgFlotationFraction = 0.0_RKIND
+            endif
+         endif
+
+
 
          block => block % next
       end do