Surface Flux Modeling for Air Quality Applications

For many gasses and aerosols, dry deposition is an important sink of atmospheric mass. Dry deposition fluxes are also important sources of pollutants to terrestrial and aquatic ecosystems. The surface fluxes of some gases, such as ammonia, mercury, and certain volatile organic compounds, can be upward into the air as well as downward to the surface and therefore should be modeled as bi-directional fluxes. Model parameterizations of dry deposition in air quality models have been represented by simple electrical resistance analogs for almost 30 years. Uncertainties in surface flux modeling in global to mesoscale models are being slowly reduced as more field measurements provide constraints on parameterizations. However, at the same time, more chemical species are being added to surface flux models as air quality models are expanded to include more complex chemistry and are being applied to a wider array of environmental issues. Since surface flux measurements of many of these chemicals are still lacking, resistances are usually parameterized using simple scaling by water or lipid solubility and reactivity. Advances in recent years have included bi-directional flux algorithms that require a shift from pre-computation of deposition velocities to fully integrated surface flux calculations within air quality models. Improved modeling of the stomatal component of chemical surface fluxes has resulted from improved evapotranspiration modeling in land surface models and closer integration between meteorology and air quality models. Satellite-derived land use characterization and vegetation products and indices are improving model representation of spatial and temporal variations in surface flux 272 processes. This review describes the current state of chemical dry deposition modeling, recent progress in bi-directional flux modeling, synergistic model development research with field measurements, and coupling with meteorological land surface models.


Introduction
Atmospheric modeling systems that simulate meteorology, chemistry, or climate include model algorithms for air-surface exchange processes.Surface fluxes of heat, moisture, and momentum are critical elements of meteorology and climate models while chemistry/air quality models also require algorithms for surface fluxes of trace chemical species.Vertical flux of any quantity in the atmospheric surface layer is well described by gradient turbulent diffusion.Momentum flux can be modeled simply by means of turbulent surface layer similarity theory since momentum goes to zero at some finite height, known as the roughness length (z o ), which depends on the characteristics of the surface and vegetation.Sensible heat flux is slightly more complex than momentum flux since thermal exchange processes must interact directly with the surface and turbulent flow cannot extent all the way to the surface because of the "no-slip" condition.Thus, for heat there is an additional consideration of laminar molecular diffusion over a very thin layer in contact with all ground and canopy surfaces.Modeling moisture fluxes adds another element of complexity because of the multiple sources of surface moisture, via evapotranspiration (ET), which is regulated by leaf stomata, soil moisture, which requires diffusion through soil, and surface wetness from rain and dew.The surface fluxes of many chemical species that are involved in air quality modeling have similarities to moisture fluxes, since they may also involve multiple pathways for exchange with the surface, including stomatal, ground, water, built surfaces, and leaf cuticles.For example, in two different field experiments where both moisture fluxes and ozone fluxes were measured by eddy covariance, one for a vineyard and cotton field in California [1] and another for a soybean field in Kentucky [2], the stomatal pathway for ozone and water vapor fluxes were shown to be highly correlated.Thus, modeling chemical dry deposition and bi-directional surface fluxes involve many of the same physical processes and can use many common algorithms with moisture flux calculations in the land surface models of meteorology and climate models, especially for the stomatal component.
While chemical constituents of the atmosphere can undergo a variety of chemical and physical transformations, dry and wet deposition are the ultimate removal mechanisms of mass emitted to the atmosphere.Surface fluxes of some chemical species, notably ammonia [3][4][5], mercury [6,7], and certain volatile organic compounds (VOCs) [8] are often bi-directional, such that the surface may be alternately a source or sink of atmospheric mass.For other chemicals, particularly biogenic VOCs such as isoprene, terpenes, and alcohols the surface emissions from biological processes are very important atmospheric sources [9].
In addition to their roles as sources and sinks of atmospheric constituents, surface fluxes also have important effects on vegetation, ecosystems, and human health.Risk assessments and setting standards for protecting vegetation and ecosystems from damage by air pollutants has usually involved concentration based indexes such as AOT40 (Accumulated Over a Threshold of 40 ppb) which is used for assessing ozone impacts.A new flux based index called AFstY (Accumulated stomatal Flux above a flux threshold Y) has been proposed to estimate O 3 dose to leaf interiors where damage occurs rather than in the air outside the leaf [10].A recent study in Europe used the Deposition of Ozone and Stomatal Exchange (DO3SE) model to compute the AFstY index to assess the O 3 impacts on key European tree species [11,12].Their results showed that the AFstY index differed significantly from the AOT40 over different regions and species.
Until recently, secondary standards for protecting ecosystems from damage caused by NOx and SOx in the U.S. have been based on ambient concentrations even though flux based standards were proposed more than 10 years ago [13].Recently, however, new SOx-NOx secondary standards are being considered to protect both terrestrial and aquatic ecosystems from acidification and nutrient enrichment.Although the U.S. Clean Air Act requires that regulations must be based on ambient concentration levels, the newly proposed standards include linkages to deposition of SOx and NOx [14] to be computed by the Community Multiscale Air Quality CMAQ model [15].Therefore, chemical surface flux modeling has become a more important process in air quality models used for regulation assessment and enforcement because environmental standards are now being promulgated that include assessments of deposition effects on ecosystems.
Since industrialization, atmospheric reactive nitrogen concentrations and deposition have increased tremendously [16][17][18].Deposition of reactive nitrogen including, ammonia, NO 2 , nitric acid, ammonium and nitrate containing aerosols, and organic nitrates has serious deleterious impacts on terrestrial and aquatic ecosystems [19][20][21].Thus, modeling anthropogenic and natural emissions, chemical transformations, wet deposition, and surface fluxes is important not only for air quality studies but also for ecosystem impact assessments.Ammonia from fertilizer and livestock operations contribute the largest portion of excess nitrogen emissions to the environment [22].Accurate accounting of fertilizer sources of ammonia requires development of new modeling techniques such as bi-directional surface fluxes, as well as integration with agricultural management models (see Section 3.3).
This article presents an overview of the current state of science and modeling for surface flux processes in air quality modeling systems.We start with a presentation of the governing principles and equations of air-surface exchange that are applicable to moisture, heat, and momentum in meteorology models and chemical surface fluxes in chemical transport models in Section 2.Then, in Section 3 we focus on modeling techniques for chemical surface fluxes, including gas and aerosol dry deposition and bi-directional fluxes.Section 4 gives an overview of recent advances in dry deposition and bi-directional modeling along with a few examples of model results.Section 5 contains a concluding discussion of limitations of current models and future research and model development directions.As developers of land surface models (LSM) [23,24], dry deposition [25], and bidirectional models [26] we often use the formulation of these models as they are applied in the Weather Research and Forecast (WRF) model [27] and the CMAQ model as examples of parameterization techniques and to illustrate some of the characteristics of model results.Note that much of the description of the dry deposition model used in CMAQ has not been previously published.

