First direct measurements of formaldehyde ﬂux via eddy covariance: implications for missing in-canopy formaldehyde sources

. We report the ﬁrst observations of formaldehyde (HCHO) ﬂux measured via eddy covariance, as well as HCHO concentrations and gradients, as observed by the Madison Fiber Laser-Induced Fluorescence Instrument during the BEACHON-ROCS 2010 campaign in a rural, Ponderosa Pine forest northwest of Colorado Springs, CO. A median noon upward ﬂux of ∼ 80 µg m − 2 h − 1 ( ∼ 24 ppt v m s − 1 ) was observed with a noon range of 37 to 131 µg m − 2 h − 1 . Enclosure experiments were performed to determine the HCHO branch (3.5 µg m − 2 h − 1 ) and soil (7.3 µg m − 2 h − 1 ) direct emission rates in the canopy. A zero-dimensional canopy box model, used to determine the apportionment of HCHO source and sink contributions to the ﬂux, underpre-dicted the observed HCHO ﬂux by a factor of 6. Simulated increases in concentrations of species similar to monoterpenes resulted in poor agreement with measurements, while simulated increases in direct HCHO emissions and/or concentrations of species similar to 2-methyl-3-buten-2-ol best improved model/measurement agreement. Given the typical diurnal variability of these BVOC emissions and direct HCHO emissions, this suggests that the source of the missing ﬂux is a process with both a strong temperature and radiation dependence.


Introduction
The oxidation of volatile organic compounds (VOCs) in the atmosphere occurs via the HO x -NO x cycle, a photochemically-driven catalytic cycling of hydrogen oxide (OH + HO 2 ) and nitrogen oxide (NO + NO 2 ) radicals. This process produces tropospheric ozone and oxidized VOCs, the latter of which may condense to form secondary organic aerosol (SOA) (Zhang et al., 2007;Jimenez et al., 2009). To accurately model both tropospheric ozone and SOA, the processes involved in VOC oxidation must be characterized. Part of the difficulty in understanding this cycle lies in the detection and quantification of all relevant species of VOCs, particularly in forest environments. Multiple studies have reported a significant discrepancy between measured and modeled OH concentrations and reactivities, suggesting errors in our understanding of the emissions or processing of VOCs Di Carlo et al., 2004;Hofzumahaus et al., 2009;Lelieveld et al., 2008;Sinha et al., 2010;Whalley et al., 2011). This discrepancy may be related to the fast in-canopy oxidation of unmeasured biogenic VOCs (BVOCs), specifically terpenes (Di Carlo et al., 2004;Goldstein et al., 2004;Holzinger et al., 2005). To confirm this, a method of determining the overall VOC oxidation rate is needed.
As a result, HCHO is an excellent tracer for overall VOC oxidation. Quantification of HCHO production in forest environments could provide a valuable constraint for the overall rate of VOC oxidation in this environment. There have been many reports of forest HCHO mixing ratios (Munger et al., 1995;Slemr et al., 1996;Lee et al., 1998;Sumner et al., 2001;Galloway et al., 2011), but a qualitative and quantitative understanding of in-canopy HCHO production is still incomplete. One recent study (Choi et al., 2010) reported a missing boundary layer HCHO production rate of as much as 1.6 ppb v h −1 , nearly double the calculated chemical production rate.
Measurements of HCHO vertical fluxes above and gradients throughout a forest canopy may yield valuable insight into production and loss of HCHO inside the canopy. Gradient measurements can give more detailed information about the sources and sinks in the canopy, while vertical flux measurements are less influenced by advection, as the area sampled by the flux is typically the area less than a kilometer upwind. HCHO fluxes have previously been estimated based on flux-gradient calculations over polar icepack (Jacobi et al., 2002;Hutterli et al., 2004), but there has been little work examining HCHO distribution in forest canopies and no reported measurements of HCHO flux by eddy covariance (EC). Of the many reported techniques to measure HCHO (Table S1: Weibring et al., 2007;Wisthaler et al., 2008;Hottle et al., 2009;McManus et al., 2010), none have reported the capability of performing the fast sampling needed for EC measurements with both the sensitivity needed to quantify small perturbations in HCHO concentration and the selectivity inherent to spectroscopic techniques.
In this work, we present HCHO gradients and EC flux observations using Fiber Laser-Induced Fluorescence (FILIF), which has the high sensitivity and high time resolution needed for EC measurements. Additionally, we discuss branch and soil enclosure experiments performed to determine HCHO emission rates. To model HCHO flux, we present a zero-dimensional box model used to apportion HCHO production and loss inside the canopy. Finally, we discuss sensitivity studies with respect to both BVOC and direct HCHO emission using the box model to ascertain their effect on measurement/model agreement.

Field campaign
All observations reported here were taken during the Biohydro-atmosphere interactions of Energy, Aerosols, Carbon, H 2 O, Organics & Nitrogen -Rocky Mountain Organic Carbon Study (BEACHON-ROCS) field campaign during 1-31 August 2010 at Manitou Experimental Forest (MEF, 39°06 02 N, 105°06 05 W, 2286 m), northwest of Colorado Springs, CO. The site has been described in detail elsewhere (Kim et al., 2010). It is located in a Central Rocky Mountains Ponderosa Pine (PPine) forest (canopy height: ∼18.5 m; leaf area index (LAI) = 1.9) with minimal undergrowth, predominately clean air masses transported from the southwest, and rare anthropogenic incursions.

