Merge #185

185: Including new vulnerability and p to theta curves  r=kmdeck a=gagnelandmanna

## Purpose 
The purpose of this PR is to change the vulnerability curve (K(P)) and water retention curve (to convert between P and theta) used in plant hydraulics, as well as to add in the infrastructure for adding new models for this/allowing for modularity.

#159 -- this will link to issue 159

Review checklist

I have:
- followed the codebase contribution guide: https://clima.github.io/ClimateMachine.jl/latest/Contributing/
- followed the style guide: https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/
- followed the documentation policy: https://github.com/CliMA/policies/wiki/Documentation-Policy
- checked that this PR does not duplicate an open PR.

In the Content, I have included 
- relevant unit tests, and integration tests, 
- appropriate docstrings on all functions, structs, and modules, and included relevant documentation.

----
- [x] I have read and checked the items on the review checklist.

Co-authored-by: Anna GL <gagnelandmanna@gmail.com>

docs/src/APIs/canopy/PlantHydraulics.md

@@ -18,12 +18,15 @@ ClimaLSM.PlantHydraulics.water_retention_curve
 ClimaLSM.PlantHydraulics.inverse_water_retention_curve
 ClimaLSM.PlantHydraulics.root_flux_per_ground_area!
 ClimaLSM.PlantHydraulics.flux
+ClimaLSM.PlantHydraulics.hydraulic_conductivity
 ```
 
 ## Plant Hydraulics Parameters
 
 ```@docs
 ClimaLSM.PlantHydraulics.PlantHydraulicsParameters
