From 6d320d2d64a38b3637bccdefa1cfdb9c9c6b0039 Mon Sep 17 00:00:00 2001 From: kmdeck Date: Wed, 5 Jun 2024 13:57:43 -0700 Subject: [PATCH] zero flux if zero lai/sai [skip ci] --- Artifacts.toml | 9 +- .../standalone/Vegetation/no_vegetation.jl | 270 +++++++++++++++++ .../standalone/Vegetation/varying_lai.jl | 270 +++++++++++++++++ .../Vegetation/varying_lai_with_stem.jl | 283 ++++++++++++++++++ src/integrated/soil_canopy_model.jl | 6 +- src/standalone/Vegetation/Canopy.jl | 11 +- src/standalone/Vegetation/PlantHydraulics.jl | 4 +- .../Vegetation/autotrophic_respiration.jl | 15 +- .../Vegetation/canopy_boundary_fluxes.jl | 155 +++++++++- .../Vegetation/canopy_parameterizations.jl | 21 +- src/standalone/Vegetation/radiation.jl | 29 +- test/standalone/Vegetation/canopy_model.jl | 261 +++++++++++++++- 12 files changed, 1267 insertions(+), 67 deletions(-) create mode 100644 experiments/standalone/Vegetation/no_vegetation.jl create mode 100644 experiments/standalone/Vegetation/varying_lai.jl create mode 100644 experiments/standalone/Vegetation/varying_lai_with_stem.jl diff --git a/Artifacts.toml b/Artifacts.toml index fdabc938f9..2793031c7f 100644 --- a/Artifacts.toml +++ b/Artifacts.toml @@ -1,9 +1,12 @@ +[ameriflux_data_US-MOz] +git-tree-sha1 = "37177736bb2a1d14416961561911e06c24942c4a" + +[era5_land_forcing_data2021] +git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a" + [soil_params_Gupta2020_2022] git-tree-sha1 = "8e28b4274b10020b6cdd54b8e7585221379d9d33" [[soil_params_Gupta2020_2022.download]] sha256 = "97dcf1158cba03b1fd397262bdfaf85a523f57038c337bcce163e32664d3616b" url = "https://caltech.box.com/shared/static/f2y23qx0lggjskftzgh7ht7fsbh36gmm.gz" - -[era5_land_forcing_data2021] -git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a" diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl new file mode 100644 index 0000000000..86f3993549 --- /dev/null +++ b/experiments/standalone/Vegetation/no_vegetation.jl @@ -0,0 +1,270 @@ +import SciMLBase +using CairoMakie +using Statistics +using Dates +using Insolation + +# Load CliMA Packages and ClimaLand Modules: + +using ClimaCore +import ClimaParams as CP +import ClimaTimeSteppers as CTS +using StaticArrays +using ClimaLand +using ClimaLand.Domains: Point +using ClimaLand.Canopy +using ClimaLand.Canopy.PlantHydraulics +import ClimaLand +import ClimaLand.Parameters as LP +const FT = Float32; +earth_param_set = LP.LandParameters(FT); +f_root_to_shoot = FT(3.5) +SAI = FT(0.0) +plant_ν = FT(2.46e-4) # kg/m^2 +n_stem = Int64(0) +n_leaf = Int64(1) +h_leaf = FT(9.5) +compartment_midpoints = [h_leaf / 2] +compartment_surfaces = [FT(0), h_leaf] +land_domain = Point(; z_sfc = FT(0.0)) +include( + joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"), +); +time_offset = 7 +lat = FT(38.7441) # degree +long = FT(-92.2000) # degree +atmos_h = FT(32) +site_ID = "US-MOz" +data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv" + +include( + joinpath( + pkgdir(ClimaLand), + "experiments/integrated/fluxnet/met_drivers_FLUXNET.jl", + ), +); + +z0_m = FT(2) +z0_b = FT(0.2) + +shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, +); +ψ_soil0 = FT(0.0) + +soil_driver = PrescribedSoil( + FT; + root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0), + ψ = t -> ψ_soil0, + α_PAR = FT(0.2), + α_NIR = FT(0.4), + T = t -> 298.0, + ϵ = FT(0.99), +); + +rt_params = TwoStreamParameters( + FT; + G_Function = ConstantGFunction(FT(0.5)), + α_PAR_leaf = FT(0.1), + α_NIR_leaf = FT(0.45), + τ_PAR_leaf = FT(0.05), + τ_NIR_leaf = FT(0.25), + Ω = FT(0.69), + λ_γ_PAR = FT(5e-7), + λ_γ_NIR = FT(1.65e-6), +) + +rt_model = TwoStreamModel{FT}(rt_params); + +cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0)) + +stomatal_model = MedlynConductanceModel{FT}(cond_params); + + +photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5)) + +photosynthesis_model = FarquharModel{FT}(photo_params); + +AR_params = AutotrophicRespirationParameters(FT) +AR_model = AutotrophicRespirationModel{FT}(AR_params); + +function fakeLAIfunction(t) + if t < 30 * 24 * 3600 + 0.0 + elseif t < (364 - 30) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 30 * 24 * 3600)), 0) + else + 0.0 + end +end + + +f_root_to_shoot = FT(3.5) +SAI = FT(0) +RAI = FT(2 * f_root_to_shoot) +ai_parameterization = + PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction), SAI, RAI) +rooting_depth = FT(1.0); + +function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) +end; + +K_sat_plant = FT(1.8e-8) +ψ63 = FT(-4 / 0.0098) +Weibull_param = FT(4) +a = FT(0.05 * 0.0098) + +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a); + +ν = FT(0.7) +S_s = FT(1e-2 * 0.0098) + +plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = ν, + S_s = S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, +); + + +plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_surfaces, + compartment_midpoints = compartment_midpoints, +); + +energy_model = ClimaLand.Canopy.BigLeafEnergyModel{FT}( + BigLeafEnergyParameters{FT}(FT(1e3)), +) + +canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = land_domain, + autotrophic_respiration = AR_model, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + energy = energy_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, +); + + +Y, p, coords = ClimaLand.initialize(canopy) +exp_tendency! = make_exp_tendency(canopy); + + +ψ_leaf_0 = FT(-2e5 / 9800) + +S_l_ini = inverse_water_retention_curve(retention_model, ψ_leaf_0, ν, S_s) + +Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) + + + +t0 = 0.0 +N_days = 10 +tf = t0 + 3600 * 24 * N_days +dt = 225.0; +evaluate!(Y.canopy.energy.T, atmos.T, t0) +set_initial_cache! = make_set_initial_cache(canopy) +set_initial_cache!(p, Y, t0); + + +n = 16 +saveat = Array(t0:(n * dt):tf) +sv = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), +) +saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); + +updateat = Array(t0:1800:tf) +updatefunc = ClimaLand.make_update_drivers(atmos, radiation) +driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) +cb = SciMLBase.CallbackSet(driver_cb, saving_cb); + +timestepper = CTS.RK4(); +ode_algo = CTS.ExplicitAlgorithm(timestepper) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), + Y, + (t0, tf), + p, +); + +sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); + +savedir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Vegetation"); +T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)] +T_atmos = [parent(sv.saveval[k].drivers.T)[1] for k in 1:length(sol.t)] +ϑ = [parent(sol.u[k].canopy.hydraulics.ϑ_l.:1)[1] for k in 1:length(sol.t)] +GPP = [ + parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] * 1e6 for + k in 1:length(sol.t) +] +resp = [ + parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] * 1e6 for + k in 1:length(sol.t) +] +SW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.SW_n)[1] for + k in 1:length(sol.t) +] +LW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.LW_n)[1] for + k in 1:length(sol.t) +] +SHF = [parent(sv.saveval[k].canopy.energy.shf)[1] for k in 1:length(sol.t)] +LHF = [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)] +RE = [ + parent(sv.saveval[k].canopy.energy.fa_energy_roots)[1] for + k in 1:length(sol.t) +] +R = [ + parent(sv.saveval[k].canopy.hydraulics.fa_roots)[1] for k in 1:length(sol.t) +] +Tr = [parent(sv.saveval[k].canopy.hydraulics.fa.:1)[1] for k in 1:length(sol.t)] +fig = Figure() +ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") +lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") +lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +axislegend(ax) +ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") +lines!(ax, sol.t ./ 24 ./ 3600, ϑ, label = "Canopy") +axislegend(ax) +ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") +lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction.(sol.t), label = "Canopy") +axislegend(ax) +save(joinpath(savedir, "no_veg_state.png"), fig) +fig2 = Figure() +ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, SW_abs, label = "SW") +lines!(ax, sol.t ./ 24 ./ 3600, LW_abs, label = "LW") +lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") +lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") +lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +axislegend(ax) +ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") +lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") + +axislegend(ax) +ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") +lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +axislegend(ax) +save(joinpath(savedir, "no_veg_fluxes.png"), fig2) diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl new file mode 100644 index 0000000000..34f61b9f2b --- /dev/null +++ b/experiments/standalone/Vegetation/varying_lai.jl @@ -0,0 +1,270 @@ +import SciMLBase +using CairoMakie +using Statistics +using Dates +using Insolation + +# Load CliMA Packages and ClimaLand Modules: + +using ClimaCore +import ClimaParams as CP +import ClimaTimeSteppers as CTS +using StaticArrays +using ClimaLand +using ClimaLand.Domains: Point +using ClimaLand.Canopy +using ClimaLand.Canopy.PlantHydraulics +import ClimaLand +import ClimaLand.Parameters as LP +const FT = Float32; +earth_param_set = LP.LandParameters(FT); +f_root_to_shoot = FT(3.5) +SAI = FT(0.0) +plant_ν = FT(2.46e-4) # kg/m^2 +n_stem = Int64(0) +n_leaf = Int64(1) +h_leaf = FT(9.5) +compartment_midpoints = [h_leaf / 2] +compartment_surfaces = [FT(0), h_leaf] +land_domain = Point(; z_sfc = FT(0.0)) +include( + joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"), +); +time_offset = 7 +lat = FT(38.7441) # degree +long = FT(-92.2000) # degree +atmos_h = FT(32) +site_ID = "US-MOz" +data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv" + +include( + joinpath( + pkgdir(ClimaLand), + "experiments/integrated/fluxnet/met_drivers_FLUXNET.jl", + ), +); + +z0_m = FT(2) +z0_b = FT(0.2) + +shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, +); +ψ_soil0 = FT(0.0) + +soil_driver = PrescribedSoil( + FT; + root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0), + ψ = t -> ψ_soil0, + α_PAR = FT(0.2), + α_NIR = FT(0.4), + T = t -> 298.0, + ϵ = FT(0.99), +); + +rt_params = TwoStreamParameters( + FT; + G_Function = ConstantGFunction(FT(0.5)), + α_PAR_leaf = FT(0.1), + α_NIR_leaf = FT(0.45), + τ_PAR_leaf = FT(0.05), + τ_NIR_leaf = FT(0.25), + Ω = FT(0.69), + λ_γ_PAR = FT(5e-7), + λ_γ_NIR = FT(1.65e-6), +) + +rt_model = TwoStreamModel{FT}(rt_params); + +cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0)) + +stomatal_model = MedlynConductanceModel{FT}(cond_params); + + +photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5)) + +photosynthesis_model = FarquharModel{FT}(photo_params); + +AR_params = AutotrophicRespirationParameters(FT) +AR_model = AutotrophicRespirationModel{FT}(AR_params); + +function fakeLAIfunction(t) + if t < 30 * 24 * 3600 + 0.0 + elseif t < (364 - 30) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 30 * 24 * 3600)), 0) + else + 0.0 + end +end + + +f_root_to_shoot = FT(3.5) +SAI = FT(0) +RAI = FT(2 * f_root_to_shoot) +ai_parameterization = + PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction), SAI, RAI) +rooting_depth = FT(1.0); + +function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) +end; + +K_sat_plant = FT(1.8e-8) +ψ63 = FT(-4 / 0.0098) +Weibull_param = FT(4) +a = FT(0.05 * 0.0098) + +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a); + +ν = FT(0.7) +S_s = FT(1e-2 * 0.0098) + +plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = ν, + S_s = S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, +); + + +plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_surfaces, + compartment_midpoints = compartment_midpoints, +); + +energy_model = ClimaLand.Canopy.BigLeafEnergyModel{FT}( + BigLeafEnergyParameters{FT}(FT(1e3)), +) + +canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = land_domain, + autotrophic_respiration = AR_model, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + energy = energy_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, +); + + +Y, p, coords = ClimaLand.initialize(canopy) +exp_tendency! = make_exp_tendency(canopy); + + +ψ_leaf_0 = FT(-2e5 / 9800) + +S_l_ini = inverse_water_retention_curve(retention_model, ψ_leaf_0, ν, S_s) + +Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) + + + +t0 = 0.0 +N_days = 180 +tf = t0 + 3600 * 24 * N_days +dt = 225.0; +evaluate!(Y.canopy.energy.T, atmos.T, t0) +set_initial_cache! = make_set_initial_cache(canopy) +set_initial_cache!(p, Y, t0); + + +n = 16 +saveat = Array(t0:(n * dt):tf) +sv = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), +) +saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); + +updateat = Array(t0:1800:tf) +updatefunc = ClimaLand.make_update_drivers(atmos, radiation) +driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) +cb = SciMLBase.CallbackSet(driver_cb, saving_cb); + +timestepper = CTS.RK4(); +ode_algo = CTS.ExplicitAlgorithm(timestepper) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), + Y, + (t0, tf), + p, +); + +sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); + +savedir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Vegetation"); +T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)] +T_atmos = [parent(sv.saveval[k].drivers.T)[1] for k in 1:length(sol.t)] +ϑ = [parent(sol.u[k].canopy.hydraulics.ϑ_l.:1)[1] for k in 1:length(sol.t)] +GPP = [ + parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] * 1e6 for + k in 1:length(sol.t) +] +resp = [ + parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] * 1e6 for + k in 1:length(sol.t) +] +SW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.SW_n)[1] for + k in 1:length(sol.t) +] +LW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.LW_n)[1] for + k in 1:length(sol.t) +] +SHF = [parent(sv.saveval[k].canopy.energy.shf)[1] for k in 1:length(sol.t)] +LHF = [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)] +RE = [ + parent(sv.saveval[k].canopy.energy.fa_energy_roots)[1] for + k in 1:length(sol.t) +] +R = [ + parent(sv.saveval[k].canopy.hydraulics.fa_roots)[1] for k in 1:length(sol.t) +] +Tr = [parent(sv.saveval[k].canopy.hydraulics.fa.:1)[1] for k in 1:length(sol.t)] +fig = Figure() +ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") +lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") +lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +axislegend(ax) +ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") +lines!(ax, sol.t ./ 24 ./ 3600, ϑ, label = "Canopy") +axislegend(ax) +ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") +lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction.(sol.t), label = "Canopy") +axislegend(ax) +save(joinpath(savedir, "varying_lai_state.png"), fig) +fig2 = Figure() +ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, SW_abs, label = "SW") +lines!(ax, sol.t ./ 24 ./ 3600, LW_abs, label = "LW") +lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") +lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") +lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +axislegend(ax) +ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") +lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") + +axislegend(ax) +ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") +lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +axislegend(ax) +save(joinpath(savedir, "varying_lai_fluxes.png"), fig2) diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl new file mode 100644 index 0000000000..2a83646828 --- /dev/null +++ b/experiments/standalone/Vegetation/varying_lai_with_stem.jl @@ -0,0 +1,283 @@ +import SciMLBase +using CairoMakie +using Statistics +using Dates +using Insolation + +# Load CliMA Packages and ClimaLand Modules: + +using ClimaCore +import ClimaParams as CP +import ClimaTimeSteppers as CTS +using StaticArrays +using ClimaLand +using ClimaLand.Domains: Point +using ClimaLand.Canopy +using ClimaLand.Canopy.PlantHydraulics +import ClimaLand +import ClimaLand.Parameters as LP +const FT = Float32; +earth_param_set = LP.LandParameters(FT); +f_root_to_shoot = FT(3.5) +plant_ν = FT(2.46e-4) # kg/m^2 +n_stem = Int64(1) +n_leaf = Int64(1) +h_leaf = FT(9.5) +h_stem = FT(9) +compartment_midpoints = [h_stem / 2, h_stem + h_leaf / 2] +compartment_surfaces = [FT(0), h_stem, h_stem + h_leaf] +land_domain = Point(; z_sfc = FT(0.0)) +include( + joinpath(pkgdir(ClimaLand), "experiments/integrated/fluxnet/data_tools.jl"), +); +time_offset = 7 +lat = FT(38.7441) # degree +long = FT(-92.2000) # degree +atmos_h = FT(32) +site_ID = "US-MOz" +data_link = "https://caltech.box.com/shared/static/7r0ci9pacsnwyo0o9c25mhhcjhsu6d72.csv" + +include( + joinpath( + pkgdir(ClimaLand), + "experiments/integrated/fluxnet/met_drivers_FLUXNET.jl", + ), +); + +z0_m = FT(2) +z0_b = FT(0.2) + +shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z0_m, + z0_b, + earth_param_set, +); +ψ_soil0 = FT(0.0) + +soil_driver = PrescribedSoil( + FT; + root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0), + ψ = t -> ψ_soil0, + α_PAR = FT(0.2), + α_NIR = FT(0.4), + T = t -> 298.0, + ϵ = FT(0.99), +); + +rt_params = TwoStreamParameters( + FT; + G_Function = ConstantGFunction(FT(0.5)), + α_PAR_leaf = FT(0.1), + α_NIR_leaf = FT(0.45), + τ_PAR_leaf = FT(0.05), + τ_NIR_leaf = FT(0.25), + Ω = FT(0.69), + λ_γ_PAR = FT(5e-7), + λ_γ_NIR = FT(1.65e-6), +) + +rt_model = TwoStreamModel{FT}(rt_params); + +cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0)) + +stomatal_model = MedlynConductanceModel{FT}(cond_params); + + +photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5)) + +photosynthesis_model = FarquharModel{FT}(photo_params); + +AR_params = AutotrophicRespirationParameters(FT) +AR_model = AutotrophicRespirationModel{FT}(AR_params); + +function fakeLAIfunction(t) + if t < 30 * 24 * 3600 + 0.0 + elseif t < (364 - 30) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 30 * 24 * 3600)), 0) + else + 0.0 + end +end + + +f_root_to_shoot = FT(3.5) +SAI = FT(1.0) +RAI = FT(3 * f_root_to_shoot) +ai_parameterization = + PrescribedSiteAreaIndex{FT}(TimeVaryingInput(fakeLAIfunction), SAI, RAI) +rooting_depth = FT(1.0); + +function root_distribution(z::T; rooting_depth = rooting_depth) where {T} + return T(1.0 / rooting_depth) * exp(z / T(rooting_depth)) +end; + +K_sat_plant = FT(1.8e-8) +ψ63 = FT(-4 / 0.0098) +Weibull_param = FT(4) +a = FT(0.05 * 0.0098) + +conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + +retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a); + +ν = FT(0.7) +S_s = FT(1e-2 * 0.0098) + +plant_hydraulics_ps = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = ν, + S_s = S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, +); + + +plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = plant_hydraulics_ps, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_surfaces, + compartment_midpoints = compartment_midpoints, +); + +energy_model = ClimaLand.Canopy.BigLeafEnergyModel{FT}( + BigLeafEnergyParameters{FT}(FT(1e3)), +) + +canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = land_domain, + autotrophic_respiration = AR_model, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + energy = energy_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, +); + + +Y, p, coords = ClimaLand.initialize(canopy) +exp_tendency! = make_exp_tendency(canopy); + + +ψ_leaf_0 = FT(-2e5 / 9800) +ψ_stem_0 = FT(-1e5 / 9800) + +S_l_ini = + inverse_water_retention_curve.( + retention_model, + [ψ_stem_0, ψ_leaf_0], + ν, + S_s, + ) + +Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1]) +Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2]) + + + +t0 = 0.0 +N_days = 10 +tf = t0 + 3600 * 24 * N_days +dt = 225.0; +evaluate!(Y.canopy.energy.T, atmos.T, t0) +set_initial_cache! = make_set_initial_cache(canopy) +set_initial_cache!(p, Y, t0); + + +n = 16 +saveat = Array(t0:(n * dt):tf) +sv = (; + t = Array{Float64}(undef, length(saveat)), + saveval = Array{NamedTuple}(undef, length(saveat)), +) +saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); + +updateat = Array(t0:1800:tf) +updatefunc = ClimaLand.make_update_drivers(atmos, radiation) +driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) +cb = SciMLBase.CallbackSet(driver_cb, saving_cb); + +timestepper = CTS.RK4(); +ode_algo = CTS.ExplicitAlgorithm(timestepper) + +prob = SciMLBase.ODEProblem( + CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), + Y, + (t0, tf), + p, +); + +sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); + +savedir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Vegetation"); +T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)] +T_atmos = [parent(sv.saveval[k].drivers.T)[1] for k in 1:length(sol.t)] +ϑ_l = [parent(sol.u[k].canopy.hydraulics.ϑ_l.:2)[1] for k in 1:length(sol.t)] +ϑ_s = [parent(sol.u[k].canopy.hydraulics.ϑ_l.:1)[1] for k in 1:length(sol.t)] +GPP = [ + parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] * 1e6 for + k in 1:length(sol.t) +] +resp = [ + parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] * 1e6 for + k in 1:length(sol.t) +] +SW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.SW_n)[1] for + k in 1:length(sol.t) +] +LW_abs = [ + parent(sv.saveval[k].canopy.radiative_transfer.LW_n)[1] for + k in 1:length(sol.t) +] +SHF = [parent(sv.saveval[k].canopy.energy.shf)[1] for k in 1:length(sol.t)] +LHF = [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)] +RE = [ + parent(sv.saveval[k].canopy.energy.fa_energy_roots)[1] for + k in 1:length(sol.t) +] +R = [ + parent(sv.saveval[k].canopy.hydraulics.fa_roots)[1] for k in 1:length(sol.t) +] +R_stem_leaf = + [parent(sv.saveval[k].canopy.hydraulics.fa.:1)[1] for k in 1:length(sol.t)] +Tr = [parent(sv.saveval[k].canopy.hydraulics.fa.:2)[1] for k in 1:length(sol.t)] +fig = Figure() +ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") +lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") +lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +axislegend(ax) +ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") +lines!(ax, sol.t ./ 24 ./ 3600, ϑ_l, label = "Leaf") +lines!(ax, sol.t ./ 24 ./ 3600, ϑ_s, label = "Stem") +axislegend(ax) +ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") +lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction.(sol.t), label = "Canopy") +axislegend(ax) +save(joinpath(savedir, "varying_lai_with_stem_state.png"), fig) +fig2 = Figure() +ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, SW_abs, label = "SW") +lines!(ax, sol.t ./ 24 ./ 3600, LW_abs, label = "LW") +lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") +lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") +lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +axislegend(ax) +ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") +lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") +lines!(ax, sol.t ./ 24 ./ 3600, R_stem_leaf, label = "R_stem_leaf") + +axislegend(ax) +ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") +lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") +lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +axislegend(ax) +save(joinpath(savedir, "varying_lai_with_stem_fluxes.png"), fig2) diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl index 754f9184b3..917cfc0b30 100644 --- a/src/integrated/soil_canopy_model.jl +++ b/src/integrated/soil_canopy_model.jl @@ -335,7 +335,7 @@ function lsm_radiant_energy_fluxes!( c = LP.light_speed(earth_param_set) h = LP.planck_constant(earth_param_set) N_a = LP.avogadro_constant(earth_param_set) - (; λ_γ_PAR, λ_γ_NIR, ϵ_canopy) = canopy_radiation.parameters + (; λ_γ_PAR, λ_γ_NIR) = canopy_radiation.parameters energy_per_photon_PAR = h * c / λ_γ_PAR energy_per_photon_NIR = h * c / λ_γ_NIR T_canopy = @@ -380,8 +380,8 @@ function lsm_radiant_energy_fluxes!( N_a * p.canopy.radiative_transfer.tpar * (1 - α_soil_PAR) - - @. LW_d_canopy = (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 # double checked + ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI + @. LW_d_canopy = ((1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4) # double checked @. LW_u_soil = ϵ_soil * _σ * p.T_ground^4 + (1 - ϵ_soil) * LW_d_canopy # double checked @. R_net_soil += ϵ_soil * LW_d_canopy - ϵ_soil * _σ * p.T_ground^4 # double checked @. LW_net_canopy = diff --git a/src/standalone/Vegetation/Canopy.jl b/src/standalone/Vegetation/Canopy.jl index 1e4baf9194..457a1b1fed 100644 --- a/src/standalone/Vegetation/Canopy.jl +++ b/src/standalone/Vegetation/Canopy.jl @@ -442,10 +442,14 @@ function ClimaLand.make_update_aux( (; g1, g0, Drel) = canopy.conductance.parameters area_index = p.canopy.hydraulics.area_index LAI = area_index.leaf + SAI = area_index.stem RAI = area_index.root (; sc, pc) = canopy.photosynthesis.parameters # update radiative transfer + @. p.canopy.radiative_transfer.ϵ = + canopy.radiative_transfer.parameters.ϵ_canopy * + (1 - exp(-(LAI + SAI))) #from CLM 5.0, Tech note 4.20 RT = canopy.radiative_transfer K = extinction_coeff.(G_Function, θs) PAR .= compute_PAR(RT, canopy.radiation, p, t) @@ -535,7 +539,7 @@ function ClimaLand.make_update_aux( conductivity_model, ψ.:($$ip1), ), - ) * (areai + areaip1) / 2 + ) * areaip1 end # We update the fa[n_stem+n_leaf] element once we have computed transpiration, below # update photosynthesis and conductance terms @@ -571,8 +575,11 @@ function ClimaLand.make_update_aux( canopy.autotrophic_respiration, Vcmax25, LAI, + SAI, RAI, - GPP, + K, + Ω, + An, Rd, β, h_canopy, diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl index 07de4f6909..e55e04b0d7 100644 --- a/src/standalone/Vegetation/PlantHydraulics.jl +++ b/src/standalone/Vegetation/PlantHydraulics.jl @@ -620,7 +620,7 @@ function root_water_flux_per_ground_area!( ) * root_distribution(root_depths[i]) * (root_depths[i + 1] - root_depths[i]) * - (area_index.root + above_ground_area_index) / 2 + above_ground_area_index else @. fa += flux( @@ -633,7 +633,7 @@ function root_water_flux_per_ground_area!( ) * root_distribution(root_depths[i]) * (model.compartment_surfaces[1] - root_depths[n_root_layers]) * - (area_index.root + above_ground_area_index) / 2 + above_ground_area_index end end end diff --git a/src/standalone/Vegetation/autotrophic_respiration.jl b/src/standalone/Vegetation/autotrophic_respiration.jl index 4245650f60..ade59018f0 100644 --- a/src/standalone/Vegetation/autotrophic_respiration.jl +++ b/src/standalone/Vegetation/autotrophic_respiration.jl @@ -57,20 +57,23 @@ function compute_autrophic_respiration( model::AutotrophicRespirationModel, Vcmax25, LAI, + SAI, RAI, - GPP, + K, + Ω, + An, Rd, β, h, ) (; ne, ηsl, σl, μr, μs, f1, f2) = model.parameters - - Nl, Nr, Ns = nitrogen_content(ne, Vcmax25, LAI, RAI, ηsl, h, σl, μr, μs) + Nl, Nr, Ns = + nitrogen_content(ne, Vcmax25, LAI, SAI, RAI, ηsl, h, σl, μr, μs) Rpm = plant_respiration_maintenance(Rd, β, Nl, Nr, Ns, f1) - Rg = plant_respiration_growth(f2, GPP, Rpm) - Ra = Rpm + Rg # Should this be a function in canopy_parameterizations.jl or is it ok here? - return Ra + Rg = plant_respiration_growth(f2, An, Rpm) + Ra = Rpm + Rg + return Ra * (1 - exp(-K * LAI * Ω)) / (K * Ω) # adjust to canopy level end Base.broadcastable(model::AutotrophicRespirationModel) = tuple(model) # this is so that @. does not broadcast on Ref(canopy.autotrophic_respiration) diff --git a/src/standalone/Vegetation/canopy_boundary_fluxes.jl b/src/standalone/Vegetation/canopy_boundary_fluxes.jl index 5199989ddf..9ee9518227 100644 --- a/src/standalone/Vegetation/canopy_boundary_fluxes.jl +++ b/src/standalone/Vegetation/canopy_boundary_fluxes.jl @@ -1,3 +1,8 @@ +using SurfaceFluxes +using Thermodynamics +using StaticArrays +import SurfaceFluxes.Parameters as SFP + import ClimaLand: surface_temperature, surface_specific_humidity, @@ -26,7 +31,7 @@ end p, t, ) where {FT} -Returns the surface resistance field of the +Returns the stomatal resistance field of the `CanopyModel` canopy. """ function ClimaLand.surface_resistance( @@ -49,6 +54,7 @@ function ClimaLand.surface_resistance( R, P_air, ) + return 1 ./ canopy_conductance # [s/m] end @@ -153,6 +159,8 @@ function canopy_boundary_fluxes!( root_water_flux = p.canopy.hydraulics.fa_roots root_energy_flux = p.canopy.energy.fa_energy_roots fa = p.canopy.hydraulics.fa + LAI = p.canopy.hydraulics.area_index.leaf + SAI = p.canopy.hydraulics.area_index.stem transpiration = p.canopy.conductance.transpiration shf = p.canopy.energy.shf lhf = p.canopy.energy.lhf @@ -160,11 +168,11 @@ function canopy_boundary_fluxes!( i_end = canopy.hydraulics.n_stem + canopy.hydraulics.n_leaf # Compute transpiration, SHF, LHF - canopy_turbulent_fluxes = turbulent_fluxes(atmos, canopy, Y, p, t) - transpiration .= canopy_turbulent_fluxes.vapor_flux - shf .= canopy_turbulent_fluxes.shf - lhf .= canopy_turbulent_fluxes.lhf - r_ae .= canopy_turbulent_fluxes.r_ae + canopy_tf = canopy_turbulent_fluxes(atmos, canopy, Y, p, t) + transpiration .= canopy_tf.vapor_flux + shf .= canopy_tf.shf + lhf .= canopy_tf.lhf + r_ae .= canopy_tf.r_ae # Transpiration is per unit ground area, not leaf area (mult by LAI) fa.:($i_end) .= PlantHydraulics.transpiration_per_ground_area( @@ -203,3 +211,138 @@ function canopy_boundary_fluxes!( ) end + +""" + function canopy_turbulent_fluxes( + atmos::PrescribedAtmosphere, + model::CanopyModel, + Y::ClimaCore.Fields.FieldVector, + p::NamedTuple, + t, + ) + +A canopy specific function for compute turbulent fluxes with the atmosphere; +returns the latent heat flux, sensible heat flux, vapor flux, and aerodynamic resistance. + + We cannot use the default version in src/shared_utilities/drivers.jl +because the canopy requires a different resistance for vapor and sensible heat +fluxes, and the resistances depend on ustar, which we must compute using +SurfaceFluxes before adjusting to account for these resistances. +""" +function canopy_turbulent_fluxes( + atmos::PrescribedAtmosphere, + model::CanopyModel, + Y::ClimaCore.Fields.FieldVector, + p::NamedTuple, + t, +) + T_sfc = ClimaLand.surface_temperature(model, Y, p, t) + ρ_sfc = ClimaLand.surface_air_density(atmos, model, Y, p, t, T_sfc) + q_sfc = ClimaLand.surface_specific_humidity(model, Y, p, T_sfc, ρ_sfc) + h_sfc = ClimaLand.surface_height(model, Y, p) + leaf_r_stomata = ClimaLand.surface_resistance(model, Y, p, t) + d_sfc = ClimaLand.displacement_height(model, Y, p) + u_air = p.drivers.u + h_air = atmos.h + return canopy_turbulent_fluxes_at_a_point.( + T_sfc, + q_sfc, + ρ_sfc, + h_sfc, + leaf_r_stomata, + d_sfc, + p.drivers.thermal_state, + u_air, + h_air, + p.canopy.hydraulics.area_index.leaf, + p.canopy.hydraulics.area_index.stem, + atmos.gustiness, + model.parameters.z_0m, + model.parameters.z_0b, + Ref(model.parameters.earth_param_set), + ) +end + +""" + function canopy_turbulent_fluxes_at_a_point( + T_sfc::FT, + q_sfc::FT, + ρ_sfc::FT, + h_sfc::FT, + leaf_r_stomata::FT, + d_sfc::FT, + ts_in, + u::FT, + h::FT, + LAI::FT, + SAI::FT, + gustiness::FT, + z_0m::FT, + z_0b::FT, + earth_param_set::EP, + ) where {FT <: AbstractFloat, EP} + +Computes the turbulent surface fluxes for the canopy at a point +and returns the fluxes in a named tuple. +""" +function canopy_turbulent_fluxes_at_a_point( + T_sfc::FT, + q_sfc::FT, + ρ_sfc::FT, + h_sfc::FT, + leaf_r_stomata::FT, + d_sfc::FT, + ts_in, + u::FT, + h::FT, + LAI::FT, + SAI::FT, + gustiness::FT, + z_0m::FT, + z_0b::FT, + earth_param_set::EP, +) where {FT <: AbstractFloat, EP} + thermo_params = LP.thermodynamic_parameters(earth_param_set) + ts_sfc = Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sfc) + + state_sfc = SurfaceFluxes.StateValues(FT(0), SVector{2, FT}(0, 0), ts_sfc) + state_in = SurfaceFluxes.StateValues( + h - d_sfc - h_sfc, + SVector{2, FT}(u, 0), + ts_in, + ) + + # State containers + sc = SurfaceFluxes.ValuesOnly( + state_in, + state_sfc, + z_0m, + z_0b, + beta = FT(1), + gustiness = gustiness, + ) + surface_flux_params = LP.surface_fluxes_parameters(earth_param_set) + conditions = SurfaceFluxes.surface_conditions( + surface_flux_params, + sc; + tol_neutral = SFP.cp_d(surface_flux_params) / 100000, + ) + _LH_v0::FT = LP.LH_v0(earth_param_set) + _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set) + + cp_m::FT = Thermodynamics.cp_m(thermo_params, ts_in) + T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in) + ΔT = T_in - T_sfc + r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc)) + ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in) + ustar::FT = conditions.ustar + r_b::FT = FT(1 / 0.01 * (ustar / 0.04)^(-1 / 2)) # CLM 5, tech note Equation 5.122 + leaf_r_b = r_b / LAI + canopy_r_b = r_b / (LAI + SAI) + E0::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch) + E = E0 * r_ae / (leaf_r_b + leaf_r_stomata + r_ae) # CLM 5, tech note Equation 5.101, and fig 5.2b, assuming all sunlit, f_wet = 0 + Ẽ = E / _ρ_liq + H = -ρ_air * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s + LH = _LH_v0 * E + return (lhf = LH, shf = H, vapor_flux = Ẽ, r_ae = r_ae) +end diff --git a/src/standalone/Vegetation/canopy_parameterizations.jl b/src/standalone/Vegetation/canopy_parameterizations.jl index ed44a8bf4e..107be105f3 100644 --- a/src/standalone/Vegetation/canopy_parameterizations.jl +++ b/src/standalone/Vegetation/canopy_parameterizations.jl @@ -935,6 +935,7 @@ end ne::FT, # Mean leaf nitrogen concentration (kg N (kg C)-1) Vcmax25::FT, # LAI::FT, # Leaf area index + SAI::FT, RAI::FT, ηsl::FT, # live stem wood coefficient (kg C m-3) h::FT, # canopy height (m) @@ -943,15 +944,13 @@ end μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1 ) where {FT} -Computes the nitrogen content of leafs (Nl), roots (Nr) and stems (Ns) -as a function of leaf area index (LAI), specific leaf density (σl), -the carbon content of roots (Rc), the carbon content of stems (Rs), -and mean leaf nitrogen concentration (nm). +Computes the nitrogen content of leafs (Nl), roots (Nr) and stems (Ns). """ function nitrogen_content( ne::FT, # Mean leaf nitrogen concentration (kg N (kg C)-1) Vcmax25::FT, # LAI::FT, # Leaf area index + SAI::FT, RAI::FT, ηsl::FT, # live stem wood coefficient (kg C m-3) h::FT, # canopy height (m) @@ -959,7 +958,7 @@ function nitrogen_content( μr::FT, # Ratio root nitrogen to top leaf nitrogen (-), typical value 1.0 μs::FT, # Ratio stem nitrogen to top leaf nitrogen (-), typical value 0.1 ) where {FT} - Sc = ηsl * h * LAI + Sc = ηsl * h * SAI Rc = σl * RAI nm = Vcmax25 / ne Nl = nm * σl * LAI @@ -998,18 +997,14 @@ end """ plant_respiration_growth( f::FT, # Factor of relative contribution - GPP::FT, # Gross primary productivity + An::FT, # Net photosynthesis Rpm::FT # Plant maintenance respiration ) where {FT} -Computes plant growth respiration as a function of gross primary productivity (GPP), +Computes plant growth respiration as a function of net photosynthesis (An), plant maintenance respiration (Rpm), and a relative contribution factor, f. """ -function plant_respiration_growth( - f2::FT, # Factor of relative contribution, usually 0.25 - GPP::FT, # Gross primary productivity - Rpm::FT, # Plant maintenance respiration -) where {FT} - Rg = f2 * (GPP - Rpm) +function plant_respiration_growth(f2::FT, An::FT, Rpm::FT) where {FT} + Rg = f2 * (An - Rpm) return Rg end diff --git a/src/standalone/Vegetation/radiation.jl b/src/standalone/Vegetation/radiation.jl index 0eeb926a7e..7b0b7b46ef 100644 --- a/src/standalone/Vegetation/radiation.jl +++ b/src/standalone/Vegetation/radiation.jl @@ -113,22 +113,6 @@ Base.@kwdef struct TwoStreamParameters{ "Leaf angle distribution function" G_Function::G end -""" - function TwoStreamParameters{FT, G}(; - ld = ConstantGFunction(FT(0.5)), - α_PAR_leaf = FT(0.3), - τ_PAR_leaf = FT(0.2), - α_NIR_leaf = FT(0.4), - τ_NIR_leaf = FT(0.25), - ϵ_canopy = FT(0.98), - Ω = FT(1), - λ_γ_PAR = FT(5e-7), - λ_γ_NIR = FT(1.65e-6), - n_layers = UInt64(20), - ) where {FT} - -A constructor supplying default values for the TwoStreamParameters struct. -""" Base.eltype(::TwoStreamParameters{FT}) where {FT} = FT @@ -192,10 +176,10 @@ Base.broadcastable(RT::AbstractRadiationModel) = tuple(RT) ClimaLand.name(model::AbstractRadiationModel) = :radiative_transfer ClimaLand.auxiliary_vars(model::Union{BeerLambertModel, TwoStreamModel}) = - (:apar, :par, :rpar, :tpar, :anir, :nir, :rnir, :tnir, :LW_n, :SW_n) + (:apar, :par, :rpar, :tpar, :anir, :nir, :rnir, :tnir, :LW_n, :SW_n, :ϵ) ClimaLand.auxiliary_types( model::Union{BeerLambertModel{FT}, TwoStreamModel{FT}}, -) where {FT} = (FT, FT, FT, FT, FT, FT, FT, FT, FT, FT) +) where {FT} = (FT, FT, FT, FT, FT, FT, FT, FT, FT, FT, FT) ClimaLand.auxiliary_domain_names(::Union{BeerLambertModel, TwoStreamModel}) = ( :surface, :surface, @@ -207,6 +191,7 @@ ClimaLand.auxiliary_domain_names(::Union{BeerLambertModel, TwoStreamModel}) = ( :surface, :surface, :surface, + :surface, ) """ @@ -243,8 +228,7 @@ function canopy_radiant_energy_fluxes!( c = FT(LP.light_speed(earth_param_set)) h = FT(LP.planck_constant(earth_param_set)) N_a = FT(LP.avogadro_constant(earth_param_set)) - (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR, ϵ_canopy) = - canopy.radiative_transfer.parameters + (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR) = canopy.radiative_transfer.parameters APAR = p.canopy.radiative_transfer.apar ANIR = p.canopy.radiative_transfer.anir energy_per_photon_PAR = h * c / λ_γ_PAR @@ -252,13 +236,14 @@ function canopy_radiant_energy_fluxes!( @. p.canopy.radiative_transfer.SW_n = (energy_per_photon_PAR * N_a * APAR) + (energy_per_photon_NIR * N_a * ANIR) - + LAI = p.canopy.hydraulics.area_index.leaf + SAI = p.canopy.hydraulics.area_index.stem + ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI # Long wave: use soil conditions from the PrescribedSoil driver T_soil::FT = s.T(t) ϵ_soil = s.ϵ _σ = FT(LP.Stefan(earth_param_set)) LW_d = p.drivers.LW_d - T_canopy = canopy_temperature(canopy.energy, canopy, Y, p, t) LW_d_canopy = @. (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 LW_u_soil = @. ϵ_soil * _σ * T_soil^4 + (1 - ϵ_soil) * LW_d_canopy diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl index bc2e8958d9..13ef24f842 100644 --- a/test/standalone/Vegetation/canopy_model.jl +++ b/test/standalone/Vegetation/canopy_model.jl @@ -232,7 +232,13 @@ for FT in (Float32, Float64) # check that this is updated correctly: # @test p.canopy.autotrophic_respiration.Ra == exp_tendency!(dY, Y, p, t0) - turb_fluxes = ClimaLand.turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) + turb_fluxes = ClimaLand.Canopy.canopy_turbulent_fluxes( + canopy.atmos, + canopy, + Y, + p, + t0, + ) @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux @test p.canopy.energy.shf == turb_fluxes.shf @@ -242,8 +248,7 @@ for FT in (Float32, Float64) h = FT(LP.planck_constant(earth_param_set)) N_a = FT(LP.avogadro_constant(earth_param_set)) _σ = FT(LP.Stefan(earth_param_set)) - (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR, ϵ_canopy) = - canopy.radiative_transfer.parameters + (; α_PAR_leaf, λ_γ_PAR, λ_γ_NIR) = canopy.radiative_transfer.parameters APAR = p.canopy.radiative_transfer.apar ANIR = p.canopy.radiative_transfer.anir energy_per_photon_PAR = h * c / λ_γ_PAR @@ -252,16 +257,17 @@ for FT in (Float32, Float64) @test p.canopy.radiative_transfer.SW_n == @. (energy_per_photon_PAR * N_a * APAR) + (energy_per_photon_NIR * N_a * ANIR) - + ϵ_canopy = p.canopy.radiative_transfer.effective_ϵ T_canopy = FT.(T_atmos(t0)) T_soil = FT.(soil_driver.T(t0)) ϵ_soil = FT.(soil_driver.ϵ) LW_d = FT.(longwave_radiation(t0)) - LW_d_canopy = (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 - LW_u_soil = ϵ_soil * _σ * T_soil^4 + (1 - ϵ_soil) * LW_d_canopy - @test Array(parent(p.canopy.radiative_transfer.LW_n))[1] ≈ - ϵ_canopy * LW_d - 2 * ϵ_canopy * _σ * T_canopy^4 + - ϵ_canopy * LW_u_soil + LW_d_canopy = @. (1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4 + LW_u_soil = @. ϵ_soil * _σ * T_soil^4 + (1 - ϵ_soil) * LW_d_canopy + @test p.canopy.radiative_transfer.LW_n == @. ( + ϵ_canopy * LW_d - 2 * ϵ_canopy * _σ * T_canopy^4 + + ϵ_canopy * LW_u_soil + ) # Penman-monteith @@ -715,7 +721,13 @@ for FT in (Float32, Float64) exp_tendency! = make_exp_tendency(canopy) dY = similar(Y) exp_tendency!(dY, Y, p, t0) - turb_fluxes = turbulent_fluxes(canopy.atmos, canopy, Y, p, t0) + turb_fluxes = ClimaLand.Canopy.canopy_turbulent_fluxes( + canopy.atmos, + canopy, + Y, + p, + t0, + ) @test p.canopy.hydraulics.fa.:1 == turb_fluxes.vapor_flux @test p.canopy.energy.lhf == turb_fluxes.lhf @test p.canopy.energy.shf == turb_fluxes.shf @@ -739,4 +751,233 @@ for FT in (Float32, Float64) ) .== FT(289), ) end + + + @testset "Zero LAI; FT = $FT" begin + domain = Point(; z_sfc = FT(0.0)) + + BeerLambertparams = BeerLambertParameters(FT) + # TwoStreamModel parameters + Ω = FT(0.69) + χl = FT(0.1) + G_Function = CLMGFunction(χl) + α_PAR_leaf = FT(0.1) + λ_γ_PAR = FT(5e-7) + λ_γ_NIR = FT(1.65e-6) + τ_PAR_leaf = FT(0.05) + α_NIR_leaf = FT(0.45) + τ_NIR_leaf = FT(0.25) + ϵ_canopy = FT(0.97) + TwoStreamparams = TwoStreamParameters( + FT; + Ω, + α_PAR_leaf, + τ_PAR_leaf, + α_NIR_leaf, + τ_NIR_leaf, + G_Function, + ) + photosynthesis_params = FarquharParameters(FT, C3()) + stomatal_g_params = MedlynConductanceParameters(FT) + + stomatal_model = MedlynConductanceModel{FT}(stomatal_g_params) + photosynthesis_model = FarquharModel{FT}(photosynthesis_params) + rt_models = ( + BeerLambertModel{FT}(BeerLambertparams), + TwoStreamModel{FT}(TwoStreamparams), + ) + energy_model = BigLeafEnergyModel{FT}(BigLeafEnergyParameters{FT}()) + earth_param_set = LP.LandParameters(FT) + LAI = FT(0.0) # m2 [leaf] m-2 [ground] + z_0m = FT(2.0) # m, Roughness length for momentum - value from tall forest ChatGPT + z_0b = FT(0.1) # m, Roughness length for scalars - value from tall forest ChatGPT + h_int = FT(30.0) # m, "where measurements would be taken at a typical flux tower of a 20m canopy" + shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}( + z_0m, + z_0b, + earth_param_set, + ) + lat = FT(0.0) # degree + long = FT(-180) # degree + + function zenith_angle( + t, + ref_time; + latitude = lat, + longitude = long, + insol_params = earth_param_set.insol_params, + ) + current_datetime = ref_time + Dates.Second(round(t)) + d, δ, η_UTC = + FT.( + Insolation.helper_instantaneous_zenith_angle( + current_datetime, + ref_time, + insol_params, + ) + ) + return Insolation.instantaneous_zenith_angle( + d, + δ, + η_UTC, + longitude, + latitude, + )[1] + end + + function shortwave_radiation( + t; + latitude = lat, + longitude = long, + insol_params = earth_param_set.insol_params, + ) + return 1000 # W/m^2 + end + + function longwave_radiation(t) + return 200 # W/m^2 + end + + u_atmos = t -> 10 #m.s-1 + + liquid_precip = (t) -> 0 # m + snow_precip = (t) -> 0 # m + T_atmos = t -> 290 # Kelvin + q_atmos = t -> 0.001 # kg/kg + P_atmos = t -> 1e5 # Pa + h_atmos = h_int # m + c_atmos = (t) -> 4.11e-4 # mol/mol + ref_time = DateTime(2005) + atmos = PrescribedAtmosphere( + TimeVaryingInput(liquid_precip), + TimeVaryingInput(snow_precip), + TimeVaryingInput(T_atmos), + TimeVaryingInput(u_atmos), + TimeVaryingInput(q_atmos), + TimeVaryingInput(P_atmos), + ref_time, + h_atmos, + earth_param_set; + c_co2 = TimeVaryingInput(c_atmos), + ) + radiation = PrescribedRadiativeFluxes( + FT, + TimeVaryingInput(shortwave_radiation), + TimeVaryingInput(longwave_radiation), + ref_time; + θs = zenith_angle, + ) + + # Plant Hydraulics + RAI = FT(0) + SAI = FT(0) + lai_fun = t -> LAI + ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}( + TimeVaryingInput(lai_fun), + SAI, + RAI, + ) + K_sat_plant = FT(1.8e-8) # m/s + ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value + Weibull_param = FT(4) # unitless, Holtzman's original c param value + a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity + plant_ν = FT(0.7) # m3/m3 + plant_S_s = FT(1e-2 * 0.0098) # m3/m3/MPa to m3/m3/m + conductivity_model = + PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param) + retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a) + root_depths = SVector{10, FT}(-(10:-1:1.0) ./ 10.0 * 2.0 .+ 0.2 / 2.0) # 1st element is the deepest root depth + function root_distribution(z::T) where {T} + return T(1.0 / 0.5) * exp(z / T(0.5)) # (1/m) + end + param_set = PlantHydraulics.PlantHydraulicsParameters(; + ai_parameterization = ai_parameterization, + ν = plant_ν, + S_s = plant_S_s, + root_distribution = root_distribution, + conductivity_model = conductivity_model, + retention_model = retention_model, + ) + Δz = FT(1.0) # height of compartments + n_stem = Int64(0) # number of stem elements + n_leaf = Int64(1) # number of leaf elements + compartment_centers = + FT.( + Vector( + range( + start = Δz / 2, + step = Δz, + stop = Δz * (n_stem + n_leaf) - (Δz / 2), + ), + ), + ) + compartment_faces = + FT.( + Vector( + range( + start = 0.0, + step = Δz, + stop = Δz * (n_stem + n_leaf), + ), + ) + ) + + ψ_soil0 = FT(0.0) + T_soil0 = FT(290) + soil_driver = PrescribedSoil( + root_depths, + (t) -> ψ_soil0, + (t) -> T_soil0, + FT(0.2), + FT(0.4), + FT(0.98), + ) + + plant_hydraulics = PlantHydraulics.PlantHydraulicsModel{FT}(; + parameters = param_set, + n_stem = n_stem, + n_leaf = n_leaf, + compartment_surfaces = compartment_faces, + compartment_midpoints = compartment_centers, + ) + autotrophic_parameters = AutotrophicRespirationParameters(FT) + autotrophic_respiration_model = + AutotrophicRespirationModel{FT}(autotrophic_parameters) + for rt_model in rt_models + canopy = ClimaLand.Canopy.CanopyModel{FT}(; + parameters = shared_params, + domain = domain, + radiative_transfer = rt_model, + photosynthesis = photosynthesis_model, + conductance = stomatal_model, + autotrophic_respiration = autotrophic_respiration_model, + energy = energy_model, + hydraulics = plant_hydraulics, + soil_driver = soil_driver, + atmos = atmos, + radiation = radiation, + ) + + Y, p, coords = ClimaLand.initialize(canopy) + Y.canopy.hydraulics .= plant_ν + Y.canopy.energy.T = FT(289) + set_initial_cache! = make_set_initial_cache(canopy) + t0 = FT(0.0) + set_initial_cache!(p, Y, t0) + + @test all(parent(p.canopy.hydraulics.fa.:1) .== FT(0)) + @test all(parent(p.canopy.energy.lhf) .== FT(0)) + @test all(parent(p.canopy.energy.shf) .== FT(0)) + @test all(parent(p.canopy.energy.fa_energy_roots) .== FT(0)) + @test all(parent(p.canopy.hydraulics.fa_roots) .== FT(0)) + @test all(parent(p.canopy.conductance.transpiration) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.LW_n) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.SW_n) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.apar) .== FT(0)) + @test all(parent(p.canopy.radiative_transfer.anir) .== FT(0)) + #@test all(parent(p.canopy.autotrophic_respiration.Ra) .== FT(0)) + end + + + end end