MTKv9 Compatibility & Species Adjustments Following GEOS-CHEM

Project.toml

@@ -4,16 +4,15 @@ authors = ["EarthSciML authors and contributors"]
 version = "0.2.1"
 
 [deps]
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
-IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
 EarthSciData = "a293c155-435f-439d-9c11-a083b6b47337"
@@ -25,14 +24,12 @@ GasChemExt = "GasChem"
 
 [compat]
 DifferentialEquations = "7"
-EarthSciData = "0.9"
-EarthSciMLBase = "0.15"
-GasChem = "0.6"
-IfElse = "0.1"
-ModelingToolkit = "8"
+EarthSciData = "0.10"
+EarthSciMLBase = "0.19"
+GasChem = "0.7"
+ModelingToolkit = "9"
 SafeTestsets = "0.1"
 StaticArrays = "1"
-Unitful = "1"
 julia = "1"
 
 [extras]

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"

docs/src/dry_deposition.md

@@ -3,41 +3,41 @@
 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
 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
-sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
+tspan = (0.0, 3600*24)
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
+sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters 
 ```
 
 ```@setup 1
 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
-sol = solve(ODEProblem(sys, u0, tspan, []),AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
+sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters 
 ```
 
 which we can plot as
@@ -49,18 +49,18 @@ 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 ```iLandUse``` is the index for the landuse, ```iSeason``` is the index for the season, ```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
+@unpack O3,z = 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)
+p1 = [z=>50]
+p2 = [z=>10]
+sol1 = solve(ODEProblem(sys, [], tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol2 = solve(ODEProblem(sys, [], 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)")
 ```

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,15 +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
-prob = ODEProblem(sys, u0, tspan, [])
+prob = ODEProblem(sys, [], tspan, []) # default initial concentration of SO₂, O₃, NO₂, H₂O₂, HNO₃, CH₂O
 sol = solve(prob,AutoTsit5(Rosenbrock23()), saveat=10.0) # default parameters
 ```
 
@@ -49,31 +48,31 @@ plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerrig
 ## Parameters
 The parameters in the model are:
 ```julia @example 1
-parameters(sys) # [cloudFrac, qrain, ρ_air, Δz]
+parameters(sys) #[ρ_air, qrain, Δz, cloudFrac]
 ```
 where ```cloudFrac``` is fraction of grid cell covered by clouds, ```qrain``` is rain mixing ratio, ```ρ_air``` is air density [kg/m3], and ```Δz``` is fall distance [m].
 
 Let's run some simulation with different value for parameter ```cloudFrac```. 
 ```@example 1
-@unpack O3 = sys
+@unpack O3,cloudFrac,qrain = sys
 
-p1 = [0.3,0.5,1.204,200]
-p2 = [0.6,0.5,1.204,200]
-sol1 = solve(ODEProblem(sys, u0, tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
-sol2 = solve(ODEProblem(sys, u0, tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
+p1 = [cloudFrac=>0.3]
+p2 = [cloudFrac=>0.6]
+sol1 = solve(ODEProblem(sys, [], tspan, p1),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol2 = solve(ODEProblem(sys, [], tspan, p2),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
 plot([sol1[O3],sol2[O3]], label = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```
 From the plot we could see that with larger cloud fraction, the wet deposition rate increases. 
 
 Let's run some simulation with different value for parameter ```qrain``` 
 ```@example 1
-p3 = [0.5,0.3,1.204,200]
-p4 = [0.5,0.6,1.204,200]
-sol3 = solve(ODEProblem(sys, u0, tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0)
-sol4 = solve(ODEProblem(sys, u0, tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0)
+p3 = [qrain=>0.3]
+p4 = [qrain=>0.6]
+sol3 = solve(ODEProblem(sys, [], tspan, p3),AutoTsit5(Rosenbrock23()), saveat=10.0)
+sol4 = solve(ODEProblem(sys, [], tspan, p4),AutoTsit5(Rosenbrock23()), saveat=10.0)
 
-plot([sol3[O3],sol4[O3]], label = ["cloudFrac=0.3" "cloudFrac=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
+plot([sol3[O3],sol4[O3]], label = ["qrain=0.3" "qrain=0.6"], title = "Change of O3 concentration due to wet deposition", xlabel="Time (second)", ylabel="concentration (ppb)")
 ```
 The graph indicates that an increase in the rain mixing ratio leads to a corresponding rise in the rate of wet deposition.
 

docs/src/wesley1989.md

@@ -39,13 +39,19 @@ a. For season index ```iSeason```, there're five seasonal categories
     4. Winter, snow on ground
     5. Transitional 
 
-b. For land use index ```iLandUse```, there're five land use categories
+b. For land use index ```iLandUse```, there're eleven land use categories
 
-    1. Evergreen–needleleaf trees
-    2. Deciduous broadleaf trees
-    3. Grass
-    4. Desert
-    5. Shrubs and interrupted woodlands 
+    1. Urban land 
+    2. agricultural land
+    3. range land
+    4. deciduous forest
+    5. coniferous forest 
+    6. mixed forest including wetland 
+    7. water, both salt and fresh 
+    8. barren land, mostly desert
+    9. nonforested wetland
+    10. mixed agricultural and range land
+    11. rocky open areas with low-growing shrubs
 
 c. For gasData ```gasData```, we include the following species: 
 

ext/EarthSciDataExt.jl

@@ -1,50 +1,58 @@
 module EarthSciDataExt
 
-using AtmosphericDeposition, EarthSciData, EarthSciMLBase, Unitful, ModelingToolkit
+using AtmosphericDeposition, EarthSciData, EarthSciMLBase, DynamicQuantities, ModelingToolkit
 
 function EarthSciMLBase.couple2(d::AtmosphericDeposition.DrydepositionGCoupler, g::EarthSciData.GEOSFPCoupler)
     d, g = d.sys, g.sys
 
     @constants(
-        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
-        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
-        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+        vK = 0.4, [description = "von Karman's constant"],
+        Cp = 1000, [unit = u"W*s/kg/K", description="specific heat at constant pressure"],
+        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
 
-    d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA)
+    # Monin-Obhukov length = -Air density * Cp * T(surface air) * Ustar^3/（vK   * g  * Sensible Heat flux）
+
+    d = param_to_var(d, :T, :z, :z₀, :u_star, :G, :ρA, :L)
+
     ConnectorSystem([
             d.T ~ g.I3₊T,
             d.z ~ g.A1₊PBLH,
             d.z₀ ~ g.A1₊Z0M,
             d.u_star ~ g.A1₊USTAR,
-            d.G ~ g. A1₊SWGDN,
-            d.ρA ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+            d.G ~ g.A1₊SWGDN,
+            d.ρA ~ g.P/(g.I3₊T*R)*MW_air,
+            d.L ~ -g.P/(g.I3₊T*R)*MW_air * Cp * g.A1₊TS * (g.A1₊USTAR)^3/(vK*gg*g.A1₊HFLUX),
         ], d, g)
 end
 
 function EarthSciMLBase.couple2(d::AtmosphericDeposition.WetdepositionCoupler, g::EarthSciData.GEOSFPCoupler)
     d, g = d.sys, g.sys
 
     @constants(
-        PaPerhPa = 100, [unit = u"Pa/hPa", description="Conversion factor from hPa to Pa"],
-        MW_air = 29, [unit = u"g/mol", description="dry air molar mass"],
-        kgperg = 1e-3, [unit = u"kg/g", description="Conversion factor from g to kg"],
+        MW_air = 0.029, [unit = u"kg/mol", description="dry air molar mass"],
         R = 8.31446261815324, [unit = u"m^3*Pa/mol/K", description="Ideal gas constant"],
+        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
 
-    d = param_to_var(d, :cloudFrac, :ρ_air)
+    # From EMEP algorithm: P = QRAIN * Vdr * ρgas => QRAIN = P / Vdr / ρgas
+    # kg*m-2*s-1/(m/s)/(kg/m3)
+
+    d = param_to_var(d, :cloudFrac, :ρ_air, :qrain)
     ConnectorSystem([
             d.cloudFrac ~ g.A3cld₊CLOUD,
-            d.ρ_air ~ g.P * PaPerhPa/(g.I3₊T*R)*kgperg*MW_air,
+            d.ρ_air ~ g.P/(g.I3₊T*R)*MW_air,
+            d.qrain ~ (g.A3mstE₊PFLCU + g.A3mstE₊PFLLSAN) / Vdr / (g.P/(g.I3₊T*R)*MW_air),
         ], d, g)
 end
 

ext/GasChemExt.jl

@@ -1,6 +1,6 @@
 module GasChemExt
 
-using AtmosphericDeposition, GasChem, EarthSciMLBase, Unitful, ModelingToolkit
+using AtmosphericDeposition, GasChem, EarthSciMLBase, DynamicQuantities, ModelingToolkit
 
 function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDeposition.DrydepositionGCoupler)
     c, d = c.sys, d.sys
@@ -9,7 +9,6 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos
         c.SO2 => d.SO2,
         c.NO2 => d.NO2,
         c.O3 => d.O3,
-        c.NO => d.NO,
         c.H2O2 => d.H2O2,
         c.CH2O => d.CH2O,
     ))
@@ -22,10 +21,8 @@ function EarthSciMLBase.couple2(c::GasChem.SuperFastCoupler, d::AtmosphericDepos
         c.SO2 => d.SO2,
         c.NO2 => d.NO2,
         c.O3 => d.O3,
-        c.CH4 => d.CH4,
-        c.CO => d.CO,
-        c.DMS => d.DMS,
-        c.ISOP => d.ISOP,
+        c.H2O2 => d.H2O2,
+        c.CH2O => d.CH2O,
     ))
 end
 

src/AtmosphericDeposition.jl

@@ -1,10 +1,13 @@
 module AtmosphericDeposition
 
-using Unitful
+using DynamicQuantities
 using StaticArrays
 using ModelingToolkit
 using EarthSciMLBase
-using IfElse
+using DataInterpolations
+using ModelingToolkit:t
+
+@register_unit ppb 1
 
 include("wesley1989.jl")
 include("dry_deposition.jl")