diff --git a/NEWS.md b/NEWS.md
index 3e24e8dce3..9f433e69a9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -10,8 +10,6 @@ main
 v0.15.0
 --------
 - Add regional simulation example PR [#757](https://github.com/CliMA/ClimaLand.jl/pull/757)
-- ![][badge-💥breaking] Extend photosynthesis mechanism parameter to support fields.
-  PR[#774](https://github.com/CliMA/ClimaLand.jl/pull/774)
 - Reduced number of dependencies, leading to faster instantiation and import time,
 - Improved compatibility of ClimaLand with older versions of packages.
   PR[#749](https://github.com/CliMA/ClimaLand.jl/pull/749)
@@ -20,6 +18,9 @@ v0.15.0
   PR[#759](https://github.com/CliMA/ClimaLand.jl/pull/759)
 - Added a regional run example.
   PR[#757](https://github.com/CliMA/ClimaLand.jl/pull/757)
+- ![breaking change][badge-💥breaking] Extend photosynthesis mechanism parameter to support fields.
+PR[#774](https://github.com/CliMA/ClimaLand.jl/pull/774)
+  - C3/C4 structs are removed. Now C3 is represented by a float of 1.0 and C4 by a float of 0.0
 
 v0.14.3
 --------
diff --git a/deprecated_features.md b/deprecated_features.md
new file mode 100644
index 0000000000..0a982943f1
--- /dev/null
+++ b/deprecated_features.md
@@ -0,0 +1,9 @@
+# Deprecated Features List
+
+This file lists deprecated features and the recommended alternative practice
+
+- ## `root_distribution`
+
+    The `root_distribution` in `PlantHydraulicsParameters` is replaced by `rooting_depth`.
+    If using a `root_distribution` function of the form P(x) = 1/d e^(z/d), then `rooting_depth`
+    is equivalent to d.
diff --git a/docs/tutorials/integrated/soil_canopy_tutorial.jl b/docs/tutorials/integrated/soil_canopy_tutorial.jl
index ce6c4acc83..69ef405b9a 100644
--- a/docs/tutorials/integrated/soil_canopy_tutorial.jl
+++ b/docs/tutorials/integrated/soil_canopy_tutorial.jl
@@ -239,9 +239,6 @@ K_sat_plant = FT(1.8e-8)
 RAI = (SAI + maxLAI) * f_root_to_shoot;
 # Note: LAIfunction was determined from data in the script we included above.
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
-function root_distribution(z::T; rooting_depth = FT(1.0)) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-end
 
 ψ63 = FT(-4 / 0.0098)
 Weibull_param = FT(4)
@@ -259,7 +256,7 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = FT(1.0),
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
index ac29d233fe..9d19a68a0e 100644
--- a/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
+++ b/docs/tutorials/standalone/Canopy/canopy_tutorial.jl
@@ -197,11 +197,6 @@ ai_parameterization =
     PrescribedSiteAreaIndex{FT}(TimeVaryingInput(LAIfunction), SAI, RAI)
 rooting_depth = FT(1.0);
 
-# Define the root distribution function p(z):
-
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth))
-end;
 
 # Create the component conductivity and retention models of the hydraulics
 # model. In ClimaLand, a Weibull parameterization is used for the conductivity as
@@ -226,7 +221,7 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = ν,
     S_s = S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 );
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index 07682b3f23..680136e3f7 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -402,7 +402,13 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = FT(0.5) # from Natan
+    rooting_depth = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rooting_depth",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -486,15 +492,11 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-    function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-        return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-    end
-
     plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
         ai_parameterization = ai_parameterization,
         ν = plant_ν,
         S_s = plant_S_s,
-        root_distribution = root_distribution,
+        rooting_depth = rooting_depth,
         conductivity_model = conductivity_model,
         retention_model = retention_model,
     )
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
index 4d7be37565..828abb472e 100644
--- a/experiments/integrated/fluxnet/ozark_pft.jl
+++ b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -190,15 +190,11 @@ photosynthesis_args =
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-end
-
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
index e41095c316..76a55ecec8 100644
--- a/experiments/integrated/fluxnet/run_fluxnet.jl
+++ b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -147,15 +147,11 @@ photosynthesis_args =
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-end
-
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
index 9f6c0589c1..39b473606f 100644
--- a/experiments/integrated/global/global_soil_canopy.jl
+++ b/experiments/integrated/global/global_soil_canopy.jl
@@ -235,15 +235,11 @@ LAIfunction = TimeVaryingInput(
 )
 ai_parameterization = Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-end
-
 plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
index c195862b7b..ecbc9d7b61 100644
--- a/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
+++ b/experiments/integrated/performance/conservation/ozark_conservation_setup.jl
@@ -135,15 +135,11 @@ photosynthesis_args =
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-end
-
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/experiments/integrated/performance/profile_allocations.jl b/experiments/integrated/performance/profile_allocations.jl
index 2ec19f9192..14a20b431e 100644
--- a/experiments/integrated/performance/profile_allocations.jl
+++ b/experiments/integrated/performance/profile_allocations.jl
@@ -271,14 +271,11 @@ photosynthesis_args =
     (; parameters = FarquharParameters(FT, is_c3; Vcmax25 = Vcmax25))
 # Set up plant hydraulics
 ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
-function root_distribution(z::T) where {T}
-    return T(1.0 / 0.5 * exp(z / 0.5)) # 1/m
-end
 plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = plant_ν,
     S_s = plant_S_s,
-    root_distribution = root_distribution,
+    rooting_depth = FT(0.5),
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 )
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index 060f0af4bf..dc965ed1f9 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -410,7 +410,13 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = FT(0.5) # from Natan
+    rooting_depth = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rooting_depth",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -496,15 +502,11 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-    function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-        return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-    end
-
     plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
         ai_parameterization = ai_parameterization,
         ν = plant_ν,
         S_s = plant_S_s,
-        root_distribution = root_distribution,
+        rooting_depth = rooting_depth,
         conductivity_model = conductivity_model,
         retention_model = retention_model,
     )
diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
index 511214f8c5..3435e0e4a5 100644
--- a/experiments/long_runs/land_region.jl
+++ b/experiments/long_runs/land_region.jl
@@ -412,7 +412,13 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
     plant_ν = FT(1.44e-4)
     plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
-    rooting_depth = FT(0.5) # from Natan
+    rooting_depth = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "rooting_depth",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
     n_stem = 0
     n_leaf = 1
     h_stem = FT(0.0)
@@ -498,15 +504,11 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 
-    function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-        return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-    end
-
     plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
         ai_parameterization = ai_parameterization,
         ν = plant_ν,
         S_s = plant_S_s,
-        root_distribution = root_distribution,
+        rooting_depth = rooting_depth,
         conductivity_model = conductivity_model,
         retention_model = retention_model,
     )
diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl
index b71240098e..ee068db32c 100644
--- a/experiments/standalone/Vegetation/no_vegetation.jl
+++ b/experiments/standalone/Vegetation/no_vegetation.jl
@@ -108,10 +108,6 @@ ai_parameterization =
     PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction), SAI, RAI)
 rooting_depth = FT(1.0);
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth))
-end;
-
 K_sat_plant = FT(1.8e-8)
 ψ63 = FT(-4 / 0.0098)
 Weibull_param = FT(4)
@@ -129,7 +125,7 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = ν,
     S_s = S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 );
diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl
index ecf9dbec3a..c8081e8f10 100644
--- a/experiments/standalone/Vegetation/varying_lai.jl
+++ b/experiments/standalone/Vegetation/varying_lai.jl
@@ -108,10 +108,6 @@ ai_parameterization =
     PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction), SAI, RAI)
 rooting_depth = FT(1.0);
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth))
-end;
-
 K_sat_plant = FT(1.8e-8)
 ψ63 = FT(-4 / 0.0098)
 Weibull_param = FT(4)
@@ -129,7 +125,7 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = ν,
     S_s = S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 );
diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
index d1068db1ec..bdb9d86ee7 100644
--- a/experiments/standalone/Vegetation/varying_lai_with_stem.jl
+++ b/experiments/standalone/Vegetation/varying_lai_with_stem.jl
@@ -107,10 +107,6 @@ ai_parameterization =
     PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction2), SAI, RAI)
 rooting_depth = FT(1.0);
 
-function root_distribution(z::T; rooting_depth = rooting_depth) where {T}
-    return T(1.0 / rooting_depth) * exp(z / T(rooting_depth))
-end;
-
 K_sat_plant = FT(1.8e-6)
 ψ63 = FT(-4 / 0.0098)
 Weibull_param = FT(4)
