2.1 System design and operation

A schematic of the system is shown in Fig. 1. Unless otherwise mentioned, all wetted surfaces are stainless steel. At both points of collection, air is drawn in through downward-facing aspirated inlets (R.M. Young, part no. 43408) to reduce thermal fractionation (Blaine et al., 2006). A nominal flow rate of 50 cm³ min⁻¹ is maintained at all times in both lines using pumps with Viton diaphragms and PTFE valves (KNF Neuberger N05-ATI) in concert with mass flow controllers (MKS M100B). Tubing on the tower (${\frac{\mathrm{3}}{\mathrm{8}}}^{\prime \prime }O.D.$ Dekabon^®) is routed into the instrument enclosure without local low spots to minimize buildup of condensed water. Within the building, water is removed from the airstream using a series of three coaxial stainless steel traps. The first is kept at 5 ^∘C with a Peltier cooler (M&C model ECP2000-SS), while the second and third traps are filled with borosilicate glass beads and immersed in a heat-transfer fluid (Syltherm XLT^®) held at −90 ^∘C using a mechanical refrigeration system (FTS VT490). Filters with a nominal pore size of 7 µm (Swagelok SS-6F-7) protect downstream equipment from particulates. Since our analyzers can only measure a single airstream, the high and low intake lines are alternately sent to analysis or waste. The intake selector valve is a pneumatically actuated, crossover ball valve with PTFE packing (Swagelok SS-43YFS2), switched every 720 s.

Figure 1Schematic diagram of the instrumentation installed at Harvard Forest.

Download

Both the O₂ and CO₂ analyzers make differential measurements of the chosen stream of ambient air relative to a cylinder of dried air (the working tank, WT). Data are recorded at 1 Hz. Periodically, the ambient air stream is diverted to waste and the WT is analyzed against one of four calibration tanks: HS, LS, MF, and LF. Details of the various tanks are given in Table 1. The frequency of the calibration runs was determined in part from pre-deployment tests in which we observed instrumental drift when running pairs of (sacrificial) tanks against one another. Other considerations include consumption of standard gases and lost time for atmospheric measurements.

Table 1Nominal composition of calibration tanks. Actual values vary slightly as cylinders are used and refilled. The O₂ values quoted here were measured in the lab of R. Keeling, while the CO₂ values were measured by the CCGG lab of NOAA ESRL GMD. Based on replicate measurements, the O₂ values are typically determined to within roughly ±1.5 per meg on the S2 scale, while the CO₂ values on the WMO X2007 scale are known to about ±0.2 µmol mol⁻¹ (Hall, 2017). The cylinders are Luxfer N265 uncoated aluminum, with 51-14B regulators from Scott Specialty Gases (now Air Liquide). Tanks and regulators are stored horizontally in an insulated box to minimize thermal and gravitational fractionation.

Download Print Version | Download XLSX

The ambient air stream (high or low) or the calibration tank is chosen for analysis with a pair of electrically actuated six-port selector valves (VICI EMT2SF6MWE and EMT2SDMWE). From these valves it passes through another mass flow controller, the NDIR CO₂ analyzer (Li-cor LI-7000), a 0.5 µm filter (Swagelok SS-6F-05), a needle-tipped metering valve (Swagelok SS-4BMG), a changeover valve assembly (four Numatics TM101V24C2 solenoid valves on a custom-built aluminum manifold), and the fuel-cell-based O₂ analyzer (Sable Systems Oxzilla II with customized stainless steel internal plumbing), and then exits through a mass flowmeter (Honeywell AWM3100V). Absolute pressures are recorded upstream of the CO₂ analyzer in three locations (Omega PX303-030A10V) and within the Oxzilla instrument (see below).

The response of the Li-cor CO₂ analyzer depends on the pressure difference between the two cells. Thus, we use a differential manometer and electronically controlled metering valve with active feedback (MKS 223BD, 247D, and 248) to ensure balanced pressures in the Li-cor analyzer.

