Response of Soil CO₂ Efflux to Shelterwood Harvesting in a Mature Temperate Pine Forest

Abstract

In forest ecosystems, soil CO₂ efflux is an important component of ecosystem respiration (R_E), which is generally driven by variability in soil temperature and soil moisture. Tree harvesting in forests can alter the soil variables and, consequently, impact soil CO₂ efflux. This study investigated the response of total soil CO₂ efflux, and its components, to a shelterwood harvesting event of a mature temperate white pine (Pinus strobus L.) forest located in Southern Ontario, Canada. The objective was to explore the response of soil CO₂ effluxes to changes in the forest microclimate, such as soil temperature and soil moisture, after shelterwood harvesting removed approximately one-third of the overstory canopy. No significant differences were found in both soil temperature and soil moisture between the pre-harvesting (2008–2011) and post-harvesting (2012–2014) periods. Despite similar soil microclimates, total soil CO₂ effluxes were significantly reduced by up to 37%. Soil CO₂ effluxes from heterotrophic sources were significantly reduced post-harvesting by approximately 27%, while no significant difference in the mineral-soil horizon sources were measured. An analysis of R_E, measured with an eddy covariance tower over the study area, showed an increase post-harvesting. However, the overall net ecosystem carbon exchange showed no significant difference between pre- and post-harvesting. This was due to an increase in the gross ecosystem productivity post-harvesting, compensating for the increased losses (i.e., increased R_E). This study highlights the complexities of soil CO₂ efflux after a disturbance, such as a harvest. The knowledge gained from this study adds to our understanding of how shelterwood harvesting may influence ecosystem carbon exchange and will be useful for forest managers focused on carbon sequestration and forest conservation.

1. Introduction

Forest ecosystems account for about 80% of the world’s biomass and are considered an important sink for atmospheric carbon dioxide (CO₂) [1]. The forest carbon cycle is strongly influenced by the environmental and edaphic factors [2,3,4,5,6,7,8,9,10]. In addition, silvicultural practices, such as forest thinning, can also have a profound impact on forest carbon dynamics by changing soil temperature, soil water content, soil organic matter, root biomass, and microbial activity [11,12,13,14,15,16]. These silvicultural practices are being adopted in North America and across the world, not only to manage forests, particularly planted stands, but also to enhance their carbon sequestration capabilities and/or to restore them to their previous native species [15].

In the Great Lakes—St. Lawrence Region of Canada, the Ministry of Natural Resources and Forestry (OMNRF) uses a number of silvicultural practices based on the characteristics of the current forest, its history and the desired forest condition. One such practice is a specific shelterwood harvest, which consists of selectively removing mature trees, while maintaining legacy trees (largest, best quality), in a series of two or more partial cuts, over the course of 10–30 years. Shelterwood harvesting encourages natural regeneration and development of seedlings in partial shade; [17,18,19]. It changes the tree spacing, number, and size distribution, which can affect the dynamics of carbon uptake through above ground productivity and carbon losses through soil respiration processes [11,20]. A more open canopy alters the radiation dynamics, interception of precipitation, wind flow patterns, water vapour deficit, soil temperature, and soil moisture [11,21,22]. Shelterwood harvesting also causes changes in root density and production altering autotrophic respiration, litter-fall input and, hence, soil microbial activity. These changes may cause variations in soil CO₂ efflux, also referred to as soil respiration [23].

Peng et al. [12] found that amongst the studies that examined the impacts of forest harvesting on soil respiration, only a few have examined selective harvesting or thinning effects, with little information on spatial and temporal patterns of soil CO₂ efflux response in these forests, particularly in conifer stands. This is important because an increasing number of temperate conifer forests in Eastern North America are being managed to enhance their carbon sequestration capabilities to restore them to their native species composition, enhance biodiversity, and conserve water resources. Therefore, there is an urgent need to enhance our understanding of the effects of silvicultural practices, in particular shelterwood harvest, on soil processes and carbon sequestration [24].

In this study, we investigated the response of soil CO₂ efflux and its components to a shelterwood harvesting event in a mature temperate white pine (Pinus strobus L.) forest, located near Lake Erie in Southern Ontario, Canada. The purpose of the study was to explore the response of soil CO₂ efflux dynamics to changes in forest microclimate, such as soil temperature and soil moisture, after a partial (i.e., selective) thinning event. Both soil temperature and soil moisture have been previously shown to strongly influence soil respiration at this site [6]. Therefore, seasonal variability of soil temperature, soil moisture and soil CO₂ efflux were measured using an automated soil CO₂ chamber system [25,26], and comparisons were made between the pre- (2008–2011) and post-harvesting (2012–2014) periods. Given that soil respiration is a major component of R_E at the site [27], we also investigated the contribution of soil CO₂ efflux after shelterwood harvesting to R_E, measured with the eddy covariance method.

2. Methods

2.1. Study Site