Governing Principles of Surface Fluxes Modeling
Turbulent fluxes of any atmospheric quantity, including heat, moisture, momentum and trace chemical species, can be computed according to Monin-Obukov similarity theory in the atmospheric surface-layer which relates fluxes at height z near the surface to the mean vertical profile: (1) where F α is the vertical flux of any quantity α, which can be measured as the mean value of the turbulent covariance of w and  wherever stationarity and horizontal homogeneity can be reasonably assumed.Also, u * is the surface friction velocity, k is the von Karman constant,  m and  h are non-dimensional profile functions derived empirically from observed data, and L is the Monin-Obukov length scale which is defined as: where <w'T o '> represents the kinematic heat flux and Τ o represents the average temperature in the surface layer.Assuming constant flux between some reference height z r within the surface layer and the surface, Equation 1 can be vertically integrated to give: where  hn is the non-dimensional temperature profile constant for neutral conditions,  h is the stability correction function for heat, and α ο is the value of α at a micro-distance above the surface where turbulence cannot penetrate.Equation (3) can be rewritten as: where R α is generally known as the aerodynamic resistance.However, Equation (4) by itself is not very useful because of the difficulty in measuring or modeling quantities at a tiny distance above the surface.We would rather replace α ο with α s (the quantity at the surface) which can be measured, estimated from surface properties, or assumed to equal zero, as for dry deposition.Again assuming that the flux is constant down to the surface, we can define two more resistances for the quasi-laminar boundary layer (R b ), which is the non-turbulent micro layer adjacent to all surfaces, and a surface resistance (R s ): where α δ represents the quantity α in the air at the surface.Combining Equations ( 4) and ( 5) yields: (6) For dry deposition α s is assumed to be zero and Equation 6 is expressed as: , , 1 For bi-directional surface flux Eq 6 can be more complicated because α s may be different for various components of the surface such as the ground and vegetation canopy as described in section 3.3.R a represents the resistance to flux by turbulent diffusion through the turbulent surface layer while R b represents molecular or Brownian diffusion across a very thin quasi-laminar boundary layer adjacent to the surface.A general form of R b for any scalar quantity was recommended by Wesely and Hicks [28] as: (8) where Sc is the Schmidt number (Sc = ν/D), with  representing the kinematic molecular viscosity, and D is the molecular diffusivity of the scalar quantity.Pr is the Prandtl number, which is the analogous quantity for heat (Pr = νν θ ), where ν θ is the molecular thermal diffusivity.The parameter B −1 is the inverse Stanton number, a dimensionless heat transfer coefficient.The value of B −1 may depend on the Reynolds number as well as the type of surface roughness.For example, Garratt and Hicks [29] suggested a value of kB −1 = 2 for fully turbulent flow over fibrous vegetative canopies.
Note that, since R a represents the aerodynamic properties of the turbulent surface layer which depends purely on meteorological parameters and surface characteristics (i.e., surface roughness) it is independent of the chemical or physical properties of the transported scalar.Thus, R a , once computed, can be used in the surface flux calculations of any scalar.R b , however, combines meteorological and fluid flow parameters (i.e., and B) with properties of the transported scalar (i.e., molecular diffusivity) and must be computed individually for each chemical species.

Chemical Surface Fluxes
Chemical surface fluxes are partially analogous to moisture fluxes, especially for the stomatal conductance, which can be adapted to any gas by scaling by the ratio of molecular diffusivity of the chemical species to the molecular diffusivity of water vapor.

