diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
index 148dda8ddff2..02fdd3c70058 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
@@ -1045,6 +1045,11 @@ subroutine advection_solver(domain, err)
       real (kind=RKIND), dimension(:), pointer :: thicknessNew
       real (kind=RKIND), dimension(:), pointer :: thickness
       real (kind=RKIND), dimension(:,:), pointer :: layerThickness
+      real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
+      real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity
+
+      integer, dimension(:), pointer :: edgeMask
+      real (kind=RKIND) :: allowableDtACFL
 
       real (kind=RKIND), dimension(:,:), pointer :: temperature
       real (kind=RKIND), dimension(:,:), pointer :: waterFrac
@@ -1094,6 +1099,18 @@ subroutine advection_solver(domain, err)
             call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
             call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
             call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
+            call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity)
+            call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
+            call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity) 
+
+            call li_layer_normal_velocity( &
+                 meshPool,               &
+                 normalVelocity,         &
+                 edgeMask,               &
+                 layerNormalVelocity,    &
+                 allowableDtACFL,        &
+                 err_tmp)
+            err = ior(err,err_tmp)
 
             call calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, err_tmp)
             err = ior(err,err_tmp)
@@ -1458,7 +1475,7 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er
       ! local variables
       !-----------------------------------------------------------------
       real (kind=RKIND), dimension(:), pointer :: thickness
-      real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, normalVelocity
+      real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, layerNormalVelocity
       integer, dimension(:,:), pointer :: cellsOnEdge
       integer, pointer :: nEdges, nVertLevels
       integer :: iEdge, cell1, cell2, k
@@ -1472,7 +1489,7 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er
       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, 'layerThicknessEdge', layerThicknessEdge)
-      call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity)
+      call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity)
 
       ! Note: SIA velocity solver uses its own local calculation of h_edge that is always 2nd order.
       ! Note: ocn_diagnostic_solve in mpas_ocn_tendency.F has 2, 3, & 4th order calculations for h_edge that can be used.
@@ -1486,7 +1503,7 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er
          cell2 = cellsOnEdge(2,iEdge)
          do k=1, nVertLevels
             ! Calculate h on edges using first order
-            VelSign = sign(1.0_RKIND, normalVelocity(k, iEdge))
+            VelSign = sign(1.0_RKIND, layerNormalVelocity(k, iEdge))
             layerThicknessEdge(k,iEdge) = max(VelSign * layerThickness(k, cell1), &
                VelSign * (-1.0_RKIND) * layerThickness(k, cell2))
             ! + velocity goes from index 1 to 2 in the cellsOnEdge array.