Management Impacts on Carbon Dynamics in a Sierra Nevada Mixed Conifer Forest

Forest ecosystems can act as sinks of carbon and thus mitigate anthropogenic carbon emissions. When forests are actively managed, treatments can alter forests carbon dynamics, reducing their sink strength and switching them from sinks to sources of carbon. These effects are generally characterized by fast temporal dynamics. Hence this study monitored for over a decade the impacts of management practices commonly used to reduce fire hazards on the carbon dynamics of mixed-conifer forests in the Sierra Nevada, California, USA. Soil CO2 efflux, carbon pools (i.e. soil carbon, litter, fine roots, tree biomass), and radial tree growth were compared among un-manipulated controls, prescribed fire, thinning, thinning followed by fire, and two clear-cut harvested sites. Soil CO2 efflux was reduced by both fire and harvesting (ca. 15%). Soil carbon content (upper 15 cm) was not significantly changed by harvest or fire treatments. Fine root biomass was reduced by clear-cut harvest (60–70%) but not by fire, and the litter layer was reduced 80% by clear-cut harvest and 40% by fire. Thinning effects on tree growth and biomass were concentrated in the first year after treatments, whereas fire effects persisted over the seven-year post-treatment period. Over this period, tree radial growth was increased (25%) by thinning and reduced (12%) by fire. After seven years, tree biomass returned to pre-treatment levels in both fire and thinning treatments; however, biomass and productivity decreased 30%-40% compared to controls when thinning was combined with fire. The clear-cut treatment had the strongest impact, reducing ecosystem carbon stocks and delaying the capacity for carbon uptake. We conclude that post-treatment carbon dynamics and ecosystem recovery time varied with intensity and type of treatments. Consequently, management practices can be selected to minimize ecosystem carbon losses while increasing future carbon uptake, resilience to high severity fire, and climate related stresses.

Introduction to minimize forest carbon losses while increasing future carbon uptake, resilience to high severity fire, and climate related stresses.

Study site
The study was conducted at Blodgett Forest (38°54 0 N, 120°39 0 W), a University of California Research Station on the western slopes of the northern Sierra Nevada approximately 21 km east of Georgetown, California. Blodgett is a mixed-conifer forest composed of sugar pine (Pinus lambertiana), ponderosa pine (Pinus ponderosa), white fir (Abies concolor), incensecedar (Calocedrus decurrens), Douglas-fir (Pseudotsuga menziesii), and California black oak (Quercus kelloggii). Elevation ranges from 1100 to 1410 m. Total annual precipitation averages about 1600 mm, falling mainly between September and May. The average minimum daily temperature in January is 0.6°C and the average maximum daily temperature in July is 28.3°C [25]. The loamy-sandy soils are underlain by Mesozoic, granitic material and are predominantly classified as the Holland and Musick series (fine-loamy, mixed, semiactive, mesic Ultic Haploxeralfs [26]).
The forest is an actively managed timberland and has been repeatedly harvested with both even-and uneven-aged methods. Wildfires have been suppressed at Blodgett in the last century, reflecting a management history common to many forests in California and elsewhere in the western United States [27]. Prior to fire suppression, the median composite fire interval at the 9-15 ha spatial scale was 4.7 years with a range of 4-28 years [28].

Forest management treatments
This study includes six treatment types: un-manipulated control (CTRL), prescribed fire (FIRE), thinning (THN), thinning plus prescribed fire (THN+FIRE), and clear-cut harvest with and without mechanical soil ripping (HARV RIP , HARV NO_RIP , respectively). The overall goal of the FIRE and THN practices were to reduce the vulnerability of these forests to highseverity fires [29]. Soil ripping is a common post-harvesting practice used in commercial forests to reduce soil compaction and prepare the soil for seedling planting.
The treatments were applied to forest stands with similar climatic conditions, as plots were established within a 10 km radius. Pre-treatment species composition was similar among stands (S1 Table). Calocedrus decurrens and Abies concolor were the most common species, accounting for about 37% (ranging 34%-40% across treatments) and 27% (ranging 24%-30% across treatments) of the trees, respectively. Pseudotsuga menziesii accounted for 13% of the trees (ranging 11%-17% across treatments) and the combination of Pinus lambertiana and Pinus ponderosa accounted for 8% of the trees (ranging 6%-13% across treatments). The FIRE, CTRL, THN and THN+FIRE treatments were part of the Fire and Fire Surrogate Study (FFS). This study began in 2000 and analyzed the ecological effects of fuel treatments on vegetation structure and other ecosystem characteristics across a network of 13 locations in the United States. Additional information about the study can be found elsewhere [30]. In addition to the FFS treatment sites, we included four sites (each between 2000 and 7000 m 2 ) that were clearcut harvested in the summer of 2010. In these sites all trees were harvested and the residual material was piled and burned. In half of the harvested sites a tractor mechanically ripped soils with a wing-tipped subsoiler. In fall 2010, one-year old ponderosa pine seedlings were planted in pairs, on a 5 x 5 m grid spacing (840 seedlings ha -1 ). For clarity, a schematic diagram showing locations and timing of treatments and measurements is included in Fig 1. We analyzed tree-ring growth and stand biomass in the CTRL, FIRE, THN and THN+FIRE treatments. Measurements included three replicates per treatment (Fig 1) and continued until 2009, seven years after the prescribed fire and mechanical treatments. Because of practical limitations and instrumentation cost, the analysis on Fs and its driving variables (i.e., soil temperature and water content, fine roots, litter and soil carbon) was conducted in one, randomly selected replicate of the CTRL and FIRE treatments (CTRL Fs and FIRE Fs ) and the newly The control, mechanical thinning, prescribed fire, thinning followed by prescribed fire treatments each had three replicates. Soil CO 2 flux and soil carbon pools were measured in the clear-cut harvest treatment area and one of the replicates of the control and prescribed fire treatments. Timing of treatments (treat.) and measurements (meas.) during a period of 17 years, from 1995 to 2012, is also shown.
installed HARV sites (Fig 1). In these four areas soil CO 2 fluxes, litter, soil carbon, and fine root biomass were measured in 2011 and 2012 (more details on timing in Fig 1).

