Merge remote-tracking branch 'donghuix/donghuix/lnd/wetland-inundatio…

…n-scheme' into next (PR #5974) In default, the infiltration rate is the same for the non-water surface and surface water storage within the grid cell. This assumption results in a significant underestimation of surface water because the infiltration in the surface water storage is not constrained. In this PR, we introduced a modified infiltration scheme to constrain the infiltration in the surface water storage. The infiltration rate in the surface water storage is smaller than the infiltration rate in the non-water surface. More details about this new scheme and corresponding evaluations can be found in the preprint: https://doi.org/10.21203/rs.3.rs-2733749/v1 The manuscript is still under review. The new infiltration scheme can be turned by adding use_modified_infil = .true. in user_nl_elm. We passed all the tests in e3sm_land_developer test suite. We also created a test for the simulation of surface water dynamics with the new infiltration scheme and calibrated parameter for NLDAS domain: surface_water_dynamics. This test was added into e3sm_land_developer test suite. [BFB] [STEALTH]
E3SM-Project · Nov 14, 2023 · ba4d10c · ba4d10c
2 parents bd765d1 + 0b7cdf4
commit ba4d10c
Show file tree

Hide file tree

Showing 15 changed files with 138 additions and 31 deletions.
diff --git a/cime_config/tests.py b/cime_config/tests.py
@@ -94,7 +94,8 @@
             "SMS.r05_r05.IELM.elm-topounit",
             "ERS.ELM_USRDAT.I1850ELM.elm-usrdat",
             "ERS.f09_f09.IELM.elm-lnd_rof_2way",
-            "ERS.r05_r05.IELM.elm-V2_ELM_MOSART_features"
+            "ERS.r05_r05.IELM.elm-V2_ELM_MOSART_features",
+            "ERS.ELM_USRDAT.IELM.elm-surface_water_dynamics"
             )
         },
 

