pass CI [skip ci]

CliMA · Sep 10, 2024 · 7cbc3fb · 7cbc3fb
1 parent 1651269
commit 7cbc3fb
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 31 deletions.
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
@@ -263,6 +263,13 @@ function make_update_boundary_fluxes(
             area_index,
             land.canopy.hydraulics.compartment_labels[1],
         )
+        # Note that we model the flux between each soil layer and the canopy as:
+        # Flux = -K_eff x [(ψ_canopy - ψ_soil)/(z_canopy - z_soil) + 1], where
+        # K_eff = K_soil K_canopy /(K_canopy + K_soil)
+
+        # 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).
 
         @. p.root_extraction =
             above_ground_area_index *

diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl
@@ -344,12 +344,12 @@ the flow path between these two points.
 
 We currently assuming an geometric
 mean for mean K_sat between the two points. Previously,
-    we used the arithmetic mean (Bonan, 2019; Zhu, 2008),
-    but then when the soil K was very low, root extraction would
-    continue. This should be modified for compartments of differing sizes.
+we used the arithmetic mean (Bonan, 2019; Zhu, 2008),
+but then when the soil K was very low, root extraction would
+continue. This should be modified for compartments of differing sizes.
 """
-function flux(z1, z2, ψ1, ψ2, K1, K2)
-    K_eff = K1 * K2 / (K1 + K2)
+function flux(z1::FT, z2::FT, ψ1::FT, ψ2::FT, K1::FT, K2::FT) where {FT}
+    K_eff = K1 * K2 / max(K1 + K2, eps(FT))
     flux = -K_eff * ((ψ2 - ψ1) / (z2 - z1) + 1)
     return flux # (m/s)
 end
@@ -582,8 +582,16 @@ end
     ) where {FT}
 
 A method which computes the water flux between the soil and the stem, via the roots,
-and multiplied by the RAI, in the case of a model running without an integrated
-soil model.
+and multiplied by the RAI, in the case of a model running without a prognostic
+soil model:
+
+Flux = -K_eff x [(ψ_stem - ψ_soil)/(z_stem - z_soil) + 1], where
+K_eff = K_soil K_stem /(K_stem + K_soil)
+
+Note that in `PrescribedSoil` mode, we compute the flux using K_soil = K_plant(ψ_soil)
+and K_stem = K_plant(ψ_stem). In `PrognosticSoil` mode, we compute the flux using
+K_soil = K_soil(ψ_soil) and K_stem = K_plant(ψ_stem). The latter is a better model, but
+our `PrescribedSoil` struct does not store K_soil, only ψ_soil.
 
 The returned flux is per unit ground area. This assumes that the stem compartment
 is the first element of `Y.canopy.hydraulics.ϑ_l`.

diff --git a/test/standalone/Vegetation/plant_hydraulics_test.jl b/test/standalone/Vegetation/plant_hydraulics_test.jl
@@ -215,7 +215,7 @@ for FT in (Float32, Float64)
             SAI,
             RAI,
         )
-        K_sat_plant = 1.8e-8 # m/s.
+        K_sat_plant = 1.0 # m/s.
         ψ63 = FT(-4 / 0.0098) # / MPa to m
         Weibull_param = FT(4) # unitless
         a = FT(0.05 * 0.0098) # 1/m
@@ -250,7 +250,7 @@ for FT in (Float32, Float64)
         )
 
         function leaf_transpiration(t)
-            T = FT(1e-8) # m/s
+            T = FT(1) # m/s
         end
 
         ψ_soil0 = FT(0.0)
@@ -284,8 +284,8 @@ for FT in (Float32, Float64)
             )
             # 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(0.0), root = RAI, stem = SAI)
-                T0A = FT(1e-8) * AI[:leaf]
+                AI = (; leaf = LAI(1.0), root = RAI, stem = SAI)
+                T0A = AI[:leaf]
                 for i in 1:(n_leaf + n_stem)
                     if i == 1
                         fa =
@@ -307,7 +307,7 @@ for FT in (Float32, Float64)
                                     vcat(root_depths, [0.0])[2:end] -
                                     vcat(root_depths, [0.0])[1:(end - 1)]
                                 ),
-                            ) * AI[:root]
+                            ) * AI[:stem]
                     else
                         fa =
                             flux(
@@ -317,16 +317,13 @@ for FT in (Float32, Float64)
                                 Y[i],
                                 PlantHydraulics.hydraulic_conductivity(
                                     conductivity_model,
-                                    ψ_soil0,
+                                    Y[i - 1],
                                 ),
                                 PlantHydraulics.hydraulic_conductivity(
                                     conductivity_model,
                                     Y[i],
                                 ),
-                            ) * (
-                                AI[plant_hydraulics.compartment_labels[1]] +
-                                AI[plant_hydraulics.compartment_labels[2]]
-                            ) / 2
+                            ) * AI[plant_hydraulics.compartment_labels[i]]
                     end
                     F[i] = fa - T0A
                 end
@@ -335,7 +332,9 @@ for FT in (Float32, Float64)
             soln = nlsolve(
                 initial_compute_exp_tendency!,
                 Vector{FT}(-0.03:0.01:0.07);
-                ftol = 1e-11,
+                ftol = eps(FT),
+                method = :newton,
+                iterations = 20,
             )
 
             S_l =
@@ -372,12 +371,8 @@ for FT in (Float32, Float64)
             canopy_exp_tendency! = make_exp_tendency(model)
             canopy_exp_tendency!(dY, Y, p, 0.0)
 
-            m = similar(dY.canopy.hydraulics.ϑ_l.:1)
-            m .= FT(0)
-            for i in 1:(n_stem + n_leaf)
-                @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0)
-            end
-            @test maximum(parent(m)) < 200 * eps(FT) # starts in equilibrium
+            tolerance = sqrt(eps(FT))
+            @test all(parent(dY.canopy.hydraulics.ϑ_l) .< tolerance) # starts in equilibrium
 
 
             # repeat using the plant hydraulics model directly
@@ -394,13 +389,12 @@ for FT in (Float32, Float64)
             standalone_exp_tendency! =
                 make_compute_exp_tendency(model.hydraulics, model)
             standalone_exp_tendency!(standalone_dY, Y, p, 0.0)
-
-            m = similar(dY.canopy.hydraulics.ϑ_l.:1)
-            m .= FT(0)
-            for i in 1:(n_stem + n_leaf)
-                @. m += sqrt(dY.canopy.hydraulics.ϑ_l.:($$i)^2.0)
-            end
-            @test maximum(parent(m)) < 10^2.1 * eps(FT) # starts in equilibrium
+            @test all(
+                parent(
+                    standalone_dY.canopy.hydraulics.ϑ_l .-
+                    dY.canopy.hydraulics.ϑ_l,
+                ) .≈ FT(0),
+            )
         end
     end