Fiber Laser-Induced Fluorescence (FILIF) of HCHO
This technique is similar to that reported by Hottle et al. (2009), the primary difference being the laser, and will only be described briefly here. The 353 nm tunable, pulsed, and narrow-bandwidth fiber laser (NovaWave Technologies, TFL Series) represents a significant improvement over previous field laser technology, as fiber lasers are inherently lighter, smaller, and more stable than traditional lasers. The ∼10 mW laser was directed into a 32-pass White-type multipass cell, and the resulting HCHO fluorescence from 390 to 500 nm was filtered using a 390 nm longpass filter then focused into a photomultiplier tube for detection. Laser power was monitored both before and after the multipass cell using photodiodes, and a fraction (∼1 mW) of the outgoing beam was directed into a cell filled with concentrated gas-phase HCHO for wavelength reference. The separation between the multipass cell mirrors was ∼25 cm. However, only ∼6 cm of each pass was through the 6 cm × 5 cm area (cell depth: ∼6 cm) through which the ambient air was flowed perpendicular to the narrow plane of the laser. The residence time of air in the cell was <25 ms in the beam volume (∼1 cm thick) at the ∼12 standard liters per minute (SLM) sampling flow rate.
Remaining volumes of the cell were purged using a zero air generator (AADCO 737-series) with a total purge flow of 500 standard cubic centimeters per minute (SCCM) regulated by a mass flow controller (MKS Instruments, M100B). Measurements were performed by dithering the laser on and off a rovibronic absorption line at 353.37 nm. The difference in fluorescence signal when the laser was centered on these two positions was proportional to the HCHO concentration. Instrument calibrations were performed weekly using a HCHO permeation tube (VICI Metronics, 100-044-2300-U45) heated to 85°C using a portable calibration gas generator (VICI Metronics, Model 120). The output of the permeation tube device as characterized by Fourier Transform Infrared (FTIR) spectroscopy was found to be 438±7 ngmin −1 ; details on this calibration can be found in the Supplement. The calibration factor varied by less than 2.5 % over the course of the campaign. Field 3σ limits of detection were typically on the order of ∼300 ppt v in 1 s, with measurement accuracies of ∼20 % limited by that of the permeation tube calibration.
Inlets for HCHO sampling with lengths of 30 to 45 m were located at heights of 25.1 m, 17.7 m, 8.5 m, and 1.6 m. Inlets consisted of ∼30 m 3/8 ID PFA Teflon tubing, short lengths of which have been found to have a negligible effect on sampling (Wert et al., 2002). To test for possible artifacts, both a 15 m and a 30 m inlet were collocated outside the instrument trailer; resulting measurements agreed within 1.5 %. Typically, ambient flow while sampling through an inlet was ∼12 SLM. Inlets were continuously purged with ambient air at ∼3 SLM when not in use. An additional scroll pump (Gast Manufacturing) with an average flow of ∼80 SLM was used to increase the flow rate of the 25.1 m inlet used for EC sampling to reduce residence time and prevent laminar flow in the inlet. The 25.1 m inlet was placed ∼0.1 m below and ∼0.5 m upwind in the primary wind direction of the center of the sonic anemometer (see Sect. 2.3).
Measurements were performed in an hourly cycle for 11-22 August. During this period, HCHO was measured from the 25.1 m inlet for the first 35 min with online and offline sampling times of 10 s and 1 s respectively at 10 Hz (for EC), following which was a 1.5 min diagnostic period. Then, each of the other three inlets was sampled sequentially with online and offline sampling times of 20 s and 10 s respectively for 7 min, following each of which was a 1.5 min diagnostic period. During 23-30 August, only EC measurements were performed with 35 min collection periods and 1.5 min diagnostic periods. For this period, as eddies with timescales on the order of 10 s contributed significantly to HCHO flux (see Sect. 2.4.2), online and offline sampling times were changed to 290 s and 5 s, respectively. This change in sampling was to test for potential EC spectral interference, which was not observed.

Other measurements
Unless otherwise noted, all other measurements used a valve switching system which changed sampling lines every 5 min and cycled through six 1/4 OD Teflon inlets mounted at 25.1 m, 17.7 m, 12.0 m, 8.5 m, 5.0 m, and 1.6 m over a 30 min period. Flow rates of ∼3.5 SLM through the sampling lines resulted in delay times between 8 to 12 s, measured by spiking a VOC pulse at each sampling inlet.
A Proton-Transfer-Reaction Mass Spectrometer (PTR-MS, Ionicon Analytik GmbH) was used for gradient measurements of selected VOCs. The instrument is based on soft chemical ionization using protonated water ions (H 3 O + ) Lindinger et al., 1998), and was operated at 2.3 mbar drift pressure and 540 V drift voltage and calibrated using two multi-component ppm v VOC standards (Karl et al., 2009).
OH, HO 2 , and RO 2 were measured using chemical ionization mass spectrometry (CIMS) as described by Tanner et al. (1997) and Hornbrook et al. (2011). The CIMS acquired measurements at ∼10 m from the tower at a height of 2.7 m with the inlet facing perpendicular to the primary wind direction. During periods with OH concentrations below the detection limit (5 × 10 5 molec cm −3 ), OH concentration was assumed to be equal to half the detection limit (2.5 × 10 5 molec cm −3 ).
Downwelling NO 2 photolysis (J NO 2 ) was measured from the top of the 30 m chemistry tower with commerciallyavailable filter radiometers (Meteorologie Consult GmbH) as described by Junkermann et al. (1989) and Volz-Thomas et al. (1996). The filtered measurement was converted to a photolysis rate by comparison with spectrally-resolved actinic flux measurements. Total J NO 2 was estimated by measurement of the ratio of upwelling to downwelling J NO 2 as measured from the tower on 10 August 2010.
OH reactivity was measured using a laser-induced pump and probe technique (Sadanaga et al., 2004) at ∼20 m from the tower and a height of ∼4 m with a 2 min sampling rate. Peroxyacetyl nitrate (PAN) was measured via Thermal Decomposition-Chemical Ionization Mass Spectrometry, as described by Zheng et al. (2011). Ozone concentrations were measured using a Model 205 Dual Beam Ozone Monitor (2B Technologies, Inc.). NO concentrations were measured using an Ecophysics CLD-88Y analyzer. NO 2 concentrations were measured using a Droplet Measurement Technologies Blue Light Converter. A LI-COR LI-7000 measured CO 2 and H 2 O concentrations at 25.1 m. LI-COR LI-190 quantum sensors measured photosynthetically active radiation (PAR) at 27.8 m and 1.8 m. Vaisala HMP35C probes measured temperature and relative humidity at 25.3 m and 7.0 m. A Vaisala PTB101B barometer measured barometric pressure. A sonic anemometer (Campbell Scientific, CSAT-3) at 25.1 m measured the three-dimensional wind vector, as well as virtual temperature, at 10 Hz.