diff --git a/components/elm/bld/ELMBuildNamelist.pm b/components/elm/bld/ELMBuildNamelist.pm
@@ -2385,6 +2385,7 @@ sub setup_logic_demand {
   $settings{'use_snicar_ad'}       = $nl_flags->{'use_snicar_ad'};
   $settings{'use_century_decomp'}  = $nl_flags->{'use_century_decomp'};
   $settings{'use_crop'}            = $nl_flags->{'use_crop'};
+  $settings{'use_modified_infil'}  = $nl_flags->{'use_modified_infil'};
 
   my $demand = $nl->get_value('clm_demand');
   if (defined($demand)) {

diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml
@@ -1993,6 +1993,15 @@ Toggle to turn on the soil erosion model.
 Toggle to turn on the soil erosion and OM pools feedback.
 </entry>
 
+<!-- ========================================================================================  -->
+<!-- Namelist options for using modifed infiltration scheme in h2osfc                          -->
+<!-- ========================================================================================  -->
+
+<entry id="use_modified_infil" type="logical" category="clm_physics"
+       group="elm_inparm" valid_values="" value=".false.">
+Toggle to use modifed infiltration scheme in h2osfc     .
+</entry>
+
 <!-- ========================================================================================  -->
 <!-- Namelist options for land river coulping                                                  -->
 <!-- ========================================================================================  -->

diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/surface_water_dynamics/shell_commands b/components/elm/cime_config/testdefs/testmods_dirs/elm/surface_water_dynamics/shell_commands
@@ -0,0 +1,10 @@
+./xmlchange DATM_MODE=CLMMOSARTTEST
+
+./xmlchange LND_DOMAIN_FILE=domain.lnd.nldas2_0224x0464_c110415.nc
+./xmlchange ATM_DOMAIN_FILE=domain.lnd.nldas2_0224x0464_c110415.nc
+./xmlchange LND_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+./xmlchange ATM_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_END -val 2000
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_START -val 2000
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_ALIGN -val 1
diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/surface_water_dynamics/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/surface_water_dynamics/user_nl_elm
@@ -0,0 +1,2 @@
+fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_nldas2_simyr2000_c181207_surfacewater.nc'
+use_modified_infil = .true.
diff --git a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 b/components/elm/src/biogeophys/CanopyHydrologyMod.F90
@@ -809,7 +809,8 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
           h2osfc       => col_ws%h2osfc       , & ! Output: [real(r8) (:)   ] surface water (mm)                                
           frac_sno     => col_ws%frac_sno     , & ! Output: [real(r8) (:)   ] fraction of ground covered by snow (0 to 1)       
           frac_sno_eff => col_ws%frac_sno_eff , & ! Output: [real(r8) (:)   ] eff. fraction of ground covered by snow (0 to 1)  
-          frac_h2osfc  => col_ws%frac_h2osfc    & ! Output: [real(r8) (:)   ] col fractional area with surface water greater than zero 
+          frac_h2osfc  => col_ws%frac_h2osfc  , & ! Output: [real(r8) (:)   ] col fractional area with surface water greater than zero 
+          frac_h2osfc_act => col_ws%frac_h2osfc_act & ! Output: [real(r8) (:)   ] col fractional area with surface water greater than zero
           )
 
        ! arbitrary lower limit on h2osfc for safer numerics...
@@ -848,6 +849,8 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
                 qflx_h2osfc2topsoi(c) = h2osfc(c)/dtime                
                 h2osfc(c)=0._r8
              endif
+
+             frac_h2osfc_act(c) = frac_h2osfc(c)
 
              if (.not. present(no_update)) then
 
@@ -869,7 +872,8 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
           else !if landunit not istsoil/istcrop, set frac_h2osfc to zero
 
              frac_h2osfc(c) = 0._r8
-
+             frac_h2osfc_act(c) = 0._r8
+
           endif
 
        end do

diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90
@@ -9,6 +9,7 @@ module SoilHydrologyMod
   use decompMod         , only : bounds_type
   use elm_varctl        , only : iulog, use_vichydro
   use elm_varctl        , only : use_lnd_rof_two_way, lnd_rof_coupling_nstep
+  use elm_varctl        , only : use_modified_infil
   use elm_varcon        , only : e_ice, denh2o, denice, rpi
   use EnergyFluxType    , only : energyflux_type
   use SoilHydrologyType , only : soilhydrology_type
@@ -111,6 +112,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
          b_infil          =>    soilhydrology_vars%b_infil_col      , & ! Input:  [real(r8) (:)   ]  VIC b infiltration parameter
          moist            =>    soilhydrology_vars%moist_col        , & ! Input:  [real(r8) (:,:) ]  soil moisture in each VIC layers (liq, mm)
          hkdepth          =>    soilhydrology_vars%hkdepth_col      , & ! Input:  [real(r8) (:)   ]  decay factor (m)
+         fover            =>    soilhydrology_vars%fover            , & ! Input:  [real(r8) (:)   ]  decay factor for saturation fraction (m)
          origflag         =>    soilhydrology_vars%origflag         , & ! Input:  logical
          fcov             =>    soilhydrology_vars%fcov_col         , & ! Output: [real(r8) (:)   ]  fractional impermeable area
          fsat             =>    soilhydrology_vars%fsat_col         , & ! Output: [real(r8) (:)   ]  fractional area with water table at surface
@@ -146,7 +148,8 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
 
       do fc = 1, num_hydrologyc
          c = filter_hydrologyc(fc)
-         fff(c) = 0.5_r8
+         g = col_pp%gridcell(c)
+         fff(c) = fover(g)
          if (zengdecker_2009_with_var_soil_thick) then
             nlevbed = nlev2bed(c)
             fff(c) = 0.5_r8 * col_pp%zi(c,nlevsoi) / min(col_pp%zi(c,nlevbed), col_pp%zi(c,nlevsoi))
@@ -261,7 +264,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
      use shr_const_mod    , only : shr_const_pi
      use elm_varpar       , only : nlayer, nlayert
      use elm_varpar       , only : nlevsoi, nlevgrnd
-     use elm_varcon       , only : denh2o, denice, roverg, wimp, pc, mu, tfrz
+     use elm_varcon       , only : denh2o, denice, roverg, wimp, mu, tfrz
      use elm_varcon       , only : pondmx, watmin
      use column_varcon    , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv
      use landunit_varcon  , only : istsoil, istcrop
@@ -317,6 +320,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
      real(r8) :: top_ice(bounds%begc:bounds%endc)           ! temporary, ice len in top VIC layers
      real(r8) :: top_icefrac                                ! temporary, ice fraction in top VIC layers
      real(r8) :: h2osoi_left_vol1                           ! temporary, available volume in the first soil layer
+     real(r8) :: pc                                         ! temporary, threhold for surface water storage to outflow
      !-----------------------------------------------------------------------
 
      associate(                                                    &
@@ -329,6 +333,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
           t_soisno             =>    col_es%t_soisno             , & ! Input:  [real(r8) (:,:) ]  soil temperature (Kelvin)
 
           frac_h2osfc          =>    col_ws%frac_h2osfc          , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by surface water (0 to 1)
+          frac_h2osfc_act      =>    col_ws%frac_h2osfc_act      , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by surface water (0 to 1) without adjustment from snow fraction
           frac_sno             =>    col_ws%frac_sno_eff         , & ! Input:  [real(r8) (:)   ]  fraction of ground covered by snow (0 to 1)
           h2osoi_ice           =>    col_ws%h2osoi_ice           , & ! Input:  [real(r8) (:,:) ]  ice lens (kg/m2)
           h2osoi_liq           =>    col_ws%h2osoi_liq           , & ! Input:  [real(r8) (:,:) ]  liquid water (kg/m2)
@@ -370,6 +375,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
           ice                  =>    soilhydrology_vars%ice_col              , & ! Input:  [real(r8) (:,:) ]  ice len in each VIC layers(ice, mm)
           i_0                  =>    soilhydrology_vars%i_0_col              , & ! Input:  [real(r8) (:)   ]  column average soil moisture in top VIC layers (mm)
           h2osfcflag           =>    soilhydrology_vars%h2osfcflag           , & ! Input:  logical
+          pc_grid              =>    soilhydrology_vars%pc                   , & ! Input:  [real(r8) (:)   ]  threshold for outflow from surface water storage
           icefrac              =>    soilhydrology_vars%icefrac_col            & ! Output: [real(r8) (:,:) ]  fraction of ice
               )
 
@@ -387,8 +393,10 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
        end do
 
        do fc = 1, num_hydrologyc
-          c = filter_hydrologyc(fc)
-          g = cgridcell(c)
+          c  = filter_hydrologyc(fc)
+          g  = cgridcell(c)
+          pc = pc_grid(g)
+
           ! partition moisture fluxes between soil and h2osfc
           if (lun_pp%itype(col_pp%landunit(c)) == istsoil .or. lun_pp%itype(col_pp%landunit(c))==istcrop) then
 
@@ -455,9 +463,24 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
                 rsurf_vic = min(qflx_in_soil(c), rsurf_vic)
                 qinmax = (1._r8 - fsat(c)) * 10._r8**(-e_ice*top_icefrac)*(qflx_in_soil(c) - rsurf_vic)
              else
-                qinmax=(1._r8 - fsat(c)) * minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
+                if ( use_modified_infil ) then
+                  qinmax=minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
+                else
+                  qinmax=(1._r8 - fsat(c)) * minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
+                end if
              end if
-
+
+             if ( use_modified_infil ) then
+                ! Assume frac_h2osfc occurs on fsat
+                if ( frac_h2osfc(c) >= fsat(c) ) then
+                   qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) -  (1.0_r8 - frac_h2osfc(c))*qinmax)
+                else
+                   qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) -  (1.0_r8 - fsat(c))*qinmax)
+                end if
+             else
+                qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) -  (1.0_r8 - frac_h2osfc(c))*qinmax)
+             end if
+
              if (use_lnd_rof_two_way) then
                 qflx_infl_excess(c) = max(0._r8,qflx_in_soil(c) -  (1.0_r8 - frac_h2osfc(c) - frac_h2orof(c))*qinmax)
              else
