Merge #238

238: Zero lai r=kmdeck a=kmdeck

## Purpose 
In the current codebase, setting the area indices of the plant to zero would result in a 0/0, because, for example, we have:
` d theta_leaf/dt = [stem_flux *LAI - canopy_transpiration]/ (LAI x Delta z)`, where canopy_transpiration involves a multiplication by LAI as well. 

Because we cannot have a changing number of prognostic variables across the domain, even when there is no vegetation present, we need to have a canopy model with some prognostic variables, parameters, etc. To turn off the canopy, setting LAI = SAI = RAI = 0 will imply zero fluxes and zero root extraction, and no change in the initial conditions of the plant prognostic variables. 

Therefore, update the codebase so that when LAI=SAI=RAI = 0, `d theta/ dt = 0` and there is no divide by zero.

Once we get closer to global runs, we should do some refactoring as it makes sense to make it so you dont need to define all these parameters where there is no vegetation. See Issue #250.

## Content
1. Change "LAI/(LAI x dz)" terms to "LAI/max(LAI x dz, eps(FT))". When LAI = 0, this gives zero flux and does not divide by zero.
2. Add consistency checks to the PlantHydraulics constructor:
- if LAI == 0, all area indices should be zero (no plant present)
- if LAI is nonzero, RAI must be nonzero (a plant cannot exist without roots)
- if n_stem == 0, the stem area index must be zero for consistency.
3. Adds unit tests to make sure consistency checks pass and that when LAI =RAI = SAI = 0, the root extraction and transpiration are both zero.

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:/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: kmdeck <[email protected]>

src/Vegetation/PlantHydraulics.jl

@@ -107,24 +107,74 @@ function PlantHydraulicsParameters(;
  )
 end
 
+"""
+ lai_consistency_check(
+ n_stem::Int64,
+ n_leaf::Int64,
+ area_index::NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}},
+ ) where {FT}
+
+Carries out consistency checks using the area indices supplied and the number of
+stem elements being modeled.
+
+Note that it is possible to have a plant with no stem compartments
+but with leaf compartments, and that a plant must have leaf compartments
+(even if LAI = 0).
 
+Specifically, this checks that:
+1. n_leaf > 0
+2. if LAI is nonzero or SAI is nonzero, RAI must be nonzero.
+3. if SAI > 0, n_stem must be > 0 for consistency. If SAI == 0, n_stem must
+be zero.
 """
- PlantHydraulicsModel{FT, PS, RE, T, B} <: AbstractPlantHydraulicsModel{FT}
+function lai_consistency_check(
+ n_stem::Int64,
+ n_leaf::Int64,
+ area_index::NamedTuple{(:root, :stem, :leaf), Tuple{FT, FT, FT}},
+) where {FT}
+ @assert n_leaf > 0
+ if area_index[:leaf] > eps(FT) || area_index[:stem] > eps(FT)
+ @assert area_index[:root] > eps(FT)
+ end
+ # If there SAI is zero, there should be no stem compartment
+ if area_index[:stem] < eps(FT)
+ @assert n_stem == FT(0)
+ else
+ # if SAI is > 0, n_stem should be > 0 for consistency
+ @assert n_stem > 0
+ end
+
+end
+
+"""
+ PlantHydraulicsModel{FT, PS, RE, T} <: AbstractPlantHydraulicsModel{FT}
 
 Defines, and constructs instances of, the PlantHydraulicsModel type, which is used
 for simulation flux of water to/from soil, along roots of different depths,
 along a stem, to a leaf, and ultimately being lost from the system by
 transpiration.
 
-This model can be used in standalone mode by prescribing the transpiration rate
-and soil matric potential at the root tips or flux in the roots, or with a
-dynamic soil model using `ClimaLSM`.
+This model can also be combined with the soil model using ClimaLSM, in which
+case the prognostic soil water content is used to determine root extraction, and
+the transpiration is also computed diagnostically. In global run with patches
+of bare soil, you can "turn off" the canopy model (to get zero root extraction, zero absorption and
+emission, zero transpiration and sensible heat flux from the canopy), by setting:
+- n_leaf = 1
+- n_stem = 0
+- LAI = SAI = RAI = 0.
+
+A plant model can have leaves but no stem, but not vice versa. If n_stem = 0, SAI must be zero.
+
+Finally, the model can be used in Canopy standalone mode by prescribing
+the soil matric potential at the root tips or flux in the roots. There is also the
+option (intendend only for debugging) to use a prescribed transpiration rate.
+
 $(DocStringExtensions.FIELDS)
 """
 struct PlantHydraulicsModel{FT, PS, RE, T} <: AbstractPlantHydraulicsModel{FT}
- "The number of stem compartments for the plant"
+ "The number of stem compartments for the plant; can be zero"
  n_stem::Int64
- "The number of leaf compartments for the plant"
+ "The number of leaf compartments for the plant; must be >=1"
  n_leaf::Int64
  "The depth of the root tips, in meters"
  root_depths::Vector{FT}
@@ -153,7 +203,8 @@ function PlantHydraulicsModel{FT}(;
  transpiration::AbstractTranspiration{FT} = DiagnosticTranspiration{FT}(),
 ) where {FT}
  args = (parameters, root_extraction, transpiration)