Eddy Covariance measurements
Eddy covariance (EC) is a widely-used micrometeorological technique for direct measurement of surface-atmosphere exchange and will be discussed here briefly; further information is available elsewhere (Baldocchi et al., 1988;Lee et al., 2004). EC uses the covariance between vertical fluctuations in wind speed, caused by atmospheric eddies, and fast variations in tracer concentration to extract the mass transport through the plane of measurement. Quantitatively, the turbulent flux of a species at a single height, assuming horizontal and vertical advection is negligible, is defined as: where w is the vertical wind speed, c is the tracer concentration, and x is the instantaneous deviation of x from the ensemble mean value (i.e. x = x −x). For this study, a sonic anemometer measured vertical wind speed, while the HCHO FILIF instrument measured tracer concentration. As eddies occur on a wide range of timescales, the averaging time to calculate the ensemble mean and fluctuating quantities can vary depending on measurement height (Berger et al., 2001). For this study, a sampling period of ∼32 min was chosen, the validity of which will be discussed in Sect.

Data reduction
Three-dimensional 10 Hz wind speeds from the sonic anemometer were rotated using the natural wind coordinate  for each 35 min flux period. Sampling periods with a friction velocity (u * ) less than 0.2 m s −1 were neglected as rotation has been shown to result in poor data quality at low wind speeds . Vertical rotation angles (e.g. tilt angles) were typically ∼ 2 ± 4°. Additionally, a delay exists between HCHO concentration and wind speed due to the residence time of the HCHO sample in the inlet tubing. A correction was determined empirically by calculating w HCHO at different time delays, or lags, to find the maximum in covariance, as shown in Fig. 1, which should be roughly equal to the residence time in the inlet tubing . In this study, an additional variable lag was present as the computers recording the sonic anemometer and HCHO data were not synchronized. This resulted in a lag time that varied considerably over the campaign. Therefore, it was necessary to divide the dataset into 4 sections with different linear trends, depending on computer resynchronization time. A sampling period was considered to have a "good" lag when the covariance was greater than 20 µg m −2 h −1 and the u * was greater than 0.3 m s −1 , and these points were used to calculate the linear trends. All sampling periods were then assumed to have a lag according to these trends. Lag times over the course of the campaign ranged from −9.4 to 1.1 s. Finally, the EC data was tested for stationarity (Foken and Wichura, 1996) by dividing each 30 min sampling period into ∼5 min periods. The average of the 5 min flux measurements for each sampling period was compared to the 30 min flux measurement for that period.
The period was considered stationary if the fluxes agreed within 30 % (Foken and Wichura, 1996). Non-stationary periods were rejected as invalid and not included in the analysis, which resulted in the removal of 48 % of daytime and 60 % of nighttime data.

Spectral analysis
To determine the validity of the remaining flux data, the cospectra of the HCHO fluxes, which may be thought of as the frequency-dependent covariance between the species, were investigated in further detail. As the cospectrum over a single period was typically quite variable, cospectra were averaged over multiple periods. Figure 2 shows the average cospectrum for HCHO and virtual temperature fluxes over daily periods during the campaign. The linear regions of each cospectrum indicative of the inertial sublayer (f >∼0.04 Hz) exhibited a lower slope than the expected value determined from a −7/3 power law . A similar effect was observed at Blodgett Forest (Farmer et al., 2006;Wolfe et al., 2009) and attributed to wake-generated turbulence present in forest canopies (Kaimal and Finnigan, 1994). Figure 2 also demonstrates that while the overall covariance The frequencies of these eddies are highly variable between different sampling periods, which makes it difficult to determine a cause. However, this variability also suggests that these negative covariance events are not likely an artifact of data collection, as this would likely result in consistent negative fluxes at a given frequency over multiple periods. As the field mission averaged cospectrum ( Fig. S1) closely resembles that of the virtual temperature flux, the negative covariance events are believed to have been due to atmospheric variability. These negative values may also be responsible for the faster drop-off of the normalized cospectrum relative to that for the temperature flux ( Fig. 3a). Spectral attenuation may be observed either when a sampling period is too short to sample low-frequency eddies, or when the sample rate is too slow to sample high-frequency eddies. The frequency-weighted cospectrum ( Fig. 3a) peaked at the frequencies contributing most to total flux. For the HCHO cospectrum, three peaks were observed, corresponding to characteristic eddy timescales of ∼0.5, 2.5, and 8 min. The 0.5 min peak corresponds to the peak in the temperature flux cospectrum as well as in the momentum flux cospectrum, which likely indicates that this is the integral time scale for turbulent transport. The 2.5 min peak is similar in timing to one observed in PAN cospectra during other forest campaigns (Turnipseed et al., 2006;Wolfe et al., 2009), which was on the same timescale of observed canopy sweep events (Holzinger et al., 2005). The 8 min peak can likely be attributed to a similar phenomenon. At frequencies greater than 0.04 Hz, the cospectrum appears to decrease more quickly than the temperature flux cospectrum. The cause of this is not understood, but similar high frequency loss has been observed in PAN cospectra (Turnipseed et al., 2006). If this loss is a result of spectral attenuation, it implies that fluxes are typically underestimated. By comparing the difference between the integrated areas of the virtual temperature and HCHO weighted cospectra, this underestimate would be on the order of ∼12 %, which was included in the error analysis as a systematic low bias.
The cospectral cumulative distribution function, or ogive ( Fig. 3b), is the cumulative contribution to the flux as a function of frequency. The HCHO ogive is significantly shifted towards lower frequencies compared to the temperature ogive indicating greater contribution to the flux by lower frequency eddies than for temperature. The lack of an asymptote toward the low frequency end of the ogive implies that the sampling period may not have been sufficient to capture all of the low frequency eddies. However, analysis during the last half of the campaign with longer sampling periods resulted in no significant gain in covariance with periods greater than 30 min.
Other potential errors in the flux measurements are discussed in the Supplement (see Sect. S2). By summing the systematic errors (response, sensor, dampening, attenuation), then propagating with this the indeterminate errors (instrument noise, lag time, calibration), we calculated the total error in the HCHO flux to be typically ∼38 %.