Dry Deposition Modeling
For chemical species which are dominated by one-way flux to the surface, the surface concentration is assumed to be zero.Thus, for such chemical scalars the gradient flux across the air-surface interface is always directed toward the surface no matter how low the air concentration.This special case is simplified as shown by Equation 7, where the flux is a simple linear combination of air concentration near the surface and the deposition velocity V d .The surface resistance R s represents all pathways for sinks of atmospheric chemical species to the surface.The major pathways for trace gasses include stomatal uptake, deposition to leaf cuticles, and deposition to ground surfaces.The relative importance of these pathways for any particular chemical species depends on their chemical and physical properties, such as solubility, reactivity, and molecular diffusivity.For example, nitric acid is both highly soluble and highly reactive such that its surface resistance is effectively zero (compared to R a and R b ) [30,31].Therefore, the deposition velocity of nitric acid can be estimated reasonably well as: 1/(R a + R b ).For other less reactive gasses the surface resistance is more important and each of the main pathways may contribute significantly to the surface flux.The simplest representation of R s for a single layer or "Big Leaf" model is represented as a combination of parallel resistances given by: where R st is stomatal resistance, R w is cuticular resistance, and R g is ground resistance.
However, since leaves are displaced above the ground there should be additional resistance for the ground pathway to account for transport through the plant canopy.Furthermore, some models account for concentration gradients within the canopy by using a multilayer approach to account for the vertical distribution of leaf surface area [32,33].Multilayer models can also calculate the light attenuation within the canopy and thereby more accurately represent the fraction of sunlit and shaded leaves at each level which has significant influence on stomatal resistance (see below the section on Scaling from Leaf to Canopy).For large scale air quality modeling, however, the complication of multi-layer models is usually not desired since the variety of vegetation contained within the large areas represented by each model grid cell is not amenable to realistic representation of canopy structure.Alternatively, a simple in-canopy resistance can be added in series with the ground surface resistance.For example, Erisman et al. [34] suggested: where LAI is the one-sided leaf area index, h c the canopy height and b an empirical constant taken as 14 m −1 .Zhang et al. [35] developed a similar parameterization for R ac related to LAI and u * that includes variations by land use category and day of the year to account for changes in canopy structure.
Another complication is that there may be an additional resistance to absorption into the mesophyll tissue inside the stomatal cavity of the leaf for some chemicals.Thus, including a mesophyll resistance R m and R ac , the total R s becomes: Thus, dry deposition velocity V d can be represented by resistances in series and parallel as shown in Figure 1.
(10) Stomatal resistance or conductance (g st = 1/R st ) controls the flux of gasses, including CO 2 and water vapor, into and out of leaves.Since stomatal resistance is the primary regulatory process for evapotranspiration, it is a key parameter for LSMs that are included in meteorology models.Thus it makes sense to use the same parameterization of the stomatal pathway for the meteorological and chemical components of air quality modeling systems [25].
The primary function of leaf stomata is to allow the intake of CO 2 for photosynthesis while minimizing the cost of water vapor loss to the atmosphere.In addition, the diffusion of water out of stomata into the atmosphere (leaf transpiration) is the major mechanism for mass flow of liquid water with any dissolved mineral nutrients from the roots to the leaves through the xylem.Stomatal opening responds to environmental conditions, including sunlight, humidity, temperature, soil moisture, CO 2 concentration, and nutrient availability.Stomatal conductance, which is directly related to the degree of stomatal opening, is usually modeled in LSM and dry deposition models following either the empirical approach of Jarvis [36] or photosynthesis-based semi-empirical models [37][38][39][40].
The Jarvis model uses a series of empirical functions as linear coefficients to the maximum stomatal conductance g smax (conductance at maximum opening), which is a characteristic of each plant reflecting the size and density of stomata.For example, Pleim and Xiu [41] used this approach in their LSM that is in the WRF model and for dry deposition in the CMAQ model to model bulk (canopy level) stomatal conductance g sb as: where the functions F 1-4 represent the fractional degree of stomatal closure caused by the environmental stresses: photosynthetically active radiation PAR, root-depth soil moisture w 2 , relative humidity at the leaf surface RH s , and air temperature in the canopy T a .Note that the scaling up from leaf to canopy appears from Equation ( 12) to be simply a linear dependence on LAI.There is, however, an additional dependence on LAI in the F 1 term that accounts for increased shading in denser canopies [42].Discussion of simple and more detailed techniques for scaling from leaf to canopy is presented below.The Jarvis-type empirical stomatal conductance model has been widely used in LSM components of meteorology [23,43,44] and dry deposition models [32,45,46].
The photosynthesis approach solves for stomatal conductance assuming that CO 2 diffusion through the stomata is equal to the net CO 2 assimilation rate for photosynthesis minus respiration.These models originated in plant biochemical studies [38,39] and are widely used in ecology and hydrology models and were later adapted for use in atmospheric models [40,47].Based on the optimality hypothesis of the study by Cowan [37] plants optimize their stomatal conductance for maximizing photosynthesis at a given amount of transpiration.Photosynthesis-based models are used in the third and fourth generations of LSMs in Global Climate Models (GCM) for modeling ET and the carbon cycle [48][49][50][51].The semi-empirical Ball-Berry model [39] computes stomatal conductance based on: where F h is an empirical function of humidity, A is the CO 2 assimilation rate, c s is the CO 2 concentration at the leaf surface, a 1 is an empirical coefficient, and g o is the minimum stomatal conductance.These models have the advantage of describing the rate of photosynthesis in a mechanistic fashion thereby providing a theoretical basis for stomatal response to environmental conditions such as temperature, solar radiation, water stress, and nutrient availability.However, there is no theoretical basis for response to air humidity or soil moisture.Since these models are combinations of theoretical biochemistry and empirical parameterizations, their predicted stomatal response to each independent environmental variable is not clearly evident.Thus the photosynthesis approach is semi-empirical, since the parameters a 1 and g o are empirically estimated for each plant type.Also, F h is an empirical function of air humidity that is similar to F 3 in Equation (12).While Jarvis-type models typically use simple functions of ambient vapor pressure deficit (D a = e sat − e a ) [43] the photosynthesis based models alternatively use F h = RH s [39] or F h = (1 + D s /D o ) −1 as suggested by Lohammer et al. [52], where D s is the vapor pressure deficit at the leaf surface (D s = e sat − e s ) and D o is an empirical constant.Note that the Lohammer humidity function results from a linear relationship between g s and ET, which has often been observed in laboratory experiments [40,53].While F h in Equation ( 13) and F 3 in Equation ( 12) seem to be equivalent, Aphalo and Jarvis [54] pointed out that since the assimilation rate has an implicit dependence on humidity, F h does not encapsulate the entire humidity dependence of these models.Pleim [55] analyzed the three humidity functions and found that for use in a Jarvis-type model F 3 = RH s is most appropriate and, therefore, is now used in the Pleim-Xiu LSM [24].
(12) (13) Recently, many LSM formulations in climate models [48,56,50] have used the photosynthesis approach to stomatal conductance.Such models have a distinct advantage over models using the Jarvis empirical approach because of the explicit dependence of stomatal conductance on CO 2 concentration.Thus, feedbacks of increasing CO 2 on ET can be included in future GCM projection simulations.While biochemical stomatal models have been developed and tested for LSM and dry deposition for use in meteorology and air quality models [33,[57][58][59] they have yet to be applied in operational mesoscale models.

Cuticle Resistance
Leaf cuticles are impregnated with waxy substances that resist water loss.Thus, evaporative fluxes from leaf surfaces (non-stomatal) occur only from external wetness from rain or dew.For some trace chemicals, however, deposition to dry leaf cuticles is a significant process for removal.Chemical deposition to wet leaves scales with Henry's law of solubility with a dependence on pH for species that readily dissociate in the aqueous phase (e.g., SO 2 and NH 3 ) [60,61].Mechanisms for deposition to dry cuticles are less well understood.Some dry deposition models simply specify constant values for cuticle resistance or even total canopy resistance for each chemical species based on empirical fits to observations without attempting to hypothesize about the mechanism [62][63][64].The problem with this approach is that it is limited to the few gasses for which field flux measurements have been made in a variety of landscapes.Other models generalize the calculation of cuticular resistance based on chemical properties such as a reactivity factor [35,60,61,65] so that estimates for R w can be made even for species for which surface fluxes have not been measured in the field.Some models also include a dependence of R w on Henry's law coefficient and/or RH even for dry cuticles particularly for highly soluble chemicals like NH 3 and SO 2 [34].Jones et al. [66] found that the R w for ammonia is a linear function of NH 3 concentration, suggesting an inhibiting effect caused by saturating the cuticular membrane.For hydrophobic organic chemical species it has been hypothesized that a small amount of mass may dissolve into the wax-like lipids in the leaf cuticle and then diffuse through the cuticle to the leaf interior.Thus, R w for chemicals that have low water solubility but high lipid solubility may be parameterized using the octanol-water partitioning coefficient [33,67].The R w values computed by such parameterizations are usually quite high resulting in very small fluxes.However, for toxic organic chemicals, often known as persistent organic pollutants (POPs), absorption into leaf cuticle of even small quantities is an important pathway for contamination of vegetation that can subsequently enter the food chain [68].
The dominant mechanism for cuticular deposition depends on the properties of chemical species.Soluble species such as SO 2 , H 2 O 2 , and NH 3 quickly dissolve into water on leaves that are wetted by rain or dew.Cuticle resistance is very low for these species when leaves are wet but resistance for dry portions of the leaves can be quite high leading to low deposition velocities in dry areas.For example, Figure 2a shows the dry deposition velocity for SO 2 from the WRF-CMAQ model system in mid-afternoon (20 UT) on July 25, 2006 on a 12 km model grid.Note that the distinct pattern of high deposition velocities closely follows the pattern of canopy water content shown in Figure 2b.Thus, it is important to account for parallel pathways to wet and dry portions of the leaves and couple the dry deposition calculation to a realistic treatment of canopy water.For the example shown in Figure 2 the dry deposition model in CMAQ is coupled to the PX LSM in WRF which integrates canopy water as a prognostic variable.In heavily vegetated areas, leaf surfaces are the most important non-stomatal pathway because of their large surface area and because concentrations are much greater in the upper canopy than near the ground.However, in more sparsely vegetated areas deposition to ground and building surfaces can be important.Many dry deposition models parameterize soil resistance similarly to other non-stomatal resistances such as cuticles.For example, Wesely [65] scales R soil like R w by Henry's law and reactivity but with a different base value.Pleim et al. [60] and Padro et al. [61] also scale ground resistance by reactivity with different magnitudes for cuticle and soil.For wet surfaces, the same Henry's law scaling can be used for any surface.For un-saturated soils, resistance can depend on soil moisture content and pH, especially for SO 2 and NH 3 , where SO 2 has greater resistance to acidic soil and NH 3 has greater resistance to alkaline soil.Erisman et al. [34] describes soil resistance values and dependencies for many chemical species.Functional dependences of SO 2 deposition on pH and relative humidity are described by Baldocchi [69].

Snow Surfaces
Observed deposition velocities to snow are generally lower than observed for soil and other surfaces.For example, Helmig et al. [70] found that deposition velocities of ozone to snow of less than 0.01 cm/s gave the best fit to observations in the Arctic.This is at the low end of the range reported in the literature for ozone which is generally less than 0.1 cm/s [71,72] with decreasing deposition as the snowpack ages.However, a few studies showed higher deposition velocities [73][74][75][76].For SO 2 , Erisman [34] suggests a surface resistance to snow of 500 s/m when the temperature is below −1 °C but decreasing to 70 m/s at +1 °C.For soluble species like SO 2 , resistance decreases markedly as the surface temperature rises above freezing and the liquid water portion of the snowpack increases.
The dry deposition model in the CMAQ model represents snow resistance R snow as parallel resistances to the ice and liquid fractions of the snow pack when there is snow and the ground temperature is above freezing [77]: where R ice is the resistance to snow or ice, R gw is the resistance to wet ground, R sndiff is resistance for diffusion through the snowpack and the liquid water fraction x m is: where x m is limited to a maximum of 0.5 and is zero when the ground temperature T g is below −1 °C.Since the value of R ice is relatively large, deposition to snow is most important for temperatures above −1 °C for soluble species.In CMAQ, R ice is parameterized for various chemical species similarly to ground and cuticle, by relative reactivity scaling.

Scaling from Leaf to Canopy
LSMs and chemical surface flux models require information on canopy structure for scaling from stomatal conductance at the leaf level to canopy conductance.While LAI accounts for the scaling of leaf area to ground area, the relationship between canopy conductance and stomatal conductance is complicated by leaf shading and different plant species in complex canopies.The simplest approach is the "big-leaf" model, such as described by Noilhan and Planton [43] and shown in Equation (12) where the effects of shading are parameterized in the radiation function.More sophisticated canopy models that are often used for field site studies and crop models [78] calculate radiation penetration through multilayer canopies with discrimination of sunlit and shaded leaves and consideration of leaf angle [79].However, these models require detailed descriptions of canopy structure and phenology that is not practical to obtain from the land use databases typically used in regional or global scale modeling systems.A study by De Pury and Farquhar [79] showed that the big-leaf model with sunlit and shaded canopy photosynthesis performs nearly as well as a multi-layer model and significantly better than the big-leaf model with a single canopy photosynthetic rate.A similar study by Zhang et al. [80] showed that the sunlit/shaded big-leaf approach also compares well to multi-layer models for representing the stomatal pathway in dry deposition models.Based on these findings a sunlit/shaded big-leaf model has been applied to dry deposition calculations in the AURAMS air quality model [46].Parameters describing the fractional coverage of vegetation and the leaf area and canopy structure are crucial for accurate land surface modeling.Typically, vegetation parameters are specified in land-use related look-up tables and plant phenological dynamics are modeled using simple timedependent functions such as a single sinusoidal curve used in the LSM by Walko et al. [81] or deep soil temperature dependent functions used by Dickinson et al. [82] and Xiu and Pleim [24].These generalized techniques, however, are not responsive to local conditions.Recently, modelers have used MODIS vegetation parameters to improve LSM performance.For example, Lawrence and Chase [83] used MODIS retrieved vegetation and land use products to modify land surface parameters and to improve plant phenology performance in the Community Land Model version 3.0 (CLM 3.0) resulting significant improvements for climate simulations.Moore et al. [84] assimilated MODIS dynamic LAI and vegetative fractional cover (FC) products into the Regional Atmospheric Modeling System (RAMS) for model simulations of East Africa.Their results show dramatic improvement in land surface temperature (LST) simulation both spatially and temporally for most of the year.Several other recent studies [84][85][86][87] have also used satellite vegetation products, such as LAI from satellite images, in land surface modeling for global circulation and mesoscale atmospheric modeling.
An alternate approach to incorporating derived vegetation products such as LAI, fraction of photosynthetic active radiation fPAR, and FC, directly into LSMs, is to use vegetation indices (VIs), such as normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), as indicators of canopy phenology and canopy structure.For example, by comparing MODIS VI data (NDVI and EVI) with four intensively measured test sites, Huete et al. [88] demonstrated that MODIS VIs represented seasonal phenology well over numerous biome types in North and South America.Zhang et al. [35] developed a method to estimate phenology events based on the curvature-change rate for MODIS time series data.In a study by Zhang et al. [89] they demonstrated vegetation phenological transition dates identified in the MODIS data are strongly correlated with MODIS land surface temperature (LST) and latitudes for all MODIS land cover types.Studies by Nagler et al. [90] and Glenn et al. [91] showed that VIs provide integrated information on leaf angles, fractional cover, canopy architecture, LAI, and chlorophyll content for vegetated areas.Therefore, assimilating MODIS VI data into a regional scale meteorology and air quality modeling system will help capture not only spatial patterns of LAI but also plant phenological and structural patterns over large regions.
Since vegetation indices are sensitive to canopy structural variations they are highly correlated with ET and canopy-level stomatal conductance.Sellers et al. [92] suggested that spectral vegetation indices obtained from satellite sensors may provide good estimates for canopy photosynthesis and conductance.Glenn et al. [93] provide a review of techniques for estimation of ET directly from VIs and numerous references to studies where various methods involving VI algorithms are applied to a wide variety of natural and cultivated ecosystems.Evaluations of these techniques often yield remarkably good agreement with surface flux measurements.However, since these are empirically derived algorithms that require calibration with ground based measurements they cannot be generalized but need multiple recalibrations for each biome.Another limitation is that VI methods for ET estimation are not responsive to short duration moisture stress (when plant physiology is not yet affected).Thus, combination of VI techniques with LSMs where soil moisture and meteorological parameters are explicitly modeled could be a way to improve canopy-level stomatal conductance for more accurate ET calculations in meteorology models and more accurate chemical fluxes in air quality models.