This study was conducted at a mature white pine forest near Turkey Point Provincial Park on the northern shore of Lake Erie in Southern Ontario, Canada. This site is part of the Turkey Point Flux Station (TPFS), or Turkey Point Observatory of the Global Water Futures (GWF) program, and global Fluxnet. The forest is dominated (>82%) by eastern white pine (Pinus strobus L.), which was planted in 1939 to stabilize local sandy soils. As of 2014, the average tree height was about 21.8 ± 1.7 m, and stand density was about 421 ± 166 stems·ha⁻¹. Other tree species include 11% balsam fir (Abies balsamea L. Mill) and species native to the hardwood forests of North American Eastern Temperate Forest Ecoregion, including 4% oak (Quercus velutina L., Q. alba L.), 2% red maple (Acer rubrum L.), and some wild black cherry trees (Prunus serotina Ehrh). The understory consists of ferns (Pteridium aquilinum L.), mosses (Polytrichum spp.), poison ivy (Rhus radicans L. ssp.) and Rubus species [28,29]. The topography at the site is fairly flat with well drained sandy soil (Brunisolic Gray Brown Luvisol, following the Canadian Soil Classification system), which is composed of ~98% sand, 1% silt, and 1% clay. Further site details are provided in Arain and Restrepo-Coupe [30] and Peichl et al. [28,29].

2.2. Shelterwood Harvest Treatment

This forested site is managed by the Ontario Ministry of Natural Resources and Forestry (OMNRF) under the shelterwood silvicultural system [17,18,19]. This system is characterized by two or more partial thinning treatments over a few decades to allow for regeneration and development of seedlings in partial shade. The forest was first thinned in 1983, where one-third of trees were randomly removed. In March of 2012, a second partial cut was conducted and approximately 30% of the dominant overstory trees were selectively removed to improve light and water availability and stimulate growth of the remaining trees. The long-term aim of the management is to convert the pine plantation back to the native stand with mixed-wood ecosystem characteristics [19]. A mechanical harvester was used to cut, de-limb and section the selected trees. To reduce soil compaction and disturbance, the harvester used pre-existing multi-use trails throughout the stand. Thinning residues (e.g., limbs, bark, crowns) were placed on the ground ahead of the harvester to limit compaction when trails were not accessible [17]. After shelterwood harvesting, the average basal area was reduced by 13% [31], while peak leaf area index was reduced by 35% [17]. It was estimated that the thinning added 43 t of foliage, 200 t of live branches, and 129 t of dead branches to the forest floor, while an additional 515 t of roots were made inactive from growth and production [17]. These components will decompose and are expected to contribute significantly to the ecosystem respiration.

2.3. Soil CO₂ Efflux Measurements

Soil CO₂ efflux was measured using an automated chamber system, as described below. The main advantage of an automated chamber system, compared to manual measurements, is the ability to take continuous long-term measurements of soil CO₂ efflux, which are useful in studying short- and long-term responses of soil CO₂ efflux to environmental variables [7]. The chamber system occupied an approximate area of 50 m × 50 m and was located about 100 m north from the eddy covariance flux tower at the site. The experimental area was chosen to complement the ecosystem carbon flux measurements by the flux tower. The individual chamber locations, although few in replicates due to logistical and budgetary constraints, were placed in locations within the footprint of the eddy covariance flux tower measurements, with the aim to provide data on the soil CO₂ efflux contribution to the overall ecosystem respiration. Chambers were placed away from tree stems, in order to avoid large roots, particularly for those chambers measuring heterotrophic respiration. Additionally, the chambers were located away from thick understory and placed on level ground. As soil respiration is comprised of contributions from the mineral-soil, the organic-litter horizon and roots with their associated mycorrhizae, some of our collars were strategically placed to elucidate the contribution from these different components to the total soil respiration.

The chamber system was established to measure total soil CO₂ efflux, including those from root exclusion and forest floor removal plots. The chamber system was installed in June 2008, with the initial installation of four chambers. One chamber was established to measure heterotrophic respiration (F_RT1) in the trenched plot whereby live tree roots were severed in the area surrounding the chamber [32]. The other three chambers were used as a control chambers and measured total soil CO₂ efflux (F_S: F_S1, F_S2, and F_S22). One of the control chambers (F_S22) was altered in May 2009 to measure the contributions from the mineral-soil horizon (F_LR) to total soil respiration, by removing the litter–fermenting–humified (LFH) layer. Additional litter, which may have settled inside this chamber, was periodically removed, including any material from the harvesting event. Additionally, in May 2009, two more control chambers (F_S3 and F_S4) were added to the chamber system to measure total soil respiration. Finally, in May 2010, two additional chambers were installed, one as a control (F_S5) and an additional chamber to measure heterotrophic respiration (F_RT2); increasing the total number of chambers in the chamber system to eight. In total, five chambers were measuring total soil CO₂ efflux, two chambers measured heterotrophic soil CO₂ efflux, and one chamber measured contributions from the mineral-soil horizon.

Soil CO₂ efflux was continuously measured at each chamber location, using automated non-steady state chambers developed by the University of British Columbia [7,25,26,33]. Each chamber consisted of a collar and a chamber lid that fit onto the collar. The collar was constructed from a polyvinyl chloride (PVC) cylinder and were inserted approximately 2–4 cm into the ground to avoid cutting roots near the soil surface. Additionally, the height of the collar was maintained at a minimum height above the ground so when the chamber lid was open, there was no reduction in air flow over the soil surface within the collar [7]. The chamber lid consisted of a transparent plastic dome fixed to a metal frame and the lid was the same diameter as the collars. A torsion spring provided enough force to close the chamber dome, while a chamber-mounted, two-way pneumatic cylinder opened the dome when compressed air was pushed through the tubing [32,33]. When the dome closed during measurements, a foam gasket sealed the dome around the collar on which the chamber rested, providing an air-tight seal. A small fan inside the chamber circulated the air, preventing stagnation zones and microclimate formation. When a chamber was not in use, the chamber lid was kept open as to allow precipitation and litter to naturally settle into the collar area. Any vegetation that grew inside the collars between visits was removed to eliminate any potential photosynthesis effect.