The response of the Oxzilla O₂ analyzer depends on both pressure and flow differences between the cells. To provide sufficient control of these parameters, we have two metering valves upstream of the Oxzilla and one downstream. We use the metering valves (as well as the active pressure-difference control upstream of the Li-cor) to tune the system so that pressures and flows are as close to equal as possible regardless of the state of the changeover valve (see below). For this tuning, and subsequent monitoring of the system, we use flow data from the Honeywell meters and pressure data recorded within the Oxzilla immediately downstream of the cells. Prior to deployment, we modified the Oxzilla's internal plumbing and electronics so that pressure transducers (Honeywell SCX15N) were teed off of each of the outlets of the fuel cells. The two transducers are read out with the single high-precision digitizer built into the Oxzilla motherboard using a multiplexer to alternately select the transducer signal sent to the digitizer.

To further improve measurement precision of O₂, we follow Stephens et al. (2007) and use a changeover valve to alternate the fuel cells into which the two gas streams flow. This allows for a difference-of-difference calculation that eliminates cell-to-cell biases in the instrument. However, the elimination of bias occurs only if that bias varies little during the two measurement periods being considered. Achieving optimal precision would seem to require high-frequency changeovers to reduce variation in the bias. However, each changeover is followed by a “dead time” when flows and instrument response need to settle. For a given dead time, more rapid changeovers mean fewer measurements are available for averaging, leading to a greater impact of random noise and thus poorer precision.

We adopt a method for quantifying this trade-off and settling on an optimal changeover frequency that was developed by Keeling et al. (2004). We determined dead time for our instrument by analyzing tank air (relative to a working tank) with the changeover running very slowly. We then measured the time after each changeover required for O₂ values to settle to within 1σ of the post-changeover equilibrium value. This dead time was 14 s for our instrument. We then determined the optimal changeover time by analyzing two tanks of air against each other with the changeover valve disabled. We processed these measurements as if the changeover valve were operating, calculating difference-of-difference values for various (hypothetical) changeover rates, imposing a 14 s dead time. This exercise indicated an optimal changeover rate of 24 s between switches.

2.2 Data reduction

2.3 System performance

While we have made great efforts to calibrate our measurements to internationally accepted CO₂ and O₂ scales using our suite of reference tanks (Table 1), we defer a detailed treatment of long-term measurement stability and accuracy for a later publication. Our focus in this paper is on the local biotic exchange ratio, as expressed in a series of short-term averages of CO₂ and O₂ covariation. Thus, our primary concern is instrumental precision and stability over periods of roughly 6 h.

To determine this instrumental precision, we use runs of calibration tanks (see Sect. 2.2.2).

There are two possible approaches: we could use the observed scatter in the three, four, or five difference-of-difference values for O₂ that are calculated during a 6 min run of a calibration tank, and likewise for the 60, 80, or 100 CO₂ measurements retained during the same run. These are the values plotted in Fig. 2. We refer to this metric as “scatter on one tank” (SOOT).

Figure 2The scatter (standard deviation from the mean) of routine CO₂ and O₂ measurements of calibration tanks. Values plotted here apply to individual CO₂ measurements (1 Hz) and single difference-of-difference O₂ values (one value every 48 s).

Download

Not surprisingly, the performance of the instrument varies over time. At times, we could link a loss of precision to problems with an individual analyzer, while on other occasions, we found water (liquid or solid) in the system to be the cause. The source of the major degradation in the precision of both analyzers that develops in late 2012 and persists until the end of our dataset remains under investigation.

An alternative approach, and the one that we adopt here, uses the observed scatter in the difference of paired standard tanks that were run immediately following one another. The calculation of the difference between the HS and LS tanks run in succession emulates the treatment of the atmospheric data by comparing one gas stream (HS in place of atmosphere) to another (LS in place of a “calibration curve”) without the involvement of the working tank. The scatter in the HS–LS difference is a conservative estimate of uncertainty in atmospheric data because it includes any variability that might be introduced by operation of the gas-selection valve, a valve that remains fixed in position during atmospheric sampling.