Aerosols
While turbulent flux through the atmospheric surface layer is similar for gases and aerosols, there are several key processes that are unique to aerosol dry deposition, including gravitational settling, Brownian diffusion, surface impaction, surface interception, and rebound.All of these processes are functions of particle diameter D p , such that the dominant process depends on position in the particle size spectrum.For instance, Brownian diffusion dominates in the nano-range (D p < 100 nm), where deposition efficiency increases with decreasing diameter, and gravitational settling and inertial impaction dominate in the super micron range, where their efficiencies increase with increasing diameter.Thus, there is a minimum of aerosol dry deposition velocity in the 0.1-1.0m range which is aptly named the accumulation range (e.g., see Figure 3).Recent reviews by Pryor et al. [94] and Petroff et al. [95] provide comprehensive literature surveys and overviews of modeling and measurements of aerosol dry deposition.Therefore, this section provides a description of the equations and techniques used for aerosol dry deposition calculations in the CMAQ model as an example of a state-of-the-science formulation in a widely used operational model.To cast aerosol deposition in a resistance and velocity form, the simultaneous turbulent and sedimentation fluxes can be combined as derived by Venkatram and Pleim [96]: where the gravitational settling velocity is: ρ is the density of the aerosol, μ is air dynamic viscosity, and the Cunningham slip correction factor is: where λ is the mean free path of air.Unlike gas deposition the quasi-laminar boundary layer resistance R b is usually the limiting resistance for aerosols because Brownian diffusion is much slower for particles than molecular diffusion is for gases.However, the effects of inertial impaction and interception by protruding micro-scale roughness elements can partially bridge the diffusion layer such that R b is inversely related to three collection efficiencies [97]: where F f is an empirical correction factor to account for increased deposition in convective conditions as suggested by Binkowski and Shankar [98]: where is the convective velocity scale.The first term in Equation ( 19) represents the collection efficiency by diffusion across the quasi-laminar boundary layer, which is similar to Equation ( 8) for gasses except that the Schmidt number for aerosols is defined as Sc = ν/D B , where D B is Brownian diffusivity which is a function of particle diameter.Slinn [97] suggested a formulation for the interception collection efficiency E in that depends on the characteristic size of microscale structures such as tiny hairs on leaves and the characteristic size of the leaves.However, since it is difficult to specify realistic estimates of these parameters over the area of typical grid cells used by air quality models (i.e.~4-20 km) on the basis of available land use data the interception efficiency is not used in the CMAQ model.
Pryor [94] summarized many formulations for the impaction efficiency E im and compares their magnitudes over a range of aerosol sizes.Slinn [99] suggested: where Stokes number is: 1 0.24 (20) Equation ( 21) was the form of E im originally used in CMAQ.However, Slinn [97] noted that Equation ( 21) was derived from empirical fits to smooth surfaces and therefore may not be appropriate for complex canopies where there is a wide range of wind speeds and length scales of collectors (leaves).Thus, for vegetated surfaces, Slinn [97] recommended a form of E im with weaker dependence on St: This form of E im was recently adopted for the CMAQ model, using the recommendation by Giorgi [100] that c = 400, so that aerosol deposition to heavily vegetated regions is better represented.Air quality models that consider aerosol size distributions either by discrete size bins (sectional models) or by log-normal modes (modal models) compute aerosol dry deposition velocities as functions of particle diameter.Thus, for sectional models Equations (16)(17)(18)(19)(20)(21)(22)(23) are evaluated using the mass mean diameter of each section.For modal models an integrated V d is computed for each mode by integrating these equations over each log-normal size distribution as described by Binkowski and Shankar, [98] and Feng [101].Figure 3 shows an example of modal integrated V d as a function of modal mass mean diameter D g and the results of the change in the impaction efficiency from Equation ( 21) that was used in CMAQv4.4 to the formulation shown by Equation (23) as used in CMAQv4.5 and later versions.This change in impaction efficiency has a profound effect on increasing the deposition velocity in the size range where the D g for the accumulation mode typically resides.

