A corrected formulation of the Multilayer Model ( MLM ) for inferring gaseous dry deposition to vegetated surfaces

The Multilayer Model (MLM) has been used for many years to infer dry deposition fluxes from measured trace species concentrations and standard meteorological measurements for national networks in the U.S., including the U.S. Environmental Protection Agency’s Clean Air Status and Trends Network (CASTNet). MLM utilizes a resistance analogy to calculate deposition velocities appropriate for whole vegetative canopies, while employing a multilayer integration to account for vertically varying meteorology, canopy morphology and radiative transfer within the canopy. However, the MLM formulation, as it was originally presented and as it has been subsequently employed, contains a non-physical representation related to the leaf-level quasi-laminar boundary layer resistance that affects the calculation of the total canopy resistance. In this note, the non-physical representation of the canopy resistance as originally formulated in MLM is discussed and a revised, physically consistent, formulation is suggested as a replacement. The revised canopy resistance formulation reduces estimates of HNO3 deposition velocities by as much as 38% during mid-day as compared to values generated by the original formulation. Inferred deposition velocities for SO2 and O3 are not significantly altered by the change in formulation (<3%). Inferred deposition loadings of oxidized and total nitrogen from CASTNet data may be reduced by 10 e20% and 5e10%, respectively, for the Eastern U. S. when employing the revised formulation of MLM as compared to the original formulation. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons. org/licenses/by/3.0/).


Introduction
The Multilayer Model (MLM; Meyers et al., 1998;Cooter and Schwede, 2000) has been used for many years to infer dry deposition fluxes from measured trace species concentrations and selected meteorological measurements for U.S. national networks (Clarke et al., 1997;Finkelstein et al., 2000;Schwede et al., 2011;USEPA, 2013). MLM has been evaluated (Meyers et al., 1998;Finkelstein et al., 2000;Sickles and Shadwick, 2002) as well as modified (Wu et al., 2003a,b) but the version as described by Cooter and Schwede (2000) is currently employed to infer dry deposition fluxes from measurements made in the CASTNET network (Schwede et al., 2001;Baumgardner et al., 2002;Sickles and Shadwick, 2007;Bowker et al., 2011;USEPA, 2013). Conceptual components of the original model have also been used to simulate the bi-directional exchange of ammonia between vegetation and the atmosphere (Wu et al., 2009). MLM utilizes the ubiquitous Ohm's Law (i.e., resistance) analogy to calculate deposition velocities appropriate for whole vegetative canopies, while employing a multilayer integration within the canopy to account for vertically varying meteorology, canopy morphology and radiative transfer. However, close inspection of the MLM formulation, as it is presented by Meyers et al. (1998) and as it is typically employed (Finkelstein et al., 2000;Finkelstein, 2001;Baumgardner et al., 2002;Sickles and Shadwick, 2007;Schwede et al., 2011), reveals a representation of the leaf-level quasilaminar boundary layer resistance, r b , that is non-physical. In this technical note, we explore the impact of this non-physical formulation and suggest a revised formulation that is more consistent with the accepted conceptual model of multilayer resistance-based dry deposition.

