A modified micrometeorological gradient method for estimating O 3 dry depositions over a forest canopy

Small pollutant concentration gradients between levels above a plant canopy result in large uncertainties in estimated air–surface exchange fluxes when using existing micrometeorological gradient methods, including the aerodynamic gradient method (AGM) and the modified Bowen ratio method (MBR). A modified micrometeorological gradient method (MGM) is proposed in this study for estimating O3 dry deposition fluxes over a forest canopy using concentration gradients between a level above and a level below the canopy top, taking advantage of relatively large gradients between these levels due to significant pollutant uptake in the top layers of the canopy. The new method is compared with the AGM and MBR methods and is also evaluated using eddy-covariance (EC) flux measurements collected at the Harvard Forest Environmental Measurement Site, Massachusetts, during 1993–2000. All three gradient methods (AGM, MBR, and MGM) produced similar diurnal cycles of O3 dry deposition velocity (Vd(O3)) to the EC measurements, with the MGM method being the closest in magnitude to the EC measurements. The multi-year average Vd(O3) differed significantly between these methods, with the AGM, MBR, and MGM method being 2.28, 1.45, and 1.18 times that of the EC, respectively. Sensitivity experiments identified several input parameters for the MGM method as first-order parameters that affect the estimated Vd(O3). A 10 % uncertainty in the wind speed attenuation coefficient or canopy displacement height can cause about 10 % uncertainty in the estimated Vd(O3). An unrealistic leaf area density vertical profile can cause an uncertainty of a factor of 2.0 in the estimated Vd(O3). Other input parameters or formulas for stability functions only caused an uncertainly of a few percent. The new method provides an alternative approach to monitoring/estimating long-term deposition fluxes of similar pollutants over tall canopies.


Introduction
Quantifying atmospheric dry and wet deposition of critical pollutants is important for assessing their life time in air and their potential impacts on various ecosystems.In chemical transport models and in monitoring networks, dry deposition is commonly estimated using the so-called inferential method, which requires a parameter -dry deposition velocity (V d ) -that is typically calculated using empirically developed dry deposition algorithms (Wesely and Hicks, 2000;Pleim and Ran, 2011).Existing dry deposition algorithms have large uncertainties, e.g., a factor of 2.0 on a long-term basis for several commonly studied species (Flechard et al., 2011;Schwede et al., 2011;Wu et al., 2011Wu et al., , 2012;;Matsuda et al., 2006).Field flux measurements are still needed to reduce these uncertainties.
Measurements of O 3 dry deposition flux mostly rely on micrometeorological methods (Wesely and Hicks, 2000).Two types of methods are commonly used: the eddycovariance (EC) technique and the flux-gradient methods.Eddy covariance is a direct measurement method that determines turbulent fluxes without application of any empirical assumption (Baldocchi et al., 1988;Stella et al., 2012).It has been extensively used to estimate turbulent fluxes of Published by Copernicus Publications on behalf of the European Geosciences Union.Z. Y. Wu et al.: A modified micrometeorological gradient method momentum, heat, and trace gases (e.g., CO 2 , H 2 O, SO 2 , O 3 ) (Baldocchi et al., 2001;Turnipseed et al., 2009;Guenther et al., 2011).However, application of EC is often limited by the difficulty of making high-quality measurements at sufficiently high frequencies (i.e., > 1 Hz) to resolve the covariance between the vertical wind velocity and scalar concentration fluctuations (Jacob, 1999).Additionally, the EC method is costly and complex to maintain.
A flux-gradient theory approach, also known as K-theory, was used as an alternative method to determine fluxes of gases that lack the fast response instrument for the EC measurement (Meyers et al., 1996;Park et al., 2014).Fluxgradient theory assumes that the turbulence flux is proportional to the product of the mean vertical concentration gradient and an eddy diffusivity (K) (Baldocchi et al., 1988).The derivation of eddy diffusivity for air pollutants currently relies on the similarity assumption which needs more verification from field measurements.Another critical aspect when employing the flux-gradient theory is to measure the concentrations of gases at different heights with sufficient accuracy and precision (Stella et al., 2012;Loubet et al., 2013).Usually measurements at two adjacent levels above a canopy are used to derive the gradient, e.g., the aerodynamic gradient method (AGM) and the modified Bowen ratio approach (MBR).Due to the small concentration gradient above the canopy and instrument measurement uncertainties, using the flux-gradient method can cause larger uncertainties in the estimated dry deposition fluxes.
Alternatively, gradients between levels above and below the canopy top are usually sufficiently large due to the significant sink in the top layers of forest canopies.Thus, if concentration gradients between levels above and below the canopy top can be used for estimating dry deposition fluxes, the uncertainties might be smaller.The present study aims to develop and evaluate such a method (hereafter referred to as MGM: the modified gradient method).It should be noted that this method is still based on the flux-gradient theory.
Long-term concurrent measurements of eddy-covariance fluxes and concentration profiles for O 3 and CO 2 have been conducted at the Harvard Forest Environmental Measurement Site (HFEMS) since 1990 (Munger et al., 1996;Urbanski et al., 2007).This data set enables us to estimate O 3 dry deposition using the existing (AGM, MBR, and EC) and newly proposed (MGM) methods and thus to evaluate the applicability and uncertainties of all the methods.The micrometeorological methods are briefly described in Sect.2, the measurement data in Sect.3, comparison results and sensitivity tests in Sect.4, and major conclusions and recommendations in Sect. 5.

