Merge pull request #11 from jialinl6/main

Update deposition
EarthSciML · Feb 16, 2024 · 49fe3ca · 49fe3ca
2 parents eb22944 + 83a45fc
commit 49fe3ca
Show file tree

Hide file tree

Showing 9 changed files with 249 additions and 113 deletions.
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -10,8 +10,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
-          - '1.6'
+          - '1.8'
+          - '1.9'
           - 'nightly'
         os:
           - ubuntu-latest

diff --git a/Project.toml b/Project.toml
@@ -4,8 +4,13 @@ authors = ["EarthSciML authors and contributors"]
 version = "0.1.0"
 
 [deps]
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+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"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -15,10 +20,13 @@ SafeTestsets = "0.0.1"
 StaticArrays = "1"
 Unitful = "1"
 julia = "1"
+EarthSciMLBase = "0.6"
+ModelingToolkit = "8"
+
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","SafeTestsets"]
+test = ["Test", "SafeTestsets"]
diff --git a/src/DepositionMTK.jl b/src/DepositionMTK.jl
@@ -2,6 +2,9 @@ module DepositionMTK
 
 using Unitful
 using StaticArrays
+using ModelingToolkit
+using EarthSciMLBase
+using IfElse
 
 include("wesley1989.jl")
 include("dry_deposition.jl")

diff --git a/src/dry_deposition.jl b/src/dry_deposition.jl
@@ -1,45 +1,35 @@
-export ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle
+export defaults, ra, mu, mfp, cc, vs, dParticle, dH2O, sc, stSmooth, stVeg, RbGas, z₀_table, A_table, α_table, γ_table, RbParticle, DryDepGas, DryDepParticle, DrydepositionG
 
-g = 9.81u"m*s^-2" # gravitational acceleration [m/s2]
-κ = 0.4 # von Karman constant
-k = 1.3806488e-23u"m^2*kg*s^-2/K" # Boltzmann constant
-M_air = 28.97e-3u"kg/mol" # molecular weight of air
-R = 8.3144621u"J/K/mol" # Gas constant
+@constants g = 9.81 [unit = u"m*s^-2", description = "gravitational acceleration"]
+@constants κ = 0.4 # von Karman constant
+@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K", description = "Boltzmann constant"]
+@constants M_air = 28.97e-3 [unit = u"kg/mol", description = "molecular weight of air"] 
+@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1", description = "Gas constant"] 
 
+@constants unit_m = 1 [unit = u"m"]
 """
 Function Ra calculates aerodynamic resistance to dry deposition 
 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]
-Based on Seinfeld and Pandis (2006) [Seinfeld and Pandis (2006)] equation 19.13 & 19.14.
+Based on Seinfeld and Pandis (2006) [Seinfeld, J.H. and Pandis, S.N. (2006) Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd Edition, John Wiley & Sons, New York.] 
+equation 19.13 & 19.14.
 """
-
 function ra(z, z₀, u_star, L)
-    if L == 0u"m"
-        ζ = 0
-        ζ₀= 0
-    else
-        ζ = z/L
-        ζ₀= z₀/L
-    end
-    print(ζ)
-    if 0 < ζ < 1
-        rₐ = 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀))
-    elseif ζ == 0
-        rₐ = 1/(κ*u_star)*log(z/z₀)
-    elseif -1 < ζ < 0
-        η₀=(1-15*ζ₀)^1/4
-        η =(1-15*ζ)^1/4
-        rₐ = 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))]
-    else
-        print("wrong rₐ")
-    end
+    ζ = IfElse.ifelse((L/unit_m == 0), 0, z/L)
+    ζ₀ = IfElse.ifelse((L/unit_m == 0), 0, z₀/L)
+    rₐ_1 = IfElse.ifelse((0 < ζ) & (ζ < 1), 1/(κ*u_star)*(log(z/z₀)+4.7*(ζ-ζ₀)), 1/(κ*u_star)*log(z/z₀))
+    η₀=(1-15*ζ₀)^1/4
+    η =(1-15*ζ)^1/4
+    rₐ = IfElse.ifelse((-1 < ζ) & (ζ < 0), 1/(κ*u_star)*[log(z/z₀)+log(((η₀^2+1)*(η₀+1)^2)/((η^2+1)*(η+1)^2))+2*(atan(η)-atan(η₀))][1], rₐ_1)
     return rₐ
 end
 