Multilayer model description and revised formulation
Within the MLM, hourly estimates of deposition velocities (V d ) for each deposited trace species are calculated from (1) and SðzÞ ¼ S orig ðzÞ ¼ LAIðzÞ r c ðzÞ (2) with 1 r c ðzÞ ¼ 1 r b ðzÞ þ r s ðzÞ þ r meso þ 2 r b ðzÞ þ r cut : Parameterizations for the various resistances in Eqs.
In this formulation, the aerodynamic resistance, R a , is in series with the whole canopy and ground surface resistances, which are represented by the parenthetical term in Eq. (1). At each level, the total canopy resistance, r c (z), is calculated as the sum of parallel resistances due to diffusion into the leaf stomata and absorption into the plant mesophyll and to adsorption onto the leaf cuticle (the factor of 2 appearing because adsorption may occur on either side of the leaf, while stomata are typically found only on the underside of leaves for hypostomatous plants). As represented in Eq. (3), the leaf boundary layer resistance r b , exists in series with each of these possible pathways for trace species deposition to a leaf. However, this representation of r b as occurring in series individually with each leaf deposition pathway is inherently nonphysical.
As a trace species molecule approaches a leaf it must first diffuse through any existing laminar boundary layer next to the leaf surface. At that juncture, the molecule can then either be deposited onto the leaf's cuticle (on the top or bottom of the leaf) or it can diffuse into a stoma and be absorbed into the mesophyll. Given this conceptual picture, it is clear that the boundary layer resistance is in series independent of the two subsequent parallel pathways for trace species deposition (i.e., onto the leaf cuticle or through the stoma into the leaf mesophyll), so that the canopy resistance should be formulated as or after rearrangement, r l ðzÞ ¼ r cut ðr s ðzÞ þ r meso Þ r cut þ 2ðr s ðzÞ þ r meso Þ : Then, r c ðzÞ ¼ r b ðzÞ þ r cut ðr s ðzÞ þ r meso Þ r cut þ 2ðr s ðzÞ þ r meso Þ ; and the revised local canopy sink becomes S rev ðzÞ ¼ LAIðzÞ r b ðzÞ þ r cut ðr s ðzÞ þ r meso Þ r cut þ 2ðr s ðzÞ þ r meso Þ À1 :

Analysis of the original and revised formulations
Consideration of limiting cases of the two formulations for r c (z) (Eqs. (3) and (7)) is instructive. First, in the case where the cuticular and mesophyllic resistances approach zero and the boundary layer resistance is much larger than the stomatal resistancer cut /0 r meso /0 r b ðzÞ >> r s ðzÞ (9) then the leaf-level resistance at height z, from Eq. (3) is given by Nomenclature h c height of the canopy (m) LAI(z) leaf area index (single-sided) of the canopy at height z (m 2 m À2 ) R a aerodynamic resistance (s m À1 ) r a,soil subcanopy aerodynamic resistance for deposition to soil (s m À1 ) r b (z) leaf boundary layer resistance at height z (s m À1 ) r c (z) total resistance of canopy elements at height z (s m À1 ) r cut leaf cuticular resistance (s m À1 ) r l (z) total leaf surface resistance at height z (s m À1 ) r meso leaf mesophyll resistance (s m À1 ) r s (z) leaf stomatal resistance at height z (s m À1 ) r soil soil uptake resistance (s m À1 ) S(z) general local sink term for trace species deposition to the canopy (s À1 ) S orig (z) local sink term of the original MLM formulation for trace species deposition to the canopy (s À1 ) S rev (z) local sink term of the revised MLM formulation for trace species deposition to the canopy (s À1 ) V d whole canopy deposition velocity (m s À1 ) z vertical coordinate from ground surface upward (m) On the other hand, the leaf-level resistance from Eq. (7) reduces to 1 r c ðzÞ ¼ 1 Clearly, the original MLM limiting case result given by Eq. (10) is inconsistent with the conceptual idea that as the surface resistances become vanishingly small the overall local canopy resistance should approach the boundary layer resistance. However, the correct limiting case result (Eq. (11)) is produced from the revised formulation of Eq. (7).
Alternatively, in the limiting case where the boundary layer resistance approaches zero - then r c (z) from Eq. (3) becomes 1 r c ðzÞ ¼ 1 r b ðzÞ þ r s ðzÞ þ r meso þ 2 r b ðzÞ þ r cut / 1 r s ðzÞ þ r meso þ 2 r cut (13) and from Eq. (7) we obtain the identical limiting result 1 r c ðzÞ ¼ r b ðzÞ þ r cut ðr s ðzÞ þ r meso Þ r cut þ 2ðr s ðzÞ þ r meso Þ À1 / r cut þ 2ðr s ðzÞ þ r meso Þ r cut ðr s ðzÞ þ r meso Þ ¼ 1 r s ðzÞ þ r meso þ 2 r cut : 4. Impact of the revised formulation on MLM whole canopy deposition velocities The revised formulation of r c (z) as presented in Eq. (7) was implemented into the version of MLM as described by Cooter and Schwede (2000). All other elements of MLM, including all resistance parameterizations, were left the same as specified in Cooter and Schwede (2000). A test meteorological data set was used to generate V d values for six monitoring sites (Table 1) representing a sampling of different vegetation species and morphology. Identical meteorological data were used for all sites to allow a direct assessment of how differing site characteristics affect V d values calculated by the new formulation. Fig. 1 presents a comparison of deposition velocities generated with the original and revised formulations for SO 2 , O 3 and HNO 3 for Oak Ridge, TN. For these three gaseous species, only the V d values for HNO 3 are significantly affected by implementation of the revised r c formulation. For SO 2 and O 3 , the boundary layer resistance is smaller than the leaf surface resistances (r s , r cut and r meso ) and thus the calculated V d values for these species are altered only slightly during the day (À3%). On the other hand, for HNO 3 with an effective Henry's law coefficient of >10 12 M atm À1 , the cuticular and mesophyllic resistances are essentially zero so that the deposition velocities for HNO 3 are reduced by as much as 38% during mid-day and 32% on average between 6:00 am and 6:00 pm LST. Similar results are obtained for the other five sites as shown in Table 1 over a variety of vegetation types and leaf area index (LAI) values. HNO 3 V d reductions across the sites range from 10% at the low LAI site Pawnee Grasslands to 32% at Oak Ridge and 27e28% at the other locations. O 3 and SO 2 values are reduced by 3% or less at all sites.
Interestingly, evaluations of MLM against measurements have not noted biases of the magnitude as are suggested here (Meyers et al., 1989;Cooter and Schwede, 2000;Baumgardner et al., 2002;Schwede et al., 2011), although Meyers et al. (1998) did report that mean biases for daytime estimates of V d for HNO 3 were positive and ranged from 0.09 cm s À1 for corn canopies to 0.47 cm s À1 for pastures. However, the reported precisions of these biases for HNO 3 were much larger, ranging from AE0.88 cm s À1 for corn to AE1.5 cm s À1 for pasture. The difficulties of obtaining accurate measurements of HNO 3 are well known (Huey et al., 1998;Kita et al., 2006), so that measurement uncertainties likely obscured the biases reported here.