+ClimaLSM.PlantHydraulics.Weibull
+ClimaLSM.PlantHydraulics.LinearRetentionCurve
 ```
 
 ## Plant Hydraulics Methods and Types

experiments/LSM/ozark/ozark.jl

@@ -97,24 +97,18 @@ photosynthesis_args = (;
 )
 # Set up plant hydraulics
 area_index = (root = RAI, stem = SAI, leaf = LAI)
-K_sat_root = FT(K_sat_plant) # m/s
-K_sat_stem = FT(K_sat_plant)
-K_sat_leaf = FT(K_sat_plant)
-K_sat = (root = K_sat_root, stem = K_sat_stem, leaf = K_sat_leaf)
 
 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{FT}(
-    area_index,
-    K_sat,
-    plant_vg_α,
-    plant_vg_n,
-    plant_vg_m,
-    plant_ν,
-    plant_S_s,
-    root_distribution,
+plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(;
+    area_index = area_index,
+    ν = plant_ν,
+    S_s = plant_S_s,
+    root_distribution = root_distribution,
+    conductivity_model = conductivity_model,
+    retention_model = retention_model,
 )
 
 plant_hydraulics_args = (
@@ -166,9 +160,7 @@ Y.soil.ϑ_l = FT(0.4)
 
 S_l_ini =
     inverse_water_retention_curve.(
-        plant_vg_α,
-        plant_vg_n,
-        plant_vg_m,
+        retention_model,
         [ψ_stem_0, ψ_leaf_0],
         plant_ν,
         plant_S_s,
@@ -344,6 +336,46 @@ plt2 = Plots.plot(
 Plots.plot(plt2, plt1, layout = (2, 1))
 Plots.savefig(joinpath(savedir, "soil_water_content.png"))
 
+dt_model = sol.t[2] - sol.t[1]
+dt_data = seconds[2] - seconds[1]
+Plots.plot(daily, cumsum(T) * dt_model, label = "Model T")
+Plots.plot!(
+    seconds ./ 24 ./ 3600,
+    cumsum(measured_T[:]) * dt_data,
+    label = "Data ET",
+)
+Plots.plot!(
+    seconds ./ 24 ./ 3600,
+    cumsum(P[:]) * dt_data * (1e3 * 24 * 3600),
+    label = "Data P",
+)
+Plots.plot!(ylabel = "∫ Water fluxes dt", xlabel = "Days", margins = 10Plots.mm)
+Plots.savefig(joinpath(savedir, "cumul_p_et.png"))
+
+# Leaf water potential data from Pallardy et al (2018)
+# Predawn Leaf Water Potential of Oak-Hickory Forest at Missouri Ozark (MOFLUX) Site: 2004-2020
+# https://doi.org/10.3334/CDIAC/ORNLSFA.004
+lwp_filename = "MOFLUX_PredawnLeafWaterPotential_2020_20210125.csv"
+lwp_artifact = ArtifactFile(
+    url = "https://caltech.box.com/shared/static/d2nbhezw1q99vslnh5qfwnqrnp3p4edo.csv",
+    filename = lwp_filename,
+)
+lwp_dataset = ArtifactWrapper(
+    @__DIR__,
+    "lwp_pallardy_etal2018",
+    ArtifactFile[lwp_artifact],
+);
+
+lwp_path = joinpath(get_data_folder(lwp_dataset), lwp_filename)
+lwp_data = readdlm(lwp_path, ',', skipstart = 1)
+# We are using 2005 data in this test, so restrict to this year
+YEAR = lwp_data[:, 1]
+DOY = lwp_data[YEAR .== 2005, 2]
+# Our t0 = Dec 31, midnight, 2004. Predawn = guess of 0600 hours
+seconds_since_t0 = FT.(DOY) * 24 .* 3600 .+ (6 * 3600)
+lwp_measured = lwp_data[YEAR .== 2005, 7] .* 1e6 # MPa to Pa
+
+
 root_stem_flux = [
     sum(sv.saveval[k].root_extraction) .* (1e3 * 3600 * 24) for
     k in 1:length(sol.t)
@@ -361,15 +393,15 @@ leaf_air_flux = [
 
 plt1 = Plots.plot(
     daily,
-    root_stem_flux,
-    label = "Soil-root-stem flux",
+    leaf_air_flux,
+    label = "Leaf-air flux",
     ylabel = "Within plant fluxes[mm/day]",
     xlim = [minimum(daily), maximum(daily)],
     size = (500, 700),
     margins = 10Plots.mm,
 )
 Plots.plot!(plt1, daily, stem_leaf_flux, label = "Stem-leaf flux")
-Plots.plot!(plt1, daily, leaf_air_flux, label = "Leaf-air flux")
+Plots.plot!(plt1, daily, root_stem_flux, label = "Soil-root-stem flux")
 
 lwp = [
     parent(sv.saveval[k].canopy.hydraulics.ψ)[2] * 9800 for k in 1:length(sol.t)
@@ -381,35 +413,24 @@ swp = [
 plt2 = Plots.plot(
     daily,
     lwp,
-    label = "Leaf",
+    label = "Model, Leaf",
     xlim = [minimum(daily), maximum(daily)],
     xlabel = "days",
     ylabel = "Water potential (Pa)",
     size = (500, 700),
     left_margin = 10Plots.mm,
 )
+Plots.plot!(plt2, daily, swp, label = "Model, Stem")
 
-Plots.plot!(plt2, daily, swp, label = "Stem")
-
-Plots.plot(plt2, plt1, layout = (2, 1))
-Plots.savefig(joinpath(savedir, "plant_hydraulics.png"))
-
-
-dt_model = sol.t[2] - sol.t[1]
-dt_data = seconds[2] - seconds[1]
-Plots.plot(daily, cumsum(T) * dt_model, label = "Model T")
-Plots.plot!(
-    seconds ./ 24 ./ 3600,
-    cumsum(measured_T[:]) * dt_data,
-    label = "Data ET",
-)
 Plots.plot!(
-    seconds ./ 24 ./ 3600,
-    cumsum(P[:]) * dt_data * (1e3 * 24 * 3600),
-    label = "Data P",
+    plt2,
+    seconds_since_t0 ./ 24 ./ 3600,
+    lwp_measured,
+    label = "Data; all species",
 )
-Plots.plot!(ylabel = "∫ Water fluxes dt", xlabel = "Days", margins = 10Plots.mm)
-Plots.savefig(joinpath(savedir, "cumul_p_et.png"))
 
 
+Plots.plot(plt2, plt1, layout = (2, 1))
+Plots.savefig(joinpath(savedir, "plant_hydraulics.png"))
+
 rm(joinpath(savedir, "Artifacts.toml"))

experiments/LSM/ozark/ozark_parameters.jl

@@ -42,12 +42,15 @@ SAI = FT(0.00242) # m2/m2
 LAI = FT(4.2) # m2/m2, from Wang et al.
 f_root_to_shoot = FT(3.5)
 RAI = (SAI + LAI) * f_root_to_shoot
-K_sat_plant = 1.8e-8 # m/s, guess, close to Natan
-plant_vg_α = FT(0.002) # 1/m, matches Natan (fitted vG to Weibull)
-plant_vg_n = FT(4.2) # unitless, matches Natan (fitted vG to Weibull)
-plant_vg_m = FT(1) - FT(1) / plant_vg_n
+K_sat_plant = 1.8e-8 # m/s
+ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value
+Weibull_param = FT(4) # unitless, Holtzman's original c param value
+a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
+conductivity_model =
+    PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
+retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)
 plant_ν = FT(0.7) # guess, m3/m3
 plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m
 rooting_depth = FT(1.0) # from Wang et al.
 z0_m = FT(2)
-z0_b = FT(0.1)
+z0_b = FT(0.2)

src/Vegetation/Canopy.jl

@@ -383,11 +383,10 @@ function ClimaLSM.make_update_aux(
         hydraulics = canopy.hydraulics
         n_stem = hydraulics.n_stem
         n_leaf = hydraulics.n_leaf
-        (; vg_α, vg_n, vg_m, S_s, ν, K_sat, area_index) = hydraulics.parameters
+        (; retention_model, conductivity_model, S_s, ν, area_index) =
+            hydraulics.parameters
         @inbounds @. ψ[1] = PlantHydraulics.water_retention_curve(
-            vg_α,
-            vg_n,
-            vg_m,
+            retention_model,
             PlantHydraulics.effective_saturation(ν, ϑ_l[1]),
             ν,
             S_s,
@@ -396,9 +395,7 @@ function ClimaLSM.make_update_aux(
         @inbounds for i in 1:(n_stem + n_leaf - 1)
 
             @. ψ[i + 1] = PlantHydraulics.water_retention_curve(
-                vg_α,
-                vg_n,
-                vg_m,
+                retention_model,
                 PlantHydraulics.effective_saturation(ν, ϑ_l[i + 1]),
                 ν,
                 S_s,
@@ -411,13 +408,14 @@ function ClimaLSM.make_update_aux(
                     hydraulics.compartment_midpoints[i + 1],
                     ψ[i],
                     ψ[i + 1],
-                    vg_α,
-                    vg_n,
-                    vg_m,
-                    ν,
-                    S_s,
-                    K_sat[hydraulics.compartment_labels[i]],
-                    K_sat[hydraulics.compartment_labels[i + 1]],
+                    PlantHydraulics.hydraulic_conductivity(
+                        conductivity_model,
+                        ψ[i],
+                    ),
+                    PlantHydraulics.hydraulic_conductivity(
+                        conductivity_model,
+                        ψ[i + 1],
+                    ),
                 ) * (
                     area_index[hydraulics.compartment_labels[i]] +
                     area_index[hydraulics.compartment_labels[i + 1]]