diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
index ac756b6d33f9..29a0784969d5 100644
--- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
+++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
@@ -20,7 +20,11 @@
 			description="The maximum allowable time step in seconds. If the allowable time step determined by the adaptive CFL calculation is longer than this, then the model will specify config_SGH_max_adaptive_timestep as the time step instead.  Defaults to 100 years (in seconds)."
 			possible_values="Any non-negative real value."
 		/>
-		<nml_option name="config_SGH_tangent_slope_calculation" type="character" default_value="from_normal_slope" units="unitless"
+                <nml_option name="config_SGH_coupling_interval" type="character" default_value="0000-00-00_00:00:01" units="unitless"
+                        description="Time interval at which the SGH model is called by MALI. If set to an interval less than the master timestep, config_SGH_coupling_interval will be superceded by the master timestep. If greater than the master timestep, config_adaptive_timestep_force_interval will be set equal to config_SGH_coupling_interval."
+                        possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS'"
+                />
+                <nml_option name="config_SGH_tangent_slope_calculation" type="character" default_value="from_normal_slope" units="unitless"
 		            description="Selection of the method for calculating the tangent component of slope at edges.
 'from_vertex_barycentric' interpolates scalar values from cell centers to vertices using the barycentric interpolation routine in operators (mpas_cells_to_points_using_baryweights) and then calculates the slope between vertices.  It works for obtuse triangles, but will not work correctly across the edges of periodic meshes.
 'from_vertex_barycentric_kiteareas' interpolates scalar values from cell centers to vertices using barycentric interpolation based on kiterea values and then calculates the slope between vertices.  It will work across the edges of periodic meshes, but will not work correctly for obtuse triangles.
@@ -180,8 +184,10 @@
                 <!-- convenience variables -->
                 <var name="effectivePressureSGH" type="real" dimensions="nCells Time" units="Pa"
                         description="effective ice pressure in subglacial hydrology system" />
-                <var name="effectivePressure" type="real" dimensions="nCells Time" units="Pa" packages="higherOrderVelocity"
-                        description="Effective pressure used in basal friction calculation.  If subglacial hydrology model is active, this will be effectivePressureSGH averaged over the subglacial hydrology model timestepping subcycles. If subglacial hydrology model is inactive, this will come from a file or a parameterization."/>
+                <var name="effectivePressure" type="real" dimensions="nCells Time" units="Pa" packages="higherOrderVelocity" default_value="-1e36"
+                        description="Effective pressure used in basal friction calculation.  If subglacial hydrology model is active, this will be effectivePressureSGH averaged over the interval defined by config_SGH_coupling_interval. If subglacial hydrology model is inactive, this will come from a file or a parameterization."/>
+                <var name="effectivePressureTimeCycleAvg" type="real" dimensions="nCells Time" units="Pa"
+                        description="effectivePressureSGH averaged over the subglacial hydrology model timestepping subcycle." /> 
                 <var name="hydropotential" type="real" dimensions="nCells Time" units="Pa"
                      description="hydropotential in subglacial hydrology system" />
                 <var name="waterFlux" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
@@ -246,7 +252,9 @@
                      description="time step length limited by pressure equation scheme in subglacial hydrology system" />
                 <var name="deltatSGH" type="real" dimensions="Time" units="s"
                      description="time step used for evolving subglacial hydrology system" />
-	     	<var name="distGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
+                <var name="deltatSGHaccumulated" type="real" dimensions="Time" units="s"
+                     description="accumulated deltatSGH throughout SGH coupling interval" />
+                <var name="distGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
 		     description="distributed discharge across the grounding line, summed from grounding line edges to adjacent ungrounded cell. Values from all edges are summed if multiple grounding line edges border a single ungrounded cell" />
 	     	<var name="totalGroundingLineDischargeCell" type="real" dimensions="Time nCells" units="m^{3} s^{-1}"
 		     description="total (channel + dist.) discharge across the grounding line, summed from grounding line edges to adjacent ungrounded cell. Values from all edges are summed if multiple grounding line edges border a single ungrounded cell" />
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 06fb2239107d..ea9508e721f2 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
@@ -111,15 +111,20 @@ subroutine li_SGH_init(domain, err)
       real (kind=RKIND), dimension(:), pointer :: bedTopography
       real (kind=RKIND), dimension(:), pointer :: iceThicknessHydro
       real (kind=RKIND), dimension(:), pointer :: hydropotential