Summary, implications and recommendations
In this work, the non-physical representation of the canopy resistance as originally formulated in MLM has been noted and a revised, physically consistent, formulation has been suggested as a replacement. The revised canopy resistance formulation reduces estimates of HNO 3 deposition velocities by as much as 38% during mid-day as compared to values generated by the original formulation. Deposition velocities for SO 2 and O 3 are only slightly altered by the change in formulation (<3%). The MLM is used by the U.S. Environmental Protection Agency's (U. S. EPA's) Clean Air Status and Trends Network (CASTNet) to derive deposition velocities for CASTNet-based deposition estimates (Cooter and Schwede, 2000;U.S. EPA, 2013). Implementation of the revised, physically consistent, formulation into MLM will result in reduced estimates of dry deposition of gaseous HNO 3 . Dry deposition of HNO 3 represents 35e40% of total oxidized nitrogen deposition and roughly onequarter of total nitrogen (i.e., oxidized þ reduced forms) deposition in the Eastern U.S. (Sickles and Shadwick, 2007). If the results of Table 1 and Fig. 1 are typical (10e30% reduction in HNO 3 deposition), estimated oxidized and total nitrogen deposition loadings would be reduced by as much as 10e20% and 5e10%, respectively, for the Eastern U. S. A re-evaluation of estimated deposition for total nitrogen over the U. S. should be undertaken (U.S. EPA,2013;Baumgardner et al., 2002;Sickles and Shadwick, 2007) to correct the historical record or at minimum to provide an indication of the sensitivity of large-scale loadings to the way in which field measurements are interpreted. Furthermore, observations of surfacee atmosphere exchange of reactive nitrogen species, particularly nitric acid, are needed to better constrain the parameterizations used to estimate total nitrogen deposition. Ideally, such efforts would leverage recent advances in in situ techniques (Neuman et al., 2000;Farmer et al., 2006;Veres et al., 2008) to directly measure nitric acid fluxes via eddy covariance over a range of vegetated environments.