update hydraulic conductivity scaling model

biogeophys/FatesPlantHydraulicsMod.F90

@@ -4596,7 +4596,7 @@ end subroutine shellGeom
 
 ! =====================================================================================
 
-function xylemtaper(pexp, dz) result(chi_tapnotap)
+function xylemtaper_Savage(pexp, dz) result(chi_tapnotap)
 
     use FatesConstantsMod, only : pi => pi_const
 
@@ -4620,7 +4620,7 @@ function xylemtaper(pexp, dz) result(chi_tapnotap)
     !
     ! !LOCAL VARIABLES:
     real(r8) :: qexp                ! total conductance exponent (as in Fig. 2b of Savage et al. (2010) minus a0 term
-    real(r8) :: lN=0.005_r8         ! petiole length[m]
+    real(r8) :: lN=0.15_r8          ! petiole length[m], updated to match the scaling for conduit width from Olson 2020 (https://nph.onlinelibrary.wiley.com/doi/10.1111/nph.16961)
     real(r8) :: n_ext=2._r8         ! number of daughter branches per parent branch, assumed constant throughout tree (self-similarity)  [-]
     real(r8) :: big_n               ! number of branching levels (allowed here to take on non-integer values): increases with tree size  [-]
     real(r8) :: r0rN                ! ratio of stem radius to terminal twig radius; r.ext0/r.extN (x-axis of Savage et al. (2010) Fig 2a)[-]
@@ -4654,8 +4654,51 @@ function xylemtaper(pexp, dz) result(chi_tapnotap)
 
     return
 
-end function xylemtaper
+end function xylemtaper_Savage
+! =====================================================================================
+! xylem hydraulic conductance  following scaling for conduit width from Olson 2020 (https://nph.onlinelibrary.wiley.com/doi/10.1111/nph.16961)
+function xylemtaper(a, L) result(ratio)
+
+  implicit none
+
+  real(r8), intent(in) :: a      ! taper exponent
+  real(r8), intent(in) :: L      ! distance from branch tip (m)
+  real(r8) :: ratio
+
+  real(r8), parameter :: x0  = 0.005_r8  ! finite cutoff (m) at petiole, to prevent singularity in function at L=0.  
+  real(r8), parameter :: eps = 1.e-8_r8
+  real(r8), parameter :: xref = 0.6_r8    ! measurement distance from branch tip (m)
+
+  real(r8) :: Leff
+  real(r8) :: expo
+  real(r8) :: denom
+
+  !------------------------------------------------------------
+  ! Enforce finite cutoff 
+  !------------------------------------------------------------
+  Leff = max(L, x0 * (1._r8 + eps))
+
+  !------------------------------------------------------------
+  ! No taper case
+  !------------------------------------------------------------
+  if (a <= 0._r8) then
+     ratio = 1._r8
+     return
+  end if
+
+  !------------------------------------------------------------
+  ! Critical case: a = 1/4
+  !------------------------------------------------------------
+  if (abs(a - 0.25_r8) < eps) then
+     ratio = (Leff - x0) / (xref * log(Leff / x0))
+  else
+     expo  = 1._r8 - 4._r8 * a
+     denom = xref**(4._r8 * a) * (Leff**expo - x0**expo)
 
+     ratio = (Leff - x0) * expo / denom
+  end if
+
+end function xylemtaper
 ! =====================================================================================
 
 subroutine Hydraulics_Tridiagonal(a, b, c, r, u, ierr)

parameter_files/fates_params_default.json

@@ -721,7 +721,7 @@
       "dims": ["fates_pft"],
       "long_name": "xylem taper exponent",
       "units": "unitless",
-      "data": [0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333, 0.333]
+      "data": [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]
     },
     "fates_hydro_pinot_node": {
       "dtype": "float",