Soil CO 2 efflux measurements
Measurements of Fs were taken using the chamber technique [31] and the soil CO 2 profile technique [32]. At each site (CTRL Fs , FIRE Fs , and HARV RIP and HARV NO_RIP ), we periodically measured Fs in multiple locations using the chamber technique and continuously monitored Fs with the profile technique. Our approach of complementing the chamber and soil CO 2 profile techniques allowed us to quantify both spatial and temporal variability in Fs [33].
a. Soil CO 2 chamber technique. We measured soil CO 2 efflux using a Li-6000, and later a Li-6400, with a 10 cm diameter soil chamber attachment (Li-Cor, Lincoln, USA) at 110 locations (total for the CTRL Fs , FIRE Fs , and HARV RIP and HARV NO_RIP sites) bi-monthly during snow free periods (June 2011 to December of 2012). In the spring of 2012 the Li-6000 and Li-6400 were compared in the field obtaining good agreement (r 2 = 0.95) and therefore, the Li-6000 data were increased for the 15% difference found on the regression slope (data not shown). Measurements across all plots were taken over a two-day period, between 0900 and 1700 hours. The order of the sites and plots measured was randomly changed during each sampling date to avoid a systematic bias based on potential temperature changes. During these manual measurements, the chamber was positioned on PVC soil collars installed 1 cm into the soil to avoid soil disturbance and to allow repeat measurements on the same locations. Soil CO 2 efflux was calculated from the change of CO 2 concentration over time and averaged for two cycles over a 10 ppm range encompassing the ambient CO 2 concentration. Soil water content (SWC, integrated over the 0-5 cm depth using a HH2 and ML2x, DeltaT devices, Cambridge, UK) and soil temperature (Ts, measured at 10 cm depth using a 6000-09TC, Licor, Lincoln, USA) were measured adjacent to each Fs collar.
In both the FIRE Fs and CTRL Fs sites, measurements were collected on 29 different locations per site distributed over 3 ha. In each the HARV RIP and HARV NO_RIP sites, we measured Fs in 20 different locations distributed over 6 ha due to the spatial arrangement of the harvest units. More information about measurement location and spatial variability can be found in Dore et al. [34].
b. Soil CO 2 profile technique. We used solid-state CO 2 sensors (CARBOCAP model GMM 220, range 0-1%, Vaisala, Helsinki, Finland) installed at 2, 8, 16, and 24 cm soil depths to measure Fs at each site (CTRL Fs , FIRE Fs , and HARV RIP and HARV NO_RIP ). Ts was measured with thermocouples at the same depth used for the CO 2 sensors, and SWC (ECH 2 O, Decagon Devices, Pullman, WA, USA) was measured at two depths, 8 and 16 cm. Site specific calibration factors of SWC sensors and their dependence on Ts were quantified in the lab using each site's mineral soil sampled over a 20 cm depth. In winter 2012 sensors were calibrated against known CO 2 concentrations in the lab, and specific parameters determined during calibration were applied to each probe at each site. CO 2 concentrations from the sensors were measured every 15 seconds and averaged over 30 minute intervals. The concentrations were corrected for soil temperature and atmospheric pressure and were used to calculate Fs (expressed in μmol m -2 s -1 ) based on Fick's first law of diffusion Fs = -Ds dC/dz [32,35]. Ds is the diffusion coefficient in the soil and C is the CO 2 concentration at depth z. Ds is calculated by the product of the CO 2 diffusion coefficient in the free air Da, and the gas tortuosity factor ξ.
The factor ξ can be quantified empirically by quantifying its relationship with soil porosity (ε), which is determined by the SWC, bulk density, and particle density for the mineral soil [35]. We determined site and depth specific ξ in the lab using intact soil samples collected adjacent to the profiles (two replicates per site at 3 depths: 0-5 cm, 10-15 and 15-20 cm; 24 samples total). A sealed 100 ml volume above a 5 cm diameter, 5 cm long aluminum cylinder containing the undisturbed soil sample was initially flushed with nitrogen. Air was able to diffuse from the bottom of the cylinder through the soil. We measured the increase of oxygen with time in the volume above the soil (using a SO-200 Apogee, Logan, USA) for a period of 3-5 minutes after a 2-minute period necessary to reach a constant rate. For each site and depth, measurements were repeated at four different water contents (n = 96 measurements). Following each diffusion measurement SWC was determined using the fresh weight of the sample and bulk density obtained from its dry weight (102°C for 24 hours and until constant weight) at the end of all measurements. The resulting equation ξ = 0.0139e 6.2889Ãε (r 2 = 0.7, p< 0.001) was common to all soils from different sites and depths (Fig 2). To calculate Fs, however, we used the relationship between the tortuosity factor ξ and SWC (to obtain an exponential trend we used 1-SWC). This relationship ξ = 0.007e 4.094Ã(1-swc) had a higher coefficient of determination (r 2 = 0.85, p< 0.001) and was able to quantify the diffusion coefficient independently from accurate site/depth specific measurements of soil bulk density and particle density for the mineral soil. The effects of these soil characteristics were included in the site/depth specific relationship between SWC and ξ.
The CO 2 profile data were adjusted to match the mean spatial values of the chamber data for 2011 and 2012 (S1 Fig). To achieve this, we used for each site and year the slope and intercept of the linear relationships between daily averages of Fs simultaneously measured by chambers and each site-profile. Thus, the adjusted values represented the treatment-scale Fs and had at the same time high temporal resolution due to continuous measurements. At each site, daily averages of the adjusted Fs were fitted with daily averages of Ts and SWC profiles for 2011 and 2012 using the equation: where S 10 is basal respiration at 10°C, Q 10 is the parameter reflecting the temperature sensitivity of Rs, and a and b are the parameters of a simplified Gompertz function [36]. The measurement depth of Ts and SWC profiles which was able to explain most of the measured seasonal  [35]. (b) Seasonal trend of the diffusion coefficient at the treatment sites calculated from the specific relationship between tortuosity factor and soil water content ξ = 0.007e 4.094 * (1-swc) and the CO 2 diffusivity in free air, which was affected by soil temperature, atmospheric pressure, and soil water content monitored at each site.
variability in Fs was used to parameterize the equation. We used Ts at 10 cm at the HARV and FIRE Fs sites and Ts at 24 cm at the CTRL Fs site. SWC at 8 cm was used at all sites. The fitted equations explained 60% to 85% of the measured Fs variability (Table 1). Annual fluxes for 2011 and 2012 were obtained by adding daily Fs values derived from Eq 1. To calculate uncertainty for annual fluxes at each site we randomly selected 1000 sets of Eq (1) parameters, with each parameter drawn from normal distributions characterized by their fitted means and standard deviation. Each set of parameters was used to compute Fs for every day of the year applying the site specific Ts and SWC. The 1000 annual values calculated from each set of parameter were used to calculate uncertainty as 95% confidence interval.
Effects of treatments on Fs were analyzed using different methods. First, we compared chamber fluxes measured on the same day at the different sites, including with this approach both the effect of treatments on basal Fs and the effect due to post-treatment differences in Ts and SWC. Second, we used a general linear model of the form Fs = f(treatment, Ts, SWC) to assess the effect of treatments on Fs correcting for the effect of Ts and SWC, thus removing across-treatment differences in environmental conditions displayed during measurement campaigns. Third, we quantified the treatment effect over the different seasons, comparing the annual sums of Fs. To separate how much of the difference in annual fluxes among sites was due to the effects of treatments on Fs basal rates and how much was due to changed microclimatic conditions, we modeled Fs at the FIRE Fs and HARV sites using the CTRL Fs site's climatic inputs, as if treatments at these sites didn't affect environmental conditions but only the Fs basal rate.

