Vertical partitioning of CO2 production within a temperate forest soil

The major driving factors of soil CO2 production – substrate supply, temperature, and water content – vary vertically within the soil profile, with the greatest temporal variations of these factors usually near the soil surface. Several studies have demonstrated that wetting and drying of the organic horizon contributes to temporal variation in summertime soil CO2 efflux in forests, but this contribution is difficult to quantify. The objectives of this study were to partition CO2 production vertically in a mixed hardwood stand of the Harvard Forest, Massachusetts, USA, and then to use that partitioning to evaluate how the relative contributions of CO2 production by genetic soil horizon vary seasonally and interannually. We measured surface CO2 efflux and vertical soil profiles of CO2 concentration, temperature, water content, and soil physical characteristics. These data were applied to a model of effective diffusivity to estimate CO2 flux at the top of each genetic soil horizon and the production within each horizon. A sensitivity analysis revealed sources of uncertainty when applying a diffusivity model to a rocky soil with large spatial heterogeneity, especially estimates of bulk density and volumetric water content and matching measurements of profiles and surface fluxes. We conservatively estimate that the O horizon contributed 40–48% of the total annual soil CO2 efflux. Although the temperature sensitivity of CO2 production varied across soil horizons, the partitioning of CO2 production by horizon did not improve the overall prediction of surface CO2 effluxes based on temperature functions. However, vertical partitioning revealed that water content covaried with CO2 production only in the O horizon. Large interannual variations in estimates of O horizon CO2 production indicate that this layer could be an important transient interannual source or sink of ecosystem C.

Most soils also experience larger variations in wetness near the surface than at depth. Inputs of carbon are also usually largest near the soil surface, although deep root inputs can be important (Nepstad et al., 1994;Jobbágy & Jackson, 2000). Hence, the major driving factors of soil CO 2 production -substrate supply, temperature, and water content -all vary vertically within the soil profile, with the greatest temporal variability usually near the soil surface.
Carbon in the organic horizons at the surface of forest soils has mean residence times ranging from a less than a year to a few decades, depending on litter quality and climate (Trumbore 2000). Significant inputs of C to the litter layer mean that decomposition in the organic horizon probably contributes significantly to the total CO 2 efflux from the soil, and the relatively short mean residence times of C in the organic horizon mean that this contribution is potentially subject to rapid change. Experiments using irrigation (Borken et al., 2003;Lee et al., 2004) and throughfall exclusion (Borken et al., 2006) have demonstrated that wetting and drying of the organic horizon accounts for a large fraction of temporal variation in summertime soil CO 2 efflux. Although the variable contribution of the organic horizon to soil CO 2 efflux is clear, it is difficult to quantify. The objectives of this study were to partition CO 2 production vertically in a mixed hardwood forest of New England, USA, and then to use that partitioning to evaluate how the relative contributions of CO 2 production by genetic soil horizon vary seasonally and interannually.

Site description
The study was conducted in a mixed deciduous forest at the Prospect Hill tract of the Harvard Forest (42132 0 N, 72111 0 W) at 340 m elevation, in Petersham, MA, USA. The present forest developed after a hurricane in 1938 and is dominated by red maple (Acer rubrum L.) and red oak (Quercus rubra L.). Black birch (Betula lenta), white pine (Pinus strobus) and hemlock (Tsuga canadensis) are also present. Before the hurricane, a similar mixed deciduous forest had developed on pastures abandoned since the late 19th century (Foster, 1992). The O horizon is 3-8 cm thick, stores about 7 kg m À2 of organic matter, and consists of Oi, Oe, and Oa horizons (Borken et al., 2006). As a result of 19th century agricultural use, the upper mineral soil is partly disturbed as indicated by varying depths of the A horizon from 2 cm down to 15 cm. The thickness of the total mineral soil varies from 40 to 65 cm and has rock contents up to 40% by volume. The soils are well-drained, fine sandy loams on glacial till, classified as the Canton series, Typic Dystrochrepts. The soils are acidic (pH of A horizon in water: 4.3; Compton & Boone, 2000). Mean annual air temperature is 8.5 1C and mean annual precipitation is 1050 mm, with precipitation evenly distributed throughout the year. Snow intermittently covers the soil from December to April. Potential evapotranspiration rates often exceed precipitation from June to August, resulting in decreasing soil moisture (Borken et al., 2006).

