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 6fbba8d5..02962898 100644
--- a/docs/src/dry_deposition.md
+++ b/docs/src/dry_deposition.md
@@ -1 +1,49 @@
-# 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 SO<sub>2</sub>, O<sub>3</sub>, NO<sub>2</sub>, NO, H<sub>2</sub>O<sub>2</sub> and CH<sub>2</sub>O change due to dry deposition. 
+
+We can create an instance of the model in the following manner:
+```julia @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.
+```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] # 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
+```julia @example 1
+using Plots
+plot(sol, xlabel="Time (second)", ylabel="concentration (ppb)", legend=:outerright)
+```
+
+## 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)")
+```
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