Lecture 5: Wet Deposition In PATMO

The PATMO paper by Danielache et al. (2023) describes PATMO as a flexible one-dimensional atmospheric chemistry code. In that framework, atmospheric composition is controlled not only by chemistry and photolysis, but also by boundary and removal processes such as surface emission, dry deposition, wet deposition, condensation, sedimentation, diffusion, and escape.

This lecture focuses on the wet-deposition part of that physical picture and connects it to the actual PATMO code used in the modern sulfur cycle case.

The original PATMO ODE code applies wet deposition like this:

! Wet Deposition
do j = 12, 2, -1
    do i = 1, chemSpeciesNumber
        dn(j,     i) = dn(j,     i) - wetdep(j, i) * n(j, i)
        dn(j - 1, i) = dn(j - 1, i) + wetdep(j, i) * n(j, i)
    end do
end do

do i = 1, chemSpeciesNumber
    dn(1, i) = dn(1, i) - wetdep(1, i) * n(1, i)
end do

In the code, wet deposition is a first-order process. For species i in layer j:

\[ \left(\frac{dn_{j,i}}{dt}\right)_{\mathrm{wet}} = -k_{\mathrm{rain},j,i} n_{j,i} \]

In PATMO notation, \(k_{\mathrm{rain},j,i}\) is stored as wetdep(j,i). The unit is \(\mathrm{s^{-1}}\). The product wetdep(j,i) * n(j,i) has units of molecules cm\(^{-3}\) s\(^{-1}\).

In tests/modern_sulfur_cycle/test.f90, PATMO computes rainout for selected sulfur species by passing an effective Henry constant:

wetdep(:,:) = 0d0

call computewetdep(1,  2.0d-2)  ! OCS
call computewetdep(11, 5.0d-2)  ! CS2
call computewetdep(28, 1.0d-1)  ! H2S
call computewetdep(18, 4.0d3)   ! SO2
call computewetdep(19, 5.0d14)  ! SO4

This calculation answers one question: how fast should species i be removed by rainout in layer j? The answer is a first-order rate constant called K_RAIN.

Read the calculation in three steps: estimate contact with water, correct for precipitation timing, then convert the result into a first-order loss rate.

First, computewetdep estimates a raw layer-dependent removal frequency \(r_{j,i}\). This is not the final rainout rate yet. It is an intermediate estimate of how quickly species i can be transferred into water in layer j:

\[ r_{j,i} = \frac{W_{\mathrm{H2O},j}/55}{N_A w_l 10^{-9} + \left(H_{\mathrm{eff},i} R T_j\right)^{-1}} \]

In this expression, \(W_{\mathrm{H2O},j}\) comes from the vertical water profile, \(H_{\mathrm{eff},i}\) is species-specific, \(T_j\) is the layer temperature, and \(R\) is the gas constant in the unit system used by the code.

Second, the code computes a correction factor \(Q_j\). This correction accounts for the fact that rainout is not simply continuous removal at the raw frequency \(r_{j,i}\):

\[ Q_j = 1 - f_z + \frac{f_z}{\gamma_j r_{j,i}} \left(1 - e^{-r_{j,i}\gamma_j}\right) \]

Third, the corrected first-order rainout rate is calculated as:

\[ K_{\mathrm{RAIN},j,i} = \frac{1 - e^{-r_{j,i}\gamma_j}}{\gamma_j Q_j} \]

This K_RAIN value is what becomes wetdep(j,i) in the ODE system. Once it is stored there, PATMO uses it exactly like any other first-order sink: K_RAIN * [species].

In the modern sulfur cycle test, wet deposition is assigned directly in tests/modern_sulfur_cycle/test.f90. The model first clears all wet-deposition rates, then assigns rainout to selected sulfur species.

wetdep(:,:) = 0d0

! calculate wet deposition
call computewetdep(1,  2.0d-2)  ! OCS
call computewetdep(11, 5.0d-2)  ! CS2
call computewetdep(28, 1.0d-1)  ! H2S
call computewetdep(18, 4.0d3)   ! SO2
! call computewetdep(26, 5d14)  ! H2SO4
call computewetdep(19, 5d14)    ! SO4

! Turco H2SO4
wetdep(12,26) = 1.77E-06
wetdep(11,26) = 3.54E-06
wetdep(10,26) = 5.31E-06
wetdep(9,26)  = 7.08E-06
wetdep(8,26)  = 8.85E-06
wetdep(7,26)  = 1.06E-05
wetdep(6,26)  = 1.24E-05
wetdep(5,26)  = 1.42E-05
wetdep(4,26)  = 1.59E-05
wetdep(3,26)  = 1.77E-05
wetdep(2,26)  = 1.95E-05
wetdep(1,26)  = 2.12E-05