@@ -128,7 +124,7 @@ plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
     ai_parameterization = ai_parameterization,
     ν = ν,
     S_s = S_s,
-    root_distribution = root_distribution,
+    rooting_depth = rooting_depth,
     conductivity_model = conductivity_model,
     retention_model = retention_model,
 );
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
index e584770469..f9458afdd1 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
@@ -154,18 +154,12 @@ function run_fluxnet(
         drivers.RAI,
     )
 
-    function root_distribution(
-        z::T;
-        rooting_depth = params.plant_hydraulics.rooting_depth,
-    ) where {T}
-        return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) # 1/m
-    end
 
     plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
         ai_parameterization = ai_parameterization,
         ν = drivers.plant_ν,
         S_s = params.plant_hydraulics.S_s,
-        root_distribution = root_distribution,
+        rooting_depth = params.plant_hydraulics.rooting_depth,
         conductivity_model = params.plant_hydraulics.conductivity_model,
         retention_model = params.plant_hydraulics.retention_model,
     )
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index 1f2b900181..085715912a 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -224,7 +224,7 @@ lsm_aux_domain_names(m::SoilCanopyModel) =
 A method which makes a function; the returned function
 updates the additional auxiliary variables for the integrated model,
 as well as updates the boundary auxiliary variables for all component
-models. 
+models.
 
 This function is called each ode function evaluation, prior to the tendency function
 evaluation.
