From 4d850a4335f8d72348bbc92dfde870f3d92119b0 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Wed, 21 Feb 2024 19:20:23 -0600 Subject: [PATCH 1/3] Change to WesleySurfaceResistance and delete commented out code --- src/dry_deposition.jl | 2 +- src/wesley1989.jl | 19 +++---------------- test/wesely1989_test.jl | 8 ++++---- 3 files changed, 8 insertions(+), 21 deletions(-) diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 9b1a0945..000256b1 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -191,7 +191,7 @@ function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason Dg = dH2O(Ts) / gasData.Dh2oPerDx # Diffusivity of gas of interest [m2/s] Sc = sc(μ, ρA, Dg) Rb = RbGas(Sc, u_star) - Rc = SurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit + Rc = WesleySurfaceResistance(gasData, G * G_unitless, (Ts * T_unitless - 273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) * Rc_unit return 1 / (Ra + Rb + Rc) end diff --git a/src/wesley1989.jl b/src/wesley1989.jl index 191cfc57..17d50344 100644 --- a/src/wesley1989.jl +++ b/src/wesley1989.jl @@ -1,4 +1,4 @@ -export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, inf, r_i, r_lu, r_ac, r_gsS, r_gsO, r_clS, r_clO, So2Data, O3Data, No2Data, NoData, Hno3Data, H2o2Data, GasData, AldData, HchoData, OpData, PaaData, OraData, Nh3Data, PanData, Hno2Data +export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, WesleySurfaceResistance, inf, r_i, r_lu, r_ac, r_gsS, r_gsO, r_clS, r_clO, So2Data, O3Data, No2Data, NoData, Hno3Data, H2o2Data, GasData, AldData, HchoData, OpData, PaaData, OraData, Nh3Data, PanData, Hno2Data const inf = 1.e25 @@ -88,14 +88,7 @@ const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid # Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew. function r_s(G, Ts, iSeason::Int, iLandUse::Int, rainOrDew::Bool) rs = 0.0 - rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * - (400.0 * 1.0 / (Ts * (40.0 - Ts)))) - # if Ts >= 39.9 || Ts <= 0.1 - # rs = inf - # else - # rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * - # (400.0 * 1.0 / (Ts * (40.0 - Ts))) - # end + rs = IfElse.ifelse((Ts >= 39.9) & (Ts <= 0.1), inf, r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) * (400.0 * 1.0 / (Ts * (40.0 - Ts)))) # Adjust for dew and rain (from "Effects of dew and rain" section). if rainOrDew rs *= 3 @@ -171,7 +164,7 @@ end # Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3). # From Wesely (1989) regarding rain and dew inputs: "A direct computation of the surface wetness would be most desirable, e.g. by estimating the amount of free surface water accumulated and then evaporated. Alternatively, surface relative humidity might be a useful index. After dewfall and rainfall events are completed, surface wetness often disappears as a result of evaporation after approximately 2 hours of good atmospheric mixing, the period of time recommended earlier (Sheih et al., 1986)". -function SurfaceResistance(gasData::GasData, G, Ts, θ, +function WesleySurfaceResistance(gasData::GasData, G, Ts, θ, iSeason, iLandUse, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) rs = r_s(G, Ts, iSeason, iLandUse, rain || dew) rmx = r_mx(gasData.Hstar, gasData.Fo) @@ -199,12 +192,6 @@ function SurfaceResistance(gasData::GasData, G, Ts, θ, rlux += correction rclx += correction rgsx += correction - # if Ts < 0.0 - # correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark - # rlux += correction - # rclx += correction - # rgsx += correction - # end r_c = 1.0 / (1.0 / (rsmx) + 1.0 / rlux + 1.0 / (rdc + rclx) + 1.0 / (rac + rgsx)) r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces. diff --git a/test/wesely1989_test.jl b/test/wesely1989_test.jl index c3995a00..2738e588 100644 --- a/test/wesely1989_test.jl +++ b/test/wesely1989_test.jl @@ -122,20 +122,20 @@ function TestWesely() for iSeason in 1:5 for ig in 1:5 G = Garr[ig] - r_c = SurfaceResistance(gasData[i], G, Ts[iSeason], θ, + r_c = WesleySurfaceResistance(gasData[i], G, Ts[iSeason], θ, iSeason, iLandUse, false, false, isSO2, isO3) if different(r_c, polData[iSeason, ig]) println(pol, iSeason, G, r_c, polData[iSeason, ig]) return false end end - r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, + r_c = WesleySurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, iSeason, iLandUse, false, true, isSO2, isO3) # dew if different(r_c, polData[iSeason, 6]) println(pol, iSeason, "dew", r_c, polData[iSeason, 6]) return false end - r_c = SurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, + r_c = WesleySurfaceResistance(gasData[i], 0.0, Ts[iSeason], θ, iSeason, iLandUse, true, false, isSO2, isO3) # rain if different(r_c, polData[iSeason, 7]) println(pol, iSeason, "rain", r_c, polData[iSeason, 7]) @@ -155,5 +155,5 @@ end @test r_lux(1.0, 1.0, 1, 1, true, false, true, false) ≈ 50 @test r_clx(1.0, 1.0, 1, 1) ≈ 9.999899563027895e24 @test r_gsx(1.0, 1.0, 1, 1) ≈ 299.9977500168749 - @test SurfaceResistance(So2Data, 1.0, 1.0, 1.0, 1, 1, true, true, true, false) ≈ 45.45454545454546 + @test WesleySurfaceResistance(So2Data, 1.0, 1.0, 1.0, 1, 1, true, true, true, false) ≈ 45.45454545454546 end From 5c7084cf3726b6d836318ce721857c708d4e8433 Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Wed, 21 Feb 2024 23:10:41 -0600 Subject: [PATCH 2/3] Add dry_deposition docs --- docs/src/dry_deposition.md | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index 6fbba8d5..a81a107a 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -1 +1,31 @@ -# Dry Deposition \ No newline at end of file +# Dry Deposition +## Model +This is an implementation of a box model used to calculate changes in gas species concentration due to dry deposition. + +## Running the model +Here's an example of how concentration of different species, such as SO2, O3, NO2, NO, H2O2 and CH2O change due to dry deposition. + +We can create an instance of the model in the following manner: +```@example 1 +using AtmosphericDeposition +using ModelingToolkit +using DifferentialEquations +using EarthSciMLBase +using Unitful + +@parameters t [unit = u"s", description="Time"] + +model = DrydepositionG(t) +``` +Before running any simulations with the model we need to convert it into a system of differential equations. +```@example 1 +sys = structural_simplify(get_mtk(model)) +tspan = (0.0, 3600*24) +u0 = [2.0,10.0,5,5,2.34,0.15] +sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) +``` +which we can plot as +```@example 1 +using Plots +plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerright) +``` \ No newline at end of file From c3f69678fa005fc57dcb8d4cef55b53115ba571f Mon Sep 17 00:00:00 2001 From: jialinl6 Date: Thu, 22 Feb 2024 12:25:42 -0600 Subject: [PATCH 3/3] Add dry deposition documentation and update Project.toml --- docs/Project.toml | 5 +++++ docs/src/dry_deposition.md | 32 +++++++++++++++++++++++++------- 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 479f84cb..9bfbf718 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,8 @@ [deps] AtmosphericDeposition = "1a52f20c-0d16-41d8-a00a-b2996d86a462" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" \ No newline at end of file diff --git a/docs/src/dry_deposition.md b/docs/src/dry_deposition.md index a81a107a..02962898 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -6,7 +6,7 @@ This is an implementation of a box model used to calculate changes in gas specie Here's an example of how concentration of different species, such as SO2, O3, NO2, NO, H2O2 and CH2O change due to dry deposition. We can create an instance of the model in the following manner: -```@example 1 +```julia @example 1 using AtmosphericDeposition using ModelingToolkit using DifferentialEquations @@ -14,18 +14,36 @@ using EarthSciMLBase using Unitful @parameters t [unit = u"s", description="Time"] - model = DrydepositionG(t) ``` Before running any simulations with the model we need to convert it into a system of differential equations. -```@example 1 +```julia @example 1 sys = structural_simplify(get_mtk(model)) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15] -sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) +u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O +sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` which we can plot as -```@example 1 +```julia @example 1 using Plots plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerright) -``` \ No newline at end of file +``` + +## Parameters +The parameters in the model are: +```julia @example 1 +parameters(sys) # [z, z₀, u_star, L, ρA, G, T, θ] +``` +where ```z``` is the top of the surface layer [m], ```z₀``` is the roughness length [m], ```u_star``` is friction velocity [m/s], and ```L``` is Monin-Obukhov length [m], ```ρA``` is air density [kg/m3], ```T``` is surface air temperature [K], ```G``` is solar irradiation [W m-2], ```Θ``` is the slope of the local terrain [radians]. + +Let's run some simulation with different value for parameter ```z```. +```julia @example 1 +@unpack O3 = sys + +p1 = [50,0.04,0.44,0,1.2,300,298,0] +p2 = [10,0.04,0.44,0,1.2,300,298,0] +sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0) +sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0) + +plot([sol1[O3],sol2[O3]], label = ["z=50m" "z=10m"], title = "Change of O3 concentration due to dry deposition", xlabel="Time (second)", ylabel="concentration (ppb)") +```