Soil respiration measurements
We have been measuring soil CO 2 efflux since 1995 at several landscape positions (Davidson et al., 1998;Savage & Davidson, 2001) within the footprint of an eddy covariance tower (Barford et al., 2001) at the Prospect Hill tract of the Harvard Forest. In this study, we use data from the 20 m Â 20 m study plot located near the base of the tower, where we also have instrumented soil pits. Six PVC rings (25 cm diameter) were installed randomly at this site in 1995. Flux measurements have been made on a weekly basis during the summer, biweekly during spring and autumn, and as conditions permit during the winter. When snow is present, we place a chamber top directly over the snow surface rather than over a ring. However, when snow cover is discontinuous we avoid measurements because of the extreme heterogeneity of the surface exchange processes.
Details of flux measurements are provided by . Briefly, a vented cover is placed over the ring and air is circulated at a rate of 0.5 L min À1 between the chamber and an infrared gas analyzer (IRGA) for a period of about 5 min. The IRGA was calibrated daily using a certified standard. The flux is calculated from the slope of the increase in CO 2 concentration of the chamber, adjusted for temperature and pressure. Davidson et al. (2002b) have evaluated the potential sources of uncertainty, error, and bias in this chamber-based system. This manual system was also in good agreement with an automated system .

Quantitative soil pits
Two quantitative soil pits were excavated at the tower site (Gaudinski et al., 2000). A grid with 10 cm Â 10 cm increments was established over the surface of each 1 m Â 1 m soil pit. As each genetic horizon was removed, the depth to the horizon below was measured at each of the grid points. For each horizon, all coarse roots and coarse rocks were collected and weighed separately. All soil material was also collected and weighed at field moist conditions. Subsamples of each horizon were retained for drying, weighing, sieving for gravel content (42 mm), and total C analysis using a CN analyzer. Bulk density (BD) of each horizon was determined by dividing the sum of the mass of soil and rocks by the volume of the horizon, based on the average depth of the horizon. Volumetric coarse rock fraction (RF) of each horizon was determined by converting the mass of the coarse rocks to volume, using an assumed particle density of 2.65 g cm À3 , and then dividing that value by the volume of the horizon.
Temperature, water, and CO 2 profiles While each soil pit was open, time domain reflectometry (TDR) probes, type T thermocouple wire for measuring temperature, and 3 mm diameter stainless steel tubes for gas sampling were inserted horizontally into the pit walls at selected depths, as described by Davidson & Trumbore (1995). In general, these probes were installed at the interfaces of each horizon and at the middle of the A and B w1 horizons. After installation, the wires and tubes were bent upwards to bring them to the surface, and the pits were backfilled. The TDR soil moisture probes (Campbell Scientific, Logan, UT, USA) were interfaced with a cable tester and a datalogger to record waveforms on an hourly basis. The dielectric constant of the soil was then converted to volumetric soil moisture based on our own laboratory calibrations using intact soil cores. Temperatures were also logged hourly.
Soil CO 2 profiles were measured on the same days as soil respiration measurements. Each tube was flushed by withdrawing and discarding 15 cm 3 of soil air, using a syringe and needle inserted into a fitting and septum at the exposed end of the tube. Duplicate samples of 10 cm 3 soil air were then withdrawn from each tube and returned to the nearby field laboratory. These samples were injected into a sample loop and then into a stream of CO 2 -free air to an IRGA, as described by Davidson & Trumbore (1995). The CO 2 concentrations were proportional to peak areas recorded by the IRGA for each injection, as determined by a standard curve using injections of known concentrations. Good agreement between consecutively sampled duplicates indicated that the volume of gas removed from the tubes did not significantly affect the CO 2 concentration at the end of the sampling tube within the soil.

