Extend plant hydraulics rooting depth to support fields (#783)

* Change PlantHydraulicsParameters to hold rooting_depth

PlantHydraulicsParameters previously held a root distribution function
that had a default value for rooting depth. This function is
always the same outside testing, so it was removed as a parameter
and placed into src. rooting_depth, an argument to the distribution function
is now a parameters that is either a Float or a field. The tests are changed to
check passing rooting_depth as a float or a field of floats that are the same everywhere.
One test in plant_hydraulics_test.jl only checks the case of when rooting_depth is a float
because changing to a field requires extensive changes to the entire test.

* Update calls to PlantHydraulicsParameters

* Change plant hydraulics integration test to use fields

The test solves a system of functions which is defined in
initial_compute_expected_tendecy!. That function is now changed
to support fields, but each equation in the system of equations still
outputs a scalar. This is done by taking a mean across the resulting field.
This is done because NLsolve does not work of equations on fields.

* Change runs to use spatially varying rooting depth

longruns/land.jl, benchmarks/land.jl and longruns/land_regional.jl
are changed to use the root depth map from the clm data

* Add message in NEWS.md

* Make suggested changes

* Edit root_distribution docstring

* Use new clm artifact in runs

* Change plant_hydraulics_test to follow previous behavior

Previously the Plant hydraulics model integration tests
were changed to work with rooting_depth as a field. In the test,
the solution of a system of equations is found using NLsolve.
The previous change redefined each equation in the system of equations
to create a field of values the equation evaluates to, and then takes
the mean of the field so the equation still results in a scalar.

To better follow the functionality of the previous test, this commit
reverts the system of equations to only work on a scalar rooting_depth.
The solution to the system of equations is returned as a vector of floats.
It is not possible to create a ClimaCore Field of vectors. Instead, the
solution is stored as a vector of fields.

To implement this, a function that returns a function representing
a system of equations is created. This results in a system of equations
for each cell. At first, each system of equations was solved once for each
variable in the vector of solutions, resulting in many repeated equations solves.
To prevent this a dictionary is used. Each variable in the solution must be
returned independently in order to create a vector of fields.

* Make changes not breaking

Add a constructor that finds rooting_depth from
root_distribution

* Remove unnecessary import

* Edit comments for hydraulics test

* Revert to breaking change

A previous commit made the change not breaking
by making constructor accept a rooting_distribtion
parameter, and then calculating rooting_depth from that.
That calculation assumes that rooting_distribuion will be in
a specific form.

* Add docstring for constructor

* Make changes non-breaking (again)

Added support for root_distribution. The PlantHydraulicsParameters
struct has root_distribtion again. Both root_distribution and
root_depth can be nothing. The constructor uses depwarn when
root_distribution is passed in. If neither root_depth or root_distribution
are passed in to the constructor, an error is thrown.

There were two uses of root_distribution. At these spots, a new
function is defined, which takes in the rooting_depth and a z as an argument.
If that rooting_depth argument is nothing, then it calls root_distribution
form the params. Otherwise it calls from src.

Removes breaking change from NEWS.md

* Add deprecated_featured.md

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
 --------

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.

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,
 )

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,
 );

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,
     )

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,
 )

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,
 )

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,
 )

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,
 )

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,
 )

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,
     )

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,
     )

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,
 );

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,
 );

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,
 );

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,
     )

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,