Soil and fine root carbon pools
Fine root biomass was measured in summer 2012. Soil samples were collected from 0-15 and 15-30 cm, with a 20 cm 2 area auger at 10 locations adjacent to randomly selected Fs collars at  To calculate the annual soil CO 2 efflux at the control, fire, and harvest rip (HARV RIP ) and not ripped (HARV NO_RIP ) sites, the daily profile soil CO 2 efflux measurements were used to fit a semi-empirical model [36]. S 10 is the basal soil CO 2 efflux at the reference temperature (10°C); Q 10 quantifies the temperature dependence of soil CO 2 efflux; the parameters a and b are the soil water content dependence of soil CO 2 efflux. The non-linear regression r 2 and 95% confidence interval (in parenthesis) of the fitted parameters are also shown. These site-specific parameters, and average daily soil temperature and soil water content were used to calculate annual soil CO 2 efflux for 2011 and 2012. Both live and dead fine roots were hand-picked and separated in < 2mm and 2-5 mm diameter classes from each soil core. Fine root dry weight was determined after drying samples at 70°C for 24 hours. The litter layer was collected inside the 10 cm diameter Fs collars (27 plots at the FIRE Fs and CTRL Fs sites, 18 plots at both HARV sites) at the end of 2012, and dried at 70°C for 4 days before dividing it into leaves and woody components. Bulk density of the mineral soil layer (0-5 cm depth) was determined with 5 cm diameter cores collected inside the Fs collars (27 plots at both the FIRE Fs and CTRL Fs sites, 18 plots at both the HARV sites) after removing the litter layer, at the end of the Fs sampling campaign. Bulk density of the deeper 5 to15 cm soil layer was determined on 8 locations (5 cm diameter, 10 cm long samples), adjacent to randomly selected Fs collars at each site. Soil carbon content was determined on 18 samples per site on a 5 cm diameter, 15 cm deep volume. Each sample was separated into two depths, 0-5 cm, and 5-15 cm. Samples were air dried, sieved (2 mm), ground, and finally their carbon content was determined using an elemental analyzer (Thermo Flash 2000, Thermo scientific, MA, USA).

