diff --git a/docs/Project.toml b/docs/Project.toml index 9bfbf718..c75f724b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,5 +4,5 @@ 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" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" 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 fc2bcfec..15fbc70a 100644 --- a/docs/src/dry_deposition.md +++ b/docs/src/dry_deposition.md @@ -3,7 +3,7 @@ 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 SO₂, O₃, NO₂, NO, H₂O₂ and CH₂O change due to dry deposition. +Here's an example of how concentration of different species, such as SO₂, O₃, NO₂, NO, H₂O₂, CH₂O and HNO₃change due to dry deposition. We can create an instance of the model in the following manner: ```julia @@ -11,16 +11,16 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = DrydepositionG(t) +model = DrydepositionG() ``` Before running any simulations with the model we need to convert it into a system of differential equations. ```julia sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O +u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -29,14 +29,14 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = DrydepositionG(t) +model = DrydepositionG() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,5,2.34,0.15] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O +u0 = [2.0,10.0,5,5,2.34,0.15,10] # initial concentrations of SO₂, O₃, NO₂, NO, H₂O₂, CH₂O, HNO₃ sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -49,16 +49,16 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig ## Parameters The parameters in the model are: ```julia -parameters(sys) # [z, z₀, u_star, L, ρA, G, T, θ] +parameters(sys) # [iLandUse, z, z₀, u_star, L, ρA, G, iSeason, 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]. +where ```iSeason``` and ```iLandUse``` are the index for the season and landuse,```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```. ```@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] +p1 = [10, 50,0.04,0.44,0,1.2,300,1,298,0] +p2 = [10, 10,0.04,0.44,0,1.2,300,1,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) diff --git a/docs/src/emep.md b/docs/src/emep.md index 4a9205c2..463b3027 100644 --- a/docs/src/emep.md +++ b/docs/src/emep.md @@ -8,17 +8,17 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = Wetdeposition(t) +model = Wetdeposition() ``` Before running any simulations with the model we need to convert it into a system of differential equations. ```julia sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,1400,275,50,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP +u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O prob = ODEProblem(sys, u0, tspan, []) sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` @@ -28,14 +28,14 @@ using AtmosphericDeposition using ModelingToolkit using DifferentialEquations using EarthSciMLBase -using Unitful +using DynamicQuantities +using ModelingToolkit:t -@parameters t [unit = u"s", description="Time"] -model = Wetdeposition(t) +model = Wetdeposition() sys = structural_simplify(model) tspan = (0.0, 3600*24) -u0 = [2.0,10.0,5,1400,275,50,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP +u0 = [2.0,10.0,5,1400,275,50,0.15,2.34,10,0.15] # initial concentration of SO₂, O₃, NO₂, CH₄, CO, DMS, ISOP, H₂O₂, HNO₃, CH₂O prob = ODEProblem(sys, u0, tspan, []) sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters ``` diff --git a/ext/EarthSciDataExt.jl b/ext/EarthSciDataExt.jl index e7bc628c..e5d19ef3 100644 --- a/ext/EarthSciDataExt.jl +++ b/ext/EarthSciDataExt.jl @@ -13,7 +13,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, gg = 9.81, [unit = u"m*s^-2", description="gravitational acceleration"], ) - # ρA in DrydepositionG(t) are in units of "kg/m3". + # ρA in DrydepositionG() are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρA = P*M/(R*T)*kgperg @@ -40,7 +40,7 @@ function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g Vdr = 5.0, [unit = u"m/s", description="droplet velocity"], ) - # ρ_air in Wetdeposition(t) are in units of "kg/m3". + # ρ_air in Wetdeposition() are in units of "kg/m3". # ρ = P*M/(R*T)= Pa*(g/mol)/(m3*Pa/mol/K*K) = g/m3 # Overall, ρ_air = P*M/(R*T)*kgperg diff --git a/ext/GasChemExt.jl b/ext/GasChemExt.jl index be976bf6..6add2acf 100644 --- a/ext/GasChemExt.jl +++ b/ext/GasChemExt.jl @@ -26,6 +26,8 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos c.CO => d.CO, c.DMS => d.DMS, c.ISOP => d.ISOP, + c.H2O2 => d.H2O2, + c.CH2O => d.CH2O, )) end diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl index 1f18c4ad..7adff712 100644 --- a/src/dry_deposition.jl +++ b/src/dry_deposition.jl @@ -260,6 +260,7 @@ function DrydepositionG(; name=:DrydepositionG) NO(t), [unit = u"ppb"], H2O2(t), [unit = u"ppb"], CH2O(t), [unit = u"ppb"], + HNO3(t), [unit = u"ppb"], ) eqs = [ @@ -269,6 +270,7 @@ function DrydepositionG(; name=:DrydepositionG) D(NO) ~ -DryDepGas(z, z₀, u_star, L, ρA, NoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO D(H2O2) ~ -DryDepGas(z, z₀, u_star, L, ρA, H2o2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * H2O2 D(CH2O) ~ -DryDepGas(z, z₀, u_star, L, ρA, HchoData, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * CH2O + D(HNO3) ~ -DryDepGas(z, z₀, u_star, L, ρA, Hno3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * HNO3 ] ODESystem(eqs, t, vars, params; name=name, diff --git a/src/wet_deposition.jl b/src/wet_deposition.jl index a19a8859..83af024e 100644 --- a/src/wet_deposition.jl +++ b/src/wet_deposition.jl @@ -72,6 +72,9 @@ function Wetdeposition(; name=:Wetdeposition) CO(t), [unit = u"ppb"], DMS(t), [unit = u"ppb"], ISOP(t), [unit = u"ppb"], + H2O2(t), [unit = u"ppb"], + HNO3(t), [unit = u"ppb"], + CH2O(t), [unit = u"ppb"], ) eqs = [ @@ -82,6 +85,9 @@ function Wetdeposition(; name=:Wetdeposition) D(CO) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CO D(DMS) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * DMS D(ISOP) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * ISOP + D(H2O2) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * H2O2 + D(HNO3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * HNO3 + D(CH2O) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * CH2O ] ODESystem(eqs, t, vars, params; name=name, diff --git a/test/connector_test.jl b/test/connector_test.jl index 264ce166..20badd7b 100644 --- a/test/connector_test.jl +++ b/test/connector_test.jl @@ -7,7 +7,9 @@ using ModelingToolkit:t composed_ode = couple(SuperFast(), FastJX(), DrydepositionG(), Wetdeposition()) combined_mtk = convert(ODESystem, composed_ode) sys = structural_simplify(combined_mtk) - @test length(unknowns(sys)) ≈ 18 + print(unknowns(sys)) + @test length(unknowns(sys)) ≈ 20 + #TODO Change 20 to 18 after the latest GasChem package include species HNO3 eqs = string(equations(sys)) wanteqs = ["Differential(t)(SuperFast₊O3(t)) ~ SuperFast₊DrydepositionG_ddt_O3ˍt(t) + SuperFast₊Wetdeposition_ddt_O3ˍt(t)"] @@ -23,7 +25,8 @@ end model = couple(SuperFast(), FastJX(), geosfp, Wetdeposition(), DrydepositionG()) sys = structural_simplify(convert(ODESystem, model)) - @test length(unknowns(sys)) ≈ 18 + @test length(unknowns(sys)) ≈ 20 + #TODO Change 20 to 18 after the latest GasChem package include species HNO3 eqs = string(equations(convert(ODESystem, model))) wanteq = "DrydepositionG₊G(t) ~ GEOSFP₊A1₊SWGDN(t)"