The chambers were controlled by a computer housed in a datalogger box, which collected and processed data, and controlled the instrument pumps and gas analyzer equipment. Cycling through the chambers, including opening and closing their lids, was done with custom MATLAB (The MathWorks Inc., Natick, MA, USA) programs. Each chamber closed individually for a one-minute interval, where the CO₂ concentration was sampled and measured by an infrared gas analyzer (model LI-840, Li-COR Inc.). Measurements were cycled through the eight collars for a total of three cycles per half-hour period. Thus, each half-hour consisted of three, minute-long CO₂ concentration measurements per chamber. One half-hourly value was produced for each collar by averaging the three measurements taken during that half-hour. These average values were used for computation and analysis. Concentrations measured during the fifteen seconds following dome closure were discarded to ensure the samples tubes were relatively free from air from the previous sampled chamber and to account for the time taken for the chamber lid to completely close.

In order to measure soil CO₂ efflux (F_S) with these chambers, the measured time rate of change in chamber headspace CO₂ mole fraction or mixing ratio are estimated [7]. The soil CO₂ efflux (µmol CO₂ m⁻² s⁻¹) from each chamber was calculated using Equation (1) as shown below:

$$F_S={\rho }_a\frac{V_e}{A}\frac{dSc}{dt}$$

(1)

where ρ_a is air density in the chamber headspace (μmol m⁻³), Ve is the effective volume of the chamber (m³), A is the area (m²) of the soil surface covered by the chamber collar, and dSc/dt is the time rate of change of CO₂ mixing ratio in the chamber headspace (μmol CO₂ mol⁻¹ dry air s⁻¹).

The mean Ve value was calculated using the Equation (2) by injecting CO₂ through the top of the chamber domes for one minute and recording the CO₂ concentration change (in ppm):

$$V_e=\frac{IRT}{PV\left(S_c-S_m\right)}$$

(2)

where S_c (μmol CO₂ μmol⁻¹ dry air s⁻¹) is the rate of CO₂ concentration increase during the calibration period, I (μmol CO₂ s⁻¹) is the rate of injection of CO₂ during the calibration period, S_m (μmol CO₂ μmol⁻¹ dry air s⁻¹) is the rate of change of CO₂ concentration, P is the atmospheric pressure (Pa), V is the volume of the chamber (m³), T is the chamber air temperature (K), and R is the universal gas constant (8.314 J μmol⁻¹ K⁻¹).

According to the work done by Jassal et al. [7], an accurate estimation of the effective chamber volume will help in acquiring good measurements of soil CO₂ efflux. Effective chamber volume varies seasonally due to snow and litterfall, and the expansion and contraction of the chambers and collars due to variation in air temperature.

Due to limited weekly field site visits, chamber volumes were measured a few times during the growing season, primarily to check changes in the litter layer within the collars and to remove vegetation. The mean value of the Ve, calculated using Equation (2), was 0.069 m³ during the snow-free season [32]. During the snow season, snow depth measurements were taken inside the chambers, when possible, and used to adjust the Ve.

Measuring soil CO₂ efflux during the winter was challenging as snowfall can completely fill the collars, heavy ice and snow on the dome can force the chamber to remain shut for extended periods of time, and cold temperatures can create leaks in sample tubes or compressed air lines [32]. Frequent site visits were made during the winter season to remove snow from the top of the dome and around the collar rim to ensure a continued tight seal for measurements.

2.4. Data Processing and Quality Control of Soil Flux Data

Soil CO₂ efflux data were quality controlled with a custom script written in MATLAB software (The MathWorks Inc., Natick, MA, USA). In the script, the slope of time versus efflux in the high-frequency (1 Hz) data for each measurement interval was estimated. This was a process by which half-hourly soil CO₂ efflux data were reviewed to identify and subsequently remove non-representative measurements. For the quality assessment, first the root-mean-square-error (RMSE) of the linear fit of chamber headspace CO₂ mole fractions versus time was calculated. Next, by applying a threshold value of 1 for the ratio of RMSE to the slope of the linear fit of dSc/dt, all values above this threshold were considered unreliable and were removed following a method similar to Jassal et al. [7]. A visual inspection was also conducted to remove any remaining erroneous data. Negative fluxes were discarded as any vegetation, whether by moss or small plants, were removed on the soil surface within the chambers.

In addition to data removal through quality control, data loss also occurred from chamber malfunctions, calibrations, the lack of compressed air available, which controls the ability of the chamber dome to close, and winter snow and ice accumulation. In 2013, hydraulic problems arose in the two-way pneumatic cylinder, resulting in approximately 80%–90% of annual data loss and data removal for most of the chambers, all occurring during the growing season. Further mechanical problems in the two-way pneumatic cylinder occurred in 2014, with about 70%–90% of the data removed. The percentage of data retained decreased over the years due to a combination of equipment age, malfunctions and removal of the chamber system (November 20, 2011 and March 27, 2012) for forest harvesting operations.

2.5. Eddy Covariance Flux and Ancillary Data

Water, carbon and energy fluxes were measured using a closed-path eddy covariance system. Flux measurements were made at 20 Hz above the canopy on top of a scaffolding tower at a height of 28 m. Above canopy fluxes of carbon and water were calculated using a methodology outlined by Baldocchi et al. [34]:

$$F=\overline{w\prime }\overline{c\prime }$$