Tree biomass and annual tree radial growth
Tree biomass and tree radial growth were measured in the CTRL, FIRE, THN, and THN+FIRE treatments. Each treatment had three replicates ranging from 12 to 28 ha (12 experimental plots total). Tree diameter at breast height (DBH), species, and status were measured in 20 circular, 0.04 ha sub-plots distributed on a systematic grid in each experimental plot (240 subplots total). Only trees with diameters larger than 5 cm were considered in this study. Measurements of the same plots were collected the year before treatments (2001) In 2010, in each of the experimental plots, 30-60 trees of the five major species (Pinus lambertiana, Pinus ponderosa, Abies concolor, Calocedrus decurrens, Pseudotsuga menziesii), distributed over five size classes (DBH < 35 cm, 35-55 cm, 55-75 cm, 75-95 cm, and DBH > 95 cm), were cored and tree-ring widths measured to 0.001mm using a sliding stage to quantify annual growth from 1995 to 2009. The chosen starting point for this analysis was 1995 to have the same length, a seven-year period, before and after treatments.
For each treatment we determined annually individual and mean tree radial growth. The individual tree radial growth is the radial growth of a hypothetical average tree. We calculated it annually for each treatment by averaging the radial growth of each size class and species (i.e. 20 increments per each treatment). The annual mean tree radial growth was calculated by averaging the annual radial growth of all trees present in each treatment. The mean tree radial growth included both the effect of treatments on the growth of each tree (a positive effect in case of reduced competition or a negative effect in case of treatment related injuries) and the effect of decreased tree density (number of surviving trees changing with treatment and year). We used for each tree the mean annual radial growth of its species, experimental plot, and size class. If a size class had no samples, its growth was estimated as the average of the two adjacent classes (same species, year and experimental plot) or as the growth of the closest class (in case of missing smallest or biggest classes). Only 6% of the trees were in size classes with no tree sampled, and of those, 67% had DBH< 35 cm.
We used the allometric equations provided in Jenkins et al. [37] [38]. Aboveground and coarse root biomass was converted to carbon assuming a carbon concentration of 48% [39]. The difference in stand biomass carbon stocks between two consecutive years was used to express the stand annual aboveground tree productivity.
In fall 2010, one-year old ponderosa pine seedlings were planted in pairs on a 5 x 5 m grid in all two HARV sites. To calculate seedling biomass and to analyze the effects of mechanical soil ripping on growth rates, diameter and height of 240 to 400 seedlings were measured in each of the HARV RIP and HARV NO_RIP treatment sites during the summer of 2012. Eightyfive seedlings were collected and dry weight of needles and wood determined. The relationship between diameter and dry weight of needles, wood, and the whole plant was determined separately. The resulting equations were b L = 1.5305 D 2.1416 (r 2 = 0.73) for needles, b s = 0.9741 D 2.2691 (r 2 = 0.87) for stems, and b P = 2.5502 D 2.1834 (r 2 = 0.82) for the whole plant, where b L, b S, and b P was the dry weight of needles, stems, and whole plant, respectively, and D was the basal diameter (in cm). The equation for the whole plants was used to calculate the carbon stored as tree biomass at the HARV sites.

Data Analysis
When Fs measured during the same measurement campaign was compared across treatment sites, effects of treatments on Fs, Ts, and SWC were calculated from the slope of the linear relationship between chamber mean at each treatment site against simultaneously measured values at the CTRL Fs site. To evaluate the treatment effects we tested if mean values of Fs, Ts, and SWC at the treatment sites were significantly different using one way repeated measures analysis of variance (ANOVA). When comparing differences of annual Fs (and diffusivity) among treatments, we used one-way ANOVA based on daily means. To compare soil carbon, litter, and fine root biomass among CTRL Fs , FIRE Fs and HARV sites we used one-way ANOVA.
To quantify the effects of the treatments on tree-ring growth, tree and stand biomass carbon stocks we calculated the effect size using the Before-After Control-Impact (BACI) approach [40]. The BACI effect size = (μca -μcb)-(μtaμtb) is the differential change in means (μ) between the control (c) and the treatment sites (t) before (b) and after (a) treatment. The effect was evaluated by testing the existence of an interaction between period and treatment using a two-way repeated measures ANOVA and using the seven year averages for pre and post-treatment periods when available. Differences were considered significant at α <0.05 level.
The Tukey post hoc test was used to make multiple comparisons among treatments if a significant difference was detected. Normality and equal variance tests were conducted and, if necessary, data were log transformed to meet these conditions. S1 Table summarizes and describes statistical tests used to compare stand characteristic, carbon pools, or fluxes among treatments.

