branching from main, not sure why previous draft pr hadnt

finished adding K(p) and p->theta and doc strings

removed some of the options, added structs

adapted the water retention functions to include vG, linear, or BC, retention params

cleaned up the file, brought minor changes to function signatures

simplified (removed model options). Completed the rhs function updates

changed param type and removed K_sat from parameter struct

Made sure psi(S_l=1)=0. Found very good paper which we should reference in our pplant paper on psi(theta) curve in wood, w complete derivation of potential energy in tree

updated plant hydraulics test

type issue

fixing broadcasting issue

updated Ozark test

debugging broadcasting issue

this works for the first iteration only, for the broadcasting error

propagating changes

runs

fixed a couple bugs, docs

add lwp data

fixed tests

added more unit tests

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