Update deposition

.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

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.4"
+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"]

src/DepositionMTK.jl

@@ -2,8 +2,12 @@ module DepositionMTK
 
 using Unitful
 using StaticArrays
+using ModelingToolkit
+using EarthSciMLBase
+using IfElse
 
 include("wesley1989.jl")
 include("dry_deposition.jl")
+include("wet_deposition.jl")
 
 end

src/dry_deposition.jl

@@ -1,45 +1,37 @@
-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"] # gravitational acceleration [m/s2]
+@constants κ = 0.4 # von Karman constant
+@constants k = 1.3806488e-23 [unit = u"m^2*kg*s^-2/K"] # Boltzmann constant
+@constants M_air = 28.97e-3 [unit = u"kg/mol"] # molecular weight of air
+@constants R = 8.3144621 [unit = u"kg*m^2*s^−2*K^-1*mol^-1"] # Gas constant
 
 """
 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.
 """
 
+@constants unit_m = 1 [unit = u"m"]
 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
 
 """
 Function mu calculates the dynamic viscosity of air [kg m-1 s-1] where T is temperature [K].
 """
+
+@constants unit_T = 1 [unit = u"K"]
+@constants unit_convert_mu = 1 [unit = u"kg/m/s"]
 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
+    #return (1.8e-5*(T/T₀)^0.85)*unit_convert_mu
 end
 
 """
@@ -66,12 +58,15 @@ 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).
 """
+# Particle diameter Dₚ greater than 20um; Stokes settling no longer applies.
+@constants unit_v = 1 [unit = u"m/s"]
 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*μ))
+    # 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
 end
 
 """
@@ -87,9 +82,12 @@ end
 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
 """
+
+@constants T_unitless = 1 [unit = u"K^-1"]
+@constants unit_dH2O = 1 [unit = u"m^2/s"] 
 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"
+    #T_unitless = T*1u"K^-1"
+    return (-2.775e-6 + 4.479e-8*T*T_unitless + 1.656e-10*(T*T_unitless)^2)*unit_dH2O
 end
 
 """
@@ -178,7 +176,7 @@ 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
@@ -195,13 +193,16 @@ irradiation [W m-2], Θ is the slope of the local terrain [radians], iSeason and
 dew and rain indicate whether there is dew or rain on the ground, and isSO2 and isO3 indicate whether the gas species of interest is either SO2 or O3, respectively. 
 Based on Seinfeld and Pandis (2006) equation 19.2.
 """
+
+@constants G_unitless = 1 [unit = u"m^2/W"]
+@constants Rc_unit = 1 [unit = u"s/m"]
 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]
     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 +219,71 @@ 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 
+

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.

src/wet_deposition.jl

@@ -0,0 +1,86 @@
+export WetDeposition, Wetdeposition, wd_defaults
+
+"""
+Calculate wet deposition based on formulas at
+www.emep.int/UniDoc/node12.html.
+Inputs are fraction of grid cell covered by clouds (cloudFrac),
+rain mixing ratio (qrain), air density (ρ_air [kg/m3]),
+and fall distance (Δz [m]).
+Outputs are wet deposition rates for PM2.5, SO2, and other gases
+(wdParticle, wdSO2, and wdOtherGas [1/s]).
+"""
+
+@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient
+@constants ρwater = 1000.0 [unit = u"kg*m^-3"]  # kg/m3
+@constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s
+
+function WetDeposition(cloudFrac, qrain, ρ_air, Δz)
+    #@constants A_wd = 5.2 [unit = u"m^3/kg/s"] # m3 kg-1 s-1; Empirical coefficient
+	E = 0.1           # size-dependent collection efficiency of aerosols by the raindrops
+	wSubSO2 = 0.15   # sub-cloud scavanging ratio
+	wSubOther = 0.5  # sub-cloud scavanging ratio
+	wInSO2 = 0.3     # in-cloud scavanging ratio
+	wInParticle = 1.0 # in-cloud scavanging ratio
+	wInOther = 1.4   # in-cloud scavanging ratio
+	# @constants ρwater = 1000.0 [unit = u"kg*m^-3"]  # kg/m3
+	# @constants Vdr = 5.0 [unit = u"m/s"] # raindrop fall speed, m/s
+
+    # precalculated constant combinations
+	AE = A_wd * E
+	wSubSO2VdrPerρwater = wSubSO2 * Vdr / ρwater
+	wSubOtherVdrPerρwater = wSubOther * Vdr / ρwater
+	wInSO2VdrPerρwater = wInSO2 * Vdr / ρwater
+	wInParticleVdrPerρwater = wInParticle * Vdr / ρwater
+	wInOtherVdrPerρwater = wInOther * Vdr / ρwater
+
+    # wdParticle (subcloud) = A * P / Vdr * E; P = QRAIN * Vdr * ρgas => wdParticle = A * QRAIN * ρgas * E
+	# wdGas (subcloud) = wSub * P / Δz / ρwater = wSub * QRAIN * Vdr * ρgas / Δz / ρwater
+	# wd (in-cloud) = wIn * P / Δz / ρwater = wIn * QRAIN * Vdr * ρgas / Δz / ρwater
+
+    wdParticle = qrain * ρ_air * (AE + cloudFrac * (wInParticleVdrPerρwater / Δz))
+	wdSO2 = (wSubSO2VdrPerρwater + cloudFrac*wSubSO2VdrPerρwater) * qrain * ρ_air / Δz
+	wdOtherGas = (wSubOtherVdrPerρwater + cloudFrac*wSubOtherVdrPerρwater) * qrain * ρ_air / Δz
+
+	return wdParticle, wdSO2, wdOtherGas
+end
+wd_defaults = [A_wd => 5.2, ρwater => 1000.0, Vdr => 5.0]
+
+# 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 concentration rate changed by wet deposition.
+Build Wetdeposition model
+# Example
+``` julia
+	@parameters t 
+	wd = Wetdeposition(t)
+```
+"""
+
+struct Wetdeposition <: EarthSciMLODESystem
+    sys::ODESystem
+    function Wetdeposition(t)
+		@parameters cloudFrac = 0.5
+		@parameters qrain = 0.5
+		@parameters ρ_air = 1.204 [unit = u"kg*m^-3"]
+		@parameters Δz = 200 [unit = u"m"]
+        @parameters t [unit = u"s"]
+
+        D = Differential(t)
+
+        @variables SO2(t) [unit = u"ppb"]
+        @variables O3(t) [unit = u"ppb"]
+
+        eqs = [
+            D(SO2) ~  -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[2] * SO2
+            D(O3) ~ -WetDeposition(cloudFrac, qrain, ρ_air, Δz)[3] * O3
+        ]
+
+        new(ODESystem(eqs, t, [SO2, O3], [cloudFrac, qrain, ρ_air, Δz]; name=:Wetdeposition))
+    end
+end