Gradient and flux profile
Daily HCHO fluxes typically showed a symmetric diurnal efflux centered at noon. The median diurnal profile of HCHO flux is shown in Fig. 4a, while the full flux time series is shown in Fig. S2 (note that positive values denote an upward flux while negative values denote a downward flux). Median noontime fluxes were ∼80 µg m −2 h −1 (∼24 ppt v m s −1 ) with maxima as high as ∼170 µg m −2 h −1 (∼50 ppt v m s −1 ). For comparison, 2-methyl-3-buten-2-ol (MBO) fluxes have been observed in PPine forests on the order of 8 to 9 mg m −2 h −1 (Baker et al., 1999;Schade et al., 2000). HCHO fluxes were also observed to have a significant dependence on both temperature and PAR (Fig. S3). Measured HCHO fluxes correspond to a median noontime net HCHO production rate of ∼3.2 ppb v h −1 below the measurement height of 25.1 m. However, the net HCHO production rate into the boundary layer from these fluxes, assuming a boundary layer height of ∼1 km, is only ∼0.079 ppb v h −1 . This is small compared to the 2 to 3 ppb v h −1 total boundary layer production rates reported by the literature Choi et al., 2010), implying that HCHO fluxes have only a small effect on boundary layer concentrations. Figure 4b shows the median diurnal HCHO concentrations for each measurement height. Nighttime hours show lower concentrations near ground level, suggesting dominance of in-canopy sinks such as deposition. The peak in concentration around 08:00 corresponds to increased wind speed and emission of precursors, followed by a sharp change in wind direction. For most of the day, a negative gradient is present, with higher concentrations near the ground. Daytime HCHO concentrations at the ground level (1.6 m) inlet were typically 15 to 20 % higher than concentrations in or above canopy. Qualitative testing of campaign-related ground equipment (e.g. tarps) at the end of the campaign suggests negligible emissions from these materials, and there was little ground level vegetation near the site. This implies a significant direct and/or photochemical ground litter source of HCHO or a significant difference in deposition loss between inside and below the canopy, with the former supported by semi-quantitative testing of the ground litter at the end of the campaign (Sect. 3.2). Canopy level enhancement of HCHO concentration was also observed in the leafy part of the canopy, likely due to either fast oxidation of emitted BVOCs or direct emission from the canopy.
The median diurnal profiles of the flux and concentration measurements do not appear to exhibit the same diurnal variation. During periods of changing wind speed and direction occurring during early morning and mid-evening, the concentration profile changes significantly while the flux profile does not. While these changes in airmass seem to affect HCHO concentration, they have little effect on the flux, implying advection is a negligible contributor to HCHO flux. J. P. DiGangi et al.: HCHO Flux via Eddy Covariance This is supported by no significant correlation between tracers of advection, such as SO 2 , CO 2 , or H 2 O, and HCHO flux. Nighttime deposition gradients are not reflected as negative fluxes, as nighttime flux observations, even those with significant turbulence (u * > 0.2 ms −1 ), are near zero. This is likely an effect of the low wind speeds in the stable nighttime boundary layer, leading to less turbulence on which the EC technique is dependent. In short, most of the expected drivers for HCHO fluxes (photochemistry, emissions, stomatal uptake and turbulence), though not those for HCHO concentrations, are linked to the solar cycle. However, we saw no evidence at this site of HCHO morning entrainment from overnight oxidative production of HCHO above the canopy, as predicted by Ganzeveld et al. (2008), in either the gradient or flux measurements.