+@constants unit_T = 1 [unit = u"K"]
+@constants unit_convert_mu = 1 [unit = u"kg/m/s"]
 """
 Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K].
 """
 function mu(T)
-    return (1.8e-5*(T/298u"K")^0.85)u"kg/m/s"
+    return (1.458*10^-6*(T/unit_T)^(3/2)/((T/unit_T)+110.4))*unit_convert_mu
 end
 
 """
@@ -61,17 +51,15 @@ function cc(Dₚ,T,P,μ)
     return 1+2*λ/Dₚ*(1.257+0.4*exp(-1.1*Dₚ/(2*λ)))
 end
 
+@constants unit_v = 1 [unit = u"m/s"]
 """
 Function vs calculates the terminal setting velocity of a
 particle where Dp is particle diameter [m], ρₚ is particle density [kg/m3], Cc is the Cunningham slip correction factor, and μ is air dynamic viscosity [kg/(s m)]. 
 From equation 9.42 in Seinfeld and Pandis (2006).
 """
 function vs(Dₚ, ρₚ, Cc, μ)
-    if Dₚ > 20.e-6u"m"
-        print("Particle diameter ", Dₚ ," [m] is greater than 20um; Stokes settling no longer applies.")
-    else
-        return Dₚ^2*ρₚ*g*Cc/(18*μ)
-    end
+    IfElse.ifelse((Dₚ > 20.e-6*unit_m), 99999999*unit_v, Dₚ^2*ρₚ*g*Cc/(18*μ))
+    # Particle diameter Dₚ greater than 20um; Stokes settling no longer applies.
 end
 
 """
@@ -83,13 +71,14 @@ function dParticle(T,P,Dₚ,Cc,μ)
     return k*T*Cc/(3*pi*μ*Dₚ)
 end
 
+@constants T_unitless = 1 [unit = u"K^-1"]
+@constants unit_dH2O = 1 [unit = u"m^2/s"]
 """
 Function dH2O calculates molecular diffusivity of water vapor in air [m2/s] where T is temperature [K]
 using a regression fit to data in Bolz and Tuve (1976) found here: http://www.cambridge.org/us/engineering/author/nellisandklein/downloads/examples/EXAMPLE_9.2-1.pdf
-"""
+""" 
 function dH2O(T)
-    T_unitless = T*1u"K^-1"
-    return (-2.775e-6 + 4.479e-8*T_unitless + 1.656e-10*T_unitless^2)u"m^2/s"
+    return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O
 end
 
 """
@@ -174,18 +163,19 @@ where Sc is the dimensionless Schmidt number, u_star is the friction velocity [m
 Dp is particle diameter [m], and iSeason and iLandUse are season and land use indexes, respectively.
 From Seinfeld and Pandis (2006) equation 19.27.
 """
-
 function RbParticle(Sc, u_star, St, Dₚ, iSeason::Int, iLandUse::Int) 
     α = α_table[iLandUse]
     γ = γ_table[iLandUse]
-    A = (A_table[iSeason,iLandUse]*10^(-3))u"m"
+    A = (A_table[iSeason,iLandUse]*10^(-3))*unit_m
     R1 = exp(-St^0.5)
     term_1 = Sc^(-γ)
     term_2 = (St/(α+St))^2
     term_3 = 1/2*(Dₚ/A)^2
     return 1/(3*u_star*(term_1+term_2+term_3)*R1)
 end
 
+@constants G_unitless = 1 [unit = u"m^2/W"]
+@constants Rc_unit = 1 [unit = u"s/m"]
 """
 Function DryDepGas calculates dry deposition velocity [m/s] for a gas species,
 where z is the height of the surface layer [m], zo is roughness length [m], u_star is friction velocity [m/s], 
@@ -198,10 +188,10 @@ Based on Seinfeld and Pandis (2006) equation 19.2.
 function DryDepGas(z, z₀, u_star, L, ρA, gasData::GasData, G, Ts, θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) 
     Ra = ra(z, z₀, u_star, L)
     μ = mu(Ts)
-    Dg = dH2O(Ts)/gasData.Dh2oPerDx #Diffusivity of gas of interest [m2/s]
+    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/u"W*m^-2", (Ts/u"K"-273), θ, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool)u"s/m"
+    Rc = SurfaceResistance(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
 
@@ -218,13 +208,70 @@ function DryDepParticle(z, z₀, u_star, L, Dp, Ts, P, ρParticle, ρA, iSeason:
     Cc = cc(Dp, Ts, P, μ)
     Vs = vs(Dp, ρParticle, Cc, μ)
     if iLandUse == 4 # dessert
-        St = stSmooth(Vs, u_star, μ, ρA)
+       St = stSmooth(Vs, u_star, μ, ρA)
     else
-        St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))u"m")
+       St = stVeg(Vs, u_star, (A_table[iSeason,iLandUse]*10^(-3))*unit_m)
     end
     D = dParticle(Ts, P, Dp, Cc, μ)
     Sc = sc(μ, ρA, D)
     Rb = RbParticle(Sc, u_star, St, Dp, iSeason, iLandUse)
     return 1/(Ra+Rb+Ra*Rb*Vs)+Vs
 end
 
