diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml
index ec1c2a015a..625a3e5ce5 100644
--- a/src/core_landice/Registry.xml
+++ b/src/core_landice/Registry.xml
@@ -941,7 +941,7 @@ is the value of that variable from the *previous* time level!
                 <var name="dynamicThickening" type="real" dimensions="nCells Time" units="m a^{-1}"
                      description="diagnostic field of dynamic thickening rate (calculated as negative of flux divergence)"
                 />
-                <var name="cellMask" type="integer" dimensions="nCells Time" units="none"
+                <var name="cellMask" type="integer" dimensions="nCells Time" units="none" time_levs="1"
                      description="bitmask indicating various properties about the ice sheet on cells.  cellMask only needs to be a restart field if config_allow_additional_advance = false (to keep the mask of initial ice extent)"
                 />
                 <var name="edgeMask" type="integer" dimensions="nEdges Time" units="none"
@@ -1137,6 +1137,10 @@ is the value of that variable from the *previous* time level!
                      units="none"  description="flag needed by external velocity solvers that indicates if the region to solve on the block's domain has changed (treated as a logical)"
                         packages="higherOrderVelocity"
                 />
+                <var name="cellIceMaskChanged"              type="integer"  dimensions="Time"
+                     units="none"  description="flag needed by external velocity solvers that indicates if the region to solve on the block's domain has changed (treated as a logical)"
+                        packages="higherOrderVelocity"
+                />
                 <var name="dirichletVelocityMask"            type="integer"     dimensions="nVertInterfaces nCells Time"
                         units="none"  time_levs="2"
                         description="mask of where Dirichlet boundary conditions should be applied to the velocity solution.  1 means apply a Dirichlet boundary condition, 0 means do not. (higher-order dycores only)"
@@ -1149,9 +1153,14 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingEdges"
                         type="integer" dimensions="nEdges Time"
                         units="unitless"
+                        time_levs="2"
                         description="edges which are floating have a value of 1.  non floating edges have a value of 0."
                         packages="higherOrderVelocity"
                 />
+                <var name="floatingEdgesChanged"              type="integer"  dimensions="Time"
+                     units="none"  description="flag needed by external velocity solvers that indicates if the floating edges ids have changed (treated as a logical)"
+                        packages="higherOrderVelocity"
+                />
         </var_struct>
 
 <!-- ================ -->
diff --git a/src/core_landice/analysis_members/mpas_li_global_stats.F b/src/core_landice/analysis_members/mpas_li_global_stats.F
index 41ddc9aebd..b367875979 100644
--- a/src/core_landice/analysis_members/mpas_li_global_stats.F
+++ b/src/core_landice/analysis_members/mpas_li_global_stats.F
@@ -283,7 +283,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
          call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal)
diff --git a/src/core_landice/analysis_members/mpas_li_regional_stats.F b/src/core_landice/analysis_members/mpas_li_regional_stats.F
index 0e0543df57..4c6df7def9 100644
--- a/src/core_landice/analysis_members/mpas_li_regional_stats.F
+++ b/src/core_landice/analysis_members/mpas_li_regional_stats.F
@@ -275,7 +275,7 @@ subroutine li_compute_regional_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
          call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
          call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied)
diff --git a/src/core_landice/mode_forward/mpas_li_advection.F b/src/core_landice/mode_forward/mpas_li_advection.F
index 7dad2e6cb3..fbd80d0afd 100644
--- a/src/core_landice/mode_forward/mpas_li_advection.F
+++ b/src/core_landice/mode_forward/mpas_li_advection.F
@@ -257,7 +257,7 @@ subroutine li_advection_thickness_tracers(&
       call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
       call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening)
       call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)
@@ -630,7 +630,7 @@ subroutine li_update_geometry(geometryPool)
 
      call mpas_pool_get_dimension(geometryPool, 'nCells', nCells)
 
-     call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+     call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
      call mpas_pool_get_array(geometryPool, 'thickness', thickness)
      call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
      call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
