diff --git a/src/clm5/biogeophys/SoilWaterMovementMod.F90 b/src/clm5/biogeophys/SoilWaterMovementMod.F90
index f86212882..f68d9a041 100644
--- a/src/clm5/biogeophys/SoilWaterMovementMod.F90
+++ b/src/clm5/biogeophys/SoilWaterMovementMod.F90
@@ -1488,20 +1488,24 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, &
 
       associate(&
          smp_l             =>    soilstate_inst%smp_l_col           , & ! Input:  [real(r8) (:,:) ]  soil matrix potential [mm]         
-         h2osoi_liq        =>    waterstate_inst%h2osoi_liq_col     , & ! Input:  [real(r8) (:,:) ]  liquid water (kg/m2)
+         h2osoi_liq        =>    waterstate_inst%h2osoi_liq_col     , & ! Output: [real(r8) (:,:) ]  liquid water (kg/m2)
+         h2osoi_ice        =>    waterstate_inst%h2osoi_ice_col     , & ! Output: [real(r8) (:,:) ]  ice lens (kg/m2)
          smpmin            =>    soilstate_inst%smpmin_col          , & ! Input:  [real(r8) (:)   ]  restriction for min of soil potential (mm)
          pfl_h2osoi_liq    =>    waterstate_inst%pfl_h2osoi_liq_col , & ! Input:  [real(r8) (:,:) ]  ParFlow soil water (mm)
          pfl_psi           =>    waterstate_inst%pfl_psi_col          & ! Input:  [real(r8) (:,:) ]  ParFlow pressure head (mm)
          )  ! end associate statement
 
 
-         ! COUP_OAS_PFL
-         ! TODO
+         ! Exchange of soil water and pressure head between ParFlow and eCLM
+         ! ParFlow has no ice phase, so that the received water it the total water content.
+         ! Keep the ice what PhaseChanged diagnosed earlier in the time step and
+         ! attribute the remainder to liquid, so that liq and ice matches the received total.
          do fc = 1, num_hydrologyc
             c = filter_hydrologyc(fc)
 
             do j = 1, nlevgrnd
-               h2osoi_liq(c,j) = pfl_h2osoi_liq(c,j)
+               h2osoi_ice(c,j) = min(h2osoi_ice(c,j), pfl_h2osoi_liq(c,j))
+               h2osoi_liq(c,j) = max(0._r8, pfl_h2osoi_liq(c,j) - h2osoi_ice(c,j))
                if (pfl_psi(c,j) <= 0) then
                   smp_l(c,j) = max(smpmin(c), pfl_psi(c,j))
                end if