Emission studies
HCHO emission rates from canopy surfaces were measured via branch and soil/litter enclosure experiments. Branch enclosures were performed using a ∼10 L Teflon chamber on a branch located 2 m above the ground. Ultra-zero air enriched with CO 2 to a final concentration of ∼410 ppm v (Scott-Marin) was flowed through the chamber at ∼6 SLM and was sampled using a 1/8" ID PTFE tube. While dry, this air was humidified by tree emission to a typical relative humidity of 20-45 %, comparable to the ambient humidity. Chamber concentration was monitored for ∼4 h. Blank experiments of the chamber without the branch were performed before and after branch sampling. Average HCHO concentration attributed to branch emission was 500±220 ppt v with an average ambient temperature of 22.3 ± 1.0°C. Total dry needle mass was measured to be ∼14.37 g, yielding an average emission rate of 15.4±6.9 ng (g dw) −1 h −1 (dw = dry weight). This is significantly lower than the 500 ng (g dw) −1 h −1 reported by Villanueva-Fierro et al. (2004) for PPine but is within the range of emissions reported for other conifers including Pinus pinea (Kesselmeier et al., 1997) and Picea abies (Cojocariu et al., 2004) (see Sect. 5.2 for discussion). As MEF had a measured specific leaf mass of 120±10 (g dw) m −2 and an LAI of 1.9 m 2 m −2 , our measurement results in an average canopy emission rate of 3.5±1.6 µg m −2 h −1 .
Soil/litter enclosure experiments were performed using a ∼22 l steel chamber, sampling at a flow rate of ∼2.5 SLM using a 1/8 ID PTFE tube. Blank experiments were performed by holding the chamber in the air to measure ambient HCHO levels, then holding the chamber firmly onto areas of ground with either undisturbed litter or soil with the surface area of litter swept away and held until the HCHO concentration equilibrated. One experiment each was performed using seemingly representative areas of ground litter and ground soil. A blank experiment was also performed by placing a clean Teflon sheet on the ground and pressing the chamber into the sheet as it was pressed into the soil, which resulted in no significant difference from the ambi-  ent blanking method. The average HCHO concentration attributed to litter & soil and bare soil were ∼900 ppt v and ∼800 ppt v respectively with an average ambient temperature of 24.0±0.2°C. Based on the ground area covered by the chamber (∼800 cm 2 ), the result is an average ground emission rate of 7.3±1.5 µg m −2 h −1 . One known interference of the HCHO instrument with these measurements is due to the significance of detection axis contamination at low flows (<8 SLM) and changing humidity conditions inside the chamber due to soil/litter moisture. Additionally, closed-chamber soil measurements have been shown to affect pressure gradients in the soil, leading to enhanced CO 2 emission (Kanemasu et al., 1974;Rayment and Jarvis, 1997;Xu et al., 2006), which may similarly affect HCHO emission. Finally, disturbance of the soil/litter and pressure gradients in the chamber itself may have also resulted in increased HCHO emission. Each of these interferences may have resulted in an overestimate of emission rate. Also, due to the heterogeneity of the ground litter, it is quite likely that these two sites do not represent the true soil/litter emission rate but provide simply a semi-qualitative estimate of soil/litter emissions.

Zero-dimensional box model
To quantify different contributions to HCHO flux, we have constructed a zero-dimensional box model to simulate HCHO flux above in the forest canopy similar to those that have been reported in the literature Choi et al., 2010). The concept behind this model is based on the need to maintain mass balance in a box vertically constrained by our HCHO flux measurement. The contribution of vertical transport (flux) to this mass balance is dependent on three processes: horizontal transport of HCHO, sources/sinks of HCHO inside the box, and changes in HCHO concentration inside the box (effectively "storing" source, sink, or transport effects). By accounting for each of these three terms, any remaining HCHO production/loss in this box must correspond to the vertical flux. While gradient data can yield vertically-resolved production/loss information, the goal is to integrate this over the entire box to determine the overall estimated HCHO flux. This mass balance is represented by the continuity equation: P and L are respectively the height-dependent chemical production and loss, E is direct emission, D is deposition, A is advection, and δF HCHO (z)/δz is the flux divergence. As the area surrounding the site was remote and reasonably homogeneous, it was assumed that horizontally-advecting airmasses were similar enough to neglect in this analysis. Solving for F HCHO (h), the modeled flux at height h, yields the following equation: where V Dep is the total deposition velocity of HCHO and the vertical dimension of the box extends from 0 m to h, the EC measurement height (25.1 m). This assumes flux at z = 0 (i.e. ground level) is zero, as soil/litter contributions are treated as direct emission. The end term corresponds to the time rate of change of the HCHO column density, referred to as storage (S). To calculate S, vertically-resolved HCHO concentrations were linearly extrapolated from the gradient data, with the concentration at heights between ground and the bottom inlet assumed to be equal to the bottom inlet concentration. For clarity, we will refer to each of the terms in Eq. (3) as the "flux contribution" for each respective process. The methods used for the determination of these different processes are outlined below. In addition to its simplicity, this model holds many advantages. Fluxes provide a convenient constraint on the vertical mixing at the measurement height, allowing this model to be independent of boundary layer height. Measurements are also available for many heights over the entire measurement volume, removing the need for concentration extrapolation. The primary disadvantage is the absence of higherorder oxidative chemistry, which may lead to significant in-canopy HCHO production from the further oxidation of VOCs formed from the oxidation of BVOCs.