Modeling effective diffusivity and CO 2 production by depth
We used the approach of modeling CO 2 production as a function of depth described by Davidson & Trumbore (1995) and Gaudinski et al. (2000), except that the more recent diffusivity model of Moldrup et al. (2000) was used as outlined here. Total porosity (cm 3 cm À3 ) was calculated from measures of BD: where BD is the measured bulk density and PD is the weighted average particle density of minerals (assumed particle density of 2.65 g cm À3 ) and organic matter (assumed particle density of 1.4 g cm À3 ): where OM is the percent organic matter, assuming that it is 58% carbon. The ratio of diffusivity in soil (D P ) to diffusivity in air (D O 5 0.139 cm 2 s À1 for CO 2 at 273 K and 1 atm) was calculated following the model of Moldrup et al. (2000): where T is the soil air temperature in kelvin, e is the ambient air-filled porosity (total porosity -volumetric water content from TDR measurements; cm 3 cm À3 ), e 100 is the air-filled porosity at À100 cm H 2 O tension (i.e. the À10 kPa treatment of an intact soil core on a pressure plate), and b is the slope of the line determined from the water retention curve: where C is the water potential and y is the volumetric water content of the pressure plate measurements. The water retention curves are from Savage & Davidson (2001) and were derived from four intact soil cores sampled at the same tower site. The curve for duplicate cores from the 0 to 10 cm depth of the mineral soil were used for the A horizon and the curve for duplicate cores from the 10-20 cm depth increment were used for B and C horizons. The e 100 term represents macroporosity, whereas the b term characterizes the pore-size distribution (Moldrop et al., 2000). This diffusivity model, like the others commonly used in the literature (e.g. Davidson & Trumbore, 1995;Rolston & Moldrup, 2002;Suwa et al., 2004), is based on studies of diffusion through porous media. In the glacial till soils of New England, the coarse RF can exceed 20%, which requires consideration of how the diffusivity models should be modified to accommodate a significant, heterogeneous, nonporous component of the soil. Moreover, the intact soil cores used for the water retention curves included gravel, but they obviously did not include large rocks. Therefore, we modified the e 100 term to adjust for the bias in the estimates of volumetric air and water contents at field capacity in cores that did not contain the average composition of coarse RF. The fraction of the volume occupied by water in small intact cores at field capacity is larger than the fraction of total soil bed volume occupied by the water at field capacity in the undisturbed soil that also contains rocks (Fig. 1). Hence, the water content at field capacity was adjusted downwards as a function of RF: where y 100 is the volumetric water content at À100 cm H 2 O tension in intact cores and RF is the volumetric RF in the quantitative soil pits. Profiles of CO 2 concentrations were fitted to an exponential model ( Fig. 2): where C z is the CO 2 concentration (percent by volume) at depth z (cm), C 1 is the fitted asymptote of CO 2 concentration at infinite depth, 'a' is a fitted parameter that characterizes the sharpness of the curve, and 0.04 is the approximate concentration of CO 2 at the soil surface. Most values for C 1 fell between 0.2 and 1.0 and most values for 'a' fell between 0.01 and 0.10. Fluxes of CO 2 were estimated at depth corresponding to the interface between two genetic horizons by applying Fick's first law: where the flux is expressed as g C m À2 h À1 , 52 700 is a units conversion factor, T is the soil air temperature in kelvin, D P is the effective diffusivity from Eqn (3), and dC/dz is the concentration gradient at depth z estimated from the first derivative of the exponential fit for the CO 2 concentration profile: Production of CO 2 in the C horizon and below was estimated as the flux of CO 2 at the depth corresponding to the top of the C horizon. Production in the B horizon was estimated as the difference between the flux at the top of the B horizon and at the top of the C horizon. Production in the A horizon was estimated as the difference between the flux at the top of the A horizon and at the top of the B horizon. Production in the O horizon was estimated as the difference between the mean flux measured by soil chambers and the flux estimated at the top of the A horizon. We did not attempt to estimate production within the O horizon directly from estimates of O horizon concentration gradients and diffusivity because both CO 2 concentrations and diffusivities are extremely difficult to measure well in the highly porous O horizon (Hirsch et al., 2004). The six chambers located at the tower site were divided into two groups of three, with the mean of the three chambers closest to each soil pit paired with the profile estimates of that pit. Efflux of CO 2 and production in each horizon were estimated for each sampling date, and these estimates were then linearly interpolated between sampling dates and were summed to obtain annual estimates.