+      real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
+      real (kind=RKIND),dimension(:), pointer :: effectivePressure
       integer, dimension(:), pointer :: cellMask
       real (kind=RKIND), pointer :: tillMax
       real (kind=RKIND), pointer :: rhoi, rhoo
       logical, pointer :: config_do_restart
       real (kind=RKIND), pointer :: config_sea_level
       integer, pointer :: config_num_halos
+      character (len=StrKIND), pointer :: config_SGH_coupling_interval
+      type (MPAS_TimeInterval_type) :: sgh_coupling_interval
+      character (len=StrKIND), pointer :: simulationStartTime
+      type (MPAS_Time_type) :: simulationStartTime_timeType
       integer :: err_tmp
 
-
       err = 0
       err_tmp = 0
 
@@ -132,8 +137,31 @@ subroutine li_SGH_init(domain, err)
          ! If SGH is not active, skip everything
          return
       endif
-
-
+      
+      ! Check SGH coupling interval is compatible with master timestep and restart interval
+      call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
+      call check_SGH_coupling_interval(meshPool, err_tmp)
+      err = ior(err, err_tmp)
+      if (err /= 0) then
+         call mpas_log_write("Error occurred in check_SGH_coupling_interval. ", MPAS_LOG_ERR)
+      endif
+     
+      ! Initiate coupling interval timer 
+      call mpas_pool_get_config(liConfigs, 'config_SGH_coupling_interval', config_SGH_coupling_interval)
+      call mpas_set_timeInterval(sgh_coupling_interval, timeString=config_SGH_coupling_interval, ierr=err_tmp)
+      err = ior(err, err_tmp)
+      call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime)
+      call mpas_set_time(simulationStartTime_timeType, dateTimeString=simulationStartTime, ierr=err_tmp)
+      err = ior(err, err_tmp)
+      call mpas_add_clock_alarm(domain%clock,'sghCouplingInterval', alarmTime=simulationStartTime_timeType, alarmTimeInterval=sgh_coupling_interval, ierr=err_tmp)
+      err = ior(err, err_tmp)
+      if (mpas_is_alarm_ringing(domain%clock, 'sghCouplingInterval', ierr=err_tmp)) then
+         err = ior(err, err_tmp)
+         call mpas_reset_clock_alarm(domain%clock, 'sghCouplingInterval', ierr=err_tmp)
+         err = ior(err, err_tmp)
+      endif
+     
+      ! begin hydro init
       call mpas_timer_start("hydro init")
 
       call mpas_log_write('Beginning subglacial hydro init.')
@@ -239,7 +267,22 @@ subroutine li_SGH_init(domain, err)
       call mpas_timer_start("halo updates")
       call mpas_dmpar_field_halo_exch(domain, 'waterPressureSmooth')
       call mpas_timer_stop("halo updates")
-     
+    
+      ! initialize effectivePressure
+      block => domain % blocklist
+      do while (associated(block))
+        
+         call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
+         call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
+         call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
+         
+         if (all(effectivePressure == -1e36_RKIND)) then
+            effectivePressure = effectivePressureSGH
+         endif
+         
+         block => block % next
+      end do
+
       ! === error check
       if (err > 0) then
           call mpas_log_write("An error has occurred in li_SGH_init.", MPAS_LOG_ERR)
@@ -316,8 +359,9 @@ subroutine li_SGH_solve(domain, err)
       real (kind=RKIND), dimension(:), pointer :: waterVelocity
       real (kind=RKIND), dimension(:), pointer :: waterVelocityCellX
       real (kind=RKIND), dimension(:), pointer :: waterVelocityCellY
-      real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
       real (kind=RKIND), dimension(:), pointer :: effectivePressure
+      real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
+      real (kind=RKIND), dimension(:), pointer :: effectivePressureTimeCycleAvg
       real (kind=RKIND), dimension(:), allocatable :: cellJunk
       integer, dimension(:), pointer :: nEdgesOnCell
       integer, dimension(:,:), pointer :: edgesOnCell