Bidirectional Surface Flux Modeling
While most gas-phase atmospheric chemical constituents predominantly deposit to the surface (e.g., O 3 , SO 2 , HNO 3 ) and some chemicals are predominantly emitted from the surface (e.g., isoprene and terpenes), a few chemicals such as ammonia, mercury, and certain oxidized VOCs exhibit distinctly bi-directional behavior where they alternately deposit and emit on diurnal, seasonal, or sometimes longer time-scales.Generally, air quality models include simultaneous emissions and dry deposition fluxes of these materials that are estimated using entirely unrelated techniques.This is clearly a less realistic treatment of the salient physical processes than a bi-directional gradient driven flux model that is responsive to changes at the model time step level.Thus, air quality modelers have recently been replacing separate representations of emissions and dry deposition with bi-directional flux calculations.
The development and application of bi-directional modeling of ammonia is particularly important since it allows for more accurate surface fluxes from fertilized fields which is a major contributor to excess nitrogen pollution that has adverse effects on human health, biodiversity, and climate change [18,22].Ammonia fluxes can vary widely in magnitude and direction on spatial scales down to agricultural fields and temporally over diurnal, multi-day, and seasonal time-scales.In areas of unfertilized natural vegetation, ammonia flux is predominately deposition although small evasive fluxes can be driven by the thermodynamics of the equilibrium between ambient NH 3 and NH 4 in the leaf apoplastic solution when ambient concentrations are low and ambient temperatures are sufficiently high.Conversely, for heavily fertilized crops such as corn, ammonia fluxes are mostly upward (23) (evasion), with the greatest fluxes immediately after fertilization and gradually decreasing over subsequent weeks and months.For minimally fertilized crops such as soybeans, fluxes can be truly bi-directional on a diurnal basis as shown in Figure 4 [26,102].Sutton et al. [3] provided an extensive review of measurements and modeling techniques for ammonia surface exchange processes.In this section we outline the model for bi-directional flux modeling in the CMAQ modeling system developed initially following Sutton et al. [103] and Nemitz et al. [104] with further refinement through comparison to field study measurements in corn and soybean [102] fields.Personne et al. [105] and Burkhardt et al. [106] developed and evaluated similar bi-directional ammonia flux models in comparison to measurements in fertilized grassland during the GRAMINAE intensive measurement campaign in spring 2000 near Braunschweig, Germany.
Although bi-directional flux modeling has many commonalities with dry deposition modeling there are also important differences that stem from the complications of compensation concentrations in surface reservoirs.Hence, the simplicity of the linear flux-concentration relationship for dry deposition, which allows V d to be computed for each chemical species upstream of the execution of the chemical transport model, is usually not feasible for bi-directional calculations because compensation points may differ for different pathways.For ammonia and mercury, compensation point concentrations are usually very different at the soil, stomata, and cuticle.
The total flux between the plant canopy and the overlying atmosphere is the sum of two bi-directional pathways, to the leaf stomata F st and the soil F g , and one uni-directional deposition pathway to the leaf cuticle F cut (Figure 5).Because each bidirectional flux pathway depends on concentration difference across resistances, an in-canopy concentration χ c is computed at the junction of the ground, stomatal, and cuticle pathways, which is often referred to as the net canopy compensation point [104].Thus the total flux is defined as: where χ a is the concentration in the air above the canopy and the component fluxes are: The resistances R a , R ac , R b , and R st are essentially the same as used for gas dry deposition with adjustments to R b and R st for molecular diffusivity.The resistances R bg , R soil , and R w , however, required specific derivation for ammonia based on comparisons to field data and previous research described in the literature.In heavily fertilized fields ammonia concentrations in the soil can be extremely large, thus the parameterizations for R bg and R soil are critical for simulation of realistic fluxes from the soil.The cuticle resistance is also crucial since a large fraction of the soil flux can be absorbed by the leaves thereby greatly reducing the flux at the top of the canopy.A study by Bash et al. [107] using a combination of in-canopy concentration measurements, above-canopy flux measurements, and an analytical model, estimated that a corn canopy was a sink for up to 70% of ammonia emitted from the soil approximately one month following fertilization.
Air concentrations of ammonia in soil pores and stomatal cavities are assumed to be in Henry's Law equilibrium with aqueous concentrations of ammonium and hydrogen ion in the soil moisture and apoplast leaf tissue, respectively.Hence, the soil and stomatal compensation points (χ g and χ st ) can be computed as [108]: where A (2.746 × 1015) and B (4507) are constants derived from the equlibria constants, T g,st is the soil and leaf temperature (K), and Γ g,st is the dimensionless NH 3 emission potential from the soil and leaf stomata where: , , 10 / , Γ , Thus, ammonia bi-directional flux models depend on methodologies for estimating stomatal and soil gamma values for all land use and vegetation types.Zhang et al. [109] proposed a model where gamma values are specified according to land use category and season based on an extensive review of measurement studies.Massad et al. [5] proposed a more dynamic methodology for estimating gamma values, also based on extensive review of measurement studies, where Γ st for unmanaged ecosystems is computed based on nitrogen input from deposition while for managed ecosystems Γ g and Γ st are functions of nitrogen fertilizer application with a decay function by time since fertilization.Cooter et al. [110] demonstrated that NH 3 flux from managed agricultural soils can be reasonably estimated by integrating the bi-directional resistance-based flux model in CMAQ with components of the Environmental Policy Integrated Climate (EPIC) model [111] to compute soil gamma values.Preliminary annual model simulations show improved estimates of NH x wet deposition over the eastern US in general (reduction in mean error), with regions of significant bias reduction as well as regions of moderate increase in bias [112].A more complete agricultural fertilizer modeling system that takes daily meteorological and nitrogen deposition data from WRF/CAMQ output and uses site, soil, crop, and management information to estimate fertilizer application method, timing, amount, and rate for specific pastures and crops on the model grid is being integrated with EPIC and CMAQ on a continental scale [113].Future scenarios, such as increased bio-fuel production or climate change that may impact crop yields and fertilizer use, can be simulated with this system to assess their impacts on air quality.
Mercury is another important pollutant that exhibits strong bi-directional surface flux behavior.Unlike ammonia, however, the main source of surface mercury that volatilizes to the air is deposition from the air.Thus, bi-directional flux modeling is necessary to account for the hop-scotching of depositing and re-emitting mercury.Bash et al. [6] provided an overview of the processes and issues of bi-directional mercury exchange between the air and soil, vegetation, and surface waters.Zhang et al. [7] provided a review of measurements and modeling techniques.Bash [107] described a bi-directional elemental mercury exchange model that has been implemented and tested in the CMAQ system.
Semivolatile organic compounds (SOCs), such as polycyclic aromatic hydrocarbons (PAHs), polychlorinated biphenyls (PCBs), and other hazardous chemicals, transfer from air to the soil by particle or gas dry deposition and by wet deposition where they can accumulate over long periods of time [114].If concentrations in soil become sufficiently high, these compounds can volatilize back to the air.In addition, residual organochlorine pesticides (OCPs), such DDT and Chlordane are still high enough in some agricultural soils to still be an atmospheric source [115].Air-surface exchange of SOCs is modeled by estimation of equilibrium partitioning with bi-directional fluxes controlled by deviations from the equilibrium state [116].