Eddy-covariance technique
EC determines the turbulent flux (F ) by calculating the covariance between the vertical wind velocity (w) and concentration of the gas (c): where the macron denotes the time average and the primes denote fluctuations from the mean (x = x (t)−x, x = mean).By convention, a positive flux is upward (emission) and a negative flux is downward (deposition).

Aerodynamic gradient method
With an assumption that turbulent transport is analogous to molecular diffusion (Baldocchi et al., 1988), the fluxgradient theory is theoretically described as follows: where K c is the eddy diffusivity for the gas, and dC/dz is the vertical concentration gradient of the gas.Two of the more popular methods for calculating K c are the aerodynamic gradient method and the modified Bowen ratio approach.
The AGM method assumes that heat and mass are transported in a similar way within a well-developed surface layer (Erisman and Draaijers, 1995).K c is related to the interstitial aerodynamic resistance (R a ) (Baldocchi, 1988) as where z 1 and z 2 indicate the heights of adjacent levels above canopy (z1>z2).Using Eqs.
(2) and (3), the deposition flux (F ) is determined as where C 1 and C 2 indicate the gas concentrations at z 1 and z 2 , respectively.R a is calculated as where κ is the von Karman's constant (0.4), u * the friction velocity u * ≡ −u w 1/2 measured at the reference height, d the zero-plane displacement height, L the Obukhov length, and h the integrated stability correction function for heat using those proposed by Businger et al. (1971) and modified by Högström (1988).

Modified Bowen ratio method
The MBR method is also based on the flux-gradient theory (Eq.2), but the eddy diffusivity (K c ) is derived from flux and gradient measurements of another scalar (e.g., sensible heat, CO 2 , H 2 O) and assumes it is equal to K c of the gas of interest.In this study, the flux and gradient measurements of CO 2 are available at the same heights as O 3 , so K c of O 3 was calculated from the CO 2 measurements as follows: where K co 2 is the eddy diffusivity of CO 2 , F co 2 is the eddycovariance flux of CO 2 , C (CO 2 ) is the concentration gradient of CO 2 over the same height interval as C (O 3 ), and z is the height interval of the concentration measurements.Using Eqs. ( 2) and ( 6), the O 3 flux (F ) is calculated as