@@ -325,6 +369,7 @@ subroutine li_SGH_solve(domain, err)
       integer, dimension(:,:), pointer :: edgeSignOnCell
       integer, dimension(:), pointer :: cellMask
       real (kind=RKIND), pointer :: deltatSGH
+      real (kind=RKIND), pointer :: deltatSGHaccumulated
       real (kind=RKIND), pointer :: masterDeltat
       real (kind=RKIND), pointer :: Cd
       real (kind=RKIND), pointer :: tillMax
@@ -332,7 +377,6 @@ subroutine li_SGH_solve(domain, err)
       integer, pointer :: nCells
       integer :: iCell, iEdge, iEdgeOnCell
       real (kind=RKIND) :: timeLeft ! subcycling time remaining until master dt is reached
-      real (kind=RKIND) :: deltatSGHaccumulated
       integer :: numSubCycles ! number of subcycles
       integer :: err_tmp
 
@@ -380,25 +424,27 @@ subroutine li_SGH_solve(domain, err)
       call mpas_timer_start("halo updates")
       call mpas_dmpar_field_halo_exch(domain, 'iceThicknessHydro')
       call mpas_timer_stop("halo updates")
-      
+     
       ! initialize while loop
       call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)  ! can get from any block
       call mpas_pool_get_array(meshPool, 'deltat', masterDeltat)
       timeLeft = masterDeltat  ! in seconds
       numSubCycles = 0
 
+      ! initialize accumulated time averages/sums throughout timecycle
       block => domain % blocklist
       do while (associated(block))
         
          call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
-         call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
-         effectivePressure = 0.0_RKIND
+         call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)
+         call mpas_pool_get_array(hydroPool, 'deltatSGHaccumulated', deltatSGHaccumulated)
+         
+         effectivePressureTimeCycleAvg = 0.0_RKIND
          deltatSGHaccumulated = 0.0_RKIND
-
+         
          block => block % next
-      end do
-
-
+      enddo
+      
       ! =============
       ! =============
       ! =============
@@ -704,10 +750,11 @@ subroutine li_SGH_solve(domain, err)
         
          call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
          call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
-         call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
+         call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)
          call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
+         call mpas_pool_get_array(hydroPool, 'deltatSGHaccumulated', deltatSGHaccumulated)
 
-         effectivePressure = (effectivePressure * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH)
+         effectivePressureTimeCycleAvg = (effectivePressureTimeCycleAvg * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH)
    
          block => block % next
       end do
@@ -722,8 +769,24 @@ subroutine li_SGH_solve(domain, err)
       ! =============
       ! =============
 
-      call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))
+      ! set up coupling interval
+      if (mpas_is_alarm_ringing(domain % clock, 'sghCouplingInterval', ierr=err_tmp)) then
+         err = ior(err, err_tmp)
+         block => domain % blocklist
+         do while (associated(block))
 
+            call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
+            call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
+            call mpas_pool_get_array(hydroPool, 'effectivePressureTimeCycleAvg', effectivePressureTimeCycleAvg)
+            
+            effectivePressure = effectivePressureTimeCycleAvg
+            block => block % next
+         end do
+         call mpas_reset_clock_alarm(domain % clock, 'sghCouplingInterval', ierr=err_tmp)
+         err = ior(err, err_tmp)
+      endif
+
+      call mpas_log_write("Hydro model subcycled $i times.", intArgs=(/numSubCycles/))
 
       ! === error check
       if (err > 0) then
@@ -2738,4 +2801,90 @@ subroutine calc_waterPressureSmooth(block,err)
       deallocate(pressureField)
    end subroutine calc_waterPressureSmooth
 