Diffusivity estimates using radon
As a check of the modeled estimates of diffusivity, radon profiles and surfaces fluxes were measured under relatively dry conditions in the late summer of 1997 and during relatively wet conditions during the spring of 2002. Radon measurements were made by direct alpha-counting of air samples as described in detail by Davidson & Trumbore (1995). In the case of soil profiles, a 60 mL sample was withdrawn by syringe from the same stainless steel tubes used for the CO 2 concentration profiles. Each syringe sample was passed through anhydrous calcium sulfate to remove water vapor and then injected into a preevacuated Lucas counting cell. Ambient air was used to fill the cell to atmospheric pressure. In the case of surface fluxes, the same soil chambers were used as in the CO 2 flux measurements, except that they were used as static chambers, rather than circulating the air. At each time interval, an 80 mL syringe sample was withdrawn from the chamber for counting in the Lucas cells. When the soils were relatively dry, the radon flux was sufficiently strong that a flux could be calculated from four samples

Soil bed space sampled by intact cores = 79%
Water content at field capacity in cores = 0.32/(0.10 + 0.32 + 0.37) = 0.405 Fig. 1 Illustration of adjustment needed for applying the water content at field capacity (y 100 ), determined in intact soil cores, to the field situation. In this example for the B horizon at our tower study site, the coarse rock fraction was 21% and the water content at -100 cm H 2 O tension in the intact cores was 0.405 cm 3 cm À3 . collected at about 20 min intervals. When the soil was wet, the intervals between sampling had to be extended to 50 or even 100 min. Because the chambers had to be kept in place over the soil for so long, which alters the concentration gradient between the soil and the overlying chamber headspace (Davidson et al., 2002b), the increase in radon concentrations within the chamber headspace was often nonlinear. In other words, the rate of increase in radon concentration decreased over time as the headspace concentration approached an asymptote. To correct for this effect, a second-order polynomial function was fit to the chamber concentration data, and the instantaneous slope of the initial increase in radon concentration at time zero was estimated from the first derivative of the polynomial function. This initial slope was combined with the chamber height to calculate a radon flux as Bq m À2 min À1 . Unlike CO 2 concentration profiles, radon concentration profiles tend to increase with depth throughout the measured profile (Fig. 2). This difference in shapes of the concentration profiles indicates that CO 2 production declined with depth as the concentration approached an asymptote, whereas there must be a strong source of radon deep in this soil, below the layer of biological activity. Radon concentrations presumably saturate at much greater depths that we were able to dig in this glacial till soil. The nearly linear profiles of radon within the depths that we sampled do not allow us to estimate the depth where radon concentrations saturate or to use an exponential fit of the radon profile to estimate diffusivity, as done by Dö rr and Mü nnich (1990). Instead, by rearranging Fick's first law (Eqn (7)), the average effective diffusivity of the profile can be estimated from the surface radon flux (F) divided by the slope of a linear fit (dC/dz) of the radon concentration gradient This estimate of D P for radon was then multiplied by 1.2 to compare with the modeled D P of CO 2 to account for the 20% higher diffusivity of CO 2 in air (D O ) compared with radon in air.

Contributions of the O horizon and mineral soil horizons
Concentrations of CO 2 within the soil generally followed the same seasonal trend as surface efflux, increasing from winter to summer (Fig. 3). The concentrations at equivalent depths within each duplicate pit at the tower site were in good agreement until 2001, when the concentrations began to diverge below 25 cm depth (Fig. 3). We observed no change in soil wetness or diffusivity that would account for accumulation of CO 2 in the deep soil of pit #2, so we speculate that roots grew near the sampling tubes in the deep soil horizons of pit #2 beginning in 2001, thus contributing to larger CO 2 concentrations there. Total CO 2 efflux varied by as much as a factor of two across years (Fig. 4, Tables 1 and 2). The O horizon contributed the largest fraction of CO 2 production compared with other horizons (Fig. 4). An unusually severe summer drought (50% of the 30-year average precipitation during the months of June-August) suppressed efflux in 1999. Estimates of CO 2 production in the O horizon became negative during 1999, indicating failure of this modeling approach under those severely dry conditions. Potential sources of error, particularly for the 1999 data, are discussed below in the discussion section on sensitivity analysis. The 1999 data have been left out of subsequent regression analyses. Based on mean annual sums, production of CO 2 in the O horizon accounted for 40-48% of the total efflux (Tables 1 and 2).