(3)

where F refers to the flux of water, carbon or energy, w’ is the covariance of vertical wind speed, measured by a sonic anemometer (CSAT3, Campbell Scientific Inc. (CSI)), and c’ is the covariance of CO₂ or water vapour concentrations, measured by an infrared gas analyzer (Li-7000, LI-COR Inc.). Half-hourly net ecosystem productivity (NEP) was calculated by adding CO₂ flux and the CO₂ storage, which is calculated as change in CO₂ concentration over time in the air column below the eddy covariance sensors [35]. Ecosystem respiration (R_E) was estimated using a non-linear logistic relationship between nocturnal half-hourly CO₂ fluxes during high turbulence conditions and soil temperature at 2 cm depth [28,29,35]. Gross ecosystem productivity (GEP) was determined by adding measured NEP to model daytime R_E. Soil temperature and soil moisture were measured at 5, 10, 20, 50, and 100 cm depths at two locations, near the chamber system location. For the analysis, soil temperature at a depth of 5 cm and soil moisture, averaged between 0 and 30 cm, were used. All flux and meteorological data were quality controlled and averaged at half-hourly intervals [35]. For further details on site’s eddy covariance tower instrumentation refer to Arain and Restrepo-Coupe [30] and Peichl et al. [28,29].

2.6. Statistical Analysis of Pre- and Post-Harvesting

To determine whether the shelterwood harvesting event can be detected within the variability of the data, a comparison was made for the pre- (2008–2011) and post-harvesting (2012–2014) periods. Half-hourly soil temperature and soil moisture values were averaged daily to allow for a comparison of adequate sized data points between the two periods [36,37]. As soil CO₂ efflux measurements had several half-hourly gaps after quality control, only available measurements taken at 11:00 for each day were compared, as the mean daily soil CO₂ efflux value was found to correspond best to values around 11:00 at this site, and in previous studies [38].

The soil temperature, soil moisture, soil CO₂ efflux, and eddy covariance flux data (all denoted by x’ in Equation 4) were normalized using the min-max normalization method where the values are rescaled to range between 0 and 1 using the following equation:

$$x\prime =\frac{x-\min \left(x\right)}{\max \left(x\right)-\min \left(x\right)}$$

(4)

where min (x) and max (x) are the minimum and maximum values of x for the complete dataset (2008–2014). This was done to normalize the variability between the various chambers and eddy covariance tower flux measurements, and allow for a more direct comparison between all. An F-test was then applied to assess whether the variances of the two periods were equal or not, followed by a t-test to determine if there was any significant difference (p-value < 0.05) in the means.

Furthermore, as these variables change seasonally, the relative contribution of component fluxes to the total soil CO₂ efflux were examined for different months of the year. Following similar work to Heinemeyer et al. [39], two-month averages of diurnal variability of soil temperature and soil CO₂ efflux rates, for each individual chamber, were constructed for March–April, July–August, August–September, and December–January. To check for changes, diurnal patterns from pre- and post-harvesting years were compared. Comparisons for the diurnal patterns from pre- and post-harvesting were made for selected two-month seasons in those years with a similar climate: when their two-month average temperature and total precipitation were similar.

We fitted a mixed effect model on our data, where natural logarithm of soil CO₂ efflux was modelled as a dependent variable of three fixed variables: soil temperature, soil moisture and shelterwood management (represented by the “thinning” categorical variable), and two random variables: “chamber number” and “day of year”. The random variables accounted for possible effects of pseudoreplication (“chamber number”) and time series correlation (“day of year”) in our dataset. The model was fitted as full (with all three fixed explanatory variables), then also without the “thinning” categorical variable, and finally without the “soil moisture” and “thinning” variables. Note that “soil temperature” and “soil moisture”, when present, were allowed to have interaction terms. Daily average values were used in these models. Akaike’s information criterion (AIC) was then used to compare the three models. This analysis was completed in R-software (R Core Team, Vienna, Austria) [40] and using the lme4 package [41].

Finally, please note that this study was not a truly replicated experiment and the results needed to be viewed and interpreted cautiously due to pseudoreplication [42].

3. Results

3.1. Seasonal Variations in Soil Temperature and Soil Moisture

During the study, the site experienced a large seasonal variation in temperature with air temperatures reaching above 30°C in the summer, and well below 0 °C in the winter. This variation was also reflected in the soil temperatures measured at 5 cm depth, but less pronounced. Daily soil temperatures fluctuated around 0 °C in the winter, increased by April, peaked in July or August, and then decreased again by October (Figure 1a). The mean soil temperature for the seven-year data period from 2008 to 2014 was 9.3 ± 7.5 °C, with the highest soil temperature of 23 °C recorded during the summer of 2013, and the lowest temperature of −2.9 °C recorded during the winter of 2012. However, in 2012, above-normal temperatures occurred in March; increasing soil temperatures well above 10 °C for approximately 20 days.

Daily soil moisture showed large variations during the seven-year study period, which is attributed to rainfall and/or snowmelt events, and influenced by the well-drained sandy soil at this site. Annually, soil moisture increased during the winter season, with maximum values observed after snowmelt (Figure 1b). Little variability in soil moisture was observed during the 2012 growing season, as very few rainfall events occurred, and this year was deemed a dry year [43]. The mean soil moisture over the seven-year study period was 0.12 ± 0.03 m³/m³, with the highest value of 0.25 m³/m³ recorded during the spring of 2013 due to a combination of snowmelt and precipitation. In each growing season, soil moisture levels in this well-drained soil declined to approximately 0.065 m³/m³.

