diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index a6cf3e947a81..3ab91cc321a6 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -268,7 +268,11 @@
                         description="If true, then the value of 'stiffnessFactor' is coupled to damage evolution, i.e. 'stiffnessFactor' = 1-damage."
 	                possible_values=".true. or .false."
 		/>
-		        <nml_option name="config_damage_gl_setting" type="character" default_value="nye" units="unitless"
+                <nml_option name="config_initialize_damage_with_stiffnessFactor" type="logical" default_value=".false." units="unitless"
+                        description="If true, then initialize damage to max((1 - stiffnessFactor), 0)"
+                        possible_values=".true. or .false."
+                />    
+                <nml_option name="config_damage_gl_setting" type="character" default_value="nye" units="unitless"
 		                description="Selection of the method for initializing damage in the first floating grid cells past the grounding line: 'nye' uses Nye's zero stress criteria (Nye, 1957, Proc. R. Soc. A, 239) for defining the damage value, 'extrapolate' extrapolates values from the floating ice downstream back to the first floating cells."
 		            possible_values="'nye', 'extrapolate'"
 		/>
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
index d6fee67ee637..f63bf787b0c5 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F
@@ -3535,6 +3535,11 @@ subroutine li_calculate_damage(domain, err)
             err = ior(err, 1)
         endif
 
+        call mpas_timer_start("halo updates")
+        call mpas_dmpar_field_halo_exch(domain, 'damage')
+        call mpas_dmpar_field_halo_exch(domain, 'damageNye')
+        call mpas_timer_stop("halo updates")
+
         if (config_print_calving_info) then
            ! End of routine accounting/reporting
            localMinInfo(1) = minval(damageSource)
@@ -3616,11 +3621,13 @@ subroutine li_finalize_damage_after_advection(domain, err)
       type (mpas_pool_type), pointer :: velocityPool
       type (mpas_pool_type), pointer :: scratchPool
       real(kind=RKIND), pointer :: config_damage_stiffness_min
+      real(kind=RKIND), pointer :: deltat
       logical, pointer :: config_damage_rheology_coupling
       logical, pointer :: config_preserve_damage
       logical, pointer :: config_print_calving_info
       character (len=StrKIND), pointer :: config_damage_gl_setting
       real (kind=RKIND), dimension(:), pointer :: damage
+      real (kind=RKIND), dimension(:), pointer :: ddamagedt
       real (kind=RKIND), dimension(:), pointer :: damageMax
       real (kind=RKIND), dimension(:), pointer :: damageNye
       real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
@@ -3653,11 +3660,12 @@ subroutine li_finalize_damage_after_advection(domain, err)
 
          ! get fields
          call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
-
+         call mpas_pool_get_array(meshPool, 'deltat', deltat)
          call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
          call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
          call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
          call mpas_pool_get_array(geometryPool, 'damage', damage)
+         call mpas_pool_get_array(geometryPool, 'ddamagedt', ddamagedt)
          call mpas_pool_get_array(geometryPool, 'damageMax', damageMax)
          call mpas_pool_get_array(geometryPool, 'damageNye', damageNye)
          call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
@@ -3739,10 +3747,18 @@ subroutine li_finalize_damage_after_advection(domain, err)
             damage = 1.0_RKIND
         end where
 
+         ! Halo update for damage before applying to stiffnessFactor
+         ! NOTE: THIS WILL NOT WORK ON MULTIPLE BLOCKS PER PROCESSOR
+         call mpas_timer_start("halo updates")
+         call mpas_dmpar_field_halo_exch(domain, 'damage')
+         call mpas_dmpar_field_halo_exch(domain, 'damageNye')
+         call mpas_timer_stop("halo updates")
+
         if (config_damage_rheology_coupling) then
             do iCell = 1, nCells
                 if (li_mask_is_floating_ice(cellMask(iCell))) then
-                    stiffnessFactor(iCell) = 1.0_RKIND - damage(iCell)
+                    stiffnessFactor(iCell) = stiffnessFactor(iCell) - &
+                                             (ddamagedt(iCell) * deltat)
                     if (stiffnessFactor(iCell) < config_damage_stiffness_min) then
                         stiffnessFactor(iCell) = config_damage_stiffness_min
                     end if
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F
index 10ee1e94b305..ca9069419764 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F
@@ -749,6 +749,7 @@ function li_core_initial_solve(domain) result(err)
       type (mpas_pool_type), pointer :: geometryPool, meshPool, velocityPool
       logical, pointer :: config_do_restart, config_write_output_on_startup, config_write_stats_on_startup
       character(len=StrKIND), pointer :: config_velocity_solver
+      logical, pointer :: config_initialize_damage_with_stiffnessFactor
       ! Variables needed for printing timestamps
       type (MPAS_Time_Type) :: currTime
       character(len=StrKIND) :: timeStamp
@@ -759,6 +760,8 @@ function li_core_initial_solve(domain) result(err)
       logical :: solveVelo
 
       integer, dimension(:), pointer :: vertexMask
+      real (kind=RKIND), dimension(:), pointer :: damage
+      real (kind=RKIND), dimension(:), pointer :: stiffnessFactor
 
 
       err = 0
@@ -770,6 +773,10 @@ function li_core_initial_solve(domain) result(err)
       call mpas_pool_get_config(liConfigs, 'config_write_output_on_startup', config_write_output_on_startup)
       call mpas_pool_get_config(liConfigs, 'config_write_stats_on_startup', config_write_stats_on_startup)
       call mpas_pool_get_config(liConfigs, 'config_velocity_solver', config_velocity_solver)
+      call mpas_pool_get_config(liConfigs, 'config_initialize_damage_with_stiffnessFactor', config_initialize_damage_with_stiffnessFactor)
+      call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
+      call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
+      call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool)
 
       currTime = mpas_get_clock_time(domain % clock, MPAS_NOW, err_tmp)
       err = ior(err, err_tmp)
@@ -777,7 +784,17 @@ function li_core_initial_solve(domain) result(err)
       err = ior(err, err_tmp)
       call mpas_log_write('Initial timestep ' // trim(timeStamp))
 
+      ! initialize damage using stiffnessFactor
+      call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
+      call mpas_pool_get_array(geometryPool, 'damage', damage)
+      if ( config_initialize_damage_with_stiffnessFactor .and. (.not. config_do_restart) ) then
+         damage(:) = 1.0_RKIND - stiffnessFactor(:)
 
+         where ( damage < 0.0_RKIND )
+            damage = 0.0_RKIND
+         end where
+
+      end if
       ! ===
       ! === Calculate Initial state
       ! ===