explicit and implicit sources

Author: kmdeck

NEWS.md

@@ -3,6 +3,7 @@ ClimaLand.jl Release Notes
 
 main
 -------
+- Add capability to step some sources implicitly PR[#1113](https://github.com/CliMA/ClimaLand.jl/pull/1113)
 
 v0.15.14
 -------

docs/list_tutorials.jl

@@ -43,5 +43,6 @@ tutorials = [
         "Intro to multi-component models" => "standalone/Usage/LSM_single_column_tutorial.jl",
         "Intro to ClimaLand Domains" => "standalone/Usage/domain_tutorial.jl",
         "Intro to forced site-level runs" => "shared_utilities/driver_tutorial.jl",
+        "Intro to implicit/explicit timestepping" => "shared_utilities/timestepping.jl",
     ],
 ]

docs/tutorials/shared_utilities/timestepping.jl

@@ -0,0 +1,74 @@
+# # Mixed implicit and explicit timestepping
+# The goal of this tutorial is to describe how
+# the timestepping of a ClimaLand model is carried out.
+# We will use forward and backward Euler as a demonstration,
+# but higher order methods are available in ClimaTimesteppers.
+
+# # Explicit vs. implicit stepping
+# Given a differential equation for a prognostic variable
+# Y
+
+# ```
+# \frac{d Y }{dt} = g(Y, t) + h(Y, t),
+# ```
+
+# an explicit (forward) Euler step would entail
+
+# ```
+# Y(t+\Delta t) = Y(t) + [g(Y(t), t) + h(Y(t),t)] \times \Delta t,
+# ```
+
+# while an implicit (backward) Euler step would entail
+
+# ```
+# Y(t+\Delta t) = Y(t) + [g(Y(t+\Delta t), t) + h(Y(t+\Delta t),t)] \times \Delta t,
+# ```
+
+# which reqires us to solve an implicit equation for  Y(t+ Δt). We
+# usually do so using Newton's method, which requires the derivative of
+# the entire right hand side with respect to the variable we are solving
+# for, `Y(t+ Δt)`. This is called the Jacobian.
+
+# Sometimes certain terms must be stepped implicit for numerical stability,
+# while others are more slowly varying or stable. In this case, a mixed
+# approach would be
+
+# ```
+# Y(t+\Delta t) = Y(t) + [g(Y(t+\Delta t), t) + h(Y(t),t)] \times \Delta t,
+# ```
+
+# assuming that `h` is the slow term, and `g` is the fast term. Note that
+# solving this implicit equation for `(Y(t+ Δt)` with Newton's method
+# would be similar to that of the fully implicit approach, but with an
+# approximated Jacobian (neglecting ∂h/∂Y).
+
+# If our timestepping scheme involves evaluating all right-hand-side
+# tendencies at the current (known) value of a prognostic variable,
+# we refer to that prognostic variable as explicit. If any of the
+# right-hand-side tendencies are evaluated at the next (unknown) value
+# of the prognostic variable, we refer to it as implicit. In the latter
+# case, the Jacobian would include a term like `∂ tendency/∂ variable`,
+# even if it is an approximate (not exact) form of the Jacobian.
+
+# # Implicit and explicit prognostic variables of the land model
+# We treat two prognostic variables of the soil model (ϑ_l,
+# ρe_int) and the canopy temperature implicitly, and 
+# the canopy water content, the soil ice content,
+# and all prognostic
+# variables of the snow model explicitly.
+
+# # Implicit vs explicit tendencies - not complete as of 4/23/25
+# Implicit:
+# - Vertical contribution of the divergence of the Darcy flux
+# - Vertical contribution of the divergence of diffusive heat flux
+# - Vertical contribution of the divergence of heat flux due to Darcy flow
+# - Canopy temperature (except root extraction of energy)
+# - SHF, LHF, evaporation, and sublimation of soil (note that these are explicit in θ_i!)
+# - Soil radiation (does not contribute to Jacobian)
+# - Subsurface runoff (this is computed in the same function the same time as surface runoff, but does not contribute to Jacobian.)
+
+# Explicit
+# - Phase changes in soil
+# - Root extraction
+# - all snow tendencies
+# - Darcy flux within the canopy

docs/tutorials/standalone/Soil/sublimation.jl

@@ -158,7 +158,7 @@ init_soil!(Y, z, soil.parameters);
 # Timestepping:
 t0 = Float64(0)
 tf = Float64(24 * 3600 * 4)
-dt = Float64(5)
+dt = Float64(10)
 
 # We also set the initial conditions of the cache here:
 set_initial_cache! = make_set_initial_cache(soil)
@@ -175,7 +175,7 @@ timestepper = CTS.ARS111()
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 3,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 )

src/integrated/soil_canopy_root_interactions.jl

@@ -125,10 +125,11 @@ in an LSM with both soil and plant hydraulic components.
 This is paired with the source term `Canopy.PrognosticSoil`:both
 are used at the same time,
 ensuring that the water flux into the roots is extracted correctly