Recent Advances in Chemical Surface Flux Modeling
In the early 1980s, two comprehensive reviews of gas and particle dry deposition were published by Hosker and Lindberg [117] and Sehmel [118].Later, Wesely and Hicks [119], reviewed the state of knowledge on dry deposition.At that time the major processes controlling the deposition of the major inorganic species such as SO 2 , O 3 were well understood because of the relative abundance of flux measurement studies of these species.For example, many flux measurements, mostly of ozone but also some of SO 2 and NO 2 , have been made by gradient and eddy correlation methods going back several decades, starting in the 1970s [120,121], with many more in the 1980s and 1990s [122][123][124][125][126]. Note that reliable eddy correlation measurements over forest canopies were not feasible until the 1990s [127][128][129].Because of the extensive measurement data in the published record the modeling of these species has been somewhat constrained.Significant improvements in the modeling gas dry deposition in recent years have resulted from advances in land surface modeling and closer coupling between air quality and meteorology models.Since stomatal resistance, aerodynamic resistance, quasilaminar boundary layer resistance, surface and canopy wetness, and snow cover are important parameters for both latent heat flux and dry deposition of many gases, air quality modeling systems that are comprised of meteorology and chemical transport models, such as the WRF-CMAQ system, can use the same parameters for both dry deposition and land surface modeling [25].
Examination of dry deposition velocities modeled by CMAQ along with concurrent parameters from the LSM in WRF demonstrates the close relationships between these processes.Figure 6 shows ozone dry deposition on a 12 km grid simulated by CMAQ for a summer afternoon (19 UT), which is within the range usually measured in the field with peak values of about 1.1 cm/s.The spatial pattern illustrates that ozone deposits primarily through the stomatal pathway since the highest values are in the most vegetated areas.It is interesting that high deposition velocities are simulated for both heavily forested areas in the East and agricultural areas in the Midwest because forests generally have higher LAI than crops but crops have much lower minimum stomatal resistance.Thus, ozone dry deposition velocity closely follows the latent heat flux (Figure 7) except where there is evaporation from rain wetted canopies.Note, however, that there is some enhancement of V d for wet canopies (e.g., in southern Quebec) because the cuticular resistance for ozone is set to a lower value than would result from Henry's law scaling (also lower than dry cuticular resistance) on the basis of a study comparing measured ozone flux and latent heat flux [2].CMAQ's dry deposition model has been evaluated through comparison to measurements over soybeans in Kentucky and mixed forest in upstate NY [130] as shown in Figure 8.While there is some lag in the morning increase of V d for the soybeans, and overprediction at night, the gradual afternoon decline which is caused by stomatal response to lower RH is well represented.The dry deposition velocity for SO 2 is far more influenced by surface wetness as shown in Figure 2a because of its high effective Henry's law solubility.The highest values of V d for SO 2 are all in the same areas where LE is high from evaporating surface wetness (see also Figure 2b).In these areas, where rain is occurring currently or recently, peak values are 3-4 cm/s, which indicates that there is very little surface resistance in these areas.Dry deposition of NO 2 and NO were studied in relationship to O 3 and ET in a series of gas exchange experiments involving tobacco and sunflower plants described by Neubert et al. [131].These experiments showed that NO 2 deposition is almost entirely stomatal with negligible mesophyll resistance, while NO had considerable mesophyll resistance that changed with light intensity.These experiments also showed that ozone deposition occurs by stomatal uptake with negligible mesophyll resistance and destruction at cuticle surfaces.A field measurement study of PAN fluxes over grassland by Doskey et al. [67] found that PAN deposition is mostly via the stomatal pathway but with a non-zero mesophyll resistance resulting in V d values 0.5-0.3 of V d for O 3 , which is in concurrence with previous studies.However, leaf level measurements reported by Sparks et al. [132] suggest unimpeded stomatal uptake of PAN, indicating insignificant mesophyll resistance and higher deposition velocities.Furthermore, recent eddy covariance flux measurements over loblolly pine by Turnipseed et al. [133] show total surface conductance about twice the stomatal conductance, which suggests significant cuticular uptake.On the basis of these measurement studies the CMAQ was revised in 2007 with higher values of surface reactivity scaling parameter for PAN resulting in maximum deposition velocities close to the values for O 3 .Recently, many more chemical species are being added to dry deposition models as measurement techniques improve in detection capability and response time and the needs of the modeling community extend to more environmental issues such as air toxics and ecosystem impacts.For example, many toxic gases have been added to the CMAQ model in recent years, thereby requiring estimates of dry deposition velocity for these chemicals.However, the lipid solubility scaling discussed in Section 3 has yet to be applied to the CMAQ dry deposition model.A recent study by Karl et al. [134] promotes the importance of the dry deposition of oxygenated VOCs, many of which are oxidation products of biogenic VOCs.They found that Methyl vinyl ketone (MVK) and Methacrolin (MAC), which account for about 80% of the initial oxidation products of isoprene, have deposition velocities of similar magnitude as O 3 .On the basis of this study CMAQ is being updated to include dry deposition of MVK and MAC.The dry deposition model in AURAMS described by Zhang et al. [46] already includes dry deposition of these species although with small values of V d than suggested by Karl et al. [134].