@@ -261,7 +261,11 @@ function make_update_boundary_fluxes(
         # Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
         # and K_canopy = K_plant(ψ_canopy). In `PrognosticSoil` mode here, we compute the flux using
         # K_soil = K_soil(ψ_soil) and K_canopy = K_plant(ψ_canopy).
-
+        # if rooting_depth param is not nothing, use root_distribution from source
+        # otherwise use root_distribution from params
+        root_likelihood(z::FT, rd::Union{FT, Nothing}) =
+            !isnothing(rd) ? root_distribution(z, rd) :
+            model.parameters.root_distribution(z)
         @. p.root_extraction =
             above_ground_area_index *
             PlantHydraulics.water_flux(
@@ -275,7 +279,10 @@ function make_update_boundary_fluxes(
                     p.canopy.hydraulics.ψ.:1,
                 ),
             ) *
-            (land.canopy.hydraulics.parameters.root_distribution(z))
+            (root_likelihood(
+                z,
+                land.canopy.hydraulics.parameters.rooting_depth,
+            ))
         @. p.root_energy_extraction =
             p.root_extraction * ClimaLand.Soil.volumetric_internal_energy_liq(
                 p.soil.T,
diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl
index 579692a322..930513a60c 100644
--- a/src/standalone/Vegetation/PlantHydraulics.jl
+++ b/src/standalone/Vegetation/PlantHydraulics.jl
@@ -37,7 +37,8 @@ export PlantHydraulicsModel,
     LinearRetentionCurve,
     Weibull,
     hydraulic_conductivity,
-    PrescribedSiteAreaIndex
+    PrescribedSiteAreaIndex,
+    root_distribution
 
 """
     AbstractPlantHydraulicsModel{FT} <: AbstractCanopyComponent{FT}
@@ -136,7 +137,8 @@ struct PlantHydraulicsParameters{
     PSAI <: PrescribedSiteAreaIndex{FT},
     CP,
     RP,
-    F <: Function,
+    RDTH <: Union{FT, ClimaCore.Fields.Field, Nothing},
+    RDST <: Union{Function, Nothing},
 }
     "The area index model for LAI, SAI, RAI"
     ai_parameterization::PSAI
@@ -148,23 +150,56 @@ struct PlantHydraulicsParameters{
     conductivity_model::CP
     "Water retention model and parameters"
     retention_model::RP
-    "Root distribution function P(z)"
-    root_distribution::F
+    "Rooting depth parameter (m) - a characteristic depth below which 1/e of the root mass lies"
+    rooting_depth::RDTH
+    "DEPRECATED: Root distribution function P(z)"
+    root_distribution::RDST
 end
 
+"""
+    PlantHydraulicsParameters(;
+        ai_parameterization::PrescribedSiteAreaIndex{FT},
+        ν::FT,
+        S_s::FT,
+        conductivity_model,
+        retention_model,
+        rooting_depth::Union{Nothing, FT, ClimaCore.Fields.Field} = nothing,
+        root_distribution::Union{Nothing, Function} = nothing, # DEPRECATED
+    )
+
+Constructor for PlantHydraulicsParameters. The root_distribution parameter is deprecated.
+If root_distribution and rooting_depth are both provided, root_distribution is ignored.
+"""
 function PlantHydraulicsParameters(;
     ai_parameterization::PrescribedSiteAreaIndex{FT},
     ν::FT,
     S_s::FT,
-    root_distribution::Function,
     conductivity_model,
     retention_model,
+    rooting_depth::Union{Nothing, FT, ClimaCore.Fields.Field} = nothing,
+    root_distribution::Union{Nothing, Function} = nothing, # DEPRECATED
 ) where {FT}
+    if !isnothing(root_distribution)
+        if !isnothing(rooting_depth)
+            Base.depwarn(
+                "root_distribution is deprecated and will be ignored when rooting_depth is provided",
+                :PlantHydraulicsParameters,
+            )
+        else
+            Base.depwarn(
+                "root_distribution keyword argument is deprecated. Use rooting_depth instead.",
+                :PlantHydraulicsParameters,
+            )
+        end
+    elseif isnothing(rooting_depth)
+        error("rooting_depth must be provided")
+    end
     return PlantHydraulicsParameters{
         FT,
         typeof(ai_parameterization),
         typeof(conductivity_model),
         typeof(retention_model),
+        typeof(rooting_depth),
         typeof(root_distribution),
     }(
         ai_parameterization,
@@ -172,6 +207,7 @@ function PlantHydraulicsParameters(;
         S_s,
         conductivity_model,
         retention_model,
+        rooting_depth,
         root_distribution,
     )
 end
@@ -358,7 +394,7 @@ mean for effective conducticity between the two layers
 
 To account for different path lengths in the two compartments Δz1 and
 Δz2, we would require the following conductance k (1/s)
-k_eff = K1/Δz1*K2/Δz2/(K1/Δz1+K2/Δz2) 
+k_eff = K1/Δz1*K2/Δz2/(K1/Δz1+K2/Δz2)
 and a water flux of
 F = -k_eff * (ψ1 +z1 - ψ2 - z2) (m/s).
 
@@ -370,6 +406,23 @@ function water_flux(z1::FT, z2::FT, ψ1::FT, ψ2::FT, K1::FT, K2::FT) where {FT}
     return flux # (m/s)
 end
 
+"""
+    root_distribution(z::FT, rooting_depth::FT)
+
+Computes value of rooting probability density function at `z`.
+
+The rooting probability density function is derived from the
+cumulative distribution function F(z) = 1 - β^(100z), which is described
+by Equation 2.23 of
+Bonan, "Climate Change and Terrestrial Ecosystem Modeling", 2019 Cambridge University Press.
+This probability distribution function is equivalent to the derivative of the
+cumulative distribution function with respect to z,
+where `rooting_depth` replaces (-1)/(100ln(β)) and z is expected to be negative.
+"""
+function root_distribution(z::FT, rooting_depth::FT) where {FT <: AbstractFloat}
+    return (1 / rooting_depth) * exp(z / rooting_depth) # 1/m
+end
+
 """
     AbstractConductivityModel{FT <: AbstractFloat}
 
@@ -621,7 +674,7 @@ function root_water_flux_per_ground_area!(
     t,
 ) where {FT}
 
-    (; conductivity_model, root_distribution) = model.parameters
+    (; conductivity_model, rooting_depth) = model.parameters
     area_index = p.canopy.hydraulics.area_index
     # We can index into a field of Tuple{FT} to extract a field of FT
     # using the following notation: field.:index
@@ -630,6 +683,11 @@ function root_water_flux_per_ground_area!(
     n_root_layers = length(root_depths)
     ψ_soil::FT = s.ψ(t)
     fa .= FT(0.0)
+    # if rooting_depth param is not nothing, use root_distribution from source
+    # otherwise use root_distribution from params
+    root_likelihood(z::FT, rd::Union{FT, Nothing}) =
+        !isnothing(rd) ? root_distribution(z, rd) :
+        model.parameters.root_distribution(z)
     @inbounds for i in 1:n_root_layers
         above_ground_area_index =
             getproperty(area_index, model.compartment_labels[1])
@@ -644,7 +702,7 @@ function root_water_flux_per_ground_area!(
                     hydraulic_conductivity(conductivity_model, ψ_soil),
                     hydraulic_conductivity(conductivity_model, ψ_base),
                 ) *
-                root_distribution(root_depths[i]) *
+                root_likelihood(root_depths[i], rooting_depth) *
                 (root_depths[i + 1] - root_depths[i]) *
                 above_ground_area_index
         else
@@ -657,7 +715,7 @@ function root_water_flux_per_ground_area!(
                     hydraulic_conductivity(conductivity_model, ψ_soil),
                     hydraulic_conductivity(conductivity_model, ψ_base),
                 ) *
-                root_distribution(root_depths[i]) *
+                root_likelihood(root_depths[i], rooting_depth) *
                 (model.compartment_surfaces[1] - root_depths[n_root_layers]) *
                 above_ground_area_index
         end
diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl
index b00bec43e6..a7b5259e0f 100644
--- a/test/standalone/Vegetation/canopy_model.jl
+++ b/test/standalone/Vegetation/canopy_model.jl
@@ -439,14 +439,11 @@ end
             PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
         retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
         root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth
-        function root_distribution(z::T) where {T}
-            return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m)
-        end
         param_set = PlantHydraulics.PlantHydraulicsParameters(;
             ai_parameterization = ai_parameterization,
             ν = plant_ν,
             S_s = plant_S_s,
-            root_distribution = root_distribution,
+            rooting_depth = FT(0.5),
             conductivity_model = conductivity_model,
             retention_model = retention_model,
         )
@@ -571,8 +568,9 @@ end
         g1_cases = (FT(790), fill(FT(790), domain.space.surface))
         Vcmax25_cases = (FT(9e-5), fill(FT(9e-5), domain.space.surface))
         mechanism_cases = (FT(1), mechanism_field)
-        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases)
-        for (g1, Vcmax25, is_c3) in zipped
+        rooting_cases = (FT(0.5), fill(FT(0.5), domain.space.surface))
+        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases, rooting_cases)
+        for (g1, Vcmax25, is_c3, rooting_depth) in zipped
             RTparams = BeerLambertParameters(FT)
             photosynthesis_params = FarquharParameters(FT, is_c3; Vcmax25)
             stomatal_g_params = MedlynConductanceParameters(FT; g1)