-from the soil.
+from the soil; treated explicitly in all prognostic variables.
 """
-struct RootExtraction{FT} <: Soil.AbstractSoilSource{FT} end
-
+@kwdef struct RootExtraction{FT} <: ClimaLand.Soil.AbstractSoilSource{FT}
+    explicit::Bool = true
+end
 """
     ClimaLand.source!(dY::ClimaCore.Fields.FieldVector,
                      src::RootExtraction,
@@ -150,5 +151,16 @@ function ClimaLand.source!(
 )
     @. dY.soil.ϑ_l += -1 * p.root_extraction
     @. dY.soil.ρe_int += -1 * p.root_energy_extraction
+
     # if flow is negative, towards soil -> soil water increases, add in sign here.
+    # Currently, these column integrals are done twice rather than add space to the cache.
+    # We can revisit this as needed.
+    ClimaCore.Operators.column_integral_definite!(
+        p.scratch1,
+        p.root_energy_extraction,
+    )
+    @. dY.soil.∫F_e_dt += p.scratch1
+
+    ClimaCore.Operators.column_integral_definite!(p.scratch1, p.root_extraction)
+    @. dY.soil.∫F_vol_liq_water_dt += p.scratch1
 end

src/integrated/soil_snow_model.jl

@@ -390,9 +390,15 @@ end
     SoilSublimationwithSnow{FT} <: AbstractSoilSource{FT}
 
 Soil Sublimation source type. Used to defined a method
-of `ClimaLand.source!` for soil sublimation with snow present.
+of `ClimaLand.source!` for soil sublimation with snow present;
+treated implicitly
+in ϑ_l, ρe_int but explicitly in θ_i.
+
 """
-struct SoilSublimationwithSnow{FT} <: ClimaLand.Soil.AbstractSoilSource{FT} end
+@kwdef struct SoilSublimationwithSnow{FT} <:
+              ClimaLand.Soil.AbstractSoilSource{FT}
+    explicit::Bool = false
+end
 
 """
      source!(dY::ClimaCore.Fields.FieldVector,
@@ -419,7 +425,10 @@ function ClimaLand.source!(
     @. dY.soil.θ_i +=
         -p.soil.turbulent_fluxes.vapor_flux_ice *
         (1 - p.snow.snow_cover_fraction) *
-        _ρ_l / _ρ_i * heaviside(z + 2 * Δz_top) # only apply to top layer, recall that z is negative
+        _ρ_l / _ρ_i * heaviside(z + 2 * Δz_top) / (2 * Δz_top) # only apply to top layer, recall that z is negative
+    @. dY.soil.∫F_vol_liq_water_dt +=
+        -p.soil.turbulent_fluxes.vapor_flux_ice *
+        (1 - p.snow.snow_cover_fraction) # The integral of the source is designed to be this
     return nothing
 end
 

src/shared_utilities/Domains.jl

@@ -935,7 +935,7 @@ function get_additional_coordinate_field_data(subsurface_space)
     z_face = ClimaCore.Fields.coordinate_field(face_space).z
     z_sfc = top_face_to_surface(z_face, surface_space)
     d = depth(subsurface_space)
-    Δz_min = minimum(Δz_top)
+    Δz_min = minimum(Δz)
     fields = (;
         z = z,
         Δz_top = Δz_top,

src/standalone/Soil/Runoff/Runoff.jl

@@ -163,13 +163,21 @@ The runoff flux is given by Equation 12 of Niu et al. (2005),
 "A simple TOPMODEL-based runoff parameterization (SIMTOP) for
 use in global climate models".
 
+This is currently treated implicitly (evaluated at the next value of
+ϑ_l) but not included in the Jacobian approximation. This is because
+`update_runoff` is called as part of the implicit tendency.
+
 $(DocStringExtensions.FIELDS)
 """
 struct TOPMODELSubsurfaceRunoff{FT} <: AbstractSoilSource{FT}
     "The subsurface runoff flux (m/s) when the depth to the water table = 1/f_over; calibrated"
     R_sb::FT
     "A calibrated parameter defining how subsurface runoff decays with depth to water table (1/m ; calibrated)"
     f_over::FT
+    "Explicit vs implicit stepping"
+    explicit::Bool
+    TOPMODELSubsurfaceRunoff{FT}(R_sb, f_over) where {FT} =
+        new{FT}(R_sb, f_over, false)
 end
 
 """
@@ -290,7 +298,7 @@ function ClimaLand.source!(
     h∇ = p.soil.h∇
     ϵ = eps(FT)
     @. dY.soil.ϑ_l -= (p.soil.R_ss / max(h∇, ϵ)) * p.soil.is_saturated # apply only to saturated layers
-    @. p.soil.∫S_θ_liq_dz -= p.soil.R_ss # the integral is designed to be this flux
+    @. dY.soil.∫F_vol_liq_water_dt -= p.soil.R_ss # the integral is designed to be this flux
 end
 
 """

src/standalone/Soil/energy_hydrology.jl

@@ -198,8 +198,8 @@ function ClimaLand.make_compute_exp_tendency(
     function compute_exp_tendency!(dY, Y, p, t)
         z = model.domain.fields.z
 
-        p.soil.∫S_θ_liq_dz .= 0
-        p.soil.∫S_ρe_int_dz .= 0
+        dY.soil.∫F_vol_liq_water_dt .= 0
+        dY.soil.∫F_e_dt .= 0
 
         # Don't update the prognostic variables we're stepping implicitly
         dY.soil.ϑ_l .= 0
@@ -218,13 +218,13 @@ function ClimaLand.make_compute_exp_tendency(
             p,
             z,
         )
-        # Source terms, also update the vertical integral of the source: p.soil.∫S_θ_liq_dz, p.soil.∫S_ρe_int_dz
-        # These change ∫S and dY by +=, which is why we zero them out above.
+        # Explicitly treated source terms
+        # These change dY by +=, which is why we ".=" them above
         for src in model.sources
-            ClimaLand.source!(dY, src, Y, p, model)
+            if src.explicit
+                ClimaLand.source!(dY, src, Y, p, model)
+            end
         end
-        dY.soil.∫F_vol_liq_water_dt .= p.soil.∫S_θ_liq_dz # source terms are stepped explicitly
-        dY.soil.∫F_e_dt .= p.soil.∫S_ρe_int_dz # source terms are stepped explicitly
     end
     return compute_exp_tendency!
 end
@@ -291,8 +291,15 @@ function ClimaLand.make_compute_imp_tendency(
                 ) * gradc2f(p.soil.ψ + z),
             )
 
-        # Don't update the prognostic variables we're stepping explicitly
         @. dY.soil.θ_i = 0
+
+        # Source terms
+        # These change dY by +=, which is why we ".=" above
+        for src in model.sources
+            if !src.explicit
+                ClimaLand.source!(dY, src, Y, p, model)
+            end
+        end
     end
     return compute_imp_tendency!
 end
@@ -493,8 +500,6 @@ of `EnergyHydrology`.
 ClimaLand.auxiliary_vars(soil::EnergyHydrology) = (
     :total_water,
     :total_energy,
-    :∫S_θ_liq_dz,
-    :∫S_ρe_int_dz,
     :K,
     :ψ,
     :θ_l,
@@ -521,8 +526,6 @@ ClimaLand.auxiliary_types(soil::EnergyHydrology{FT}) where {FT} = (
     FT,
     FT,
     FT,
-    FT,
-    FT,
     boundary_var_types(
         soil,
         soil.boundary_conditions.top,
@@ -536,8 +539,6 @@ ClimaLand.auxiliary_types(soil::EnergyHydrology{FT}) where {FT} = (
 )
 
 ClimaLand.auxiliary_domain_names(soil::EnergyHydrology) = (
-    :surface,
-    :surface,
     :surface,
     :surface,
     :subsurface,
@@ -647,9 +648,12 @@ end
 """
     PhaseChange{FT} <: AbstractSoilSource{FT}
 
-PhaseChange source type.
+PhaseChange source type; treated explicitly in all
+prognostic variables.
 """
-struct PhaseChange{FT} <: AbstractSoilSource{FT} end
+@kwdef struct PhaseChange{FT} <: AbstractSoilSource{FT}
+    explicit::Bool = true
+end
 
 """
      source!(dY::ClimaCore.Fields.FieldVector,
@@ -659,7 +663,8 @@ struct PhaseChange{FT} <: AbstractSoilSource{FT} end
              model
              )
 
-Computes the source terms for phase change.
+Computes the source terms for phase change
+explicitly in time.
 
 """
 function ClimaLand.source!(
@@ -728,7 +733,7 @@ function ClimaLand.source!(
             _grav,
         )
 
-    # These source/sink terms conserve total water, so we do not include them in ∫S_θ_liq_dz because their
+    # These source/sink terms conserve total water, so we do not include them when checking conservation because their
     # contribution is zero by design.
 end
 
@@ -737,10 +742,12 @@ end
     SoilSublimation{FT} <: AbstractSoilSource{FT}
 
 Soil Sublimation source type. Used to defined a method
-of `ClimaLand.source!` for soil sublimation.
+of `ClimaLand.source!` for soil sublimation; treated implicitly
+in ϑ_l, ρe_int but explicitly in θ_i.
 """
-struct SoilSublimation{FT} <: AbstractSoilSource{FT} end
-
+@kwdef struct SoilSublimation{FT} <: AbstractSoilSource{FT}
+    explicit::Bool = false
+end
 """
      source!(dY::ClimaCore.Fields.FieldVector,
              src::SoilSublimation{FT},
@@ -767,7 +774,7 @@ function ClimaLand.source!(
     @. dY.soil.θ_i +=
         -p.soil.turbulent_fluxes.vapor_flux_ice * _ρ_l / _ρ_i *
         heaviside(z + 2 * Δz_top) / (2 * Δz_top) # only apply to top layer, recall that z is negative
-    @. p.soil.∫S_θ_liq_dz += -p.soil.turbulent_fluxes.vapor_flux_ice # The integral of the source is designed to be this
+    @. dY.soil.∫F_vol_liq_water_dt += -p.soil.turbulent_fluxes.vapor_flux_ice # The integral of the source is designed to be this
 end
 
 ## The functions below are required to be defined