Diurnal pattern of soil temperatures over the two-month period for each season of the study showed a typical decrease overnight, then an increase coinciding with an increase in daylight, with peak values in late -afternoon (Figure 2a–d). The diurnal pattern of mean soil moisture over the July to August and August to September months showed increased overnight values that decreased as the day progressed (Figure 2f,g). For the March–April and December–January months, mean diurnal values of soil moisture stayed relatively flat (Figure 2e,h). Soil moisture for the winter was higher than the growing season, due to the lack of evapotranspiration and an increase in precipitation events.

3.2. Seasonal Variability of Soil CO₂ Efflux

From 2008 to 2014, the soil CO₂ efflux rates showed seasonal patterns following seasonal soil temperature, where soil CO₂ effluxes were low during the winter and high during the summer. There was a large amount of variability in soil CO₂ effluxes during the growing season and much less during the winter months. Soil CO₂ efflux from the control chambers (Figure 1c,d) showed the highest daily average values at 3.1 ± 2.0, 2.5 ± 1.8, 3.2 ± 2.2, 3.1 ± 2.1, and 3.3 ± 2.4 μmol of C m² day⁻¹ for F_S1, F_S2, F_S3, F_S4, and F_S5, respectively. Both F_S1 and F_S2 showed high CO₂ effluxes in 2009, compared with the following years. Soil CO₂ effluxes for F_S4 were higher in 2010 and 2014 compared with the other control chambers. Equipment malfunctions limited the number of observations for 2013 and 2014, with chamber data only from F_S3 and F_S5 in 2013, and F_S4 in 2014, passing quality control. No observations during the 2011/12 winter season were recorded due to the removal of the automatic chamber system for shelterwood harvesting.

Before removal of the litter-layer, F_LR produced similar soil CO₂ effluxes to both F_S1 and F_S2 (Figure 1e). After litter-layer removal, soil CO₂ effluxes dropped and remained lower than the control chambers, producing daily average value of 2.2 ± 1.5 μmol of C m² day⁻¹ over six years.

The two chambers where soil was trenched to exclude tree roots also produced lower values of soil CO₂ efflux compared to the control chambers, with daily average values of 2.0 ± 1.5 and 2.2 ± 1.5 μmol of C m² day⁻¹ for F_RT1 and F_RT2, respectively. These averaged daily soil CO₂ efflux values were similar to those from F_LR for most years, except 2009 and 2014.

Mean averaged half-hourly soil CO₂ efflux rates, compared to soil temperature and soil moisture, showed diurnal fluctuations that mimic the diurnal pattern of soil temperature, while responding to increases in soil moisture (Figure 3). An example from June 2009 showed that soil CO₂ effluxes from the control (F_S), heterotrophic (F_RT) and mineral-soil horizon (F_LR) increased with an increase in both soil temperature and soil moisture. However, in this example, June 17 showed an instance where the mineral-soil component decreased at the onset of a rainfall event. Conversely, and less common, instances of a decrease in the mineral-soil component occurred in the following years.

In January and December 2013, large pulses of soil CO₂ effluxes were observed from F_LR, F_RT1, F_RT2, and F_S2 chambers, often exceeding values recorded during growing seasons (Figure 1c–e). Based on observed meteorological data, these soil CO₂ efflux pulses seem to have corresponded to an increase in soil temperature, from temperatures below freezing to a few degrees above 0 °C, and an increase in soil moisture (from 0.14 to 0.20 m³/m³), both of which may have stimulated microbial activity. These rare observations could also be due to the melting of ice or dense snow cover on top of the soil within the chambers, which was blocking the release of soil CO₂ accumulated in the near-surface soil layer.

Two-month diurnal average of soil CO₂ effluxes for each season showed differences in the range of magnitude. Although soil CO₂ efflux had more variability than soil temperature and soil moisture, for each season, the soil CO₂ efflux increased during the daylight, peaking around 21:00, about 3–5 hours after soil temperature peaks. Afterwards, emissions decreased overnight, with the lowest usually observed in the early morning hours. Soil CO₂ effluxes at approximately 11:00 were found to be close to the mean daily average for total soil CO₂ efflux and its components.

3.3. Comparison of Soil Temperature and Soil Moisture Pre- and Post-Harvesting