Modified gradient method
The newly proposed MGM method is also based on the fluxgradient theory (Eq.2).It is noted that the flux-gradient theory has long been questioned within the plant canopy due to infrequent but predominant large eddies within the canopy (Wilson, 1989;Raupach, 1989).For example, Bache (1986) suggested that the flux-gradient theory was a reasonable assumption for estimating wind profiles in the upper portion of the canopy but failed to reproduce the secondary wind maximum that was often observed within the trunk space of forests.It should also be noted that most O 3 uptake occurs in the upper layers of the canopy where most canopy leaves grow.Within these upper layers the vertical length scales of turbulence are probably smaller than the distance associated with changes in concentration and wind speed gradients (Baldocchi, 1988).Thus, the flux-gradient theory is likely applicable for estimating the vertical flux distribution of air pollutants within a plant canopy as has been used in previous studies (e.g., Baldocchi, 1988;Bash et al., 2010;Wolfe and Thornton, 2011).
Applying the flux-gradient theory within the canopy, a height-dependent flux (F (z)) can then be calculated as where z ≤ h, and K c (z) is the vertical eddy diffusivity.Based on Eq. ( 8), the O 3 flux at canopy top (F (h)) is defined as where C h and C 3 are the concentrations at canopy top (h) and the height of z 3 (z 3 <h), respectively.R a (h : According to the aerodynamic gradient method (Eq.4), the O 3 flux above the canopy can be calculated from the concentration gradient between the reference height z 1 and the canopy top h (z 1 > h) as follows: Based on the assumption of a constant flux layer above the canopy, the O 3 flux above the canopy calculated in Eq. ( 11) should be equal to the O 3 flux at the canopy top derived from Eq. ( 9).Using Eqs. ( 9) and ( 11), we can derive that R a (z 1 : h) is calculated using Eq. ( 5).R a (h : z 3 ) is integrated vertically between the two heights within the canopy using Eq. ( 10).
K c (z) is assumed to equal 0.8K m (z), which is the withincanopy eddy diffusivity for momentum transfer (Halldin and Lindroth, 1986).As described in Baldocchi (1988) where a(z) is the leaf area density at height z, and u(z) is the horizontal wind speed within the canopy.Similar to Baldocchi (1988), K m (z) is assumed to be constant below crown closure (about 0.7h) and equal to K m at 0.7h.Thus we suggest here that the level of concentration measurement below canopy (z 3 ) should not be lower that the crown closure of canopy.
The effective drag coefficient (C m (z)) is assumed to be constant with height (see Thom, 1975) following Baldocchi (1988): where LAI is the canopy leaf area index, u m the mean wind speed within canopy, and u(z 1 ) the wind speed at the reference height z 1 .The bulk canopy drag coefficient (C am ) is computed as The mean within-canopy wind speed (u m ) is calculated as The within-canopy wind speed profile (u(z)) follows Cionco (1972): where u h is wind speed at the canopy top, and α is wind speed attenuation coefficient.The above-canopy logarithmic wind profile is used to scale the wind speed measured at the reference height z 1 to the canopy height h: where z 0 is the roughness length for momentum, and m is the integrated stability correction function for momentum as proposed by Businger et al. (1971) and modified by Högström (1988).
Assuming a zero concentration on the absorbing surface, the dry deposition velocity (V d ) of O 3 can be determined as where C(z 1 ) is the O 3 concentration measured at the reference height z 1 .
3 Field measurements used in this study

Site description
The HFEMS (42.54 • N, 72.18 • W) is located in central Massachusetts at an elevation of 340 m above sea level.The forest is 80 years old on average and consists of red maple (Acer rubrum) and red oak (Quercus rubra) with scattered stands of Eastern hemlock (Tsuga canadensis), red pine (Pinus resinosa), and white pine (Pinus strobus).The canopy height near the observation tower is up to 23 m with a peak LAI of ∼ 5.0 m 2 m −2 during summer.The nearest sources of significant pollution are a secondary road about 2 km to the west of the site and a main highway about 5 km to the north.A permanent 30 m Rohn 25G tower has been utilized at HFEMS to measure eddy-covariance fluxes of sensible heat, H 2 O, momentum, CO 2 , and O 3 along with vertical profiles of CO 2 and O 3 since 1990 (Fig. 1).Eddy-covariance fluxes were measured at a height of 29 m above the ground.For the profile measurements air was continuously sampled from heights of 29, 24.1, 18.3, 12.7, 7.5, 4.5, 0.8, and 0.3 m a.g.l. to determine the concentrations of CO 2 and O 3 .In this study, the upper three levels were used to derive the gradients.Details on the site and the instrumental methods can be found in Munger et al. (1996).Data used in this study are available online at http://atmos.seas.harvard.edu/lab/data/nigec-data.html.Zhao et al. (2011) retrieved the vertical profile of leaf area density at Harvard Forest from a ground-based lidar scanning.Two tree species groups (i.e., hardwood and conifer) were chosen.According to the species composition around the measurement tower, the average leaf area density used in this study was calculated as 75 % of that of hardwood and The monthly averaged LAI at HFEMS was derived from the ground-based measurements for most years between 1998 and 2013 using the LICOR LAI-2000 system at 30-40 plots around the tower (Urbanski et al., 2007).As the measurements during January and February were not available, these values were obtained based on extrapolation (Fig. 2).The roughness length (z 0 ) and displacement height (d) were calculated as a function of canopy height (h) and LAI, following Meyers et al. (1998) (see Fig. 2):

Data selection
A total of 10 252 hourly measuring points, recorded at HFEMS during 1993-2000, were screened to eliminate the influence of periods associated with instrumental and measurement problems and violation of the use of the fluxgradient theory.
In order to reduce random measurement error in the concentration gradient, O 3 concentrations below 1 ppbv were rejected, resulting in approximately 0.1 % of the data being omitted.In addition, periods with [O 3 ] < [NO y ] (1.9 %) were excluded to avoid periods when O 3 chemical reactions may have exceeded O 3 deposition (Munger et al., 1996).Wind speed below 1.0 m s −1 (1.2 %) and drag coefficient below 0.02 (6.6 %) were removed because of probable invalid fluxgradient relationships (Feliciano et al., 2001).Outliers in the data (2.9 %) were removed, omitting any deposition velocity exceeding the maximum achievable deposition velocity ) by more than a factor of 1.5 (Matsuda et al., 2006).Periods with counter-gradient profiles (69.8 %) which represent a downward flux (from EC measurements) while with a negative gradient (upper level minus lower level) or vice versa were rejected (Park et al., 2014).The counter-gradient transport should be mainly due to the non-local nature of turbulent transport within canopies.Large sweep-ejection air motions, associated with coherent structures that can deeply penetrate into the canopy, are believed to be largely responsible for the exchange of momentum, heat, and mass between air above-and within-canopy (e.g, Shaw et al., 1983;Thomas and Foken, 2007).A total of 74.0 % of the data was omitted in the following analysis.This percentage value is slightly smaller than the sum of those from all the criteria due to the overlap of some data points between the criteria.
Figure 3 shows the mean diurnal cycles of O 3 concentration at different heights derived from the original data set and from the data after selection.The O 3 concentration typically increased during the early morning to reach a daily maximum of over 40 ppbv in the early afternoon and then decreased to ∼ 30 ppbv at night.As shown in Fig. 3a, the gradient between the two heights above the canopy (i.e., 29 and 24.1 m) was only about 0.4 ppbv on average, smaller than that between the levels above the canopy (24.1 m) and inside the canopy (18.3 m) (∼ 0.8 ppbv).The gradients were relatively small during the morning (e.g., 0.1 ppbv at 11:00 LST) compared to the other periods of the day.In the morning, the most effective turbulent exchange between the air aboveand within-canopy would substantially reduce the gradients (Sörgel et al., 2011).It is worth mentioning that many earlier studies suggested that the effects of chemistry on O 3 flux divergence in the near surface were generally small, likely because the chemical reactions for O 3 have larger timescales than the turbulent transport (e.g., Gao et al., 1991;De Arellano and Duynkerke, 1992;Duyzer et al., 1997;Padro et al., 1998;Stella et al., 2012).After screening the data with the criteria, the gradients among these three levels were signif- icantly larger, reaching up to 1.0 and 1.6 ppbv, respectively (see Fig. 3b).