- @assert n_leaf != 0
+ area_index = parameters.area_index
+ lai_consistency_check(n_stem, n_leaf, area_index)
  @assert (n_leaf + n_stem) == length(compartment_midpoints)
  @assert (n_leaf + n_stem) + 1 == length(compartment_surfaces)
  for i in 1:length(compartment_midpoints)
@@ -408,7 +459,17 @@ end
 A function which creates the compute_exp_tendency! function for the PlantHydraulicsModel.
 The compute_exp_tendency! function must comply with a rhs function of OrdinaryDiffEq.jl.
 
-Below, `fa` denotes a flux multiplied by the relevant cross section (per unit area ground).
+Below, `fa` denotes a flux multiplied by the relevant cross section 
+(per unit area ground, or area index, AI). The tendency for the 
+ith compartment can be written then as:
+∂ϑ[i]/∂t = 1/(AI*dz)[fa[i]-fa[i+1]).
+
+Note that if the area_index is zero because no plant is present, 
+AIdz is zero, and the fluxes `fa` appearing in the numerator are 
+zero because they are scaled by AI.
+
+To prevent dividing by zero, we change AI/(AI x dz)" to
+"AI/max(AI x dz, eps(FT))"
 """
 function make_compute_exp_tendency(model::PlantHydraulicsModel, _)
  function compute_exp_tendency!(dY, Y, p, t)
@@ -417,13 +478,17 @@ function make_compute_exp_tendency(model::PlantHydraulicsModel, _)
  n_leaf = model.n_leaf
  fa = p.canopy.hydraulics.fa
  fa_roots = p.canopy.hydraulics.fa_roots
-
+ FT = eltype(t)
  @inbounds for i in 1:(n_stem + n_leaf)
- AIdz =
+ # To prevent dividing by zero, change AI/(AI x dz)" to
+ # "AI/max(AI x dz, eps(FT))"
+ AIdz = max(
  area_index[model.compartment_labels[i]] * (
  model.compartment_surfaces[i + 1] -
  model.compartment_surfaces[i]
- )
+ ),
+ eps(FT),
+ )
  if i == 1
  # All fluxes `fa` are per unit area of ground
  root_flux_per_ground_area!(
@@ -488,34 +553,40 @@ function root_flux_per_ground_area!(
  ψ_base = p.canopy.hydraulics.ψ[1]
  n_root_layers = length(model.root_depths)
  ψ_soil::FT = re.ψ_soil(t)
- @inbounds for i in 1:n_root_layers
- if i != n_root_layers
- @. fa +=
- flux(
- model.root_depths[i],
- model.compartment_midpoints[1],
- ψ_soil,
- ψ_base,
- hydraulic_conductivity(conductivity_model, ψ_soil),
- hydraulic_conductivity(conductivity_model, ψ_base),
- ) *
- root_distribution(model.root_depths[i]) *
- (model.root_depths[i + 1] - model.root_depths[i]) *
- area_index[:root]
- else
- @. fa +=
- flux(
- model.root_depths[i],
- model.compartment_midpoints[1],
- ψ_soil,
- ψ_base,
- hydraulic_conductivity(conductivity_model, ψ_soil),
- hydraulic_conductivity(conductivity_model, ψ_base),
- ) *
- root_distribution(model.root_depths[i]) *
- (FT(0) - model.root_depths[n_root_layers]) *
- area_index[:root]
-
+ if area_index[:root] < eps(FT)
+ fa .= FT(0)
+ else
+ @inbounds for i in 1:n_root_layers
+ if i != n_root_layers
+ @. fa +=
+ flux(
+ model.root_depths[i],
+ model.compartment_midpoints[1],
+ ψ_soil,
+ ψ_base,
+ hydraulic_conductivity(conductivity_model, ψ_soil),
+ hydraulic_conductivity(conductivity_model, ψ_base),
+ ) *
+ root_distribution(model.root_depths[i]) *
+ (model.root_depths[i + 1] - model.root_depths[i]) *
+ area_index[:root]
+ else
+ @. fa +=
+ flux(
+ model.root_depths[i],
+ model.compartment_midpoints[1],
+ ψ_soil,
+ ψ_base,
+ hydraulic_conductivity(conductivity_model, ψ_soil),
+ hydraulic_conductivity(conductivity_model, ψ_base),
+ ) *
+ root_distribution(model.root_depths[i]) *
+ (FT(0) - model.root_depths[n_root_layers]) *
+ (
+ area_index[:root] +
+ area_index[model.compartment_labels[1]]
+ ) / 2
+ end
  end
  end
 end

src/soil_plant_hydrology_model.jl

@@ -133,7 +133,10 @@ function make_interactions_update_aux(
  z = ClimaCore.Fields.coordinate_field(land.soil.domain.space).z
  (; area_index, conductivity_model) = land.canopy.hydraulics.parameters
  @. p.root_extraction =
- area_index[:root] *
+ (
+ area_index[:root] +
+ area_index[land.canopy.hydraulics.compartment_labels[1]]
+ ) / 2 *
  PlantHydraulics.flux(
  z,
  land.canopy.hydraulics.compartment_midpoints[1],

test/Vegetation/canopy_model.jl

@@ -104,7 +104,7 @@ include(joinpath(pkgdir(ClimaLSM), "parameters", "create_parameters.jl"))
 
  # Plant Hydraulics
  RAI = FT(1)
- SAI = FT(1)
+ SAI = FT(0)
  area_index = (root = RAI, stem = SAI, leaf = LAI)
  K_sat_plant = FT(1.8e-8) # m/s
  ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value