@@ -682,14 +680,11 @@ end
             retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
             root_depths =
                 SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth
-            function root_distribution(z::T) where {T}
-                return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m)
-            end
             param_set = PlantHydraulics.PlantHydraulicsParameters(;
                 ai_parameterization = ai_parameterization,
                 ν = plant_ν,
                 S_s = plant_S_s,
-                root_distribution = root_distribution,
+                rooting_depth = rooting_depth,
                 conductivity_model = conductivity_model,
                 retention_model = retention_model,
             )
@@ -837,8 +832,9 @@ end
         g1_cases = (FT(790), fill(FT(790), domain.space.surface))
         Vcmax25_cases = (FT(9e-5), fill(FT(9e-5), domain.space.surface))
         mechanism_cases = (FT(1), mechanism_field)
-        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases)
-        for (g1, Vcmax25, is_c3) in zipped
+        rooting_cases = (FT(0.5), fill(FT(0.5), domain.space.surface))
+        zipped = zip(g1_cases, Vcmax25_cases, mechanism_cases, rooting_cases)
+        for (g1, Vcmax25, is_c3, rooting_depth) in zipped
             BeerLambertparams = BeerLambertParameters(FT)
             # TwoStreamModel parameters
             Ω = FT(0.69)
@@ -971,14 +967,11 @@ end
             retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
             root_depths =
                 SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth
-            function root_distribution(z::T) where {T}
-                return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m)
-            end
             param_set = PlantHydraulics.PlantHydraulicsParameters(;
                 ai_parameterization = ai_parameterization,
                 ν = plant_ν,
                 S_s = plant_S_s,
-                root_distribution = root_distribution,
+                rooting_depth = rooting_depth,
                 conductivity_model = conductivity_model,
                 retention_model = retention_model,
             )
diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl
index 481b78770f..2930059065 100644
--- a/test/standalone/Vegetation/plant_hydraulics_test.jl
+++ b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -226,9 +226,6 @@ for FT in (Float32, Float64)
             PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
         retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
         root_depths = FT.(-Array{FT}(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth
-        function root_distribution(z::T) where {T}
-            return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m)
-        end
         compartment_midpoints = Vector{FT}(
             range(
                 start = Δz / 2,
@@ -241,15 +238,6 @@ for FT in (Float32, Float64)
             range(start = 0.0, step = Δz, stop = Δz * (n_stem + n_leaf)),
         )
 
-        param_set = PlantHydraulics.PlantHydraulicsParameters(;
-            ai_parameterization = ai_parameterization,
-            ν = plant_ν,
-            S_s = plant_S_s,
-            root_distribution = root_distribution,
-            conductivity_model = conductivity_model,
-            retention_model = retention_model,
-        )
-
         function leaf_transpiration(t)
             T = FT(1e-8) # m/s
         end
@@ -259,18 +247,28 @@ for FT in (Float32, Float64)
 
         soil_driver = PrescribedSoil(FT)
 
-        plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(;
-            parameters = param_set,
-            transpiration = transpiration,
-            n_stem = n_stem,
-            n_leaf = n_leaf,
-            compartment_surfaces = compartment_surfaces,
-            compartment_midpoints = compartment_midpoints,
-        )
         autotrophic_parameters = AutotrophicRespirationParameters(FT)
         autotrophic_respiration_model =
             AutotrophicRespirationModel{FT}(autotrophic_parameters)
         for domain in domains
+            param_set = PlantHydraulics.PlantHydraulicsParameters(;
+                ai_parameterization = ai_parameterization,
+                ν = plant_ν,
+                S_s = plant_S_s,
+                rooting_depth = fill(FT(0.5), domain.space.surface),
+                conductivity_model = conductivity_model,
+                retention_model = retention_model,
+            )
+
+            plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(;
+                parameters = param_set,
+                transpiration = transpiration,
+                n_stem = n_stem,
+                n_leaf = n_leaf,
+                compartment_surfaces = compartment_surfaces,
+                compartment_midpoints = compartment_midpoints,
+            )
+
             model = ClimaLand.Canopy.CanopyModel{FT}(;
                 parameters = shared_params,
                 domain = domain,
@@ -283,70 +281,110 @@ for FT in (Float32, Float64)
                 atmos = atmos,
                 radiation = radiation,
             )
-            # Set system to hydrostatic equilibrium state by setting fluxes to zero, and setting LHS of both ODEs to 0
-            function initial_compute_exp_tendency!(F, Y)
-                AI = (; leaf = LAI(1.0), root = RAI, stem = SAI)
-                T0A = FT(1e-8) * AI[:leaf]
-                for i in 1:(n_leaf + n_stem)
-                    if i == 1
-                        fa =
-                            sum(
-                                water_flux.(
-                                    root_depths,
+            # Create a system of equations for a given rooting_depth for each cell
+            function create_exp_tendency_func(rooting_depth)
+                # Solve for hydrostatic equilibrium
+                return function initial_compute_exp_tendency!(F, Y)
+                    AI = (; leaf = LAI(1.0), root = RAI, stem = SAI)
+                    T0A = FT(1e-8) * AI[:leaf]
+                    for i in 1:(n_leaf + n_stem)
+                        if i == 1
+                            fa =
+                                sum(
+                                    water_flux.(
+                                        root_depths,
+                                        plant_hydraulics.compartment_midpoints[i],
+                                        ψ_soil0,
+                                        Y[i],
+                                        PlantHydraulics.hydraulic_conductivity(
+                                            conductivity_model,
+                                            ψ_soil0,
+                                        ),
+                                        PlantHydraulics.hydraulic_conductivity(
+                                            conductivity_model,
+                                            Y[i],
+                                        ),
+                                    ) .*
+                                    ClimaLand.Canopy.PlantHydraulics.root_distribution.(
+                                        root_depths,
+                                        rooting_depth,
+                                    ) .* (
+                                        vcat(root_depths, [0.0])[2:end] -
+                                        vcat(root_depths, [0.0])[1:(end - 1)]
+                                    ),
+                                ) * AI[:stem]
+                        else
+                            fa =
+                                water_flux(
+                                    plant_hydraulics.compartment_midpoints[i - 1],
                                     plant_hydraulics.compartment_midpoints[i],
-                                    ψ_soil0,
+                                    Y[i - 1],
                                     Y[i],
                                     PlantHydraulics.hydraulic_conductivity(
                                         conductivity_model,
-                                        ψ_soil0,
+                                        Y[i - 1],
                                     ),
                                     PlantHydraulics.hydraulic_conductivity(
                                         conductivity_model,
                                         Y[i],
                                     ),
-                                ) .* root_distribution.(root_depths) .* (
-                                    vcat(root_depths, [0.0])[2:end] -
-                                    vcat(root_depths, [0.0])[1:(end - 1)]
-                                ),
-                            ) * AI[:stem]
-                    else
-                        fa =
-                            water_flux(
-                                plant_hydraulics.compartment_midpoints[i - 1],
-                                plant_hydraulics.compartment_midpoints[i],
-                                Y[i - 1],
-                                Y[i],
-                                PlantHydraulics.hydraulic_conductivity(
-                                    conductivity_model,
-                                    Y[i - 1],
-                                ),
-                                PlantHydraulics.hydraulic_conductivity(
-                                    conductivity_model,
-                                    Y[i],
-                                ),
-                            ) * AI[plant_hydraulics.compartment_labels[i]]
+                                ) * AI[plant_hydraulics.compartment_labels[i]]
+                        end
+                        F[i] = fa - T0A
                     end
-                    F[i] = fa - T0A
                 end
             end
+            #=======================
+            nlsolve does not solve systems of equations that return fields, so each cell
+            must have its own system of equations that is solved for. ClimaCore fields
+            cannot hold vectors, so multiple fields must be used to hold the vector solution
+            to the system of equations. Here we solve for the steady state of the hydraulics
+            system. Then, the solution is used to check that evaluating the
+            tendecy of the model also results in a steady state. This check is repeated using
+            the plant hydraulics model directly.
+            =======================#
+            # dict to prevent recalculation when returning different index of solution
+            solutions_for_rooting_depth = Dict{FT, Vector{FT}}()
+            # solves system of equations at each cell and returns specific index of solution
+            function solve_sys_for_i(rooting_depth::FT, i) where {FT}
+                if haskey(solutions_for_rooting_depth, rooting_depth)
+                    return get(
+                        solutions_for_rooting_depth,
+                        rooting_depth,
+                        FT(i),
+                    )[i]
+                else
+                    soln = nlsolve(
+                        create_exp_tendency_func(rooting_depth),
+                        Vector{FT}(-0.03:0.01:0.07);
+                        ftol = eps(FT),
+                        method = :newton,
+                        iterations = 20,
+                    )
+                    return get!(
+                        solutions_for_rooting_depth,
+                        rooting_depth,
+                        soln.zero,
+                    )[i]
+                end
+            end
+            # instead of a field of vectors, we have a vector of fields
+            soln_field_vector = map(1:(n_stem + n_leaf)) do i
+                solve_sys_for_i.(plant_hydraulics.parameters.rooting_depth, i)
+            end
 
-            soln = nlsolve(
-                initial_compute_exp_tendency!,
-                Vector{FT}(-0.03:0.01:0.07);
-                ftol = eps(FT),
-                method = :newton,
-                iterations = 20,
-            )
-
-            S_l =
+            S_l = map(soln_field_vector) do soln_i_field
                 inverse_water_retention_curve.(
                     retention_model,
-                    soln.zero,
+                    soln_i_field,
                     plant_ν,
                     plant_S_s,
                 )
+            end
 
-            ϑ_l_0 = augmented_liquid_fraction.(plant_ν, S_l)
+            ϑ_l_0 = map(S_l) do S_l_i
+                augmented_liquid_fraction.(plant_ν, S_l_i)
+            end
 
             Y, p, coords = initialize(model)
             if typeof(domain) <: ClimaLand.Domains.Point
@@ -516,17 +554,14 @@ for FT in (Float32, Float64)
             PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
         retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
         root_depths = [FT(0.0)]# 1st element is the deepest root depth
-        function root_distribution(z::T) where {T}
-            return T(0) # (1/m)
-        end
         compartment_midpoints = [h_canopy]
         compartment_surfaces = [FT(0.0), h_canopy]
-
+        # set rooting_depth param to largest possible value to test no roots
         param_set = PlantHydraulics.PlantHydraulicsParameters(;
             ai_parameterization = ai_parameterization,
             ν = plant_ν,
             S_s = plant_S_s,
-            root_distribution = root_distribution,
+            rooting_depth = maxintfloat(FT),
             conductivity_model = conductivity_model,
             retention_model = retention_model,
         )