Comparison of V d (O 3 ) by the eddy-covariance and gradient methods
O 3 dry deposition velocity (V d (O 3 )) measured by the eddycovariance technique at Harvard Forest typically ranged from 0.14 to 0.53 cm s −1 with a median value of 0.30 cm s −1 during the study period (Table 1).Since the screened deposition velocities still include certain outlying data, the mean value was calculated using data between the 10th and 90th percentiles in order to reduce the influence of the outlying data.Following this approach, the mean V d (O 3 ) by the EC technique was 0.34 cm s −1 , which was significantly smaller than the means calculated by the gradient methods (Table 1).The ratios of the mean V d (O 3 ) calculated by the modified gradient, modified Bowen ratio, and aerodynamic gradient methods to the mean of the EC technique were 1.18, 1.45, and 2.28, respectively.Previous studies on the inter-comparisons of these methods for O 3 are few and the results varied.Muller et al. (2009) found that the mean V d (O 3 ) by the AGM method was 1.60-3.47times that of the EC technique at a grassland in southern Scotland.Loubet et al. (2013) showed that the AGM method gave 40 % larger V d (O 3 ) than the EC technique over a mature maize field in Paris.Keronen et al. (2003) found that V d (O 3 ) by the AGM and EC methods generally agreed well at a Nordic pine forest, as did Stella et al. (2012) over bare soil in Paris.Droppo (1985) found close V d (O 3 ) values with the MBR and EC methods at a northeastern US grassland site.
Figure 4 shows the diurnal cycles of V d (O 3 ) by the EC and gradient methods.Although the trends were similar, the MBR and AGM V d (O 3 ) were consistently larger than the EC V d (O 3 ).The EC V d (O 3 ) was about 0.2 cm s −1 on average during the night and reached a daily maximum of 0.5 cm s −1 around noon.The V d (O 3 ) by the MBR and AGM methods reached around 0.8 and 1.3 cm s −1 during the daytime, respectively, and remained about 0.4 cm s −1 during the night.The MGM V d (O 3 ) agreed well with the EC V d (O 3 ) during the daytime but was slightly larger at night.This discrepancy has been identified in previous studies (Keronen et al., 2003;Stella et al., 2012) and could be due to the fact that nocturnal conditions affect both EC and gradient measurements.The EC technique has been found to underestimate flux during calm nighttime periods at Harvard Forest (Goulden et al., 1996).The stability correction functions used in the gradient methods (AGM and MGM) are subject to large uncertainties under stable conditions (Högström, 1988).
The very large differences in V d (O 3 ) between the AGM and EC methods are likely caused by a combination of various factors.As can be seen from Eq. ( 4), any underestimation in the calculation of aerodynamic resistance (R a ) would directly result in the overestimation of V d .Uncertainties in R a from using different formulas are generally on the order of 30 % over a whole canopy (Zhang et al., 2003).In the case of by the eddy covariance (EC) and three gradient methods (MGM: the modified gradient method; MBR: the modified Bowen ratio method; AGM: the aerodynamic gradient method).In each box, the central mark is the median, and the edges of the box are the 10th and 90th percentiles.Note that the average is the arithmetic mean of the data between the 10th and 90th percentiles.
Eq. ( 4), uncertainties can be larger than 30 % if other uncertainties from the related parameters are larger.The potential underestimation of R a (Eq.4) also explains the small overestimation in V d from the MGM method, in which the same R a formula is used, although it plays a secondary role.Measurement uncertainties in the concentration gradients could also cause large discrepancies between the AGM and EC methods, especially under small gradient conditions.This is supported by the finding that the MBR method also overestimated V d when compared to the EC measurements.
As shown in Fig. 5, the EC V d (O 3 ) exhibited a significant seasonal pattern with peak values in summer (∼ 0.5 cm s −1 ) and small values in winter (0.15-0.28 cm s −1 ).Both the MGM and MBR methods captured this seasonal cycle, but the MGM method produced a higher V d (O 3 ) than the EC technique during winter (December-February) and the MBR method produced a significant overestimation in summer (June-September).The monthly AGM V d (O 3 ) was consistently larger than the EC V d (O 3 ) and exhibited a less clear seasonal pattern with alternating increases and decreases in the V d (O 3 ).