diff --git a/src/core_landice/mode_forward/mpas_li_calving.F b/src/core_landice/mode_forward/mpas_li_calving.F
index 91146b7c05..6dff0a04dc 100644
--- a/src/core_landice/mode_forward/mpas_li_calving.F
+++ b/src/core_landice/mode_forward/mpas_li_calving.F
@@ -402,7 +402,7 @@ subroutine li_restore_calving_front(domain, err)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'restoreThickness', restoreThickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          ! get required fields from the thermal pool
          call mpas_pool_get_array(thermalPool, 'temperature', temperature)
@@ -870,7 +870,7 @@ subroutine floating_calving(domain, calvingFraction, err)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          calvingThickness = 0.0_RKIND
 
@@ -929,7 +929,7 @@ subroutine remove_small_islands(meshPool, geometryPool)
       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       do iCell = 1, nCellsSolve
@@ -1039,7 +1039,7 @@ subroutine topographic_calving(domain, calvingFraction, err)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
 
          calvingThickness = 0.0_RKIND
@@ -1152,7 +1152,7 @@ subroutine eigencalving(domain, err)
          call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
          call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'eigencalvingParameter', eigencalvingParameter)
@@ -1336,7 +1336,7 @@ subroutine distribute_calving_flux(meshPool, geometryPool, scratchPool, calvingF
       ! get fields
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       call mpas_pool_get_array(meshPool, 'deltat', deltat)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'requiredCalvingVolumeRate', requiredCalvingVolumeRate)
@@ -1566,7 +1566,7 @@ subroutine calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       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, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       calvingFrontMask = 0 !initialize
@@ -1656,7 +1656,7 @@ subroutine remove_icebergs(domain)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
          call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
          call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve)
@@ -1710,7 +1710,7 @@ subroutine remove_icebergs(domain)
             call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
             call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
             call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
-            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
             call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
             call mpas_pool_get_array(geometryPool, 'contiguousDynamicIceOld', contiguousDynamicIceOld)
             call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
@@ -1765,7 +1765,7 @@ subroutine remove_icebergs(domain)
       block => domain % blocklist
       do while (associated(block))
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
diff --git a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
index ab47f0d226..a8f2664137 100644
--- a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
+++ b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
@@ -258,10 +258,9 @@ subroutine li_calculate_apparent_diffusivity(meshPool, velocityPool, scratchPool
       call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel=1)
       call mpas_pool_get_array(geometryPool, 'apparentDiffusivity', apparentDiffusivity)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'uReconstructZonal', uReconstructAxis1)
       call mpas_pool_get_array(velocityPool, 'uReconstructMeridional', uReconstructAxis2)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
       call mpas_pool_get_field(scratchPool, 'workCell', slopeReconstructXField)
       call mpas_allocate_scratch_field(slopeReconstructXField, .true.)
       call mpas_pool_get_field(scratchPool, 'workCell2', slopeReconstructYField)
diff --git a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
index 95f6716fa8..0b99e9d843 100644
--- a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
+++ b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
@@ -247,7 +247,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
 
          ! get fields from the scratch pool
@@ -853,7 +853,7 @@ subroutine calc_iceshelf_draft_info(domain, GLdepth, CFdepth, err)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
          call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
          call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
diff --git a/src/core_landice/mode_forward/mpas_li_sia.F b/src/core_landice/mode_forward/mpas_li_sia.F
index b22001c4f2..d400c12e35 100644
--- a/src/core_landice/mode_forward/mpas_li_sia.F
+++ b/src/core_landice/mode_forward/mpas_li_sia.F
@@ -316,7 +316,7 @@ subroutine li_sia_solve(meshPool, geometryPool, thermalPool, velocityPool, err)
       call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'slopeEdge', slopeEdge)
       call mpas_pool_get_array(geometryPool, 'normalSlopeEdge', normalSlopeEdge)
       call mpas_pool_get_array(thermalPool, 'temperature', temperature)