Effects of temporal variation of temperature and water content
Production estimates within each horizon were correlated to temperatures measured within each horizon, except in the C horizon, where the production signal was too weak to permit statistically significant correlation with temperature (Fig. 5). The temperature sensitivity, expressed as Q 10 values, increased from the O horizon (Q 10 5 2.3; 95% confidence interval: 1.7-3.1) to the B w1 horizon (Q 10 5 4.4; 95% confidence interval: 3.4-5.7). The sums of the production predicted in each horizon based on these temperature functions accounted for 61% of the variance in measured surface CO 2 efflux. A temperature function correlating surface CO 2 efflux and temperature at only 10 cm depth (efflux 5 0.090e (0.987(Temp-10)/10) ; Q 10 5 2.7; 95% confidence interval: 2.3-3.1) also accounted for 61% of the variation of CO 2 efflux and predicted nearly identical fluxes as the sum of temperature functions across horizons. Hence, it appears that no predictive capacity of surface CO 2 efflux based on soil temperatures was gained by partitioning production by horizon as a function of temperature in All values are in units of g C m À2 yr À1 . Means, standard deviations, and coefficients of variation (CV) are given across years. All values are in units of g C m À2 yr À1 . Means, standard deviations, and coefficients of variation (CV) are given across years. each respective horizon. Soil temperatures at depth exhibited less seasonal amplitude and a lag of about 10-20 days in minima and maxima compared with temperatures near the soil surface (Fig. 6). However, the combination of only modest temperature lags among depths and the relatively small contributions of CO 2 production in B and C horizons (Fig. 5) resulted in little effect of partitioning the temperature function by depth on predictions of surface fluxes. Soil water content was significantly correlated with soil CO 2 efflux and O horizon production during the months of July-September, but was not correlated with the flux calculated at the top of the A horizon or for production within the A horizon (Fig. 7). Although water content is an input to the estimate of diffusivity within the A horizon (Eqn (3)), other factors, such as the CO 2 concentration gradient, apparently cancelled this effect. Production of CO 2 within the A horizon did not appear to respond significantly to variation in soil water content, indicating that most of the temporal variation of CO 2 production because of wetting and drying occurred in the O horizon. Adding a linear function of volumetric water content at 6 cm depth (the TDR probe at the O/A interface) to the empirical temperature model (CO 2 efflux 5 0.308e (1.149(Temp-10)/10) 1.000y 6 cm ) increased the R 2 value from 0.61 to 0.69, which is consistent with previous regression analyses of temporal variation of soil respiration at this study site (Savage & Davidson, 2001). Similarly, adding a linear water function to the temperature model of O horizon CO 2 production (O horizon production 5 0.547e (1.002(Temp-10)/10) 0.285y 6 cm ) also improved the R 2 value from 0.22 to 0.42. A soil moisture power function (O horizon production 5 0.553e (1.147(Temp-10)/10) y 6 cm 2.087 ) further increases the R 2 value to 0.49. The increase in the exponential parameter for apparent temperature sensitivity from 0.833 to 1.147 implies an increase in the apparent Q 10 value from 2.3 to 3.2 when the moisture power function term is added.

Comparison of modeled diffusivity of CO 2 and measured diffusivity of radon
Modeled diffusivity was well correlated with radonbased measured diffusivity, but with an offset and with a slope of less than unity (Fig. 8). This is true regardless of whether the comparison is made with modeled diffusivity of the A horizon only or with an average of all horizons, although the slope is somewhat lower and the Y-intercept somewhat higher for the regression with the modeled A horizon diffusivities. The radon-based diffusivity estimate is an integrated average of the entire profile. On average, the modeled diffusivities are about 36% higher than the diffusivities based on radon measurements. Radon-based estimates could be subject to errors in radon surface flux measurements Day of year Fig. 6 Seasonal variation in soil temperature at five depths from hourly measurements averaged across 8 years (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003). and in measuring and fitting the radon concentration gradients. Assuming that the radon-based estimates are correct and that the modeled diffusivities in the mineral soil are too high, then the model would overestimate CO 2 flux at the top of the A horizon and thus underestimate production of CO 2 in the O horizon. If modeled diffusivity in the A horizon is overestimated by 36%, then the O horizon production would constitute 54-65% of the total CO 2 surface efflux. Hence, the model estimates of O horizon production reported above, at 40-48% of surface CO 2 efflux, may be conservative.