Sensitivity of V d (O 3 ) by the modified gradient method to the key parameters/formulas
As shown in Sect.4.1, the MGM method performed better than the MBR and AGM methods.This improvement is mainly attributable to smaller errors in the O 3 concentration gradients.However, the MGM method increased the com- b Vertical profile of leaf area density from Meyers et al. (1998) as shown in Fig. 7. c D74: Dyer (1974); P70: Paulson (1970); W70: Webb (1970).d The arithmetical mean of the data between the 10th and 90th percentiles.[1993][1994][1995][1996][1997][1998][1999][2000] as measured by the eddy covariance (EC) and three gradient methods (MGM: the modified gradient method; MBR: the modified Bowen ratio method; AGM: the aerodynamic gradient method).Note that the average is the arithmetical mean of data between 10th and 90th percentiles.
plexity of the algorithm and added more model parameters, which may in turn increase the uncertainty of the estimated V d (O 3 ).
To test the sensitivity of the MGM V d (O 3 ) values to the key parameters/formulas, calculations were conducted by changing the parameters/formulas within a reasonable range.For some single-value parameters (i.e., roughness length, displacement height, wind speed attenuation coefficient, and leaf area index), sensitivity tests were conducted by increasing or decreasing the value by 10 %.
As shown in Fig. 6 and Table 2, the MGM V d (O 3 ) was highly sensitive to the changes in wind speed attenuation coefficient and displacement height.Higher wind speed attenuation coefficients could result in lower within-canopy wind speed (Eq.17) and thus a lower eddy exchange coefficient and V d (O 3 ) (Table 2).Based on a least-square fitting of within-canopy wind profiles measured at Harvard Forest for noon periods in summer, the attenuation coefficient was estimated to be ∼ 10.6.Cionco (1972) suggested that the attenuation coefficient varies with leaf area.Therefore, the application of this value throughout the whole year could produce an uncertainty in the estimated V d (O 3 ).
The MGM V d (O 3 ) increased when the displacement height increased or vice versa (Fig. 6, Table 2).Sakai et al. (2001) calculated the displacement height at Harvard Forest using noon-period measurements and indicated the ratio of displacement height to canopy height was 0.77 in summer  (Dyer, 1974) ψ h (Paulson, 1970) ψ h (Webb, 1970) (Dyer, 1974) ψ m (Paulson, 1970) ψ m (Webb, 1970)  with foliated canopy and 0.6 in winter with leafless canopy.In this study, we estimated a similar value in summer (0.79) and a slightly higher value in winter (0.66) using the method proposed by Meyers et al. (1998) (Fig. 2).The overestimation of the displacement height could partly explain the overestimation of V d (O 3 ) by the MGM method during December to February (Fig. 5).
Figure 6 shows that the MGM V d (O 3 ) was less sensitive to the changes in roughness length and leaf area index.The relative differences in the estimated V d (O 3 ) were less than 2 % when the roughness length and leaf area index varied by 10 % (Table 2).Harvard Forest Profile 1 (Meyers et al., 1998) Profile 2 (Meyers et al., 1998) Profile 3 (Meyers et al., 1998) Figure 7. Vertical profiles of leaf area density in Harvard Forest and those used in sensitivity experiments.Meyers et al. (1998) provided three typical types of leaf area density profiles, which are significantly different in shape from the profile in Harvard Forest used in this study (see Fig. 7).We conducted sensitivity experiments by replacing the Harvard Forest profile with those of Meyers et al. (1998) to assess the impact of the vertical profile of leaf area density on the determination of V d (O 3 ).As shown in Fig. 6 and Table 2, the vertical profile of leaf area density impacted the estimated V d (O 3 ) greatly, with a relative difference in V d (O 3 ) of greater than 50 %.The profile with higher leaf density in the upper canopy (profile 3) resulted in a higher V d (O 3 ) while the profile with abundant understory plants (profile 1) led to a lower V d (O 3 ).
In this study, the stability correction functions proposed by Businger et al. (1971) and modified by Högström (1988) were used, but several others exist, such as those by Dyer (1974), Paulson (1970), andWebb (1970).Figure 6 indicated that uncertainties in the stability correction functions for heat ( h ) and momentum ( m ) had little impact on the MGM V d (O 3 ) values.The relative difference of V d was less than 4 % for different h and less than 1 % for different m .Stella et al. (2012) found that the variation of V d (O 3 ) on different h was roughly 10 % on average when using the AGM method.h influences the estimation of V d due to its impact on the calculation of turbulent transfer above the canopy.Since the MGM method considered both the aboveand within-canopy turbulent transfer, the MGM V d (O 3 ) values were thus less sensitive to the choice of h .

