diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 4d9e549003..a3853d70e4 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -40,10 +40,6 @@ steps: command: "julia --color=yes --project=experiments experiments/standalone/Soil/evaporation.jl" artifact_paths: "experiments/standalone/Soil/evaporation*png" - - label: "ozark_test_hydrology_only" - command: "julia --color=yes --project=experiments experiments/integrated/ozark/hydrology_only/ozark.jl" - artifact_paths: "experiments/integrated/ozark/hydrology_only/*png" - - label: "ozark_test" command: "julia --color=yes --project=experiments experiments/integrated/ozark/ozark.jl" artifact_paths: "experiments/integrated/ozark/*png" diff --git a/experiments/integrated/ozark/hydrology_only/ozark.jl b/experiments/integrated/ozark/hydrology_only/ozark.jl deleted file mode 100644 index 13a7168a19..0000000000 --- a/experiments/integrated/ozark/hydrology_only/ozark.jl +++ /dev/null @@ -1,481 +0,0 @@ -import SciMLBase -import ClimaTimeSteppers as CTS -using ClimaCore -import CLIMAParameters as CP -using Plots -using Statistics -using Dates -using Insolation - -using ClimaLSM -using ClimaLSM.Domains: LSMSingleColumnDomain -using ClimaLSM.Soil -using ClimaLSM.Canopy -using ClimaLSM.Canopy.PlantHydraulics -import ClimaLSM -import ClimaLSM.Parameters as LSMP -include(joinpath(pkgdir(ClimaLSM), "parameters", "create_parameters.jl")) -const FT = Float64 -earth_param_set = create_lsm_parameters(FT) -climalsm_dir = pkgdir(ClimaLSM) -# This reads in the data from the flux tower site and creates -# the atmospheric and radiative driver structs for the model -include( - joinpath( - climalsm_dir, - "experiments/integrated/ozark/ozark_met_drivers_FLUXNET.jl", - ), -) -include(joinpath(climalsm_dir, "experiments/integrated/ozark/ozark_domain.jl")) -include( - joinpath(climalsm_dir, "experiments/integrated/ozark/ozark_parameters.jl"), -) -include( - joinpath( - climalsm_dir, - "experiments/integrated/ozark/hydrology_only/ozark_simulation.jl", - ), -) - -# Now we set up the model. For the soil model, we pick -# a model type and model args: -soil_domain = land_domain.subsurface -soil_ps = Soil.RichardsParameters{FT, vanGenuchten{FT}}( - soil_ν, - vanGenuchten(; α = soil_vg_α, n = soil_vg_n), - soil_K_sat, - soil_S_s, - θ_r, -) - -soil_args = (domain = soil_domain, parameters = soil_ps) -soil_model_type = Soil.RichardsModel{FT} - -# Now we set up the canopy model, which we set up by component: -# Component Types -canopy_component_types = (; - radiative_transfer = Canopy.BeerLambertModel{FT}, - photosynthesis = Canopy.FarquharModel{FT}, - conductance = Canopy.MedlynConductanceModel{FT}, - hydraulics = Canopy.PlantHydraulicsModel{FT}, -) -# Individual Component arguments -# Set up radiative transfer -radiative_transfer_args = (; - parameters = BeerLambertParameters{FT}(; - Ω = Ω, - ld = ld, - α_PAR_leaf = α_PAR_leaf, - α_NIR_leaf = α_NIR_leaf, - λ_γ_PAR = λ_γ_PAR, - λ_γ_NIR = λ_γ_NIR, - ) -) -# Set up conductance -conductance_args = (; - parameters = MedlynConductanceParameters{FT}(; - g1 = g1, - Drel = Drel, - g0 = g0, - ) -) -# Set up photosynthesis -photosynthesis_args = (; - parameters = FarquharParameters{FT}( - Canopy.C3(); - oi = oi, - ϕ = ϕ, - θj = θj, - f = f, - sc = sc, - pc = pc, - Vcmax25 = Vcmax25, - Γstar25 = Γstar25, - Kc25 = Kc25, - Ko25 = Ko25, - To = To, - ΔHkc = ΔHkc, - ΔHko = ΔHko, - ΔHVcmax = ΔHVcmax, - ΔHΓstar = ΔHΓstar, - ΔHJmax = ΔHJmax, - ΔHRd = ΔHRd, - ) -) -# Set up plant hydraulics -ai_parameterization = PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI) - -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(; - ai_parameterization = ai_parameterization, - ν = plant_ν, - S_s = plant_S_s, - root_distribution = root_distribution, - conductivity_model = conductivity_model, - retention_model = retention_model, -) - -plant_hydraulics_args = ( - parameters = plant_hydraulics_ps, - n_stem = n_stem, - n_leaf = n_leaf, - compartment_midpoints = compartment_midpoints, - compartment_surfaces = compartment_surfaces, -) - -# Canopy component args -canopy_component_args = (; - radiative_transfer = radiative_transfer_args, - photosynthesis = photosynthesis_args, - conductance = conductance_args, - hydraulics = plant_hydraulics_args, -) -# Other info needed -shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( - z0_m, - z0_b, - earth_param_set, -) - -canopy_model_args = (; parameters = shared_params, domain = land_domain.surface) - -# Integrated plant hydraulics and soil model -# If you wish to pass in a prescribed transpiration function (for -# debugging purposes), you would pass it as another element -# of this tuple. The default if it is not passed is -# a diagnostic (big leaf-computed) transpiration. -soil_driver = PrognosticSoil(; soil_α_PAR = soil_α_PAR, soil_α_NIR = soil_α_NIR) -land_input = (atmos = atmos, soil_driver = soil_driver, radiation = radiation) -land = SoilPlantHydrologyModel{FT}(; - land_args = land_input, - soil_model_type = soil_model_type, - soil_args = soil_args, - canopy_component_types = canopy_component_types, - canopy_component_args = canopy_component_args, - canopy_model_args = canopy_model_args, -) -Y, p, cds = initialize(land) -exp_tendency! = make_exp_tendency(land) -imp_tendency! = make_imp_tendency(land) - -#Initial conditions -Y.soil.ϑ_l = SWC[Int(round(t0 / 1800))] # Get soil water content at t0 -ψ_stem_0 = FT(-1e5 / 9800) -ψ_leaf_0 = FT(-2e5 / 9800) - -S_l_ini = - inverse_water_retention_curve.( - retention_model, - [ψ_stem_0, ψ_leaf_0], - plant_ν, - plant_S_s, - ) - -for i in 1:2 - Y.canopy.hydraulics.ϑ_l.:($i) .= - augmented_liquid_fraction.(plant_ν, S_l_ini[i]) -end - -set_initial_aux_state! = make_set_initial_aux_state(land) -set_initial_aux_state!(p, Y, t0); - -# Set up timestepper and jacobian for soil -update_jacobian! = make_update_jacobian(land.soil) - -ode_algo = CTS.IMEXAlgorithm( - timestepper, - CTS.NewtonsMethod( - max_iters = max_iterations, - update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), - convergence_checker = conv_checker, - ), -) - -W = RichardsTridiagonalW(Y) -jac_kwargs = (; jac_prototype = W, Wfact = update_jacobian!) - -# Simulation -sv = (; - t = Array{FT}(undef, length(saveat)), - saveval = Array{NamedTuple}(undef, length(saveat)), -) -cb = ClimaLSM.NonInterpSavingCallback(sv, saveat) - -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLSM.dss!, - ), - Y, - (t0, tf), - p, -) -sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - callback = cb, - adaptive = false, - saveat = saveat, -) - -# Plotting -daily = sol.t ./ 3600 ./ 24 -savedir = joinpath(climalsm_dir, "experiments/integrated/ozark/hydrology_only/") -model_GPP = [ - parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] for - k in 1:length(sv.saveval) -] - -plt1 = Plots.plot(size = (500, 700)) -Plots.plot!( - plt1, - daily, - model_GPP .* 1e6, - label = "Model", - xlim = [minimum(daily), maximum(daily)], - xlabel = "days", - ylabel = "GPP [μmol/mol]", -) -Plots.plot!( - plt1, - seconds ./ 3600 ./ 24, - GPP .* 1e6, - label = "Data", - margins = 10Plots.mm, - lalpha = 0.3, -) -plt2 = Plots.plot( - seconds ./ 3600 ./ 24, - FT.(SW_IN), - label = "Data", - xlim = [minimum(daily), maximum(daily)], - ylabel = " SW radiation (W/m^2)", - size = (500, 700), -) -Plots.plot(plt1, plt2, layout = (2, 1)) -Plots.savefig(joinpath(savedir, "GPP.png")) - - -T = - [ - parent(sv.saveval[k].canopy.conductance.transpiration)[1] for - k in 1:length(sv.saveval) - ] .* (1e3 * 24 * 3600) -measured_T = LE ./ (LSMP.LH_v0(earth_param_set) * 1000) .* (1e3 * 24 * 3600) -plt1 = Plots.plot(size = (500, 700)) -Plots.plot!( - plt1, - daily, - T, - label = "Modeled T", - xlim = [minimum(daily), maximum(daily)], - xlabel = "days", - ylabel = "Vapor Flux [mm/day]", -) -Plots.plot!( - plt1, - seconds ./ 3600 ./ 24, - measured_T, - label = "Data (vapor flux)", - margins = 10Plots.mm, - lalpha = 0.3, -) - -plt2 = Plots.plot( - seconds ./ 3600 ./ 24, - VPD, - label = "Measured", - ylabel = "VPD (Pa)", - xlim = [minimum(daily), maximum(daily)], - size = (500, 700), - margins = 10Plots.mm, -) -Plots.plot(plt2, plt1, layout = (2, 1)) -Plots.savefig(joinpath(savedir, "ET.png")) -plt1 = Plots.plot(size = (500, 700)) - -β = [parent(sv.saveval[k].canopy.hydraulics.β)[1] for k in 1:length(sol.t)] -plt1 = Plots.plot( - daily, - β, - label = "Model", - xlim = [minimum(daily), maximum(daily)], - xlabel = "days", - ylabel = "Moisture Stress", - size = (500, 700), -) -g_stomata = - [parent(sv.saveval[k].canopy.conductance.gs)[1] for k in 1:length(sol.t)] -plt2 = Plots.plot( - daily, - g_stomata, - label = "Model", - xlim = [minimum(daily), maximum(daily)], - xlabel = "days", - ylabel = "Stomatal conductance", - ylim = [0, 0.3], - size = (500, 700), -) -Plots.plot(plt2, plt1, layout = (2, 1)) -Plots.savefig(joinpath(savedir, "stomatal_conductance.png")) - - -plt1 = Plots.plot(size = (500, 700)) -Plots.plot!( - plt1, - daily, - [parent(sol.u[k].soil.ϑ_l)[end] for k in 1:1:length(sol.t)], - label = "10cm", - xtickfontsize = 5, - ytickfontsize = 5, - xlim = [minimum(daily), maximum(daily)], - ylim = [0.05, 0.5], - xlabel = "Days", - ylabel = "SWC [m/m]", - margins = 10Plots.mm, - color = "blue", -) - -plot!( - plt1, - daily, - [parent(sol.u[k].soil.ϑ_l)[end - 1] for k in 1:1:length(sol.t)], - label = "20cm", - color = "red", -) - -plot!( - plt1, - daily, - [parent(sol.u[k].soil.ϑ_l)[end - 2] for k in 1:1:length(sol.t)], - label = "30cm", - color = "purple", -) - -Plots.plot!(plt1, seconds ./ 3600 ./ 24, SWC, label = "Data") -plt2 = Plots.plot( - seconds ./ 3600 ./ 24, - P .* (-1e3 * 24 * 3600), - label = "Data", - ylabel = "Precipitation [mm/day]", - xlim = [minimum(daily), maximum(daily)], - size = (500, 700), -) -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] -# Find which index in the data our simulation starts at: -idx = argmin(abs.(seconds .- sol.t[1])) -Plots.plot( - daily, - cumsum(T) * dt_model .+ cumsum(measured_T[:])[idx] * dt_data, - 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) -] - -stem_leaf_flux = [ - parent(sv.saveval[k].canopy.hydraulics.fa)[1] .* (1e3 * 3600 * 24) for - k in 1:length(sol.t) -] -leaf_air_flux = [ - parent(sv.saveval[k].canopy.hydraulics.fa)[2] .* (1e3 * 3600 * 24) for - k in 1:length(sol.t) -] - - -plt1 = Plots.plot( - daily, - 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, root_stem_flux, label = "Soil-root-stem flux") - -lwp = [ - parent(sv.saveval[k].canopy.hydraulics.ψ)[2] * 9800 for k in 1:length(sol.t) -] -swp = [ - parent(sv.saveval[k].canopy.hydraulics.ψ)[1] * 9800 for k in 1:length(sol.t) -] -ψ_soil = [ - sum( - parent(sv.saveval[k].soil.ψ) .* - root_distribution.(parent(cds.subsurface.z)), - ) / sum(root_distribution.(parent(cds.subsurface.z))) * 9800 for - k in 1:length(sol.t) -] -plt2 = Plots.plot( - daily, - lwp, - 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, ψ_soil, label = "Model, Mean soil") -Plots.plot!( - plt2, - seconds_since_t0 ./ 24 ./ 3600, - lwp_measured, - label = "Data; all species", -) - - -Plots.plot(plt2, plt1, layout = (2, 1)) -Plots.savefig(joinpath(savedir, "plant_hydraulics.png")) - -rm(joinpath(savedir, "Artifacts.toml")) diff --git a/experiments/integrated/ozark/hydrology_only/ozark_simulation.jl b/experiments/integrated/ozark/hydrology_only/ozark_simulation.jl deleted file mode 100644 index 5634520b41..0000000000 --- a/experiments/integrated/ozark/hydrology_only/ozark_simulation.jl +++ /dev/null @@ -1,10 +0,0 @@ -t0 = FT(140 * 3600 * 24) -N_days = 140 -tf = t0 + FT(3600 * 24 * N_days) -dt = FT(225) -n = 16 -saveat = Array(t0:(n * dt):tf) -timestepper = CTS.ARS222() -norm_condition = CTS.MaximumError(FT(1e-8)) -conv_checker = CTS.ConvergenceChecker(; norm_condition) -max_iterations = 20 diff --git a/src/ClimaLSM.jl b/src/ClimaLSM.jl index 6bd9861f42..28daaa670f 100644 --- a/src/ClimaLSM.jl +++ b/src/ClimaLSM.jl @@ -292,7 +292,6 @@ import .Canopy.PlantHydraulics: root_flux_per_ground_area! ### Concrete types of AbstractLandModels ### and associated methods include("integrated/soil_energy_hydrology_biogeochemistry.jl") -include("integrated/soil_plant_hydrology_model.jl") include("integrated/pond_soil_model.jl") include("integrated/soil_canopy_model.jl") diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl index 504692474c..64b7099fd7 100644 --- a/src/integrated/soil_canopy_model.jl +++ b/src/integrated/soil_canopy_model.jl @@ -324,6 +324,51 @@ function net_radiation_at_ground( ) end +""" + PlantHydraulics.root_flux_per_ground_area!( + fa::ClimaCore.Fields.Field, + s::PrognosticSoil, + model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, + Y::ClimaCore.Fields.FieldVector, + p::NamedTuple, + t::FT, + ) where {FT} + +An extension of the `PlantHydraulics.root_flux_per_ground_area!` function, + which returns the +net flux of water between the +roots and the soil, per unit ground area, +when both soil and plant +hydraulics are modeled prognostically. This is for use in an LSM. + +It is computed by summing the flux of water per ground area between +roots and soil at each soil layer. +""" +function PlantHydraulics.root_flux_per_ground_area!( + fa::ClimaCore.Fields.Field, + s::PrognosticSoil, + model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, + Y::ClimaCore.Fields.FieldVector, + p::NamedTuple, + t::FT, +) where {FT} + fa .= sum(p.root_extraction) +end + +""" + RootExtraction{FT} <: Soil.AbstractSoilSource{FT} + +Concrete type of Soil.AbstractSoilSource, used for dispatch +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. +""" +struct RootExtraction{FT} <: Soil.AbstractSoilSource{FT} end + + """ ClimaLSM.source!(dY::ClimaCore.Fields.FieldVector, src::RootExtraction, diff --git a/src/integrated/soil_plant_hydrology_model.jl b/src/integrated/soil_plant_hydrology_model.jl deleted file mode 100644 index 888c8fdb8d..0000000000 --- a/src/integrated/soil_plant_hydrology_model.jl +++ /dev/null @@ -1,238 +0,0 @@ -export SoilPlantHydrologyModel -""" - struct SoilPlantHydrologyModel{ - FT, - SM <: Soil.AbstractSoilModel{FT}, - VM <: Canopy.CanopyModel{FT}, - } <: AbstractLandModel{FT} - soil::SM - canopy::VM - end - -A concrete type of land model used for simulating systems with a -canopy and a soil hydrology component. - -This is primarily for testing purposes. -$(DocStringExtensions.FIELDS) -""" -struct SoilPlantHydrologyModel{ - FT, - SM <: Soil.AbstractSoilModel{FT}, - VM <: Canopy.CanopyModel{FT}, -} <: AbstractLandModel{FT} - "The soil model to be used" - soil::SM - "The canopy model to be used" - canopy::VM -end - - -""" - SoilPlantHydrologyModel{FT}(; - land_args::NamedTuple = (;), - soil_model_type::Type{SM}, - soil_args::NamedTuple = (;), - canopy_component_types::NamedTuple = (;), - canopy_component_args::NamedTuple = (;), - canopy_model_args::NamedTuple = (;), - ) where {FT, SM <: Soil.AbstractSoilModel{FT}} - -A constructor for the `SoilPlantHydrologyModel`, which takes in the concrete model -type and required arguments for each component, constructs those models, -and constructs the `SoilPlantHydrologyModel` from them. - -In this model, the soil_driver, of type PrognosticSoil, must be defined by the -user and passed in within the land_args tuple. The soil_driver is used to -provide the canopy model with the soil albedo, and for the full soil + canopy -model, the soil model contains the albedo, and the constructor for the model is -able to define and pass the soil_driver, but here for the soil model, the soil -albedo is not defined, and therefore the user must manually define the -soil driver and provide it in land_args. - - -Each component model is constructed with everything it needs to be stepped -forward in time, including boundary conditions, source terms, and interaction -terms. - -Runoff is not currently supported for this configuration, as this model is for testing only; - please see the documentation for the -`SoilCanopyModel`, `EnergyHydrology`, or `RichardsModel` for models that support runoff. -""" -function SoilPlantHydrologyModel{FT}(; - land_args::NamedTuple = (;), - soil_model_type::Type{SM}, - soil_args::NamedTuple = (;), - canopy_component_types::NamedTuple = (;), - canopy_component_args::NamedTuple = (;), - canopy_model_args::NamedTuple = (;), -) where {FT, SM <: Soil.AbstractSoilModel{FT}} - - # These may be passed in, or set, depending on use scenario. - (; soil_driver, atmos, radiation) = land_args - precipitation = atmos.liquid_precip - # These should always be set by the constructor. - sources = (RootExtraction{FT}(),) - - boundary_conditions = (; - top = (water = FluxBC((p, t) -> eltype(t)(precipitation(t))),), - bottom = (water = Soil.FreeDrainage(),), - ) - - soil = soil_model_type(; - boundary_conditions = boundary_conditions, - sources = sources, - soil_args..., - ) - - if :transpiration in propertynames(land_args) - transpiration = land_args.transpiration - else - transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}() - end - - canopy = Canopy.CanopyModel{FT}(; - radiative_transfer = canopy_component_types.radiative_transfer( - canopy_component_args.radiative_transfer..., - ), - photosynthesis = canopy_component_types.photosynthesis( - canopy_component_args.photosynthesis..., - ), - conductance = canopy_component_types.conductance( - canopy_component_args.conductance..., - ), - hydraulics = canopy_component_types.hydraulics(; - transpiration = transpiration, - canopy_component_args.hydraulics..., - ), - soil_driver = soil_driver, - atmos = atmos, - radiation = radiation, - canopy_model_args..., - ) - - return SoilPlantHydrologyModel{FT, typeof(soil), typeof(canopy)}( - soil, - canopy, - ) -end - -interaction_vars(m::SoilPlantHydrologyModel) = (:root_extraction,) - -interaction_types(m::SoilPlantHydrologyModel{FT}) where {FT} = (FT,) - -interaction_domain_names(m::SoilPlantHydrologyModel) = (:subsurface,) - -""" - make_interactions_update_aux( - land::SoilPlantHydrologyModel{FT, SM, RM}, - ) where {FT, SM <: Soil.RichardsModel{FT}, RM <: Canopy.CanopyModel{FT}} - -A method which makes a function; the returned function -updates the auxiliary variable `p.root_extraction`, which -is needed for the soil model and -for the canopy model. - -This function is called each ode function evaluation. - -Root extraction is in units of 1/s and is equivalent to: -RAI * flux per cross section of roots * root distribution (z). -""" -function make_interactions_update_aux( - land::SoilPlantHydrologyModel{FT, SM, RM}, -) where {FT, SM <: Soil.RichardsModel{FT}, RM <: Canopy.CanopyModel{FT}} - function update_aux!(p, Y, t) - z = ClimaCore.Fields.coordinate_field(land.soil.domain.space).z - (; conductivity_model) = land.canopy.hydraulics.parameters - area_index = p.canopy.hydraulics.area_index - @. p.root_extraction = - ( - area_index.root + getproperty( - area_index, - land.canopy.hydraulics.compartment_labels[1], - ) - ) / 2 * - PlantHydraulics.flux( - z, - land.canopy.hydraulics.compartment_midpoints[1], - p.soil.ψ, - p.canopy.hydraulics.ψ.:1, - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - p.soil.ψ, - ), - PlantHydraulics.hydraulic_conductivity( - conductivity_model, - p.canopy.hydraulics.ψ.:1, - ), - ) * - (land.canopy.hydraulics.parameters.root_distribution(z)) - end - return update_aux! -end - -""" - PlantHydraulics.root_flux_per_ground_area!( - fa::ClimaCore.Fields.Field, - s::PrognosticSoil, - model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, - Y::ClimaCore.Fields.FieldVector, - p::NamedTuple, - t::FT, - ) where {FT} - -An extension of the `PlantHydraulics.root_flux_per_ground_area!` function, - which returns the -net flux of water between the -roots and the soil, per unit ground area, -when both soil and plant -hydraulics are modeled prognostically. This is for use in an LSM. - -It is computed by summing the flux of water per ground area between -roots and soil at each soil layer. -""" -function PlantHydraulics.root_flux_per_ground_area!( - fa::ClimaCore.Fields.Field, - s::PrognosticSoil, - model::Canopy.PlantHydraulics.PlantHydraulicsModel{FT}, - Y::ClimaCore.Fields.FieldVector, - p::NamedTuple, - t::FT, -) where {FT} - fa .= sum(p.root_extraction) -end - -""" - RootExtraction{FT} <: Soil.AbstractSoilSource{FT} - -Concrete type of Soil.AbstractSoilSource, used for dispatch -in an LSM with both soil and plant hydraulic components. - -This is paired with the source term `Canopy.PlantHydraulics.PrognosticSoilPressure`:both -are used at the same time, -ensuring that the water flux into the roots is extracted correctly -from the soil. -""" -struct RootExtraction{FT} <: Soil.AbstractSoilSource{FT} end - -""" - ClimaLSM.source!(dY::ClimaCore.Fields.FieldVector, - src::RootExtraction, - Y::ClimaCore.Fields.FieldVector, - p::NamedTuple - model::RichardsModel) - -An extension of the `ClimaLSM.source!` function, - which computes source terms for the -soil model; this method returns the water loss or gain due -to roots when a plant hydraulic prognostic model is included. -""" -function ClimaLSM.source!( - dY::ClimaCore.Fields.FieldVector, - src::RootExtraction, - Y::ClimaCore.Fields.FieldVector, - p::NamedTuple, - model::RichardsModel, -) - @. dY.soil.ϑ_l += -1 * p.root_extraction - # if flow is negative, towards soil -> soil water increases, add in sign here. -end