Sensitivity analysis
In addition to comparing modeled and measured diffusivities within the soil profile, we conducted a sensitivity analysis of the factors that affect modeled diffusivity in the mineral soil, and the resulting estimate of CO 2 production in the O horizon. This modeling approach requires the input of four temporal variables (Fig. 9a) and six fixed parameters (Fig. 9b). Using pit #1 of the tower site as a test case, we systematically varied each parameter or variable by a factor of 0.5, 0.75, 0.9, 1.0, 1.1, 1.25, and 1.5, while keeping all others constant. The mean annual production of the O horizon was used as the output variable to test the sensitivity of variation in each parameter and variable. The estimate of O horizon production was most sensitive to variation in mineral soil BD and TDR measurements of volumetric water content (Fig. 9), both of which directly affect the estimate of air-filled porosity (e) in Eqn (3). The sensi-tivity was asymmetric, because of the nonlinearity of the diffusivity function (Eqn (3)). The estimates of diffusivity in the mineral soil and, hence, production of CO 2 in the O horizon, are especially sensitive to water content at the low end of the range of observed water contents (Fig. 9a).
The O horizon production was also linearly sensitive to estimates of C 1 and a (Fig. 9a), which are derived from curve fitting of the CO 2 concentration profiles (Eqn (6)). A 10% variation in estimates of C 1 and a each caused about a 12% change in the estimated O horizon production. This sensitivity analysis indicates that O horizon production was less sensitive to all other parameters and variables (Fig. 9), although hidden effects of RF estimates are also discussed below.

Importance of the organic horizon of the forest soil
The estimate of 300-400 g C m À2 yr À1 mean annual production in the O horizon is about 40-48% of mean annual soil CO 2 efflux. This estimate is slightly higher than that of Bowden et al. (1993), who used a litter removal experiment to estimate that the litter layer contributed 37% of total soil CO 2 efflux at a different study area within the Harvard Forest property.  Fig. 9 Analysis of sensitivity of annual O horizon CO 2 production estimates to variation in each input variable for temporally varying factors (upper panel: volumetric water content, C 1 and a from Eqns (6) and (8), and soil temperature) and fixed parameters (lower panel: y 100 and RF from Eqn (5), BD from Eqn (1), carbon content from Eqn (2), horizon depth from Eqns (6) and (8), and the fitted 'b' term of the water retention curve from Eqn (4)).
Buchmann (2000) also used O horizon removal experiments to estimate that 30-40% of the soil CO 2 efflux was generated in the O horizon of several spruce forests of Germany. Based on lower diffusivities in the mineral soil measured by the radon method compared with the model in the present study, our estimate of the importance of the O horizon is likely to be conservative, and its contribution could be as high as 50-60% of total soil CO 2 production. Our annual estimates of soil respiration are within the range of 60-80% of total ecosystem respiration measured by eddy covariance (Goulden et al., 1996;Barford et al., 2001). Our O horizon production estimate is about twice the rate of mean annual litterfall C (Davidson et al., 2002a). Leaching of DOC from the O horizon to the mineral soils, estimated at 22 g C m À2 yr À1 (Currie et al., 1996), is a relatively small component of the O horizon C budget. Therefore, assuming that the O horizon is at steady state over the long term, no more than half of the O horizon CO 2 production, on average, can come from decomposition of organic matter. Root and mycorrhizal respiration and decomposition of root tissue, root exudates, and mycorrhizal C must comprise the other half (or more) of the CO 2 production in the O horizon.
O horizon production of CO 2 was more variable among years (Tables 1 and 2) than was litterfall (170-270 g C m À2 yr À1 ; Davidson et al., 2002a), indicating that the O horizon could serve as a significant transient source or sink of C on an interannual basis. The variation in the estimate of CO 2 production in the O horizon is larger than interannual variation in net ecosystem exchange of C at the tower (Barford et al., 2001). A transient sink in the O horizon because of incomplete decomposition of C inputs would be most likely in years with below average summertime precipitation. Reduced summertime soil respiration because of drought is a common observation in the literature (Davidson et al., 1998;Savage & Davidson, 2001;Xu & Qi, 2001;Borken et al., 2002Borken et al., , 2006Irvine & Law, 2002;Rey et al., 2002;Hanson et al., 2003;Janssens & Pilegaard, 2003;Subke et al., 2003;Curiel Yuste et al., 2004), but the soil horizons most affected have not been clear, and the evidence has been conflicting at times (cf. Buchmann, 2000;Subke et al., 2003). A transient O horizon C sink because of incomplete decomposition of leaf and root litter inputs is consistent with eddy covariance estimates of relatively high net ecosystem exchange of C during summer drought years at the Harvard Forest (Barford et al., 2001) and lower ecosystem respiration during the dry season in a tropical forest (Saleska et al., 2003). A 5-year throughfall displacement experiment caused an increase of 220 g C m À2 in the litter layer stock of C in a temperate forest of Tennessee, USA (Johnson et al., 2002).
Higher Q 10 values for CO 2 production in the mineral soil than in the O horizon might reflect greater relative importance of seasonality of root growth and root inputs to the mineral soil. In another study at the Harvard Forest, Boone et al. (1998) calculated a Q 10 of 4.7 for root-related CO 2 production (seasonal growth and respiration). On the other hand, lower Q 10 values in the O horizon may simply reflect water limitation in that horizon in the late summer, when the soil is warm throughout the profile, but only the O horizon is sufficiently dry to experience suppressed CO 2 production.