diff --git a/src/core_landice/mode_forward/mpas_li_statistics.F b/src/core_landice/mode_forward/mpas_li_statistics.F
index 47a7be44e8..99de7c4939 100644
--- a/src/core_landice/mode_forward/mpas_li_statistics.F
+++ b/src/core_landice/mode_forward/mpas_li_statistics.F
@@ -235,7 +235,7 @@ subroutine li_compute_statistics(domain, itimestep)
         ! Geometry arrays
         call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
         call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
-        call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+        call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
         call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
         call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel = 1)
         call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel = 1)
diff --git a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
index 9e2e4f79ca..596119d01c 100644
--- a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
+++ b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
@@ -163,7 +163,7 @@ subroutine li_SGH_init(domain, err)
          waterPressure = max(0.0_RKIND, waterPressure)
          waterPressure = min(waterPressure, rhoi * gravity * thickness)
          ! set pressure correctly under floating ice and open ocean
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          where ( (li_mask_is_floating_ice(cellMask)) .or. &
                  ( (.not. li_mask_is_ice(cellMask)) .and. (bedTopography < config_sea_level) ) )
@@ -379,7 +379,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
          call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
          call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          WtillOld = Wtill
          Wtill = Wtill + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - Cd)
@@ -420,7 +420,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity)
          call mpas_pool_get_array(hydroPool, 'waterVelocityCellX', waterVelocityCellX)
          call mpas_pool_get_array(hydroPool, 'waterVelocityCellY', waterVelocityCellY)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          allocate(cellJunk(nCells+1))
          call mpas_reconstruct(meshPool, waterVelocity, waterVelocityCellX, waterVelocityCellY, &
               cellJunk, cellJunk, cellJunk)
@@ -476,7 +476,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          divergence(:) = 0.0_RKIND  ! zero div before starting
          ! loop over locally owned cells
@@ -551,7 +551,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'divergence', divergence)
          call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
          call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          waterThicknessOld = waterThickness
          waterThickness = waterThicknessOld + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - divergence  &
@@ -755,7 +755,7 @@ subroutine calc_edge_quantities(block, err)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal)
       call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal)
       call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
       call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeTangent', hydropotentialBaseSlopeTangent)
@@ -1217,7 +1217,7 @@ subroutine calc_pressure(block, err)
       call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed)
       call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:))
@@ -1338,7 +1338,7 @@ subroutine calc_pressure_diag_vars(block, err)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase)
       call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
       call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
       effectivePressure = rhoi * gravity * thickness - waterPressure
          ! < this should evalute to 0 for floating ice if Pw set correctly there.
diff --git a/src/core_landice/mode_forward/mpas_li_thermal.F b/src/core_landice/mode_forward/mpas_li_thermal.F
index 5bceb3c9d2..296d5070c6 100644
--- a/src/core_landice/mode_forward/mpas_li_thermal.F
+++ b/src/core_landice/mode_forward/mpas_li_thermal.F
@@ -820,7 +820,7 @@ subroutine li_thermal_solver(domain, err)
          ! get fields from the geometry pool
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
          call mpas_pool_get_array(geometryPool, 'basalWaterThickness', basalWaterThickness)
diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F
index e1e543e929..3bd6541abe 100644
--- a/src/core_landice/mode_forward/mpas_li_velocity.F
+++ b/src/core_landice/mode_forward/mpas_li_velocity.F
@@ -269,8 +269,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       integer, pointer :: nEdges
       integer, pointer :: nVertInterfaces
       integer, pointer :: nCells
-      integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask, vertexMaskOld
+      integer, dimension(:), pointer :: edgeMask, cellMask, cellMaskOld, vertexMask, vertexMaskOld
       integer, dimension(:,:), pointer :: dirichletVelocityMaskOld, dirichletVelocityMaskNew
+      integer, dimension(:), pointer :: floatingEdgesOld, floatingEdgesNew
       real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, normalVelocityInitial
       real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, &
          uReconstructZonal, uReconstructMeridional
@@ -288,7 +289,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       real (kind=RKIND), dimension(:), pointer :: upperSurface
       integer, dimension(:), pointer :: indexToEdgeID
       integer, pointer :: anyDynamicVertexMaskChanged