+defaults = [g => 9.81, κ => 0.4, k => 1.3806488e-23, M_air => 28.97e-3, R => 8.3144621, unit_T => 1, unit_convert_mu => 1, T_unitless => 1, unit_dH2O => 1, unit_m => 1, G_unitless => 1, Rc_unit => 1, unit_v => 1]
+
+# Add unit "ppb" to Unitful 
+module MyUnits
+using Unitful
+@unit ppb "ppb" Number 1 / 1000000000 false
+end
+Unitful.register(MyUnits)
+
+"""
+Description: This is a box model used to calculate the gas species concentration rate changed by dry deposition.
+Build Drydeposition model (gas)
+# Example
+``` julia
+	@parameters t 
+	d = DrydepositionG(t)
+```
+"""
+struct DrydepositionG <: EarthSciMLODESystem
+    sys::ODESystem
+    function DrydepositionG(t)
+        iSeason = 1
+        iLandUse = 10
+        rain = false
+        dew = false
+		@parameters z = 50 [unit = u"m"]
+		@parameters z₀ = 0.04 [unit = u"m"]
+		@parameters u_star = 0.44 [unit = u"m/s"]
+		@parameters L = 0 [unit = u"m"]
+        @parameters ρA = 1.2 [unit = u"kg*m^-3"]
+        @parameters G = 300 [unit = u"W*m^-2"]
+        @parameters T = 298 [unit = u"K"]
+        @parameters θ = 0
+        @parameters t [unit = u"s"]
+
+        D = Differential(t)
+
+        @variables SO2(t) [unit = u"ppb"]
+        @variables O3(t) [unit = u"ppb"]
+        @variables NO2(t) [unit = u"ppb"]
+        @variables NO(t) [unit = u"ppb"]
+        @variables H2O2(t) [unit = u"ppb"]
+        @variables CH2O(t) [unit = u"ppb"]
+
+        eqs = [
+            D(SO2) ~  -DryDepGas(z, z₀, u_star, L, ρA, So2Data, G, T, θ, iSeason, iLandUse, rain, dew, true, false) / z * SO2
+            D(O3) ~ -DryDepGas(z, z₀, u_star, L, ρA, O3Data, G, T, θ, iSeason, iLandUse, rain, dew, false, true) / z * O3
+            D(NO2) ~ -DryDepGas(z, z₀, u_star, L, ρA, No2Data, G, T, θ, iSeason, iLandUse, rain, dew, false, false) / z * NO2
+            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
+        ]
+
+        new(ODESystem(eqs, t, [SO2, O3, NO2, NO, H2O2, CH2O], [z, z₀, u_star, L, ρA, G, T, θ]; name=:DrydepositionG))
+    end
+end 
+
diff --git a/src/wesley1989.jl b/src/wesley1989.jl
@@ -88,12 +88,14 @@
 	# 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
-		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))))
+		# 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
 		# Adjust for dew and rain (from "Effects of dew and rain" section).
 		if rainOrDew
 			rs *= 3
@@ -112,7 +114,7 @@
 	end
 
 	#Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]).
-	function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat
+	function r_smx(r_s, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat
 		return r_s * Dh2oPerDx + r_mx
 	end
 
@@ -193,12 +195,16 @@
 		rac = r_ac[iSeason, iLandUse]
 
 		# Correction for cold temperatures from page 4 column 1.
-		if Ts < 0.0
-			correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark
-			rlux += correction
-			rclx += correction
-			rgsx += correction
-		end
+		correction = IfElse.ifelse((Ts < 0.0), 1000.0 * exp(-Ts-4), 0) # [s m-1] #mark
+		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.