Mean values of soil temperature over the pre- and post-harvesting time periods were 9.4 ± 7.5 °C and 9.1 ± 7.5 °C, respectively, while the mean values of soil moisture were 0.11 ± 0.03 m³/m³ and 0.12 ± 0.03 m³/m³, respectively. A t-test with normalized daily data did not show a significant difference for both soil variables (

Table 1. Mean soil temperature, mean soil moisture and soil CO₂ efflux, from total soil CO₂ efflux (F_S1, F_S2, F_S3, F_S4 and F_S5), heterotrophic respiration (F_RT1 and F_RT2) and mineral-soil horizon respiration (F_LR), for pre- (2008–2011) and post-harvesting (2012–2014) periods, along with the percent difference. * denotes significant difference (α = 0.05) of the normalized daily means for soil temperature and soil moisture and for those soil CO₂ efflux measurements taken at 11:00 between pre- and post-harvesting periods.

Environmental Variable	Pre-Thinning (no. of observations)	Post-Thinning (no. of observations)	% Difference	t-statistic	p-value
Soil Temperature	9.4 ± 7.5 °C (n = 1274)	9.1 ± 7.5 °C (n = 1095)	−3%	−0.90	0.37
Soil Moisture	0.11 ± 0.03 m3/m3 (n = 1274)	0.12 ± 0.03 m3/m3 (n = 1095)	3%	2.48	0.01
	Mean Soil CO2 Efflux (μmol of C m2·s−1)
FS1	6.09 ± 3.66 (n = 805)	4.01 ± 3.23 (n = 226)	−34%	−4.69*	0.00
FS2	4.64 ± 3.20 (n = 919)	4.03 ± 3.51 (n = 261)	−13%	−3.48*	0.00
FS3	5.83 ± 3.90 (n = 654)	5.71 ± 4.03 (n = 253)	−2%	−1.35	0.18
FS4	5.91 ± 3.78 (n = 578)	4.90 ± 3.70 (n = 256)	−17%	−4.30*	0.00
FS5	6.73 ± 4.68 (n = 409)	4.68 ± 4.09 (n = 282)	−30%	−6.02*	0.00
FLR	4.06 ± 2.76 (n = 882)	3.30 ± 2.77 (n = 313)	−19%	69.2*	0.00
FRT1	3.84 ± 2.80 (n = 627)	2.58 ± 2.63 (n = 214)	−33%	−7.44*	0.00
FRT2	4.44 ± 2.79 (n = 283)	3.65 ± 2.84 (n = 191)	−18%	−4.81*	0.00

). While soil temperature did not change much between the pre- and post-harvesting periods, an increase of 3% in soil moisture was observed, even with the data from a 2012 drought year included in the post-harvesting period. A comparison of the median and range of pre- and post-harvesting data confirmed that soil temperature was similar and soil moisture showed a slight increase after harvesting (Figure 4). As soil temperature at this site was correlated to air temperature, the median and range of pre- and post-harvesting air temperature was also compared and determined that the measurements were similar between the two periods.

3.4. Comparison of Soil CO₂ Effluxes Pre- and Post-Harvesting

For March to April, representing spring, a comparison was made for years 2010 and 2012, which had similar observed climate. For the control (Figure 2i,m) and the litter-layer removed chambers (Figure 2q) minimal change was observed as the magnitude and diurnal pattern were similar. One trenched chamber (the other trenched chamber not being installed until May 2010), showed an increase with the average diurnal soil CO₂ efflux increasing from 0.86 ± 0.07 to 1.09 ± 0.11 μmol of C m² s⁻¹ (Figure 2u).

In July to August, a comparison of 2011 to 2012 showed a decrease in the average diurnal soil CO₂ efflux for the control chambers, but an increase in the standard deviation (Figure 2j,n). The largest change was measured for F_S5, with a decrease in the average by 1.42 μmol of C m² s⁻¹ and an increase in the standard deviation of 0.13 μmol of C m² s⁻¹. F_LR had an increase in both the mean and standard deviation, (3.09 ± 0.07 μmol of C m² s⁻¹ to 3.45 ± 0.18 μmol of C m² s⁻¹; Figure 2r), while both F_RT1 and F_RT2 decreased in the average diurnal soil CO₂ efflux, and an increase in the standard deviation (3.07 ± 0.09 μmol of C m² s⁻¹ to 2.72 ± 0.23 μmol of C m² s⁻¹ and 4.02 ± 0.18 μmol of C m² s⁻¹ to 2.93 ± 0.21 μmol of C m² s⁻¹ for F_RT1 and F_RT2, respectively; Figure 2v).

Comparing the two-month diurnal average of August to September, between 2009 and 2012, the four available control chambers all had a decrease in soil CO₂ effluxes (Figure 2k,o). The largest difference was seen in F_S4, with a daily average increase of 1.24 μmol of C m² s⁻¹ and an increase in the daily average standard deviation of 0.02 μmol of C m² s⁻¹. F_LR had a slight increase in both the daily average and standard deviation post-harvesting (Figure 2s). The daily average soil CO₂ efflux for F_RT1 decreased by 0.29 μmol of C m² s⁻¹ between 2009 and 2012, while the standard deviation increased by 0.10 μmol of C m² s⁻¹ (Figure 2w). Despite similar soil temperature and soil moisture conditions, soil CO₂ effluxes decreased post-harvesting during the growing season, except for F_LR, which showed an increase.

For the December–January months, a comparison of 2010–2013 two-month diurnal averages showed an increase in soil CO₂ effluxes for all five control chambers, with the largest average increase of 1.25 ± 0.33 μmol of C m² s⁻¹ seen in F_S2 (Figure 2l,p). F_LR had an increase in the daily average soil CO₂ efflux, with an average of 0.84 μmol of C m² s⁻¹ in 2013, compared to an average of 0.31 μmol of C m² s⁻¹ in 2010 (Figure 2t). Both F_RT1 and F_RT2 had a higher soil CO₂ efflux in 2013, compared with 2010 (Figure 2x). F_RT2 increased from an average of 0.57 ± 0.06 μmol of C m² s⁻¹ in 2010, to an average of 0.82 ± 0.12 μmol of C m² s⁻¹ in 2013. F_RT1 had the largest difference between the two years, with an average of 0.37 ± 0.05 μmol of C m² s⁻¹ in 2010, to an average of 2.58 ± 0.52 μmol of C m² s⁻¹ in 2013. The difference was attributed to a large amount of CO₂ released at the end of December 2013, which persisted for approximately five days and was comparable to values observed during the growing seasons. Note that the number of observations recorded during the winter months were much less than during the growing season, so these recordings in December 2013 have a stronger influence on the winter average.