@@ -472,10 +495,22 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
              !5. surface runoff from h2osfc
              if (h2osfcflag==1) then
                 ! calculate runoff from h2osfc  -------------------------------------
-                if (frac_h2osfc(c) <= pc) then
-                   frac_infclust=0.0_r8
+                if (use_modified_infil) then
+                  if (frac_h2osfc_act(c) <= pc .and. frac_h2osfc(c) <= pc) then 
+                     frac_infclust=0.0_r8
+                  else
+                      if (frac_h2osfc(c) <= pc) then
+                        frac_infclust=(frac_h2osfc_act(c)-pc)**mu
+                      else
+                        frac_infclust=(frac_h2osfc(c)-pc)**mu
+                      endif
+                  endif
                 else
-                   frac_infclust=(frac_h2osfc(c)-pc)**mu
+                  if (frac_h2osfc(c) <= pc) then
+                    frac_infclust=0.0_r8
+                  else
+                    frac_infclust=(frac_h2osfc(c)-pc)**mu
+                  endif
                 endif
              endif
 
@@ -514,7 +549,17 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
                 h2osfc(c) = 0.0
                 qflx_h2osfc_drain(c)= 0._r8
              else
-                qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax,h2osfc(c)/dtime)
+                if ( use_modified_infil ) then
+                   ! Assume frac_h2osfc occurs on top of fsat
+                   if (frac_h2osfc(c) <= fsat(c)) then
+                     qflx_h2osfc_drain(c)=0._r8
+                   else
+                     qflx_h2osfc_drain(c)=min((frac_h2osfc(c)-fsat(c))*qinmax,h2osfc(c)/dtime)
+                   endif
+                else
+                   ! Original scheme
+                   qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax,h2osfc(c)/dtime)
+                end if
              endif
 
              if(h2osfcflag==0) then