Chemical production
HCHO chemical production is predicted from the first-order oxidation of different VOCs by the following equation: where α i,HCHO is the yield of HCHO and k VOC i ·Ox is the rate constant for the respective VOC and oxidant (Ox). Table S2 shows a full list of modeled reactions with yields and rates (Atkinson and Arey, 2003;Hasson et al., 2004;Atkinson et al., 2006;Lee et al., 2006;Carrasco et al., 2007;Jenkin et al., 2007;Dillon and Crowley, 2008). Isoprene and its oxidation products were neglected in this analysis, due to the low reported concentrations (0.1 to 0.3 ppb v ) of isoprene at this site (Kim et al., 2010) and the short daytime lifetime of HCHO (midday: 1 to 5 h), which likely limits the impact of HCHO advected from upwind production sources. As a result, the total PTR-MS signal at m/z = 69 was considered to be MBO. Monoterpene (MT) speciation was determined by previous observations at this site (Kim et al., 2010), where α-pinene, β-pinene, and 3-carene were found to be 22 %, 26 %, and 21 % of total MT, respectively. The remaining MT (31 %) were assumed to have a reaction rate and HCHO yield equal to the average of the other three. HCHO production from the CH 3 O 2 radical was calculated from methane and peroxyacetyl (PA) radical concentrations, where PA concentrations were calculated using the steady state model presented by LaFranchi et al. (2009) in a method similar to that used by Choi et al. (2010) (see Sect. S3). The oxidants used were OH and ozone. Nighttime oxidation by NO 3 was neglected due to low NO x concentrations at this site. Ozone gradients were available during the measurement period while OH gradients were not. As a result, OH concentration was assumed constant throughout the canopy. This assumption was validated by a series of vertical gradient studies later in the campaign.

Chemical loss
Chemical destruction of HCHO can proceed via reaction with OH or photolysis. Loss due to OH was calculated with the rate constant described by Atkinson et al. (2006) and assuming the OH concentrations were equal at all heights. Typical midday HCHO lifetime with respect to OH was ∼13 h. Photolysis rates for HCHO were determined by weighting the measured downwelling J NO 2 values by the ratio of clear-sky HCHO and NO 2 photolysis rates estimated using the Tropospheric Ultraviolet and Visible (TUV) Radiation Model (http://cprm.acd.ucar.edu/Models/TUV). To account for light extinction in the canopy, the photolysis rates were weighted by the leaf area distribution function (LADF) using a modified Weibull distribution (Teske and Thistle, 2004), for which parameters were determined by destructive harvesting measurements of PPine at a similar PPine forest (Wolfe and Thornton, 2011). The extinction ratio was then calculated by: where SZA is solar zenith angle calculated from the TUV model and k rad = 0.75, an empirical parameter to scale the ground level extinction to be ∼25 % at noon to match the measured photosynthetically active radiation (PAR) profile. Integrating these photolysis rates over the entire canopy, this yielded a typical noon lifetime of HCHO due to photolysis of ∼3.5 h above the canopy and ∼14 h near the ground. Actual loss of HCHO to photolysis was determined by calculating the height-dependent loss using the HCHO gradient, then integrating to calculate the overall HCHO photolysis loss.

Direct emission
Emission flux contributions were extrapolated from the chamber experiments using a simple exponential model (E HCHO = A · exp(βT )), where T is temperature in°C and β = 0.07°C −1 is an empirical constant found for HCHO by Villanueva-Fierro et al. (2004). The PPine emissions were weighted by a factor of (0.85 × P AR/P AR 0 + 15), where PAR 0 is the average, clear-sky, noontime measured PAR, thereby fixing nighttime emissions to 15 % of daytime emissions as observed by (Villanueva-Fierro et al., 2004). The pre-exponential factors (A) determined for both soil and branch emissions from the emission rates found in the experiments (Sect. 3.2) were 1.52 and 0.74 µg m −2 h −1 , respectively.

Dry deposition
Total dry deposition was estimated using a resistance model similar to that used for PAN deposition in previous flux budget studies (Turnipseed et al., 2006;Wolfe et al., 2009). The resistance model calculates the total deposition resistance (R Dep ) as the sum of resistances from separate physical processes (Wesely, 1989;Wesely and Hicks, 2000): R a and R b were calculated using standard literature methods (see Sect. S4) (Monteith, 1965;Wesely, 1989;Hummelshoj, 1995, 1997;Massman, 1998). R c is the surface resistance, or resistance to actual uptake or loss on the leaf, and consists of two parallel terms, stomatal (R ST ) and nonstomatal (R NS ) resistance. As stomatal uptake is negligible at night, R NS was estimated from the nighttime HCHO deposition velocity. At night, the lack of thermal turbulence leads to very small fluxes. Therefore, we can estimate the nighttime HCHO deposition rate by using Eq. (3), setting F HCHO to zero, and solving for deposition: our estimates for nighttime chemical production/loss and emissions. Dividing the values of the non-stomatal deposition flux contribution during the relatively constant nighttime hours (23:00 to 04:00) by the average canopy [HCHO] resulted in an average nighttime deposition velocity of 0.18 ± 0.08 cm s −1 . R NS was then calculated by inverting the following equation: where R a, night and R b,night are R a and R b averaged over the relatively constant nighttime hours. It should be noted that this represents the total non-stomatal deposition velocity, to which both cuticular and soil/ground uptake contribute, but are mathematically inseparable by this method. Similarly, it was necessary to assume that the R b,soil is equal to the calculated R b for a pine needle. Literature values using the boundary layer budget method report HCHO nighttime deposition velocity as ranging from 0.65 to 0.84 cm s −1 Choi et al., 2010). The discrepancy between this work and the literature likely lies in the different assumptions on which either model is based. The boundary layer method assumes similarity between HCHO and ozone deposition and usually depends on literature estimates of ozone deposition. This method also assumes that deposition is the only nighttime loss process and there are no production processes. Finally, the boundary layer method is based on a single measurement and assumes a continuous concentration throughout the boundary layer. The gradient method used in this work makes no assumptions on the HCHO profile, as it is measured directly, and does not depend on literature ozone deposition. The gradient method also estimates nighttime production and loss via the model terms. However, the gradient method still has limitations in that it is much more dependent on direct emission measurements/estimates and assumes the canopy gradient is well represented by the available measurements (in this case, four heights).
R ST was calculated by the following equation (Wesely, 1989).
Mesophyll resistance (R m ) is the resistance to absorption into the plant mesophyll once inside the stomata, which is negligible for HCHO due to its large Henry's law constant (Wesely, 1989;Zhang et al., 2002). R ST,H 2 O was calculated using the Penman-Monteith equation (Monteith, 1965;Monteith and Unsworth, 1990). The resulting average daily minimum R c was ∼180 s m −1 . An alternative method used for estimating R c was the parameterization described by Wesely (1989) for an autumn coniferous forest, which yielded a comparable daily minimum average of ∼226 s m −1 . This latter method was not used in the final model, as the measurement-based method was considered more accurate. The daytime-maximum median V Dep determined by this method was 0.39±0.11 cm s −1 , and had a diurnal profile peaking at 09:00, then gradually decreasing until a sharp decrease at dusk. Similar to the nighttime deposition velocity, this daytime deposition velocity is considerably smaller than the literature value of 1.5 cm s −1 (Krinke and Wahner, 1999). These discrepancies may partly result from the lower LAI and less underbrush at the BEACHON site compared to the literature sites. Additionally, the deposition term is highly dependent on litter emission, which makes it very sensitive to the temperature-dependent method we use to extrapolate litter emission rates. However, the method used in this work is also not dependent on measured ozone deposition velocities, which may be influenced by chemistry as well as deposition (Kurpius and Goldstein, 2003). Direct comparison of deposition using the ozone similarity method described in the literature (e.g. the boundary-layer budget approach)  was not possible for this dataset, as nighttime concentrations did not exhibit clear first-order decay.