Conclusions and recommendations
A modified micrometeorological gradient method was developed to quantify O 3 dry depositions over a forest canopy, making use of concentration gradients between levels above and below the canopy top.The MGM method produced V d (O 3 ) values close to the eddy-covariance measurements at Harvard Forest during daytime but slightly overestimated the measurements at night.The modified gradient method seemed to be an improvement compared to the two existing flux-gradient methods (AGM and MBR) in terms of predicting the long-term mean, diurnal, and seasonal cycles of V d (O 3 ).Sensitivity tests showed that the model parameters for MGM including the wind speed attenuation coefficient, canopy displacement height and vertical distribution of leaf density were first-order parameters affecting the V d (O 3 ) estimates.Model results were less sensitive to roughness length, leaf area index, and stability function for heat and momentum.
The newly developed MGM method has the potential to be applied routinely to monitor/estimate long-term deposition fluxes of O 3 and other similar pollutants over tall canopies.The within-canopy measurement height should be close to but not lower than the canopy closure height where most of the flux exchange occurs.Key model parameters mentioned above need to be characterized as accurately as possible.For example, seasonal profiles of the vertical distribution of leaf area density, canopy displacement height, and vertical wind profile related parameters are needed.

Figure 1 .
Figure 1.Schematic of flux and concentration gradient measurements at Harvard Forest Environmental Measurement Site.

