eCLM-ParFlow: Couple surface water and route all layer-1 water fluxes to ParFlow

src/clm5/biogeophys/SnowHydrologyMod.F90

@@ -287,6 +287,10 @@ subroutine SnowWater(bounds, &
     real(r8) :: vol_ice(bounds%begc:bounds%endc,-nlevsno+1:0)      ! partial volume of ice lens in layer
     real(r8) :: eff_porosity(bounds%begc:bounds%endc,-nlevsno+1:0) ! effective porosity = porosity - vol_ice
     real(r8) :: mss_liqice(bounds%begc:bounds%endc,-nlevsno+1:0)   ! mass of liquid+ice in a layer
+#ifdef COUP_OAS_PFL
+    real(r8) :: snowliq_to_soil(bounds%begc:bounds%endc)           ! liquid from the bottom snow layer to soil layer [mm]
+                                                                   ! <0: snow liquid deficit refilled from soil
+#endif
     !-----------------------------------------------------------------------
 
     associate( &
@@ -333,6 +337,9 @@ subroutine SnowWater(bounds, &
     do fc = 1,num_snowc
        c = filter_snowc(fc)
        l=col%landunit(c)
+#ifdef COUP_OAS_PFL
+       snowliq_to_soil(c) = 0._r8
+#endif
 
        wgdif = h2osoi_ice(c,snl(c)+1) &
             + frac_sno_eff(c) * (qflx_dew_snow(c) - qflx_sub_snow(c)) * dtime
@@ -352,6 +359,11 @@ subroutine SnowWater(bounds, &
              if (wgdif >= 0._r8) exit
              h2osoi_liq(c,j) = 0._r8
              h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + wgdif
+#ifdef COUP_OAS_PFL
+             ! Deficit moved from the bottom snow layer into the soil, which is
+             ! overwritten by ParFlow. Take it from ParFlow instead.
+             if (j == 0) snowliq_to_soil(c) = wgdif
+#endif
           enddo
        end if
     end do
@@ -548,6 +560,9 @@ subroutine SnowWater(bounds, &
 
        qflx_top_soil(c) = (qout(c) / dtime) &
             + (1.0_r8 - frac_sno_eff(c)) * qflx_rain_grnd(c)
+#ifdef COUP_OAS_PFL
+       qflx_top_soil(c) = qflx_top_soil(c) + snowliq_to_soil(c) / dtime
+#endif
        int_snow(c) = int_snow(c) + frac_sno_eff(c) &
                      * (qflx_dew_snow(c) + qflx_dew_grnd(c) + qflx_rain_grnd(c)) * dtime
     end do
@@ -831,6 +846,9 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, &
          snw_rds          => waterstate_inst%snw_rds_col         , & ! Output: [real(r8) (:,:) ] effective snow grain radius (col,lyr) [microns, m^-6]
 
          qflx_sl_top_soil => waterflux_inst%qflx_sl_top_soil_col , & ! Output: [real(r8) (:)   ] liquid water + ice from layer above soil to top soil layer or sent to qflx_qrgwl (mm H2O/s)
+#ifdef COUP_OAS_PFL
+         qflx_pfl_top     => waterflux_inst%qflx_pfl_top_col     , & ! Output: [real(r8) (:)   ] water added to soil layer 1 after ParFlow exchange (mm H2O/s)
+#endif
 
          snl              => col%snl                             , & ! Output: [integer  (:)   ] number of snow layers
          dz               => col%dz                              , & ! Output: [real(r8) (:,:) ] layer depth (m)
@@ -874,7 +892,16 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, &
           ! use 0.01 to avoid runaway ice buildup
           if (h2osoi_ice(c,j) <= .01_r8) then
              if (ltype(l) == istsoil .or. urbpoi(l) .or. ltype(l) == istcrop) then
+#ifdef COUP_OAS_PFL
+                if (j == 0 .and. col%hydrologically_active(c)) then
+                   ! Pass the snow water to ParFlow. The ice phase is kept in eCLM.
+                   qflx_pfl_top(c) = qflx_pfl_top(c) + (h2osoi_liq(c,j) + h2osoi_ice(c,j))/dtime
+                else
+                   h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
+                end if
+#else
                 h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) + h2osoi_liq(c,j)
+#endif
                 h2osoi_ice(c,j+1) = h2osoi_ice(c,j+1) + h2osoi_ice(c,j)
 
                 if (j == 0) then
@@ -996,7 +1023,16 @@ subroutine CombineSnowLayers(bounds, num_snowc, filter_snowc, &
              ! this is where water is transfered from layer 0 (snow) to layer 1 (soil)
              if (ltype(l) == istsoil .or. urbpoi(l) .or. ltype(l) == istcrop) then
                 h2osoi_liq(c,0) = 0.0_r8
+#ifdef COUP_OAS_PFL
+                if (col%hydrologically_active(c)) then
+                   ! Pass the snow water to ParFlow
+                   qflx_pfl_top(c) = qflx_pfl_top(c) + zwliq(c)/dtime
+                else
+                   h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c)
+                end if
+#else
                 h2osoi_liq(c,1) = h2osoi_liq(c,1) + zwliq(c)
+#endif
              end if
              if (ltype(l) == istwet) then
                 h2osoi_liq(c,0) = 0.0_r8

src/clm5/biogeophys/SoilHydrologyMod.F90

@@ -450,6 +450,14 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
                 qflx_evap(c)=qflx_ev_soil(c)
              endif
 
+#ifdef COUP_OAS_PFL
+             ! All surface input minus evaporation (and surface runoff) from soil and
+             ! surface water is passed to ParFlow, without the eCLM infiltration capacity.
+             qflx_infl(c) = qflx_top_soil(c) - qflx_surf(c) &
+                  - (1.0_r8 - fsno - frac_h2osfc(c))*qflx_evap(c) &
+                  - frac_h2osfc(c)*qflx_ev_h2osfc(c)
+             qflx_h2osfc_surf(c) = 0._r8
+#else
              !1. partition surface inputs between soil and h2osfc
              qflx_in_soil(c) = (1._r8 - frac_h2osfc(c)) * (qflx_top_soil(c)  - qflx_surf(c))
              qflx_in_h2osfc(c) = frac_h2osfc(c) * (qflx_top_soil(c)  - qflx_surf(c))          
@@ -545,6 +553,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
              !7. remove drainage from h2osfc and add to qflx_infl
              h2osfc(c) = h2osfc(c) - qflx_h2osfc_drain(c) * dtime
              qflx_infl(c) = qflx_infl(c) + qflx_h2osfc_drain(c)
+#endif
           else
              ! non-vegetated landunits (i.e. urban) use original CLM4 code
              if (snl(c) >= 0) then
@@ -2279,6 +2288,9 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, &
      ! !LOCAL VARIABLES:
      integer  :: c,j,fc,i                                ! indices
      real(r8) :: dtime                                   ! land model time step (sec)
+#ifdef COUP_OAS_PFL
+     real(r8) :: h2osoi_ice1_old                         ! ice of soil layer 1 before condensation/sublimation (kg/m2)
+#endif
      !-----------------------------------------------------------------------
 
      associate(                                                            & 
@@ -2288,6 +2300,9 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, &
           frac_h2osfc        =>    waterstate_inst%frac_h2osfc_col       , & ! Input:  [real(r8) (:)   ]                                                    
           qflx_dew_grnd      =>    waterflux_inst%qflx_dew_grnd_col      , & ! Input:  [real(r8) (:)   ]  ground surface dew formation (mm H2O /s) [+]      
           qflx_dew_snow      =>    waterflux_inst%qflx_dew_snow_col      , & ! Input:  [real(r8) (:)   ]  surface dew added to snow pack (mm H2O /s) [+]    
+#ifdef COUP_OAS_PFL
+          qflx_pfl_top       =>    waterflux_inst%qflx_pfl_top_col       , & ! Output: [real(r8) (:)   ]  water added to soil layer 1 after ParFlow exchange (mm H2O/s)
+#endif
           qflx_sub_snow      =>    waterflux_inst%qflx_sub_snow_col       & ! Output: [real(r8) (:)   ]  sublimation rate from snow pack (mm H2O /s) [+]   
           )
 
@@ -2303,14 +2318,24 @@ subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, &
           if (snl(c)+1 >= 1) then
 
              ! make consistent with how evap_grnd removed in infiltration
+#ifdef COUP_OAS_PFL
+             ! Pass the dew to ParFlow.
+             qflx_pfl_top(c) = qflx_pfl_top(c) + (1._r8 - frac_h2osfc(c))*qflx_dew_grnd(c)
+             h2osoi_ice1_old = h2osoi_ice(c,1)
+#else
              h2osoi_liq(c,1) = h2osoi_liq(c,1) + (1._r8 - frac_h2osfc(c))*qflx_dew_grnd(c) * dtime
+#endif
              h2osoi_ice(c,1) = h2osoi_ice(c,1) + (1._r8 - frac_h2osfc(c))*qflx_dew_snow(c) * dtime
              if (qflx_sub_snow(c)*dtime > h2osoi_ice(c,1)) then
                 qflx_sub_snow(c) = h2osoi_ice(c,1)/dtime
                 h2osoi_ice(c,1) = 0._r8
              else
                 h2osoi_ice(c,1) = h2osoi_ice(c,1) - (1._r8 - frac_h2osfc(c)) * qflx_sub_snow(c) * dtime
              end if
+#ifdef COUP_OAS_PFL
+             ! Pass the ice change to ParFlow in addition. SPo: check for eff. poro. coupling
+             qflx_pfl_top(c) = qflx_pfl_top(c) + (h2osoi_ice(c,1) - h2osoi_ice1_old) / dtime
+#endif
           end if
 
        end do
@@ -2369,10 +2394,13 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
 
      associate(                                                            &
           dz                 =>    col%dz                                , & ! Input:  [real(r8) (:,:) ] layer depth (m)
+          qflx_pfl_top       =>    waterflux_inst%qflx_pfl_top_col       , & ! Input:  [real(r8) (:)   ] water added to soil layer 1 after ParFlow exchange (mm H2O/s)
+          qflx_h2osfc_to_ice =>    waterflux_inst%qflx_h2osfc_to_ice_col , & ! Input:  [real(r8) (:)   ] surface water converted to ice (mm H2O /s)
           qflx_snwcp_liq     =>    waterflux_inst%qflx_snwcp_liq_col     , & ! excess rainfall due to snow capping (mm H2O /s) [+]
           qflx_drain         =>    waterflux_inst%qflx_drain_col         , & ! sub-surface runoff (mm H2O /s)
           qflx_drain_perched =>    waterflux_inst%qflx_drain_perched_col , & ! perched wt sub-surface runoff (mm H2O /s)         
           qflx_qrgwl         =>    waterflux_inst%qflx_qrgwl_col         , & ! qflx_surf at glaciers, wetlands, lakes (mm H2O /s)
+          qflx_surf          =>    waterflux_inst%qflx_surf_col          , & ! surface runoff (mm H2O /s)
           qflx_rsub_sat      =>    waterflux_inst%qflx_rsub_sat_col      , & ! soil saturation excess [mm h2o/s]
           qflx_infl          =>    waterflux_inst%qflx_infl_col          , & ! infiltration (mm H2O /s)
           qflx_rootsoi       =>    waterflux_inst%qflx_rootsoi_col       , & ! vegetation/soil water exchange (mm H2O/s) (+ = to atm) 
@@ -2395,11 +2423,20 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
                    qflx_parflow(c,j) = -qflx_rootsoi(c,j) !mm/s
                end if
             end do
+
+            ! Water added to the top soil layer by eCLM after the ParFlow state was applied
+            ! (condensation/sublimation in RenewCondensation, snow water from CombineSnowLayers)
+            qflx_parflow(c,1) = qflx_parflow(c,1) + qflx_pfl_top(c)
+            ! Remove what froze into ice from surface water
+            qflx_parflow(c,1) = qflx_parflow(c,1) - qflx_h2osfc_to_ice(c)
+            ! Excess rainfall due to snow capping pass to ParFlow instead of river routing
+            qflx_parflow(c,1) = qflx_parflow(c,1) + qflx_snwcp_liq(c)
+
             ! Compute subsurface run-off (mm/s)
             qflx_drain(c) = -sum(qflx_parflow(c,:))
             qflx_drain_perched(c) = 0._r8  
             qflx_rsub_sat(c)      = 0._r8
-            qflx_qrgwl(c)         = qflx_snwcp_liq(c)   ! Set imbalance for snow capping
+            qflx_qrgwl(c)         = 0._r8
             ! Convert eCLM fluxes (mm/s) to ParFlow fluxes (1/hr):
             !         1/hr            =           [mm/s]          *   [s/hr]   *  [m/mm]  *         [1/m]
             qflx_parflow(c,1:nlevsoi) = qflx_parflow(c,1:nlevsoi) * sec_per_hr * m_per_mm * (1._r8/dz(c,1:nlevsoi))
@@ -2412,6 +2449,17 @@ subroutine ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
                qflx_drain(c) = 0._r8
                ! This must be done for roofs and impervious road (walls will be zero)
                qflx_qrgwl(c) = qflx_snwcp_liq(c)
+
+               ! Instead of leaving as river runoff, impervious urban runoff enters
+               ! ParFlow through layer 1. It is zeroed here so it is not counted twice.
+               qflx_parflow(c,1:nlevsoi) = 0._r8
+               qflx_parflow(c,1)         = qflx_surf(c) + qflx_qrgwl(c)
+               qflx_surf(c)              = 0._r8
+               qflx_qrgwl(c)             = 0._r8
+
+               qflx_drain(c) = -sum(qflx_parflow(c,:))
+
+               qflx_parflow(c,1:nlevsoi) = qflx_parflow(c,1:nlevsoi) * sec_per_hr * m_per_mm * (1._r8/dz(c,1:nlevsoi))
             end if
          end do
      end associate

src/clm5/biogeophys/SoilWaterMovementMod.F90

@@ -1503,6 +1503,7 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, &
          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)
+         qflx_pfl_top      =>    waterflux_inst%qflx_pfl_top_col    , & ! Output: [real(r8) (:)   ]  water added to soil layer 1 after exchange (mm H2O/s)
          pfl_psi           =>    waterstate_inst%pfl_psi_col          & ! Input:  [real(r8) (:,:) ]  ParFlow pressure head (mm)
          )  ! end associate statement
 
@@ -1517,11 +1518,13 @@ subroutine soilwater_parflow(bounds, num_hydrologyc, &
             do j = 1, nlevgrnd
                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
+               ! Saturated/ponded cells (pfl_psi > 0) have zero matric potential
+               smp_l(c,j) = max(smpmin(c), min(0._r8, pfl_psi(c,j)))
             end do
 
+            ! Reset flux after passing ParFlow state to eCLM
+            qflx_pfl_top(c) = 0._r8
+
             ! ParFlow replaces the eCLM soil water solver. Recompute hk_l from the
             ! ParFlow water content using the same formulation as compute_hydraulic_properties.
             nlayers = nbedrock(c)

src/clm5/biogeophys/WaterfluxType.F90

@@ -74,6 +74,7 @@ module WaterfluxType
      real(r8), pointer :: qflx_rootsoi_col         (:,:) ! col root and soil water exchange [mm H2O/s] [+ into root]
 #ifdef COUP_OAS_PFL
      real(r8), pointer :: qflx_parflow_col         (:,:) ! col source/sink flux per soil layer sent to ParFlow [1/hr] [- out from root]
+     real(r8), pointer :: qflx_pfl_top_col         (:)   ! col water added to soil layer 1 by eCLM within time step, sent to ParFlow (mm H2O/s)
 #endif
      real(r8), pointer :: qflx_infl_col            (:)   ! col infiltration (mm H2O /s)
      real(r8), pointer :: qflx_surf_col            (:)   ! col surface runoff (mm H2O /s)
@@ -219,6 +220,7 @@ subroutine InitAllocate(this, bounds)
     allocate(this%qflx_rootsoi_col         (begc:endc,1:nlevsoi))    ; this%qflx_rootsoi_col         (:,:) = nan
 #ifdef COUP_OAS_PFL    
     allocate(this%qflx_parflow_col         (begc:endc,1:nlevsoi))    ; this%qflx_parflow_col         (:,:) = nan
+    allocate(this%qflx_pfl_top_col         (begc:endc))              ; this%qflx_pfl_top_col         (:)   = nan
 #endif  
     allocate(this%qflx_infl_col            (begc:endc))              ; this%qflx_infl_col            (:)   = nan
     allocate(this%qflx_surf_col            (begc:endc))              ; this%qflx_surf_col            (:)   = nan

src/clm5/main/clm_driver.F90

@@ -1224,6 +1224,9 @@ subroutine clm_drv_init(bounds, &
     use WaterStateType     , only : waterstate_type
     use WaterFluxType      , only : waterflux_type
     use EnergyFluxType     , only : energyflux_type
+#ifdef COUP_OAS_PFL
+    use landunit_varcon    , only : istsoil, istcrop
+#endif
     !
     ! !ARGUMENTS:
     type(bounds_type)     , intent(in)    :: bounds  
@@ -1257,6 +1260,7 @@ subroutine clm_drv_init(bounds, &
 #ifdef COUP_OAS_PFL
          pfl_psi            => waterstate_inst%pfl_psi_col               , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
          pfl_h2osoi_liq     => waterstate_inst%pfl_h2osoi_liq_col        , & ! Input:  [real(r8) (:,:) ]  COUP_OAS_PFL
+         h2osfc             => waterstate_inst%h2osfc_col                , & ! Output: [real(r8) (:)   ]  surface water (mm)
 #endif
          elai               => canopystate_inst%elai_patch               , & ! Input:  [real(r8) (:)   ]  one-sided leaf area index with burying by snow    
          esai               => canopystate_inst%esai_patch               , & ! Input:  [real(r8) (:)   ]  one-sided stem area index with burying by snow    
@@ -1313,6 +1317,11 @@ subroutine clm_drv_init(bounds, &
           g = col%gridcell(c)
           pfl_psi(c,:) = atm2lnd_inst%pfl_psi_grc(g,:)
           pfl_h2osoi_liq(c,:) = atm2lnd_inst%pfl_h2osoi_liq_grc(g,:)
+          ! ParFlow's positive pressure head of the top cell is the ponding depth.
+          ! Surface water is diagnosed from it to sync eCLM and ParFlow state.
+          if (lun%itype(col%landunit(c)) == istsoil .or. lun%itype(col%landunit(c)) == istcrop) then
+            h2osfc(c) = max(0._r8, pfl_psi(c,1))
+          end if
         end if
       end do
 #endif