Merge pull request #758 from CliMA/kd/root_extraction_nans

extraction fix
CliMA · Sep 12, 2024 · 8f582be · 8f582be
2 parents 6627c51 + 8208c89
commit 8f582be
Show file tree

Hide file tree

Showing 7 changed files with 85 additions and 52 deletions.
diff --git a/docs/src/APIs/canopy/PlantHydraulics.md b/docs/src/APIs/canopy/PlantHydraulics.md
@@ -17,7 +17,7 @@ ClimaLand.PlantHydraulics.augmented_liquid_fraction
 ClimaLand.PlantHydraulics.water_retention_curve
 ClimaLand.PlantHydraulics.inverse_water_retention_curve
 ClimaLand.PlantHydraulics.root_water_flux_per_ground_area!
-ClimaLand.PlantHydraulics.flux
+ClimaLand.PlantHydraulics.water_flux
 ClimaLand.PlantHydraulics.hydraulic_conductivity
 ```
 

diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
@@ -28,7 +28,8 @@ import ClimaAnalysis
 import ClimaAnalysis.Visualize as viz
 import ClimaUtilities
 
-import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
@@ -47,7 +48,7 @@ using Dates
 import NCDatasets
 
 const FT = Float64;
-
+time_interpolation_method = LinearInterpolation(PeriodicCalendar())
 regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
@@ -83,6 +84,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     snow_precip = TimeVaryingInput(
@@ -92,6 +94,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     u_atmos = TimeVaryingInput(
@@ -100,20 +103,23 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "q",
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sp",
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
 
     T_atmos = TimeVaryingInput(
@@ -122,6 +128,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     h_atmos = FT(10)
 
@@ -145,6 +152,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
     LW_d = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -153,6 +161,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
 
     function zenith_angle(
@@ -489,6 +498,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         file_reader_kwargs = (;
             preprocess_func = (data) -> data > 0.05 ? data : 0.0,
         ),
+        method = time_interpolation_method,
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
@@ -613,6 +623,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         start_date;
         output_writer = nc_writer,
         output_vars = :short,
+        average_period = :monthly,
     )
 
     diagnostic_handler =
@@ -627,7 +638,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 14
+    tf = 60 * 60.0 * 24 * 365
     Δt = 900.0
     nelements = (101, 15)
     if greet

diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
@@ -28,7 +28,8 @@ using ClimaAnalysis
 import ClimaAnalysis.Visualize as viz
 using ClimaUtilities
 
-import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput
+import ClimaUtilities.TimeVaryingInputs:
+    TimeVaryingInput, LinearInterpolation, PeriodicCalendar
 import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaUtilities.Regridders: InterpolationsRegridder
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
@@ -45,7 +46,7 @@ using Dates
 import NCDatasets
 
 const FT = Float64;
-
+time_interpolation_method = LinearInterpolation(PeriodicCalendar())
 regridder_type = :InterpolationsRegridder
 context = ClimaComms.context()
 device = ClimaComms.device()
@@ -81,6 +82,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     snow_precip = TimeVaryingInput(
@@ -90,6 +92,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,),
+        method = time_interpolation_method,
     )
 
     u_atmos = TimeVaryingInput(
@@ -98,20 +101,23 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     q_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"),
         "q",
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     P_atmos = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
         "sp",
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
 
     T_atmos = TimeVaryingInput(
@@ -120,6 +126,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         surface_space;
         reference_date = start_date,
         regridder_type,
+        method = time_interpolation_method,
     )
     h_atmos = FT(10)
 
@@ -143,6 +150,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
     LW_d = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"),
@@ -151,6 +159,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         reference_date = start_date,
         regridder_type,
         file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,),
+        method = time_interpolation_method,
     )
 
     function zenith_angle(
@@ -430,7 +439,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 340
+    tf = 60 * 60.0 * 24 * 365
     Δt = 900.0
     nelements = (101, 15)
     if greet

diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
@@ -263,18 +263,22 @@ 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 *
-            PlantHydraulics.flux(
+            PlantHydraulics.water_flux(
                 z,
                 land.canopy.hydraulics.compartment_midpoints[1],
                 p.soil.ψ,
                 p.canopy.hydraulics.ψ.:1,
-                PlantHydraulics.hydraulic_conductivity(
-                    conductivity_model,
-                    p.soil.ψ,
-                ),
+                p.soil.K,
                 PlantHydraulics.hydraulic_conductivity(
                     conductivity_model,
                     p.canopy.hydraulics.ψ.:1,

diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl
@@ -524,7 +524,7 @@ function ClimaLand.make_update_aux(
             # Compute the flux*area between the current compartment `i`
             # and the compartment above.
             @. fa.:($$i) =
-                PlantHydraulics.flux(
+                PlantHydraulics.water_flux(
                     hydraulics.compartment_midpoints[i],
                     hydraulics.compartment_midpoints[ip1],
                     ψ.:($$i),

diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl
@@ -23,7 +23,7 @@ import ClimaLand:
     name
 export PlantHydraulicsModel,
     AbstractPlantHydraulicsModel,
-    flux,
+    water_flux,
     effective_saturation,
     augmented_liquid_fraction,
     water_retention_curve,
@@ -328,7 +328,7 @@ end
 
 
 """
-    flux(
+    water_flux(
         z1,
         z2,
         ψ1,
@@ -337,18 +337,25 @@ end
         K2,
     ) where {FT}
 
-Computes the water flux given the absolute potential (pressure/(ρg))
- at the center of the two compartments z1 and z2,
-and the conductivity along
-the flow path between these two points.
+Computes the water flux given the absolute potential ψ (pressure/(ρg))
+ and the conductivity K (m/s) at the center of the two layers
+with midpoints z1 and z2.
 
-We currently assuming an arithmetic
-mean for mean K_sat between the two points (Bonan, 2019; Zhu, 2008)
-to take into account the change in K_sat halfway
-between z1 and z2; this is incorrect for compartments of differing sizes.
+We currently assuming a harmonic
+mean for effective conducticity between the two layers
+(see CLM Technical Documentation).
+
+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) 
+and a water flux of
+F = -k_eff * (ψ1 +z1 - ψ2 - z2) (m/s).
+
+This currently assumes the path lengths are equal.
 """
-function flux(z1, z2, ψ1, ψ2, K1, K2)
-    flux = -(K1 + K2) / 2 * ((ψ2 - ψ1) / (z2 - z1) + 1)
+function water_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
 
@@ -580,8 +587,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`.
@@ -610,7 +625,7 @@ function root_water_flux_per_ground_area!(
 
         if i != n_root_layers
             @. fa +=
-                flux(
+                water_flux(
                     root_depths[i],
                     model.compartment_midpoints[1],
                     ψ_soil,
@@ -623,7 +638,7 @@ function root_water_flux_per_ground_area!(
                 above_ground_area_index
         else
             @. fa +=
-                flux(
+                water_flux(
                     root_depths[i],
                     model.compartment_midpoints[1],
                     ψ_soil,