Figure 2 .
Figure 2. Monthly variation of leaf area index (LAI), the displacement height (d) to canopy height (h) ratio, and the roughness length (z 0 ) to canopy height ratio at Harvard Forest.

Figure 3 .
Figure 3. Mean diurnal cycles of O 3 concentration at heights of 29, 24.1, and 18.3 m a.g.l. at Harvard Forest during 1993-2000.Panel (a) was derived from the original data and (b) was from the data after selection.
Figure 4. (a) The box plot of hourly V d (O 3 ) and (b) diurnal average cycles of V d (O 3 ) at Harvard Forest during 1993-2000 as measuredby the eddy covariance (EC) and three gradient methods (MGM: the modified gradient method; MBR: the modified Bowen ratio method; AGM: the aerodynamic gradient method).In each box, the central mark is the median, and the edges of the box are the 10th and 90th percentiles.Note that the average is the arithmetic mean of the data between the 10th and 90th percentiles.

Figure 6 .
Figure 6.Diurnal average cycles of V d (O 3 ) over Harvard Forest during 1993-2000 by the modified gradient method (MGM) with different parameter/formula changes and compared with that by the eddy-covariance (EC) technique: (a) roughness length, (b) displacement height, (c) wind speed attenuation coefficient, (d) leaf area index, (e) vertical profile of leaf area density, (f) stability correction functions for heat, and (g) stability correction functions for momentum.

Table 1 .
Statistics of hourly V d (O 3 ) (cm s −1 ) at Harvard Forest during 1993-2000 as measured by the eddy covariance (EC) and three gradient methods (MGM: the modified gradient method; MBR: the modified Bowen ratio method; AGM: the aerodynamic gradient method).

Table 2 .
Relative difference between V d (O 3 ) determined by the modified gradient method with different parameters/formulas (%) a .