Conclusions
Modeling of surface fluxes for air quality modeling involves micrometeorology, plant physiology and biology, surface chemistry, soil hydraulics and physics, agricultural management, spatial descriptions of geophysical characteristics, and most importantly, surface flux measurements for a variety of chemical substances over a variety of landscapes.Uncertainties in parameterizations, incomplete knowledge of processes, limitations to theory, and paucity of detailed data and measurements, combine to limit the realism and accuracy of surface flux models used in large-scale (mesoscale to global) air quality modeling systems.Progress is being made on many of these fronts but significant uncertainties remain.Even for ozone dry deposition, which is the most widely measured substance, current model results still show substantial errors compared to observations (Figure 8).Some of the errors in chemical fluxes have their source in common with the meteorological fluxes (heat, moisture, and momentum), such as the aerodynamic resistance, while other errors are due to incomplete understanding of chemical interactions with surfaces and leaf tissue.Further errors can be attributed to insufficient detail and accuracy in the spatial and temporal representations of landscape, land use, and vegetation in meso to global scale grid models.
Surface flux modeling relies on flux-profile relationships in accordance with Monin-Obukov similarity theory (MOST) in the atmospheric surface layer [135] (Equation 1).The MOST relationships are defined by dimensionless universal profile functions that have been empirically determined from the 1968 KANSAS experiment [136].Thus, the validity of MOST is limited to flat terrain with uniform landscape and landcover [137].Not only are MOST-based model calculations less accurate in non-ideal conditions, but also flux measurements are generally unreliable in complex landscapes.For substances that are mainly limited by surface resistances (e.g., O 3 , NO 2 , PAN), the limitation of MOST and resulting inaccuracies in estimates of R a in non-ideal conditions such as sloping terrain or patchy land-use will have relatively small effects on V d calculations [119].However, for species not limited by surface uptake, such as HNO 3 or other highly soluble species such as H 2 O 2 or SO 2 when surfaces are wet, V d estimates should be considered especially suspect in hilly or patchy landscapes.
There are distinct advantages to an integrated or coupled approach to LSM and chemical flux modeling.As meteorology and air quality models become more integrated such as, WRF-Chem [138], WRF-CMAQ [139], COSMO-ART [140], GATOR [141], to name a few, there is more opportunity for integrated surface flux algorithms that share key parameters like R a , R st , and R b between meteorological and chemical flux models.Such integrated approaches not only ensure consistency but also potentially benefit the chemical flux calculations by using a more realistic stomatal conductance estimates from the ET modeling in the meteorological LSM.For example, in the WRF-CMAQ system the latent heat flux is optimized by an indirect soil moisture data assimilation scheme that effectively tunes the stomatal conductance through its dependence on soil moisture [41].Thus, the resulting stomatal conductance is more constrained and generally more accurate than estimates by stand-alone dry deposition models.In this way, chemical flux models and air quality models automatically benefit as land surface models are improved.
Advancements in land surface and chemical surface flux modeling will depend on improved descriptions of land use and vegetation characteristics.For example, high resolution land use and land cover data from the NLCD, based on 30 m resolution LANDSAT thematic mapper, have recently been incorporated into the WRF-CMAQ system for land surface and chemical flux modeling [142].Furthermore, a more detailed vegetation dataset is being developed by combining the high spatial resolution land use data from NLCD with county-level tree species data from the Forest Inventory Analysis (FIA) and crop data from the National Agricultural Statistical Service (NASS) [113].By combining with ecoregion information this dataset will be able to better define key parameters such as minimum stomatal conductance, LAI, and canopy height across continental gridded modeling domains.Ultimately, the same data set with high spatial and vegetation-type resolution will be used for chemical surface fluxes, meteorological land surface modeling, and biogenic emissions.
Chemical, meteorological, and hydrological surface fluxes are key interface processes that link the atmosphere to the terrestrial, aquatic, and marine components of the Earth system.The dry deposition components of air quality models are expanding in scope to include many more chemical species that are important not only for atmospheric constituents but also as inputs to terrestrial and aquatic ecosystems.In addition, there is growing interest in the deposition of toxic contaminants that enter the food chain via air-surface exchange and subsequently impact human health.There has also been an extension to modeling substances that have surface compensation points (e.g., ammonia and mercury) necessitating the development of bi-directional flux capability.The new bi-directional systems often require adjunct land models such as agricultural management models (e.g., EPIC) that are used for ammonia from fertilizer.All of these extensions to the traditional dry deposition models urgently need evaluation and constraint from field and laboratory measurements.

Figure 1 .
Figure 1.Schematic of dry deposition model for soil and leaf.

5 M
od a l M as s M e a n D iam e te r ( m )

Figure 4 .
Figure 4. Average ammonia surface flux measured in a soybean field in Warsaw, North Carolina and modeled by the bi-directional flux scheme in the CMAQ model averaged over 2 months in summer 2002.

Figure 5 .
Figure 5. Schematic of bi-direction flux model for soil and leaf.

Figure 6 .
Figure 6.Ozone dry deposition velocity (cm/s) from a 12 km grid resolution CMAQ simulation on 25 July 2006 at 19 UT.