+!***********************************************************************
+!
+!  routine check_SGH_coupling_interval
+!
+!> \brief   Perform various checks on the SGH coupling interval setting
+!> \author  Matt Hoffman, modified by Alex Hager
+!> \date    Sep 2024
+!> \details
+!>  This routine checks that the SGH coupling interval is an even multiple
+!>  of the adaptive timestep force inverval and divides evenly into the
+!>  restart interval. 
+!
+!-----------------------------------------------------------------------
+    subroutine check_SGH_coupling_interval(meshPool, err)
+
+      use mpas_timekeeping
+      use mpas_stream_manager
+      use mpas_derived_types, only : MPAS_STREAM_PROPERTY_RECORD_INTV
+
+      type (mpas_pool_type), intent(in) :: meshPool  !< mesh information
+      integer, intent(out) :: err
+
+      ! local variables
+      character (len=StrKIND), pointer :: config_SGH_coupling_interval
+      logical, pointer :: config_adaptive_timestep
+      character (len=StrKIND), pointer :: config_adaptive_timestep_force_interval, config_dt
+      type (MPAS_TimeInterval_Type) :: coupling_interval, force_interval, dt_interval, zero_interval
+      type (MPAS_Time_Type) :: start_time
+      character (len=StrKIND), pointer :: simulationStartTime
+      integer (kind=I8KIND) :: n_intervals
+      type (MPAS_TimeInterval_Type) :: remainder
+      integer :: err_tmp
+
+      err = 0
+      err_tmp = 0
+
+      call mpas_log_write("")
+      call mpas_log_write("-- Checking consistency of config_SGH_coupling_interval and other settings --")
+
+      ! define zero interval for comparing against below
+      call mpas_set_timeInterval(zero_interval, dt = 0.0_RKIND)
+      ! get start time as a reference time
+      call mpas_pool_get_array(meshPool, 'simulationStartTime', simulationStartTime)
+      call mpas_set_time(start_time, dateTimeString=simulationStartTime)
+      ! define SGH coupling interval as a timeInterval type
+      call mpas_pool_get_config(liConfigs, 'config_SGH_coupling_interval', config_SGH_coupling_interval)
+      call mpas_set_timeInterval(coupling_interval, timeString=config_SGH_coupling_interval, ierr=err_tmp)
+      err = ior(err, err_tmp)
+
+      call mpas_pool_get_config(liConfigs, "config_adaptive_timestep", config_adaptive_timestep)
+      if (config_adaptive_timestep) then
+         ! for adaptive dt, check that config_adaptive_timestep_force_interval divides evenly into config_SGH_coupling_interval
+         call mpas_pool_get_config(liConfigs, "config_adaptive_timestep_force_interval", config_adaptive_timestep_force_interval)
+         call mpas_set_timeInterval(force_interval, timeString=config_adaptive_timestep_force_interval, ierr=err_tmp)
+         err = ior(err, err_tmp)
+         call mpas_interval_division(start_time, coupling_interval, force_interval, n_intervals, remainder)
+         if (remainder .EQ. zero_interval) then
+            call mpas_log_write("config_adaptive_timestep_force_interval divides into config_SGH_coupling_interval $i times " // &
+                 "with no remainder - check passes", intArgs=(/int(n_intervals)/))
+         else
+            call mpas_log_write("config_adaptive_timestep_force_interval divides into config_SGH_coupling_interval $i times " // &
+                 "with nonzero remainder", MPAS_LOG_ERR, intArgs=(/int(n_intervals)/))
+            err = ior(err, 1)
+         endif
+      else
+        ! For fixed dt, check that dt divides evenly into config_SGH_coupling_interval
+         call mpas_pool_get_config(liConfigs, "config_dt", config_dt)
+         call mpas_set_timeInterval(dt_interval, timeString=config_dt, ierr=err_tmp)
+         err = ior(err, err_tmp)
+         call mpas_interval_division(start_time, coupling_interval, dt_interval, n_intervals, remainder)
+         if (remainder .EQ. zero_interval) then
+            call mpas_log_write("config_dt divides into config_SGH_coupling_interval $i times " // &
+                 "with no remainder - check passes", intArgs=(/int(n_intervals)/))
+         else
+            call mpas_log_write("config_dt divides into config_SGH_coupling_interval $i times " // &
+                 "with nonzero remainder", MPAS_LOG_ERR, intArgs=(/int(n_intervals)/))
+            err = ior(err, 1)
+         endif
+      endif
+
+      call mpas_log_write("")
+
+   !--------------------------------------------------------------------
+   end subroutine check_SGH_coupling_interval
+
+
 end module li_subglacial_hydro