We find the scatter in the HS–LS differences for 12 d periods is around ±2 ppm equivalent for O₂ and ±0.2 µmol mol⁻¹ for CO₂ until the middle of 2008. From 2012 onward, the scatter was ±9 ppm equivalent and ±0.6 µmol mol⁻¹.

Each HS and LS O₂ value was typically the average from three difference-of-difference pairs, while the CO₂ values were usually an average of 60 measurements. In contrast, a 12 min atmospheric measurement on either the high or low intake is an average of six O₂ difference-of-difference pairs and 120 CO₂ measurements (Sect. 2.2.1). Thus, the atmospheric measurements should have uncertainties given by

Based on these equations, we believe that the atmospheric measurements carry uncertainties of ${\mathit{\sigma }}_{{{\mathrm{O}}_{\mathrm{2}}}_{\mathrm{atm}}}=\mathrm{1}$ ppm equivalent and ${\mathit{\sigma }}_{{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{atm}}}=\mathrm{0.1}$ µmol mol⁻¹ prior to mid-2010. For measurements beginning in 2012, the uncertainties were ${\mathit{\sigma }}_{{{\mathrm{O}}_{\mathrm{2}}}_{\mathrm{atm}}}=\mathrm{4.5}$ ppm equivalent and ${\mathit{\sigma }}_{{{\mathrm{CO}}_{\mathrm{2}}}_{\mathrm{atm}}}=\mathrm{0.3}$ µmol mol⁻¹.

Our choice to work backwards from the observed scatter in HS–LS pairs may not be the optimal technique for extracting uncertainties, but we are confident it yields O₂ and CO₂ errors that bear the correct relative size. Thus, we use the value listed above when calculating slopes in our Deming regressions (Sect. 3).

4.1 The 6 h provenance of atmospheric variability

We use a Lagrangian transport model to estimate the surface regions to which air parcels are exposed during the 6 h prior to their arrival at the EMS tower. Specifically, we created back-trajectories for a range of dates and times, day and night, across seasons and years using the HYSPLIT Lagrangian transport model (NOAA, 2018) forced with NAM12 meteorology. We generated many sets of six 6 h back-trajectories originating hourly during the daytime or nighttime intervals to determine the provenance of the air we analyze at Harvard Forest. An example of one such set of trajectories is shown in Fig. 7.

Figure 7An example of 6 h back-trajectories for Harvard Forest on 2 February 2015. Trajectories start hourly between 22:00 and 04:00 EST, corresponding to the nighttime interval for slope plots. State lines for New Hampshire, Vermont, and Massachusetts are in gray. On this particular night, the source region was fairly constant, so changes in O₂ and CO₂ are relatively likely to be primarily due to local influences.

Download

In general, the 6 h provenance ranges from 50 to 200 km, with little spatial variation over the 6 h interval. Given the location of our study site, a ∼100 km provenance, primarily from the west, means the air we analyze has traveled over a mostly forested landscape, occasionally passing small urban centers. This analysis is consistent with the work of Gerbig et al. (2006), who found a 10-fold (or greater) falloff in the influence of fluxes more than ∼100 km from the study site.

4.2 “Mixing” of slopes

To test the applicability of an end-member approach, we considered the simplest model with two contributors: generic fossil fuel, with an exchange ratio of −1.4, and biotic activity, with an exchange ratio of −1.05. These characterized the stoichiometry of surface fluxes used in a one-dimensional box model, in which parcels of air travel forward in time in 12 min steps, passing over a fictitious landscape with forested, suburban, and urban components. In 30 of these time steps the parcels exchange CO₂ and O₂ with the underlying sources and sinks, changing the composition of the parcel according to

where C is the mixing ratio of O₂ or CO₂ in micromoles per mole, F is the flux from the landscape underneath the box in micromoles per square meter per second, t is the time step in seconds, h is the height of the box in meters, and the factors of 22.4 (the molar volume) and 1000 convert moles to liters and liters to cubic meters, respectively. For h we use one-half of the planetary boundary layer (PBL) height. We spin up the model for 24 h and then record O₂ and CO₂ concentrations of consecutive parcels as they “arrive” in their final time steps over the same 6 h day and night intervals used for the observations.