Comparison of the soil CO₂ effluxes measured pre- and post-harvesting showed that there was a reduction from each of the chambers after the shelterwood harvest event (Table 1, Figure 5). Reduction in total soil CO₂ efflux ranged from 1 to 37%, with the highest difference from F_S1, and minimal change in F_S3 measurements. The chamber with the litter-layer removed (F_LR) had a 2% reduction, while both measurements from the trenched chambers (F_RT1 and F_RT2) decreased by approximately 27%. A significant difference for all but two of the chambers were shown for the normalized data (Table 1). Both F_S3 and F_LR showed no significant difference between pre- and post-harvesting periods (Table 1). Overall, we observed a decrease in emissions after the harvesting, with a reduction observed in control chambers and those with severed roots (Figure 5). In Figure 5, with the year 2013, soil CO₂ efflux values appear to be much higher, however, this was due to large data gaps where most of the data in that year came from the growing season and was biased to higher values.

Results from mixed effects model analysis corroborated above analysis in that soil CO₂ efflux was positively related to both soil temperature and soil moisture. Including the “thinning” variable in the model improved the model fit, suggesting that there was a statistical difference in fluxes and the relationship between pre- and post-harvesting soil CO₂ effluxes. Furthermore, differences in the estimated model coefficients also suggest differences in the response of soil CO₂ emissions to changes in soil temperature and soil moisture, with the largest variability detected in the “soil moisture” variable (see Appendix A for details).

3.5. Comparison of Soil CO₂ Effluxes to Ecosystem Respiration

To assess the impact of shelterwood harvesting on overall carbon budget of the forest, an analysis of net ecosystem productivity (NEP) using the eddy covariance flux data from the nearby flux tower was conducted following methods outlined by Trant (2014). Using data from 2008 to 2014, as compared to the 2003–2012 for the analysis by Trant (2014), showed that thinning, following harvesting, did not significantly impact the response of NEP (Figure 6). Post-harvesting NEP fluxes were within the range of interannual variability over the study period (2008–2014). Mean pre-harvesting NEP was 237 ± 81 g C m⁻² year⁻¹ as compared to post-harvesting NEP of 185 ± 97 g C m⁻² year⁻¹. Since NEP is a difference of two large fluxes, Gross Ecosystem Productivity (GEP) and ecosystem respiration (R_E), the differences between these components were also considered. Mean annual post-harvesting (2012 to 2014) GEP was 1519 ± 78 g C m⁻² year⁻¹ as compared to pre-harvesting (2008 to 2011) GEP of 1479 ± 74 g C m⁻² year⁻¹. Median and range of R_E fluxes pre-harvesting (1233 ± 152 g C m⁻² year⁻¹) and post-harvesting (1326 ± 55 g C m⁻² year⁻¹) were also similar (Figure 4d). Based on the results, there was no significant impact on the flux measurements following harvesting, however, both the GEP and R_E fluxes increased after the shelterwood harvesting, while NEP fluxes decreased (less overall carbon sink). The observed reduction in soil CO₂ efflux rates did not reflect a substantial reduction in the ecosystem respiration at our site.

4. Discussion

Soil temperature and soil moisture have been shown to be a major influence on the soil CO₂ efflux dynamics [44,45,46]. Results from this study were in agreement with these studies as soil CO₂ efflux followed the seasonal pattern of soil temperature, with the lowest CO₂ effluxes observed during the winter, and highest during the growing season. Observed half-hourly soil CO₂ efflux rates showed diurnal fluctuations that mimic the diurnal pattern of soil temperature, while responding to increases in soil moisture. Pulses of soil CO₂ efflux were observed during the winter of 2013/2014, with measured values similar to those observed during the growing season. These pulses seem to correspond to an increase in soil temperature and melting of the snow cover. The formation of dense snow and/or ice lenses in a snowpack can impede the upward diffusion of soil CO₂ efflux [45,47]. Monson et al. [45] observed a similar increase in soil CO₂ efflux due to meltwater, from the above snow cover layer, and an increase in soil temperature, which increased the activity of soil microbes during spring for a subalpine ecosystem. DeForest et al. [48] observed pulses in soil respiration following full leaf senesce that corresponded with an input of fresh litter decay. This pattern did not occur annually, and may be attributed to a warmer season. In relation to the possibility of warming soils, these pulses could have implications on the cycling of nutrients released from the decomposition of the litter-layer [48]. Improvement of collecting robust and accurate soil CO₂ efflux data, during the winter season, would help in the understanding the response of soil CO₂ efflux to winter environmental factors, such as snow, ice accumulation, and freeze/thaw events.

Forest management practices, such as shelterwood harvesting, produces disturbances that alter the physical structure, microclimate, physiochemical properties, and biological activity of forest soils and vegetation, affecting soil respiration [12,46,49,50,51,52]. Mixed results from the literature [23,49,50,51,52,53,54,55,56,57] highlights the need for comprehensive and long-term studies, to explore the impact of shelterwood harvesting, which is being increasingly adopted to manage forest ecosystems, on soil CO₂ efflux and the overall carbon budget. This study is an important step forward in this direction, where soil CO₂ efflux measurements were conducted, over the seven-year period (2008–2014), in a managed temperate forest.