Model results and discussion
Modeled fluxes were calculated using data from 13-21 August. The major canopy-integrated HCHO production and loss terms for the base (unaltered) version of the model are shown in Fig. 5, while values for all terms are shown in Table S3. The dominant production terms are direct emission  from both PPine and ground litter, OH oxidation of MBO, CH 4 , and acetaldehyde, and chemical destruction of PA radicals. MT oxidation and ozonolysis in general contribute minimally to the HCHO production. The total production diurnal cycle is similar in form to the radiative diurnal cycle, reflecting the production dependence on temperature and ambient radiation. HCHO loss was dominated by dry deposition, as expected for an in-canopy airmass. The total loss diurnal cycle therefore mostly reflects the diurnal cycle in stomatal uptake. As shown in Fig. 6, the base model underpredicts the noontime HCHO fluxes by a factor of 6 during the day. Modeled nighttime fluxes agree much better with observations, but this is expected as we have constrained nighttime deposition via an assumption of no flux at night.

General sensitivity analyses
We performed a sensitivity analyses on a number of input parameters to determine what model conditions resulted in the best agreement with measurements. Many of these made little or no significant difference in model-measurement agreement. For example, we assumed that OH mixing ratios below the CIMS detection limit were equal to half of the detection limit value, or 2.5 × 10 5 molec cm −3 . To test this, we performed model calculations with OH mixing ratios ranging from zero to the CIMS detection limit (5 × 10 5 molec cm −3 ). This effect was found to be negligible (<5 %) on the order of the missing HCHO flux.
We also tested the effect of separating HCHO PPine emission and stomatal deposition. Strictly, PPine direct emission and stomatal deposition are not independent processes and are related by the HCHO compensation point, the ambient HCHO concentration above which stomatal deposition is expected and below which stomatal emission is expected. As the compensation point can vary by tree species and environment (Seco et al., 2007(Seco et al., , 2008, it was not possible to treat this explicitly at this site. However, as an upper limit to the error this assumption could contribute to the missing flux, we can neglect stomatal deposition by assuming that we are strictly in an emission-only regime. This results in only a ∼10 % reduction in noontime missing HCHO flux. In an extreme case, we can also assume that our measured soil/litter emission rate also represents the sum of both soil/litter emission and deposition. This was simulated by also neglecting the non-stomatal deposition component (therefore also neglected cuticular deposition) and resulted in only a ∼25 % reduction in noontime missing HCHO flux. As a result, this cannot explain the majority of the missing HCHO flux.

Pine emission sensitivity (E350)
In an attempt to explain this missing flux, we scaled the modeled PPine emission rate to reach the best match to the measured flux. We achieved the best match at a PPine emission rate of 350 ng (g dw) −1 h −1 (E350). The diurnal cycle of this case matches the measured flux quite well, though the model we used for PPine emission was directly dependent on temperature and PAR. The emission rate used in E350 is more than an order of magnitude greater than the emission rate predicted by our branch enclosure studies. It is comparable to the 500 ± 400 ng (g dw) −1 h −1 measured by Villanueva-Fierro et al. (2004). However, the formaldehyde rates observed by Villanueva-Fierro et al. (2004) are consistently an order of magnitude higher than those reported for similar tree species by other investigators (Kesselmeier et al., 1997;Cojocariu et al., 2004). The cause of this discrepancy is unclear, but the climate of the area studied by Villanueva-Fierro et al. (2004) was different, and there may have been differences in other factors such as stress conditions. The formaldehyde quantification technique used for all of these studies (collection and storage on DNPH cartridges followed by analysis with high pressure liquid chromatography) has potentially large errors associated with background subtraction and differences due to analytical and enclosure techniques may have contributed to the discrepancy in reported rates. It should also be noted that emission rate is dependent on the β value used in the exponential model, as described in Sect. 4.3. However, with no method of separating these quantities, we continued to use the value found in Villanueva-Fierro et al. (2004) throughout this analysis.

MBO sensitivity (VOC-I)
In another case, we simulated an increase in MBO concentrations, using it as a proxy for a precursor with both a temperature and PAR dependent emission profile. The best match to measured flux was an increase by a factor of 10 (VOC-I). This implies that HCHO production could be significantly impacted by either contributions from higher-order oxidation products of MBO or oxidation of an unmeasured BVOC/combination of BVOCs with a similar temperature/PAR-dependent emission profile. As MBO emission is both a temperature and PAR dependent process, the VOC-I and E350 model cases demonstrate that the HCHO flux corresponds to a temperature/PAR dependent emission profile. However, if these unmeasured BVOCs are assumed to have an OH reactivity similar to MBO, they would contribute 9× the OH reactivity of MBO (median noontime MBO concentration: ∼1.1 ppb v ; median noontime MBO contribution to OH reactivity: ∼1.3 s −1 ), which would be on the order of ∼12 s −1 . As the measured median noontime OH reactivity is on the order of 6 to 7 s −1 in the canopy during the campaign, this suggests that the unmeasured BVOC does not have a similar OH reactivity to MBO. Therefore, in order to form HCHO inside the canopy faster than vertical transport out of the canopy, the primary oxidation pathway of this unmeasured BVOC would need to be through a species other than OH (e.g. ozone).

Monoterpene sensitivity (VOC-II)
As the missing BVOCs thought to cause the OH reactivity gap have been attributed to terpenes, a final sensitivity analysis was simulating an increase in MT concentrations by a factor of 10 (VOC-II). As MT concentrations are highest at night, but oxidation is highest during the day, the result was an increase essentially independent of the time of day. This does not match the observed HCHO flux diurnal cycle, suggesting that measurements of MT are unlikely to be under-predicted. In Fig. 6, values are shown for VOC-II while using the same dry deposition rates as the base case model. When dry deposition was calculated the same as for the other cases, the net effect was an inverse diurnal cycle as HCHO production from MT is greatest at night, a poor match to the measured fluxes. Additionally, the model was no longer able to predict zero flux at night, as the nighttime deposition velocity reached the aerodynamic limit due to a significantly decreased non-stomatal resistance. This further supports that species with a temperature-dependent, PARindependent emission profile, as with MT at this site, are unlikely to be the source of the missing flux.

Conclusions
In this work, we demonstrate the first published measurements using the FILIF technique and the first published measurements of HCHO flux by eddy covariance. The ability to use this emerging class of fiber laser technology now allows for more complex spectroscopic techniques to be used in field conditions, which was previously quite difficult due to the sensitivities of traditional lasers. These advantages allow the FILIF technique to be one of the fastest and most sensitive methods for HCHO detection, with laboratory limits of detection (3σ ) as low as ∼25 ppt v in 1 s. HCHO fluxes were found to have a median diurnal cycle quite similar to that of PAR, with a median midday maximum of ∼80 µg m −2 h −1 (∼24 ppt v m s −1 ). Strong HCHO gradients were observed at night implying deposition. Moderate inverted gradients were observed during the day with higher concentrations near the ground during midday, implying ground litter emission. These gradients were also observed in the canopy during mid/late afternoon, implying PPine emission and/or fast, in-canopy, photochemical production. Branch and soil chamber experiments confirmed HCHO emission of 3.5±1.6 µg m −2 h −1 from PPine and 7.3±1.5 µg m −2 h −1 from soil and ground litter. While typical midday canopy HCHO net production rates are ∼3.2 ppb v h −1 , this corresponds to only 0.079 ppb v h −1 over the entire boundary layer, insignificant with respect to the HCHO budget. Additionally, these measured emissions, along with the gradient profiles, clarify the need to account for not only HCHO emissions from the canopy and undergrowth in forests, but the soil and ground litter as well.
A zero-dimensional box model of the forest canopy using first-order chemical production of HCHO was shown to under-predict HCHO fluxes by a factor of 6. A sensitivity analysis showed that the model would agree with measurements by increasing either the PPine emission rate or the concentration of a species with an emission and reactivity profile similar to MBO. This suggests that the missing HCHO flux is caused by a process that is dependent on temperature and PAR. The disagreement of the measured flux with the VOC-II case, with simulated increased MT, further supports this argument, as the dominant MT at this site have a primarily temperature-dependent emission profile (Kim et al., 2010). Potential explanations for this discrepancy are higher HCHO emission, production by fast, higher-order chemistry of MBO oxidation products, or the processing of unmeasured BVOCs with emission profiles also dependent on tempera-ture and PAR by oxidants other than OH. A model including explicit chemistry of these oxidation products would distinguish between the latter two of these possibilities. The lack of agreement between both non-stomatal HCHO deposition and emission rates between this work and the literature also highlights a need to parameterize HCHO compensation points, emission and deposition rates for trees and soil as functions of temperature, radiation, and humidity.
These measurements provide a constraint on the oxidation in a forest canopy of unmeasured BVOCs, which have been attributed as a cause of the model/measurement mismatch in OH reactivity and concentrations. To conclusively determine this effect, it will be necessary to determine the amount of missing flux that is not due to either higher order chemistry or direct emission. Calculations of OH reactivity compared to measurements have shown that the missing flux cannot solely result from oxidation of missing VOC by OH. Additionally, the minimal emissions of sesquiterpenes at this site (Kim et al., 2010) and the expected OH reactivities suggest that VOC oxidation cannot explain the entire missing flux. As a result, direct emission must be the cause of at least a portion of the missing flux, and this study does not remove the possibility that it may be entirely due to this effect. Future investigations into not only HCHO emission rates from the canopy, but also the soil and ground litter, will be crucial to correctly apportioning HCHO flux.