diff --git a/components/elm/src/biogeophys/SoilHydrologyType.F90 b/components/elm/src/biogeophys/SoilHydrologyType.F90
@@ -65,7 +65,9 @@ Module SoilHydrologyType
      real(r8), pointer :: max_infil_col     (:)     => null()! col VIC maximum infiltration rate calculated in VIC
      real(r8), pointer :: i_0_col           (:)     => null()! col VIC average saturation in top soil layers
      real(r8), pointer :: ice_col           (:,:)   => null()! col VIC soil ice (kg/m2) for VIC soil layers
-
+     real(r8), pointer :: fover             (:)     => null()! decay factor for surface runoff
+     real(r8), pointer :: pc                (:)     => null()! surface water threshold probability
+
    contains
 
      procedure, public  :: Init
@@ -154,6 +156,9 @@ subroutine InitAllocate(this, bounds)
     allocate(this%max_infil_col     (begc:endc))                 ; this%max_infil_col     (:)     = spval
     allocate(this%i_0_col           (begc:endc))                 ; this%i_0_col           (:)     = spval
     allocate(this%ice_col           (begc:endc,nlayert))         ; this%ice_col           (:,:)   = spval
+
+    allocate(this%fover             (begg:endg))                 ; this%fover             (:)     = spval
+    allocate(this%pc                (begg:endg))                 ; this%pc                (:)     = spval
 
   end subroutine InitAllocate
 
@@ -231,8 +236,8 @@ subroutine InitCold(this, bounds)
     use elm_varctl      , only : fsurdat, iulog, use_vichydro, use_var_soil_thick
     use elm_varpar      , only : nlevsoi, nlevgrnd, nlevsno, nlevlak, nlevurb
     use elm_varcon      , only : denice, denh2o, sb, bdsno
-    use elm_varcon      , only : h2osno_max, zlnd, tfrz, spval, pc
-    use elm_varcon      , only : nlvic, dzvic, pc, mu, grlnd
+    use elm_varcon      , only : h2osno_max, zlnd, tfrz, spval
+    use elm_varcon      , only : nlvic, dzvic, grlnd
     use landunit_varcon , only : istice, istwet, istsoil, istdlak, istcrop, istice_mec
     use column_varcon   , only : icol_shadewall, icol_road_perv
     use column_varcon   , only : icol_road_imperv, icol_roof, icol_sunwall