+      integer, pointer :: cellIceMaskChanged
       integer, pointer :: dirichletMaskChanged
+      integer, pointer :: floatingEdgesChanged
       ! truly local variables
       integer :: cell1, cell2
       integer :: cell3, cell4, thisCell
@@ -301,7 +304,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       real (kind=RKIND) :: maxThicknessOnProc, maxThicknessAllProcs
       real (kind=RKIND) :: xVelEdge, yVelEdge
       integer :: blockDynamicVertexMaskChanged, procDynamicVertexMaskChanged
+      integer :: blockCellIceMaskChanged, procCellIceMaskChanged
       integer :: blockDirichletMaskChanged, procDirichletMaskChanged
+      integer :: blockFloatingEdgesChanged, procFloatingEdgesChanged
 
       call mpas_timer_start("velocity solve")
 
@@ -386,7 +391,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          ! Don't bother checking this with SIA because it requires an extra global reduce
          maxThicknessOnProc = 0.0_RKIND  ! initialize to
          procDynamicVertexMaskChanged = 0
+         procCellIceMaskChanged = 0
          procDirichletMaskChanged = 0
+         procFloatingEdgesChanged = 0
 
          block => domain % blocklist
          do while (associated(block))
@@ -398,7 +405,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             maxThicknessOnProc = max(maxThicknessOnProc, maxval(thickness))
 
             ! The interface expects an array where 1's are floating edges and 0's are non-floating edges.
-            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
             call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
             floatingEdges = li_mask_is_floating_ice_int(edgeMask)
             call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1)
@@ -408,7 +415,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! TODO:  There may be some aspects of the mask that are ok change for external dycores,
             !        but for now just check the whole thing.
             call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMaskOld, timeLevel=2)
-            if ( sum(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld)) /= 0 ) then
+            if ( sum(abs(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld))) /= 0 ) then
                 blockDynamicVertexMaskChanged = 1
             else
                 blockDynamicVertexMaskChanged = 0
@@ -417,11 +424,26 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! Determine if any blocks on this processor had a change to the vertex mask
             procDynamicVertexMaskChanged = max(procDynamicVertexMaskChanged, blockDynamicVertexMaskChanged)
             !print *,'procVertexMaskChanged', procVertexMaskChanged
+            
+            ! Determine if the cell mask changed during this time step for this block (needed for external dycores)
+            ! TODO:  There may be some aspects of the mask that are ok change for external dycores,
+            !        but for now just check the whole thing.
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMaskOld, timeLevel=2)
+            if (( sum(abs(li_mask_is_dynamic_ice_int(cellMask) - li_mask_is_dynamic_ice_int(cellMaskOld))) /= 0 ) .or. &
+                ( sum(abs(li_mask_is_ice_int(cellMask) - li_mask_is_ice_int(cellMaskOld))) /= 0 )) then
+                blockCellIceMaskChanged = 1
+            else
+                blockCellIceMaskChanged = 0
+            endif
+      
+            ! Determine if any blocks on this processor had a change to the cell mask
+            procCellIceMaskChanged = max(procCellIceMaskChanged, blockCellIceMaskChanged)
 
-            ! Also check to see if the Dirichlet b.c. mask has changed
+            ! Check to see if the Dirichlet b.c. mask has changed
             call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskOld, timeLevel=2)
             call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskNew, timeLevel=1)
-            if ( sum(dirichletVelocityMaskNew - dirichletVelocityMaskOld) /= 0 ) then
+            if ( sum(abs(dirichletVelocityMaskNew - dirichletVelocityMaskOld)) /= 0 ) then
                 blockDirichletMaskChanged = 1
             else
                 blockDirichletMaskChanged = 0
@@ -429,6 +451,17 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! Determine if any blocks on this processor had a change to the vertex mask
             procDirichletMaskChanged = max(procDirichletMaskChanged, blockDirichletMaskChanged)
 