Comparison of soil temperature and soil moisture, between pre- and post-harvesting periods, showed no significant difference, although there was a small (0.3 ± 7.5 °C) decrease in soil temperature, while soil moisture showed a small increase (0.01 ± 0.03 m³/m³) following post-harvesting. The slight increase in soil moisture post-harvesting may reflect the reduced root density and subsequent reduced water uptake. It may also reflect increased insulation from the litter left behind by the harvesting event, which could help reduce soil evaporation and hence water loss. Despite similar soil microclimates, total soil CO₂ effluxes were significantly reduced by up to 37%. Soil CO₂ effluxes from heterotrophic sources were significantly reduced post-harvesting by approximately 27%, while no significant difference in the mineral-soil horizon sources were measured. Although root respiration may decrease with a decrease in root biomass post-harvesting, root respiration from the remaining trees may also increase as a result of increased photosynthetic rate and growth of new roots due to a decrease in competition for water and nutrients among remaining trees [49]. In the thinned stands, dead roots from harvested trees can also contribute more soil CO₂ efflux after the treatment [49]. In contrast, Shabaga et al. [50] found that partial harvesting, in northern mixed deciduous forests, increased soil respiration due to an increase in soil temperatures, and enhanced decomposition rates of harvested residue, such as leaves/needles and woody debris left on ground, and the decomposition of dead roots. They also reported that, after thinning, forest understory regrowth contributed to a reduction in soil temperatures, but influenced patterns of soil CO₂ efflux by increasing summer respiration and Q₁₀ values, likely due to increased root respiration. Wang et al. [51] found no change in soil respiration between control and thinned larch plantation sites, however, the distribution in different components were altered, as both heterotrophic and autotrophic respiration were higher at the thinned site, whereas litter decomposition was lower. Wang et al. [51] attributed the higher levels of heterotrophic respiration from the mineral-soil to a change in soil temperature, something not found in this study.

As suggested by Sullivan et al. [52], a possible reason for the lower soil CO₂ efflux can be explained by a decline in live tree root respiration being greater than a possible increase in soil CO₂ efflux, resulting from an increase in root production by surviving trees, and a possible increase in heterotrophic respiration associated with the decomposition of newly dead roots and understory plant litter. The shelterwood harvest may have altered the soil environment, by changing amounts and sources of belowground carbon for microbial metabolism, thus reducing soil CO₂ efflux [49]. Since this study only examined the short-term effects of shelterwood harvesting, the increase in quantity of high-quality carbon available for microbial decomposition, usually associated with this harvesting, may not yet be available for decomposition [52]. Tang et al. [49] found that soil respiration decreased by 13% post-thinning at their site, for data normalized with soil temperature and soil moisture. They note that due to natural climate variability and the effect of thinning on soil temperature and moisture, no significant trends were found in actual non-normalized soil respiration post-thinning. They suggested that the net effects of thinning on soil respiration are the result of complex and numerous interacting factors. Tang et al. [49] argued that apart from thinning impacts, that changes soil temperature, soil moisture and biotic factors, soil CO₂ efflux in a thinned forest may also be influenced by climate variability (such as drought in 2012 for this site).

An analysis of the eddy covariance flux data, specifically net ecosystem productivity (NEP), gross ecosystem productivity (GEP), and ecosystem respiration (R_E) from the nearby flux tower, showed that harvesting had no significant impact. Since there was a reduction in soil CO₂ efflux post-harvesting while both GEP and R_E increased, it appears that the surge in R_E was primarily caused by the increase in autotrophic respiration due to the enhancement in forest growth. Indeed, personal communications with above-ground growth researchers at our site (personal communication, Shawn McKenzie) have shown that post-harvesting, above-ground tree biomass increased with significant increases in diameter at breast-height and basal area detected. This suggests that despite the decrease in soil CO₂ emissions, an increase in above-ground R_E may have compensated for any reductions in soil CO₂ efflux. Finally, an increase in R_E can also be attributed to a disparity between flux footprints for the eddy covariance fluxes, which encompass the whole stand, as compared to a localized area where the soil CO₂ efflux chambers were installed. Soil CO₂ efflux chambers often overestimate soil respiration as compared to R_E measured by the eddy covariance system [58]. This study highlights the complexities and challenges associated with soil CO₂ efflux studies and the need for comprehensive and long-term field studies to better understand the response of these fluxes to a change in the future climate.

5. Conclusions

Silvicultural practices, such as shelterwood harvesting, can have a profound impact on forest carbon dynamics by changing soil temperature, and soil water content, soil organic matter, root biomass, and microbial activity. These silvicultural practices are being adopted in North America and across the world, not only to manage forests, particularly planted stands, but also to enhance their carbon sequestration capabilities, and/or to restore them to their previous native species, such as mixed woodland in the Great Lakes region. By comparing soil CO₂ effluxes between pre- and post-harvesting years, this study showed that soil temperature and soil moisture did not show any significant changes post-harvesting compared to pre-harvesting. Despite similar climate conditions, soil CO₂ effluxes post-harvesting were lower than pre-harvesting. Further collection of soil CO₂ efflux measurements as the forest re-establishes are needed to determine the full effect of the shelterwood system on soil respiration and, particularly, the impact to heterotrophic and autotrophic respiration following post-harvesting. Results presented in this study are useful for forest management practitioners, specifically those focused on carbon sequestration and forest conservation. The knowledge gained from this research will help in developing polices for better management of forest ecosystems in Eastern Canada, in addition to assisting in planning and developing realistic strategies to offset fossil fuel CO₂ emissions to improve environmental quality.