The CO₂ surface fluxes vary with time of day and day of year and are based on observations at Harvard Forest (Munger and Wofsy, 2017) and in urban (Bergeron and Strachan, 2011; Velasco and Roth, 2010; Ward et al., 2015) and suburban (Bergeron and Strachan, 2011) settings. The O₂ surface fluxes were derived from the CO₂ fluxes by assuming exchange ratios of 1.4 for fossil fuel combustion and 1.05 for biospheric activity. The PBL values are climatological hourly average values, calculated separately for summer and winter from the NAM12 dataset (NOAA, 2018). Within the STILT model (Gerbig et al., 2006), only air within the lower half of the PBL is influenced directly by surface fluxes. Air above this height is closely matched to the composition of the free troposphere. Because we use the STILT model to infer the provenance of air we analyze, we define our box height as half the height of the PBL so that the two models (STILT and our 1-D box model) are conceptually consistent. We performed sensitivity studies, varying the landscape composition (ranging from 100 % forest to 100 % urban), surface flux magnitude (values characteristic of winter and summer seasons), and PBL height (25 %–75 % of the NAM12 values). In all cases, the CO₂ and O₂ values of the parcels entering the model were held constant, ensuring that any variability predicted by the model was a result of the fluxes from the landscape over which the parcels traveled in the previous 6 h and/or changes in PBL height.

We emphasize that the purpose of the 1-D model is to explore the validity of the end-member mixing approach. In particular, we ask whether a mix of sources will relate simply to a mix of slopes, despite the fluxes and the box height having their own diel cycles.

In all but a few special circumstances, the model predicts slopes that are completely consistent with the specified admixture of the source fluxes over which the air parcels have traveled. Exceptions arise when the covariation in the boundary layer height and diel cycles in the fluxes conspire to bias slopes by a few percent. Results of various sensitivity tests and a discussion of these exceptions can be found in Table 5.1 of Conley (2018), along with a more complete description of the model and its forcing fluxes.

The 1-D model results give us confidence that slopes close to −1.0 are a closer representation of the LBER, while more negative slopes have ever larger contributions from fossil fuel combustion. Though reassuring, we acknowledge the model results are not definitive due to the many simplifications involved. More sophisticated studies are needed to give a truly robust confirmation of the end-member mixing framework.

Aside from the slopes, the full range of O₂ and CO₂ variability predicted by the model is in the general range of our observations. Because of the crude nature of our model, this result is somewhat surprising. Our values for the box height and the molar volume are almost certainly wrong, so the data–model agreement may be coincidental.

5.1 Information from subsets of data

Considering particular subsets of the data (Table 2) may provide insight into the processes determining observed slope values. Thus, we categorize our 6 h slopes by day and night, high and low, winter (December, January, February) and summer (June, July, August).

We begin with data aggregated in broad categories. Each group's average slope, calculated using the 3σ cut described above, is shown in Fig. 8. These groupings convey some unequivocal messages: summer slopes are less negative than winter slopes. Day slopes are less negative than night slopes. Low-intake slopes are less negative than high-intake slopes. Daytime data from the low intake are least negative of all.

Figure 8Average slopes for various subsets of the data. Error bars show the standard error on the mean. Subsets within each group (i.e., sharing a color) are nonoverlapping.

Download

We also consider narrower categories, each with a unique set of slopes. The average slopes for each category are shown in Fig. 9. For these averages, we impose a 2.5σ iterative cut because a 3σ cut fails to reject several unphysical fits (with slopes >2.0) in the winter subsets, badly biasing the mean.

Figure 9Average slopes for various nonoverlapping subsets of the data. Error bars show the standard error on the mean. Winter data are shown in blue, summer in red. Categories are ordered by increasingly negative mean values for convenience.

Download

The most striking feature of this plot is the large scatter in the winter data. Even with the 2.5σ iterative cut, the slopes within a category range widely and are difficult to interpret. In contrast, the summer data comprise much tighter groups and have well-defined mean values, showing the same day–night and high–low relationships as the full dataset.