+            ! Check to see if the Floating Edges Ids have changed
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesOld, timeLevel=2)
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesNew, timeLevel=1)
+            if ( sum(abs(floatingEdgesNew - floatingEdgesOld)) /= 0 ) then
+                blockFloatingEdgesChanged = 1
+            else
+                blockFloatingEdgesChanged = 0
+            endif
+            ! Determine if any blocks on this processor had a change to the vertex mask
+            procFloatingEdgesChanged = max(procFloatingEdgesChanged, blockFloatingEdgesChanged)
+
             ! Set beta for solver to use
             ! (this could potentially be applied to an SIA solver, so it is calculated in this module)
             call calculate_beta(block, err_tmp)
@@ -452,10 +485,16 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged)
          call mpas_dmpar_max_int(domain % dminfo, procDynamicVertexMaskChanged, anyDynamicVertexMaskChanged)
          !print *,'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged
+         ! Do the same for the cell mask
+         call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged)
+         call mpas_dmpar_max_int(domain % dminfo, procCellIceMaskChanged, cellIceMaskChanged)
          ! Do the same for the Dirichlet b.c. mask
          call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged)
          call mpas_dmpar_max_int(domain % dminfo, procDirichletMaskChanged, dirichletMaskChanged)
-         !print *,'dirichletMaskChanged', dirichletMaskChanged
+         ! Do the same for the Floating edges
+         call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged)
+         call mpas_dmpar_max_int(domain % dminfo, procFloatingEdgesChanged, floatingEdgesChanged)
+         !print *,'floatingEdgesChanged', floatingEdgesChanged
 
       endif
 
@@ -479,7 +518,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
          call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
 
 
@@ -844,7 +883,7 @@ subroutine calculate_beta(block, err)
 
       call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve)
       call mpas_pool_get_array(velocityPool, 'beta', beta)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
       betaSolve = beta
 
@@ -1028,7 +1067,7 @@ subroutine remove_fast_ice(meshPool, uReconstructX, uReconstructY, normalVelocit
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
       call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F
index 53a8e32e17..f1fa8bd529 100644
--- a/src/core_landice/mode_forward/mpas_li_velocity_external.F
+++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F
@@ -452,7 +452,9 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       logical, pointer :: config_output_external_velocity_solver_data
       real (kind=RKIND), pointer :: config_ice_density
       integer, pointer :: anyDynamicVertexMaskChanged
+      integer, pointer :: cellIceMaskChanged
       integer, pointer :: dirichletMaskChanged
+      integer, pointer :: floatingEdgesChanged
       integer, pointer :: nEdges
       integer, pointer :: nVertLevels
       integer, pointer :: timestepNumber
@@ -486,7 +488,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
       call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
       call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
 
@@ -501,9 +503,11 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       call mpas_pool_get_array(velocityPool, 'uReconstructZ', uReconstructZ)
       call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve)
       call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged)
+      call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged)
       call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged)
       call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1)
-      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
+      call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged)
       call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
 
       ! Hydro variables
@@ -521,7 +525,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       ! Note these functions will always be called on the first solve because we
       ! initialize vertexMask to garbage which sets anyDynamicVertexMaskChanged to 1.
       if ((anyDynamicVertexMaskChanged == 1) .or. (config_always_compute_fem_grid) .or. &
-          (dirichletMaskChanged == 1) ) then
+          (dirichletMaskChanged == 1) .or. (floatingEdgesChanged == 1) .or. (cellIceMaskChanged == 1)) then
          call mpas_log_write("Generating new external velocity solver FEM grid.", flushNow=.true.)
          call generate_fem_grid(config_velocity_solver, vertexMask, cellMask, dirichletVelocityMask, &
               floatingEdges, layerThicknessFractions, lowerSurface, thickness, err)
@@ -805,13 +809,13 @@ subroutine li_velocity_external_write_albany_mesh(domain)
       call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
       call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
       call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
 
       ! Velocity variables
       call mpas_pool_get_array(velocityPool, 'beta', beta)
-      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)