Uncertainties of estimates of CO 2 production by horizon
The accuracy of the estimate of CO 2 production as a function of depth using this profile modeling approach is difficult to evaluate. When the diffusivity of the A horizon is overestimated, then the flux at the top of the A horizon is overestimated, which results in an underestimate of the production in the O horizon. Using data at the same site from only 1996, and using higher estimates of diffusivity generated by the Millington & Quirk (1961) model, Gaudinski et al. (2000 estimated that O horizon production was 190 g C m À2 yr À1 . We estimate twice as much mean annual O horizon production here using the Moldrup et al. (2000) model of mineral soil diffusivity. If we used our lower radonbased estimates of mineral soil diffusivity, the estimate of CO 2 production in the O horizon would be even higher. We suspect that both models may overestimate diffusivity in this mineral soil for two reasons. First, the radon measurements indicate lower diffusivities. Second, we encountered estimates of CO 2 flux at the top of the A horizon that exceeded mean chamber flux measurements during the dry summer conditions of 1999.
Several factors could explain this 1999 inconsistency, including: (1) an unknown source of error in our chamber flux measurements that was unique to 1999; (2) a greater importance of advection processes under dry conditions, which contribute to CO 2 efflux that is inadequately sampled by chamber measurements (Risk et al., 2002b); (3) inaccurate volumetric water content estimates, especially at the extreme dry end of the TDR calibration curve; (4) high sensitivity of the diffusivity model to small differences in water content, especially when the soil is very dry (Fig. 9a); (5) CO 2 concentration profiles in dry soils that conform more to a linear function than an exponential function of increasing concentration with depth; and (6) a systematic overestimation of mineral soil diffusivity of unknown origin from our application of the Moldrup et al. (2000) model.
Although the sensitivity analysis (Fig. 9b) showed that O horizon production was not strongly sensitive to the value of the RF parameter used to adjust water content at field capacity in Eqn (5), the presence of rocks also affects other measured parameters and may thus contribute indirectly, but systematically, to uncertainties in the model output. Rocks are included in our field and laboratory measurements of BD, and the resulting BD estimates are inserted directly into the model parameterization. Because we measured BD rather than modeled BD as a function of RF, the RF does not affect BD in the sensitivity analysis. Errors in sampling rocks would, however, clearly affect the field estimate of BD, which is one of the parameters to which the model is very sensitive (Fig. 9b), because it determines the calculation of total porosity (Eqn (1)). Hence, our estimates of porosity and diffusivity are also affected by rocks to the extent that we measured them adequately in 1 m Â 1 m soil pits. The same could be true for coarse roots, although they are usually less abundant than rocks in this glacial till soil. Some methods for measuring BD exclude rocks and some databases report only rock-free BD, which would cause an error in the estimate of total porosity in Eqn (1). The transition from O to A horizons is often difficult to identify unequivocally, and the presence of rocks makes it even more difficult to judge the depth of each horizon, which also influences the BD estimate. Finally, diffusivity models are based on studies of porous media, and the effects of the presence of large rocks interspersed in the matrix of small soil particles and pore spaces have not been addressed theoretically. The models address tortuosity of pore spaces among soil particles ranging in size from clay to sand and among soil aggregates, but they do not specifically address the effects of rocks on the tortuosity of diffusional pathways within soil profiles. Hence, our application of the model of Moldrup et al. (2000) could include some unrecognized biases because of the absence of consideration of the role of rocks in the underlying model structure.
Mismatching the profile modeling with the appropriate surface chamber flux measurements is also a potential source of error in our method of estimating O horizon CO 2 production. For preliminary data analysis during the first several years of our study, we used the mean of all six chambers at the tower site for comparing with profile estimates of CO 2 flux from both soil pits. However, when we noticed that the CO 2 concentrations at 26 and 55 cm depths in pit #2 began to diverge from concentrations at equivalent depths in pit #1 in 2001 (Fig. 3), we examined the fluxes from individual chambers to see if there was a similar divergence in surface CO 2 efflux. The three chambers closest to pit #2 did, indeed, yield higher fluxes than the other three chambers from 2001 onward (Fig. 3a). Therefore, we split the surface flux measurements into two means of three chambers each, and paired these with their respective nearby profile measurements. If we had not done so, the estimated flux at the top of the A horizon in pit #2 would have exceeded the overall site mean of the chamber flux measurements from 2001 onward.
Finally, it should be noted that this profile modeling approach assumes that the CO 2 concentrations are temporarily at steady state at the moment of measurement. We know that this assumption cannot always be true, especially following wetting events. Rapid changes of CO 2 concentrations at one or more depths could cause the concentration profile to diverge from a typical shape fitted by our exponential function (Eqn (6); Fig. 2). Data were not used when the exponential curve fit was not statistically significant at Po0.05, which may have eliminated some dates when CO 2 profiles were not near steady state. However, it is still possible that production estimates by horizon could be underestimated or overestimated if CO 2 concentrations were undergoing rapid changes at the time of sampling. Jassal et al. (2004Jassal et al. ( , 2005 demonstrated that this source of error would be brief in the coarse texture forest soil that they studied because diffusion of CO 2 through the soil profile was rapid relative to changes in rates of CO 2 production. Furthermore, because the sensitivity analysis in our study indicated that estimates of O horizon production of CO 2 were only moderately affected by uncertainties in the 'a' and 'C 1 ' parameters of the CO 2 concentration profile fits (Fig. 9a), departure from steady-state conditions of the CO 2 profile probably caused only modest errors.

Conclusions
Large uncertainties remain regarding the use of diffusivity models in soils, especially where rock content and spatial heterogeneity are large. The uncertainties of this measurement approach do not permit us to assign a specific error term to the estimates of production in the O horizon. However, our modeled estimate that the O horizon contributes at least 40% of the total annual soil CO 2 efflux at this site is probably conservative.
Despite uncertainty in the accuracy of the relative contributions of each horizon on an annual basis, vertical partitioning provided insight into the depth dependence of seasonal patterns of CO 2 production. The temperature sensitivity of CO 2 production varied across soil horizons, although the partitioning of CO 2 production by horizon did not improve the overall prediction of surface CO 2 effluxes based on temperature functions alone. On the other hand, partitioning of fluxes by horizon did demonstrate that the O horizon is a source of large seasonal and interannual variations in CO 2 production, and that water content covaried with CO 2 production only in the O horizon. A logical next step to this modeling approach would be to include estimates of accumulation and depletion of readily available carbon substrates by horizon, although this would require temporal estimates of carbon inputs that are difficult to obtain.
Although the O horizon of a mature or advanced successional forest may be at steady state over the long term of several years, the large interannual variation in estimates of O horizon CO 2 production demonstrated in this study indicate that this layer could be a transient interannual source or sink of ecosystem carbon. This transient source or sink of carbon in the O horizon is of similar magnitude as eddy covariance measurements of interannual variation in net ecosystem exchange of carbon.