Explaining these patterns is challenging. Possible factors determining the observed slope include the magnitudes of the fluxes, the dynamics of the atmosphere, and variation in the biospheric exchange ratio over time.

5.2 Temporal variability in slopes

Our data are sufficiently abundant that we can examine the values of slopes as a function of time on both interannual and seasonal timescales. Figure 10 shows that one year looks much like the next, but stacked years reveal a seasonal variation in the slopes.

Figure 10(a) Slopes (averaged weekly) for a climatological year. (b) Slopes (averaged monthly) for the entire period of study. The gray line in (a) is a locally weighted least-squares fit with 46 % smoothing. It is simply intended as an objective guide to the eye.

Download

The seasonal cycle in slopes may be due to seasonal changes in LBER itself. Trees synthesize many different compounds throughout the year (Seibt et al., 2004), each with their own oxidative ratios (Bloom et al., 1989). However, direct measurements of leaves and tree rings indicate oxidative ratios in these tissues appear invariant through the year (Gallagher et al., 2017).

Conversely, the total amount of biological activity unquestionably varies seasonally, while the fossil fuel contribution remains roughly constant. As mentioned above, earlier work at Harvard Forest using carbon monoxide and acetylene (Potosnak et al., 1999) shows clearly that fossil fuel combustion is a driver of CO₂ variability in the winter, but is overwhelmed by biotic signals during the summer. Thus, the simple seasonal shift in the biotic–fossil fuel balance is a likely explanation for some, if not all, of the seasonal cycle in slopes.

Spectral analysis (using Lomb–Scargle periodograms to accommodate our irregularly spaced data; Press et al., 2007) of various subsets of the data shows some evidence of the seasonal cycle mentioned above. It appears to be mostly in the night–high-intake subset. There is also a suggestion of multi-year variation (0.6 cycles per year) in the day–low dataset and 3-week periodicity in the day–high and night–low data. Identical spectral analyses of environmental variables measured at Harvard Forest (temperature, net ecosystem exchange, ecosystem respiration, photosynthetically active radiation, and u^*) all show highly significant seasonal and diel cycles, as expected.

The 3-week cycle in the day–high and night–low records is perplexing. No such cycle appears in the environmental variables listed above. Thus, our working hypothesis is that the 3-week cycles are an experimental artifact resulting from the timing of runs of calibration tanks and swapping of working tanks.

5.3 Selecting periods of heightened biogenic signal

Given the competition between biogenic and anthropogenic influences in our observations, we searched for periods when the biogenic contribution might be most dominant. We selected times when the fluxes in the immediate neighborhood were likely to dominate. Our assumption was that background variation (i.e., compositional variability due to influences beyond the ∼100 km provenance of our measurements) is likely to carry signals of continental-scale gradients with more substantial fossil fuel influences embedded therein.

One approach we considered was based on u^*. We hypothesized that slopes calculated when u^* was low were more likely to be biogenic since the air at our intakes was exchanging minimally with the overlying boundary layer. However, a scatter plot of slopes as a function of u^* showed no correlation. Similarly, slope values were insensitive to discrete cuts on average u^* (<40 cm s⁻¹), maximum u^* (<50 cm s⁻¹), and the scatter in u^* (${\mathit{\sigma }}_{u^{*}}<\mathrm{10}$ cm s⁻¹).

We also explored the possibility of a link between slopes and the ranges of CO₂ and O₂ during the 6 h periods. Again, no correlation was apparent, whether with discrete cuts during summer and winter periods separately, or with scatter plots of slope vs. range.

We tested the hypothesis that summer periods of extended high barometric pressure with languid anticyclonic circulation might be dominated by local fluxes. Again, we observed no correlation between slope and barometric pressure, whether with discrete cuts or with scatter plots.

Finally, we looked for correlations between slopes and wind speed, direction, and variability. We found none, even when winds were from the relatively polluted sector to the southwest.

In short, relatively simple indicators of reduced advective transport show no correlation with shallower slopes. This insensitivity of our result adds support for our model of a geographically modest footprint comprised of a landscape in which biogenic fluxes are much stronger than anthropogenic ones.