Effect of treatments on microclimate
Both fire and tree harvesting reduced tree density of the stands (p < 0.001; Table 1 and S1  Table). Consequently the amount of energy reaching the ground increased, as the direct comparison of Ts simultaneously measured at the different sites showed. Soil temperature was 30% higher in the HARV site and 16% higher in the FIRE Fs site (p < 0.001; Fig 3b, S1 Table) compared to the CTRL Fs site. SWC was 29% higher at the HARV site and 11% higher at the FIRE Fs site compared to the CTRL Fs site (p < 0.001 ; Fig 3c, S1 Table). Thus, differences in both Ts and SWC with the CTRL Fs site were highest at the HARV sites where all vegetation was removed. Seasonal variations of Ts and SWC were highest at the HARV sites and lowest at the CTRL Fs site.

Effect of treatments on soil CO 2 efflux
Soil CO 2 efflux (measured using the chamber technique) was reduced by treatment (p < 0.001; S1 Table). Comparing Fs measured at the different sites on the same measurement date, the FIRE Fs and HARV sites (average of HARV RIP and HARV NO_RIP ) were 14% and 18% lower, respectively, compared to the CTRL Fs site (Fig 3, S1 Table). Soil CO 2 efflux was not significantly different between HARV RIP and HARV NO_RIP sites (p > 0.05).
CO 2 concentrations at HARV sites were higher and showed greater seasonal variation (Fig  4) than at the FIRE Fs and CTRL Fs sites, especially at deeper depths (16 and 24 cm). Concentrations of CO 2 in a soil layer are the result of production and the capacity of CO 2 to exit the soil (i.e. diffusion). The diffusion coefficient was highest at the CTRL Fs site and lowest at the HARV RIP site (p < 0.001; Fig 2b, S1 Table).
FIRE Fs and HARV sites had reduced annual Fs compared to the CTRL Fs site (p < 0.001; Table 1 and S1 Table). The annual emissions of treated sites decreased despite having higher Ts and SWC (Table 1 and S1 Table). Seasonally, sites were more similar during winter and spring, when differences between treated and control site were on average 0.5 μmolm -2 s -1 (Fig 5a) compared to the summer, when difference were 2.5 μmolm -2 s -1 . In July-August, when soil water availability decreased (Fig 5c), Fs decreased at all sites. Decrease in Fs started 2-5 weeks earlier at the treated sites than at the CTRL Fs site for both years (Fig 5a).

Effect of treatments on soil carbon, fine roots, and litter
Soil carbon (up to a depth of 15 cm) was different among treatment sites when expressed as % content, however not when expressed as g C m -2 , due to the different soil bulk densities of the sites (p = 0.03; Table 2 and S1 Table). Fine root biomass was similar at the FIRE Fs site compared to the CTRL Fs site (583 and 517 g C m -2 respectively), but was reduced by 60-70% (to 143 g C m -2 at the HARV NO_RIP and 191 g C m -2 at the HARV RIP site) two years after the harvest. The litter layer was reduced at all treatment sites compared to the CTRL Fs site, ranging from -40% in the FIRE Fs (p < 0.008) and HARV NO_RIP sites to -80% in the HARV RIP site (p < 0.001; Table 2 and S1 Table). If we consider carbon stored as mineral soil, fine roots, and litter layer, three years after the second fire (initial 2002, re-burn in 2009, Fig 1), the FIRE Fs site had 19% less carbon than the CTRL Fs site, and two years after the clear-cut the HARV site had 26% less carbon than the CTRL Fs site.