@@ -542,18 +547,36 @@ subroutine InitCold(this, bounds)
        fdrain(:) = 2.5_r8
     end if
     call ncd_pio_closefile(ncid)
+
+    call getfil (fsurdat, locfn, 0)
+    call ncd_pio_openfile (ncid, locfn, 0)
+    call ncd_io(ncid=ncid, varname='fover', flag='read', data=this%fover, dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) then
+      this%fover(:) = 0.5_r8
+    end if
+    call ncd_pio_closefile(ncid)
+
+    call getfil (fsurdat, locfn, 0)
+    call ncd_pio_openfile (ncid, locfn, 0)
+    call ncd_io(ncid=ncid, varname='pc', flag='read', data=this%pc, dim1name=grlnd, readvar=readvar)
+    if (.not. readvar) then
+      this%pc(:) = 0.4
+    end if
+    call ncd_pio_closefile(ncid)
 
     associate(micro_sigma => col_pp%micro_sigma)
       do c = bounds%begc, bounds%endc
-
+
+         g = col_pp%gridcell(c)
+
          ! determine h2osfc threshold ("fill & spill" concept)
          ! set to zero for no h2osfc (w/frac_infclust =large)
-
+         
          this%h2osfc_thresh_col(c) = 0._r8
          if (micro_sigma(c) > 1.e-6_r8 .and. (this%h2osfcflag /= 0)) then
             d = 0.0
             do p = 1,4
-               fd   = 0.5*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0)))) - pc
+               fd   = 0.5*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0)))) - this%pc(g)
                dfdd = exp(-d**2/(2.0*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0*shr_const_pi))
                d    = d - fd/dfdd
             enddo
@@ -569,12 +592,8 @@ subroutine InitCold(this, bounds)
          endif
 
          ! set decay factor
-         if (use_lnd_rof_two_way) then
-            g = col_pp%gridcell(c)
-            this%hkdepth_col(c) = 1._r8/fdrain(g)
-         else
-            this%hkdepth_col(c) = 1._r8/2.5_r8
-         endif
+         g = col_pp%gridcell(c)
+         this%hkdepth_col(c) = 1._r8/fdrain(g)
 
       end do
     end associate

diff --git a/components/elm/src/biogeophys/SoilStateType.F90 b/components/elm/src/biogeophys/SoilStateType.F90
@@ -14,7 +14,7 @@ module SoilStateType
   use landunit_varcon , only : istice, istdlak, istwet, istsoil, istcrop, istice_mec
   use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv 
   use elm_varcon      , only : zsoi, dzsoi, zisoi, spval, namet, grlnd
-  use elm_varcon      , only : secspday, pc, mu, denh2o, denice, grlnd
+  use elm_varcon      , only : secspday, mu, denh2o, denice, grlnd
   use landunit_varcon , only : istice, istdlak, istwet, istsoil, istcrop, istice_mec
   use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
   use elm_varctl      , only : use_cn, use_lch4,use_dynroot, use_fates

diff --git a/components/elm/src/biogeophys/WaterStateType.F90 b/components/elm/src/biogeophys/WaterStateType.F90
@@ -362,7 +362,7 @@ subroutine InitCold(this, bounds, &
     use column_varcon   , only : icol_shadewall, icol_road_perv
     use column_varcon   , only : icol_road_imperv, icol_roof, icol_sunwall
     use elm_varcon      , only : denice, denh2o, spval, sb, bdsno 
-    use elm_varcon      , only : h2osno_max, zlnd, tfrz, spval, pc
+    use elm_varcon      , only : h2osno_max, zlnd, tfrz, spval
     use elm_varctl      , only : fsurdat, iulog
     use spmdMod         , only : masterproc
     use abortutils      , only : endrun
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,2 @@
		fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_nldas2_simyr2000_c181207_surfacewater.nc'
		use_modified_infil = .true.