Effects of treatments on forest biomass and tree growth
Pre-treatment biomass and tree density did not differ significantly among sites (p >0.05; Table 1 and S1   Over seven years after the treatments, individual tree annual radial growth decreased after FIRE (-10%; p < 0.001), increased 25% after THN (p < 0.001), and was not significantly different after THN+FIRE (p = 0.51 Fig 7a, Table 3 and S1 Table). For each species (Pinus lambertiana, Pinus ponderosa, Abies concolor, Calocedrus decurrens, Pseudotsuga menziesii) FIRE decreased the annual radial growth by 7-20%, whereas THN increased annual radial growth by 13-30%. The THN+FIRE treatment had a mixed effect, increasing annual radial growth in some species (Abies concolor + 33%) but decreasing it in others (Pseudotsuga menziesi-13%; Table 3 and S1 Table). Effect of treatments on individual tree annual radial growth depended on tree size. Annual radial growth increased most in smaller diameter trees (<35cm diameter) for all treatments, with the highest increase at the THN treatment (52% increase compared to a 3% decrease at the FIRE treatment, Table 3).
When averaging the tree annual radial growth over the stand, because of the lower posttreatment tree density (Fig 7b), the before-after analysis showed mean annual radial growth decreased 27% after FIRE (p = 0.03) and was not affected by THN (p = 0.4) and THN+FIRE treatments (p = 0.42; Table 3 and S1 Table). All active treatments reduced stand aboveground productivity (expressed as difference between biomass of two consecutive years), however the THN+FIRE treatment had the lowest post-treatment average productivity of 3 t C ha -1 (2003-2009), approximately 50% lower than before treatment (p < 0.001; Fig 7c, Table 3).
The CTRL treatment was still accumulating carbon during the observed 14-year period. However, radial growth decreased between the periods 1995-2001 and 2003-2009 from 3.04 (±0.14 SD) to 2.65 (±0.31 SD) mm yr -1 . From 2001 to 2009, carbon stored as tree biomass  increased from 205 (±37 SD) to 240 (±39 SD) t C ha -1 at the CTRL treatment, was reduced ca. 10-15% by FIRE (p = 0.005) or THN (p = 0.012), and 36% by the THN+FIRE treatment (p <0.001; Table 3 and S1 Table). All active treatments decreased tree density, however the greatest decrease occurred in the THN+FIRE (54% effect size, p <0.001; Fig 6, Table 3 and S1 Table), and was higher after THN than after FIRE (43% decrease compared to 32% respectively). In the FIRE treatment the number of live trees decreased substantially (> 100 trees   Table 3 and S1 Table).
In the first year after treatment, the FIRE treatment did not lose any carbon stored as tree biomass compared to losses of 29 (±0.10 SD) t C ha -1 in the THN treatment and the 52 (±0.11 SD) t C ha -1 at the THN+FIRE treatment (Table 3 and S1 Table). However, between the first and seventh post-treatment years, the FIRE and the THN+FIRE treatments both accumulated a total of ca. 5 (±3 SD) t C ha -1 compared to 30 (±3 SD) t C ha -1 accumulated by the THN treatment, and 27 (±3 SD) t C ha -1 at the CTRL. Both the FIRE and the THN treatments recovered to pre-treatment biomass carbon stocks levels by the seventh post-treatment year (p = 0.13; Table 3 and S1 Table).
At the HARV sites (Table 4), soil ripping did not benefit the seedlings 20 months after planting. No differences were found between HARV RIP and HARV NO_RIP in survivorship, height, and diameter and thus biomass carbon stocks, still~1 g C m -2 .

Discussion
Prescribed fire, thinning, and clear-cut harvest all affected measured pools of carbon in mixed conifer forests in the north-central Sierra Nevada. These treatments also affected ecosystem carbon fluxes, decreasing carbon sequestered annually as tree biomass, or released from the soil. Prescribed fire and thinning affected biomass carbon stocks and growth differently. Over seven post-treatment years, the THN treatment showed a stronger sink strength and reached higher biomass levels compared to treatments that include fire. Our results confirmed that biomass carbon stocks and growth depend on the intensity and type of treatment, and that it is possible to design forest restoration treatments that achieve both increased forest resiliency and carbon storage.

Treatments and soil CO 2 efflux
The effect on Fs was in part due to the reduction of tree cover, which affected the main drivers of ecosystem processes: energy and water (Table 1 and S1 Table). In addition, forest management treatments altered basal Fs (Table 1 and S1 Table), which is related to biotic factors such as soil microbes, root biomass, and soil organic carbon content. The lower basal respiration at the HARV site can be explained by the lack of autotrophic CO 2 production in the soil.
When we remove the effect of environmental conditions by modeling a hypothetical Fs at the FIRE Fs and HARV sites using the CTRL Fs site's climatic inputs (Fig 5d), we obtained annual Fs lower (circa 300 g C m -2 year -1 ) than the annual current rates. Thus effects of forest management treatments and effects on microclimate altered soil basal Fs in opposite directions: the reduced basal Fs lowered the annual Fs, but the increase in soil temperature and SWC increased Fs. Because at our treated sites Fs was decreased by treatment (Table 1 and S1 Table), Forest Management Effects on Carbon Dynamics the effect of basal Fs was stronger than the effect of microclimate. Complex effects of disturbances on Fs were also found in other studies, where factors acted independently on autotrophic and heterotrophic processes and compensatory mechanisms caused disturbances effects in part to cancel out [7,41].
In several studies analyzing the effects of fire and thinning on Fs, researchers also found an increase of soil temperature and soil water content after forest management [17,[41][42][43]. However, effects on Fs varied. Burning increased and thinning reduced Fs in forests similar to those studies here [42]. Previous work in young ponderosa pine forests demonstrated that Fs generally decreased only when thinning and burning were combined [41,44]. Misson et al. [7] and Concilio et al. [45] concluded that the response to forest management treatments is proportional to their intensity. Disagreement in results could be due in part to the difficulty to detect treatment effects. In a previous study conducted at our research site we determined that disturbance increased spatial variability of Fs [34], mirroring post-disturbance increases in spatial variability of Ts and SWC and reducing the ability to accurately quantify Fs and to detect differences between treatments [34].
Differences in Fs among sites were not a result of large differences in maximum and minimum Fs rates, but to the timing of their fluctuations. For example, we measured a decline in Fs during summer drought, as observed in several other studies [7,[46][47][48], and treatments changed the onset and duration of these periods (Fig 5a). During the low SWC summer periods, the CTRL Fs site showed higher Fs than the FIRE Fs and HARV sites (Fig 5a), even if SWC was lower at this site (Fig 5b and 5c). This could be due to the difficulties in measuring SWC at the appropriate depth, so that the top 5 cm measurement does not represent the layer effectively controlling Fs. Summer SWC integrated over the root depth could be higher at CTRL Fs than at the treated sites and should be used in future research. An alternative explanation is the prevailing contribution of the autotrophic component at this site, together with an independent control of Ts and SWC over this component. Tang and others [41] found that the heterotrophic component of Fs was more susceptible to seasonal drought than the autotrophic component in a nearby forest.
Our values of Q 10 , ranging from 3.8 (± 0.18 95% CI) to 2.2 (± 0.43 95% CI), showed a decreasing trend going from the CTRL Fs , to the FIRE Fs and the HARV sites (Table 1), such that temperature effects were reduced where seasonal temperature variations were higher. Our values matched the values reported by Kobziar et al. [44] in a similar forest located about 100 km south from Blodgett, ranging from 1.9 to 5.38. A lower Q 10 at the HARV sites, where autotrophic respiration is lowest, agrees with the lower Q 10 for heterotrophic decomposition than for autotrophic respiration reported in Ngao et al. [36]. Our higher Q 10 at the CTRL Fs with older trees concurs with the higher Q 10 for the older stand in a chronosequence of coastal Douglas-fir stands [49]. In our sites Q 10 was higher at the CTRL Fs , where basal respiration was highest, a finding also reported in Davidson et al [46]. Finally Q 10 decreased at the HARV sites, where temperature was highest, showing a dependence on soil temperature also reported in Tjoelker et al. [50].
We measured seasonal peak daily Fs of 6-8 μmol m -2 s -1 and annual (control treatments) values ranging between 1200 and 1300 g C m -2 yr -1 . Our values are in the range of fluxes measured in previous studies in similar forested ecosystems [25,[41][42]44]. The relatively high annual Fs (for example compared to the range of values for temperate coniferous forests in Subke et al. [51]) was in part explained by a Fs activity that persisted late in the fall. Soil temperature was > 10°C until November, and warm and wet conditions are favorable to both heterotrophic and autotrophic processes [46]. Our high winter Fs values (2-4 μmol m -2 s -1 ) are supported by Xu et al. [25] winter Fs values in a ponderosa pine plantation adjacent to Blodgett Forest (3 μmol m -2 s -1 ).
The diffusion coefficient of CO 2 was lower at the HARV and FIRE Fs sites compared to the CTRL Fs site (Fig 2, S1 Table). Differences were not due to modifications in soil physical properties because the relationship between tortuosity and soil water status were not different (Fig  2a). However SWC was higher at these sites compared to the CTRL Fs site. Thus soil water content had a prominent role in preventing CO 2 to exit from the soil surface [46] and could be the cause of the higher CO 2 concentration recorded at the HARV site (Fig 4). Higher CO 2 concentration could also be a result of higher CO 2 production, due to the higher temperature and SWC favorable to decomposing processes, or more abundant substrate for decomposition (however we did not find this in terms of soil carbon, litter, or fine roots).

Fuel Treatment Impacts to Biomass and Tree Growth
After seven years, prescribed fire reduced tree density by 32% and resulted in a 28% reduction in tree biomass productivity and a 10% reduction in individual tree annual radial growth. These reductions were from a combination of direct mortality from fire and secondary mortality [10,52]. Secondary mortality had a stronger influence on reducing post-fire stand biomass than increased growth rates of the surviving trees, which was also found in other research [17].
Trees that survived the prescribed fire were not able to make efficient use of the increased pool of inorganic nitrogen after fire consumed a large amount of forest floor [53]. A previous related study assessed the impacts of fire and thinning on tree vigor and vulnerability [54] and found negative impacts for fire and positive for thinning. The authors proposed that the negative fire impacts could be minimized by delaying prescribed fire 2-5 years following thinning/ mastication, or removing whole-trees instead of mastication which left all activity fuels on site [54].
Over the seven post-disturbance years observed in this study, thinning had a smaller impact on ecosystem carbon dynamics than treatments that included fire. Post-thinning tree density was relatively stable compared to the FIRE treatment. The reduced competition for resources in the THN treatment had a positive impact on tree radial growth on both individual tree and mean tree growth of the treatment. Accumulated carbon was similar than the CTRL treatment both of which were much larger than the FIRE treatment. Nevertheless, if we limited our inventories to the first post-disturbance year, we would have found the opposite conclusion with a stronger impact of thinning than fire (THN reducing biomass carbon stocks 18% and tree density 41%, compared to FIRE reducing biomass carbon stocks 4% and the number of trees by 14%). Fire induced mortality persisted almost a decade, whereas when fuel and biomass were treated mechanically the impacts were concentrated in the first year post treatment.
When analyzing effects of thinning and fire on forests carbon dynamics it is important to include not only the carbon stored in trees, but also the nature and fate of the carbon removed during management operations. With a prescribed fire, most of the carbon removed is released during combustion directly to the atmosphere [8,10]. Alternatively wood partially burned and standing dead (snags) can take decades to decompose. When biomass is removed by logging, only carbon stored in residual matter that is burned, or in wood used for energy production, is immediately released to the atmosphere. Carbon used as paper and shipping material is stored for longer periods (1-6 years), and for the longest period when wood is used as solid and wood composite products, particularly if used in home construction [55]. Litter, coarse woody debris, and understory vegetation may be only minimally affected by thinning [3,11,53]. However carbon is also emitted for transportation of wood products to final destinations and by logging equipment during thinning. At Blodgett forest, thinning did not significantly affect the litter layer or coarse woody debris, whereas both were reduced by prescribed fire [10,53]. Wood logged by thinning at this forest was sold (with a positive revenue) mostly as construction lumber [55]. Fifty percent of wood removed had a life span of 50 years, and only circa 15% was released to the atmosphere in the first year [10]. In this forest, compared to prescribed fire, thinning caused smaller carbon losses during treatment, kept carbon stored for a longer period, did not remove carbon stored at the soil surface, and increased tree and stand growth.
Despite the decreased productivity and the differences in effects between FIRE and THN treatments, after seven years both of these treatments reached pre-treatment biomass carbon stocks levels, and still maintained a lower fire hazards because of the reduction or removal of surface and ladder fuels [6,20,31]. For the THN+FIRE treatment, where 54% of trees were removed or killed, we did not observe any recovery during the first 7 post treatment years. The THN+FIRE treatment displayed the beneficial effect of reduced competition among trees on radial tree growth (also observed at the THN treatment) for two species (A. concolor and C. decurrens). However, the positive effects on these species were attenuated with the negative effect caused on other species (especially P. menziesii).
In comparison with the partial disturbance due to fire or thinning, the clear-cut harvest had the strongest effect on carbon pools and fluxes. Two years after harvesting, carbon stored above ground and as tree biomass was reduced 80% and stand productivity was close to zero. At the same time soil fluxes were reduced by a similar amount as the fire site (ca. 15%). Soil ripping did not have a positive initial effect on seedling growth but it increased the impacts of harvesting on carbon dynamics, as indicated by the larger decrease, compared to the HARV NO_RIP site, of soil carbon and litter.
The intensity of the treatments impacted ecosystem carbon dynamics. A large decrease in tree density (55% reduction) reduced stand biomass carbon stocks and productivity for at least seven years in the thinning followed by fire treatment. When tree density was decreased to a smaller extent (40% reduction) by thinning alone, tree growth was higher than at the untreated stands, and biomass carbon stocks quickly recovered. The intensity of the treatment was not the only determinant factor, the type of disturbance was also important. Fire and thinning, applied singularly or combined, had different impacts on stand productivity and on the time needed to recover to pre-treatment carbon levels. The prescribed fire removed only a relatively small number of trees, with no change in biomass carbon stocks in the first post-treatment year, but tree injuries caused secondary tree mortality and low growth. Therefore, long-term studies are needed to assess the effects of management treatments on carbon dynamics, because studies limited to the first post-disturbance year will not capture the multiple and complex interactions involved. Finally, these results provide insights for managers pursuing forests carbon sequestration to select treatment types and intensities that minimize carbon losses, speed recovery time, and increases post-disturbance carbon uptake.  Table. Statistical analysis. a) Description of the statistical analysis of soil carbon pools and CO 2 efflux (Fs) in a mixed conifer forest subject to fire (FIRE), clear cut harvesting with (RIP) and without (NO_RIP) soil ripping, and an undisturbed forest (CTRL) at Blodgett Forest Research Station in the Sierra Nevada. 1, 2 compares soil CO 2 efflux measured using the chamber technique (chamber). 1 Test for total effect on Fs, where the general linear model 2 characterizes interactions among Fs, soil temperature, and soil water content. b) Statistical analysis of stand characteristics, biomass, growth, and productivity of a mixed conifer forest subject to fire (FIRE), thinning (THN), thinning plus burning (THN+FIRE) and an undisturbed forest (CTRL) at Blodgett Forest Research Station in the north-central Sierra Nevada. (DOCX)