Harvesting Intensity and Aridity Are More Important Than Climate Change in Affecting Future Carbon Stocks of Douglas-Fir Forests

Robinson, Alyssa J.; Defrenne, Camille E.; Roach, W. Jean; Dymond, Caren C.; Pickles, Brian J.; Simard, Suzanne W.

doi:10.3389/ffgc.2022.934067

ORIGINAL RESEARCH article

Front. For. Glob. Change, 11 August 2022

Sec. Forest Management

Volume 5 - 2022 | https://doi.org/10.3389/ffgc.2022.934067

This article is part of the Research Topic Retaining Ecosystem Legacies in Forest Management: Effects on Forest Structures and Functions View all 9 articles

Harvesting Intensity and Aridity Are More Important Than Climate Change in Affecting Future Carbon Stocks of Douglas-Fir Forests

$\r\nAlyssa J. Robinson*$ Alyssa J. Robinson^1*

¹Department of Forest and Conservation Sciences, Faculty of Forestry, The University of British Columbia, Vancouver, BC, Canada
²College of Forest Resources and Environmental Science, Michigan Technological University, Houghton, MI, United States
³Forest Carbon and Climate Services Branch, Government of British Columbia, Victoria, BC, Canada
⁴School of Biological Sciences, University of Reading, Reading, United Kingdom

Improved forest management may offer climate mitigation needed to hold warming to below 2°C. However, uncertainties persist about the effects of harvesting intensity on forest carbon sequestration, especially when considering interactions with regional climate and climate change. Here, we investigated the combined effects of harvesting intensity, climatic aridity, and climate change on carbon stocks in Douglas-fir [Pseudotsuga menziesii Mirb. (Franco)] stands. We used the Carbon Budget Model of the Canadian Forest Sector to simulate the harvest and regrowth of seven Douglas-fir stand types covering a 900 km-long climate gradient across British Columbia, Canada. In particular, we simulated stand growth under three regimes (+17%, −17% and historical growth increment) and used three temperature regimes [historical, representative concentration pathways (RCP) 2.6 and RCP 8.5]. Increasing harvesting intensity led to significant losses in total ecosystem carbon stocks 50 years post-harvest. Specifically, forests that underwent clearcutting were projected to stock about 36% less carbon by 2,069 than forests that were left untouched. Belowground carbon stocks 50 years into the future were less sensitive to harvesting intensity than aboveground carbon stocks and carbon losses were greater in arid interior Douglas-fir forests than in humid, more productive forests. In addition, growth multipliers and decay due to the RCP’s had little effect on total ecosystem carbon, but aboveground carbon declined by 7% (95% confidence interval [−10.98, −1.81]) in the high emissions (RCP8.5) scenario. We call attention to the implementation of low intensity harvesting systems to preserve aboveground forest carbon stocks until we have a more complete understanding of the impacts of climate change on British Columbia’s forests.

Introduction

Forests have natural potential in mitigating climate change (Chen et al., 2000; IPCC, 2006; Pan et al., 2011; Shukla et al., 2019) as trees and understory plants capture CO₂ from the atmosphere through photosynthetic processes and store carbon (C) both above- and belowground (Streck and Scholz, 2006). Globally, forests have sequestered the equivalent of 20% of fossil fuel emissions over the past three decades (Le Quéré et al., 2018) and they store approximately 70% of aboveground plant C and 70% of terrestrial soil C (Sedjo, 1993; Dixon et al., 1994; Goodale et al., 2002; Bar-On et al., 2018).

The impact of climate change on forest C stocks is not fully understood and projections vary widely among climate change scenarios. These scenarios (i.e., Representative Concentration Pathways, RCP’s) are used by the Intergovernmental Panel on Climate Change (IPCC) to describe potential climate futures. Each RCP represents the degree to which the atmosphere traps C and is based on the volume of C released by human activity and natural processes such as wildfires and volcanic activity. Four RCP scenarios have been the focus of much research over the last 8 years (IPCC, 2014): (i) the RCP 2.6 predicts a global temperature rise of ≤1.5°C (compared to the pre-industrial baseline) by 2085–2100 and requires the implementation of technological and behavioral climate change mitigation strategies, (ii) the RCP 4.5 predicts that a global temperature rise of >2.0°C by 2085–2100 is likely, (iii) the RCP 6.0 predicts a global temperature rise of 2.0–4.0°C by 2085–2100 and (iv) the RCP 8.5 predicts a global temperature rise of 4.0°C by 2085–2100 and involves no climate change mitigation efforts. Note that more recently, the IPCC has introduced new climate change scenarios called Shared Socioeconomic Pathways (SSP’s) which project CO₂ concentrations based on implemented climate policies (IPCC, 2021). Globally, forest C sequestration potential may increase as a result of enhanced vegetative Gross Primary Productivity (GPP) due to CO₂ fertilization (Wenzel et al., 2016) and increased nitrogen deposition (Magnani et al., 2007; de Vries et al., 2008; Sutton et al., 2008), and as a result of the expansion of productive forested zones (Metsaranta et al., 2011) such as the Interior Cedar Hemlock and Interior Douglas-Fir zones of British Columbia (BC), Canada (Wang et al., 2012). However, a changing climate may also lead to conditions favorable to insect (Carroll et al., 2004; Kurz et al., 2008a) and pathogen outbreaks (Woods et al., 2005; Sturrock et al., 2011; Weed et al., 2013), which may lead to increased forest respiration, in turn, negatively impacting forest C stocks (Kurz et al., 2008b; Dymond et al., 2010). Furthermore, as many regions of the world become more drought prone, trees are expected to experience more rapid early growth resulting in shorter lifespans, potentially creating net carbon sources in the long term (Stovall et al., 2019; Brienen et al., 2020; Appiah Mensah et al., 2021; Socha et al., 2021).

The world’s forests may assist in providing one quarter of the climate mitigation needed to stabilize warming to below 2°C by 2030 (Griscom et al., 2017), and it is therefore crucial to understand how to promote forest C sequestration in a warmer future (Temperli et al., 2012; Hof et al., 2017; Dobor et al., 2020; Dymond et al., 2020). Generally, increasing forest C sequestration can be accomplished following three approaches: (i) improved forest stewardship, (ii) reduced deforestation through increased conservation and preservation, and (iii) increased reforestation and afforestation (Field and Mach, 2017). Since the demand for forest products is increasing steadily (as seen by a 12% increase in global roundwood production between the years 2000 and 2020; Food and Agriculture Organization [FAO], 2020), the need for improved forest stewardship is becoming evident.

Long-term research shows that harvesting typically produces forests that do not fully regain their total ecosystem C stock compared to their original state (Chatterjee et al., 2009; Keith et al., 2014). For example, old growth forests of southern Australia (i.e., 120–250 years old, but which may contain trees aged 400–500 years old; Wood et al., 2010) have been shown to store 55% more carbon above- and belowground than forests logged every 50 years (Keith et al., 2014). In addition, forests that are partially harvested, mixed, or uneven-aged have greater C storage capabilities compared to clearcut forests or uniform plantations (Zhou et al., 2013; Forrester and Bauhus, 2016). In a meta-analysis of 81 studies, Zhou et al. (2013) showed that lower harvesting intensities, such as large patch retention, removed approximately a third of aboveground biomass C, whereas higher harvesting intensities such as single tree retention removed close to half of aboveground biomass C compared to unharvested forests. Lower intensity harvesting has additional co-benefits to biodiversity by preserving individuals that effectively maintain biological legacies and provide habitat for resident species (Franklin and Donato, 2020). The effects of harvesting on belowground C stocks are not yet fully understood but some studies have reported no significant change with harvesting (Hume et al., 2018). Yet, the authors showed that forest floor C stocks of stands harvested at higher intensity were smaller than those of stands harvested at lower intensity. Additionally, regional climate has been shown to interact with harvesting intensity to affect C stocks. For example, in a study investigating 1-year post harvest effects of regional climate and harvesting intensity, Simard et al. (2020) found that clearcutting reduced total ecosystem C stocks of arid Douglas-fir forests by 60%, whereas it reduced C stocks of humid forests by 50% (with untouched forest as baseline).

Forest C models have been previously used to assess the impacts of harvesting, climate change, and/or regional climate on forest C stocks. For example, using a C flux model (Biome-BGC) projecting total ecosystem C from 2006 to 2089, Boisvenue and Running (2010) showed that a drier climate projection (precipitation projections from the GFDL-CM2.0 model) led to a decrease in C accumulation whereas a wetter projection [precipitation projections from the CGM3.1 (T63) model] led to an increase in C accumulation in high-elevation forests of the western United States. Using LANDIS-II (a mechanistic model), Dymond et al. (2016) highlighted that within temperate coniferous forests of BC, which are expected to be warmer in the future, projected regional climate change resulted in increased C sinks, but only within previously cold and wet regions. In addition, a simulation using the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3; Kurz et al., 2009) looking at the effect of harvesting intensity on forest C stocks showed that partial harvesting resulted in increased total ecosystem C over a 240-year simulation, whereas clearcutting resulted in decreased C in red spruce forests of Canada (Taylor et al., 2008). Altogether, climate change, harvesting intensity, and local climate all have an impact on forest C stocks, however, their combined effect is unclear due to multiple interactive effects, and hence requires further study.

The province of BC covers 95 million ha with multiple, high relief mountain ranges and variable ecological conditions. Due to this variety, the impact of climate change is expected to be diverse across BC, yet, overall, summers in southern BC are expected to experience declines in precipitation of up to 50% and increasing temperatures between +2.5 and 7.5°C (Spittlehouse, 2008). Increased aridity could result in increases in the frequency and intensity of wildfires, droughts, and insect outbreaks (Haughian et al., 2012; Spittlehouse and Dymond, 2022). Conversely, central and coastal BC are expected to experience increased winter precipitation of up to 40% (Spittlehouse, 2008). Such an increase in precipitation is conducive to fungal tree pathogens, flooding, and slope failures (Haughian et al., 2012). Climate change in BC will not only impact the disturbance regimes of landscapes, but risks change to individual tree growth and mortality as well. For example, Douglas fir [Pseudotsuga menziesii Mirb. (Franco)], a common softwood species in BC, could experience growth changes of ±17% in site index (Coops et al., 2010; Weiskittel et al., 2012).

Here, we investigated the combined effect of harvesting intensity, climatic aridity, and climate change on C stocks in 2069 in Douglas-fir stands of BC. We used CBM-CFS3 to simulate the harvest and regrowth of seven Douglas-fir stand types covering a 900 km-long climate gradient across BC, Canada; these stands were previously described by Roach et al. (2021) and empirical 1-year harvesting impacts by Simard et al. (2020). Growth was simulated under three regimes: +17%, −17% and historical growth increment over 50 years. Decay is temperature dependent and therefore we instituted three temperature regimes: historical, RCP 2.6, and RCP 8.5. We analyzed the 50-year C stocks in three summary pools: total ecosystem, aboveground, and belowground. We hypothesized that: (i) increasing harvesting intensity will result in lower ecosystem C stocks in 2069, (ii) increased temperatures will result in lower belowground C stocks across all sites regardless of harvesting intensity, (iii) increasing or decreasing growth rates will result in higher or lower aboveground carbon stocks, respectively, and (iv) more arid sites will store less C.

Materials and Methods

Study Areas

This study took place on seven Douglas-fir dominated locations covering a 900 km-long climate gradient in BC (Table 1 and Figure 1). Six locations [John Prince Research Forest (JPRF), Venables, Alex Fraser Research Forest (AFRF), Jaffray, Redfish Creek, and Two-Bit Creek] were in the Montane Cordillera ecozone of interior BC (54.65°N, 124.43°W to 49.21°N, 115.37°W) and the last location, Malcolm Knapp Research Forest (MKRF), was in the Pacific Maritime ecoregion of BC near Maple Ridge (49.32°N, 122.54°W). The locations span the following four biogeoclimatic ecosystem classification (BEC) zones: Interior Douglas-fir (IDF), Interior Cedar-Hemlock (ICH), Coastal Western Hemlock (CWH), and Sub-Boreal Spruce (SBS). For each location, we obtained climate normal data (1981–2010 averages) from ClimateNA (v 5.50; Wang et al., 2016) based on latitude, longitude, and elevation and characterized the regional climate using mean annual temperature (MAT), mean annual precipitation (MAP) and an index of climatic aridity; annual heat: moisture index (AHM; [(MAT + 10)/(MAP/1,000)],°C mm^–1).

TABLE 1

Table 1. Site characteristics of the seven forest locations across British Columbia, Canada.

FIGURE 1

Figure 1. Map of the seven research locations across British Columbia, Canada used in this study.

Each location encompassed three 20-ha replicate sites which were divided into five 4-ha treatment blocks with harvesting intensity treatments randomly assigned to the blocks, for a total of 15 blocks per location (randomized block design; see Supplementary Table 1). Exceptions included Two-bit which only had two replications (10 blocks) for all the treatments and Venables, which had two replications for the uncut control for a total of 14 blocks. The harvesting treatments included, in decreasing order of intensity: (i) clearcut (no tree retention), (ii) single tree retention (retention of 25 large Douglas-fir stems ha^–1 at 25 m spacing), (iii) small patch retention (30 m wide unconnected patches containing 30% of the stand area), (iv) large patch retention (retention of 60% of the stand area, thinned from below), and (v) no harvest (control). Forests were harvested in the winter of 2018/19 and data (tree species, size, age, and health) were collected in the summer prior to harvesting according to the Canadian National Forest Inventory protocols (see Measurement and Sampling Methods in Simard et al., 2020).

Modeling Approach

Variable Density Yield Projection

We used the Variable Density Yield Projection (VDYP) program (ver. 7; BC MFLNRO, 2019) to develop timber growth and yield curves for our stands that were eventually used in CBM-CFS3. This step was necessary as it created accurate growth and yield curves that properly reflect the field data. Since the VDYP program reflects the growth of unmanaged stands, the model is unable to provide results on the responses of tree growth to forest management, changes in species composition, or succession (BC MFLNRO, 2019). Although this may be considered a limitation, we applied management schemes in CBM-CFS3 to conduct our investigation.

VDYP functions by inputting stand data, including BEC zone, ecozone, tree species presence and percent cover, stand age and basal area, tree density as well as potential height. To determine potential height from field data, dominant and codominant trees without signs of disturbance (i.e., intact stems, no signs of insect/pathogen disturbance, etc.) were selected and their heights averaged (BC MFLNRO, 2019).

Data collected from the seven forest locations were used to parameterize VDYP for each replicate site and to create one growth and yield table per block (reflecting the growth and yield of all tree species present within a block). Stockable area was set to 95%, additional attributes were left as pre-computed values and the minimum diameter at breast height was set to 12.5 cm for pine, and 17.5 cm for all other species. The model was run from age 0 to age 250 in increments of 5 years.

Carbon Budget Model of the Canadian Forest Sector

CBM-CFS3 was developed in 1989 to simulate and track C dynamics at the levels of the stand and landscape. The model is used for reporting to the IPCC (Kurz et al., 2008a) and for investigating the effects of disturbance and management on C dynamics (Kull et al., 2019). For more details on CBM-CFS3, see Kurz and Apps (1999) and Kull et al. (2019). CBM-CFS3 has a few limitations, some of which are valuable to understand regarding our study. The model is unable to represent the processes which cause climate change effects on forest productivity. Specifically, for decay processes only, CBM-CFS3 only reflects the effects of changes in temperature but not precipitation. Hence, the effects of a changing climate will not be detected in live tree C pools, but will be detected in dead organic matter pools (Kull et al., 2019). Consequently, we have applied a growth multiplier to reflect the potential increase or decrease in live tree growth due to climate change. The multiplier (±17% in site index by age 50; Coops et al., 2010; Weiskittel et al., 2012) was applied gradually over the 50 years by slowly increasing the percent increase or decrease in growth (Table 2).

TABLE 2

Table 2. Growth multipliers applied in CBM-CFS3 where the first row represents gradual increased growth to reach +17% site index by year 50 and the bottom row represents gradual decreased growth to reach −17% site index by year 50.

We used the stand level project creator to run a simulation for harvesting type at each replicate site, and one stand was created for each individual block. Administrative boundaries and ecological boundaries were set to BC and Montane Cordillera, respectively, for all sites, except MKRF, which was set to BC and Pacific Maritime. Species composition was set to the leading species, usually Douglas-fir, and growth curves for each of the tree species were drawn from VDYP. We only used a single species since the growth curve produced by VDYP already factored in the presence of multiple species at a single site. In addition, stand age in 2019 was based on field data (e.g., MKRF was 68 years-old at the time of harvesting), stand area was set to 1 ha and historic disturbance to wildfire. For future disturbance, harvesting type specific to each block was implemented in the first year of the simulation, and planting of Douglas-fir was implemented in year two. To make use of the growth multiplier to reflect increased or decreased growth of forests due to climate change, a disturbance type with no impact on C levels had to be used. We implemented natural succession as this disturbance type. Every 5 years, a natural succession disturbance event was applied to the remaining stand and to the planted Douglas-fir (example growth curves can be found in Supplementary Figure 1).

The climate data editor was used to reflect the specific climate of the given stand’s location being modeled. For the baseline condition, time step zero (which represents past climate) and time step one (which represents future climate) temperatures were both set to the historical 50-year climate averages from Climate NA (Wang et al., 2016; See Supplementary Table 2 for details). In scenarios where alternative temperature regimes were implemented, time step 0 was unchanged, but time step one was set to the 50-year average temperature based on results from Climate BC v6.21 CanESM2 projection scheme from 2019 to 2069 (Wang et al., 2016). Simulations were run for 50 years, and results were collected using the results explorer.

Data Analyses

Statistical analyses were conducted in R version 4.1.1 (R Core Team, 2021) and results were considered statistically significant at P < 0.05. We used linear mixed-effects models to test the interacting effects of harvesting intensity, climatic aridity and climate change scenarios on projected Douglas-fir forest C stocks while accounting for the nested structure of the data (sites nested within locations). Aridity was a continuous variable, harvesting intensity was a categorical variable with five levels: control (baseline), clearcut, single tree retention, small patch retention and large patch retention, and climate change scenarios was a categorical variable with three levels: no-emission baseline, low emissions RCP 2.6 and high emissions RCP 8.5. We fitted a model for each C pool (ecosystem, above- and belowground) and growth multiplier (none, +17%, −17%) separately. Harvesting intensity, aridity and climate change scenarios were added as interacting fixed effects, stand age was added as a covariate and site (three replicate sites per location) and location (seven forest locations) were added as nested random effects. Aridity was scaled and centered using the scale() function to remove the correlation between main effects and their interactions, which leads to inflated variance inflation factor values (Schielzeth, 2010; Francoeur, 2013; Harrison et al., 2018; scaled and centered estimates can be found in Supplementary Table 3, summary statistics can be found in Supplementary Table 4).

We fitted models with lmer() (package lme4; Bates et al., 2014) and used step() (package lmerTest; Kuznetsova et al., 2017) to perform a backward elimination of non-significant effects. The best-fit models were as follows:

(i) Ecosystem C stocks ∼ harvesting intensity × aridity + stand age + (1| location/site)

(ii) Aboveground C stocks ∼ harvesting intensity × aridity + climate change scenarios + stand age + (1| location/site)

(iii) Belowground C stocks ∼ harvesting intensity × aridity + stand age + (1| location/site)

We checked model adequacy with check_model() (package performance; Lüdecke et al., 2020) and response variables were log-transformed when needed to meet model assumptions. We checked model fit using r2_nakagawa() (package performance), tested model significance with a likelihood ratio test [anova(), package stats], significance of fixed effects with a type II Wald χ² test [Anova(), package car; Fox and Weisberg, 2019] and plotted standardized coefficients using plot_model() (package sjPlot; Lüdecke, 2021).

Results

Increased harvesting intensity significantly reduced C stocks of Douglas-fir forests 50 years post-harvest (Figures 2A,B and Supplementary Table 5). In 2069, forests that underwent clearcutting (highest harvesting intensity) in 2019 are projected to stock about 36% less ecosystem C than forests that were left untouched, whereas forests where small and large tree patches (lowest harvesting intensities) were retained are projected to stock about 30 and 28% less C than uncut forests, respectively. The net C losses due to the single tree treatment were smaller than expected (5% less than the control) because C stocks in humid forests where single trees were retained were almost as large as those in the humid forests left untouched. Lastly, consistent among the three C pools, younger forests stored more C than older forests (which stored 1% less C per every year older).

FIGURE 2

Figure 2. Response of model projected Douglas-fir forest carbon stocks (2069) to harvesting intensity, climatic aridity (1981–2010 averages), and climate change (low-emissions RCP 2.6 and high-emissions RCP 8.5). Carbon stocks were projected using the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3; Kurz et al., 2009). Panels (A,C,E) present standardized coefficients of linear mixed-effects models for ecosystem, above- and belowground carbon stocks, respectively. Standardized coefficients illustrate the effect of the predictors on the response variable in terms of their standardized effect sizes. Circles indicate average estimates and lines indicate 95% confidence intervals (see Table 3). For harvesting treatments, each coefficient (green color) represents the difference between the uncut control and a given treatment and for climate change, each coefficient represents the difference between no emission and a given scenario. ***P < 0.001; **P < 0.01; *P < 0.05. Panels (B,D,F) present the response of ecosystem, above- and belowground carbon stocks to harvesting intensity, climatic aridity, and climate change, respectively. The aridity index was calculated as [mean annual temperature + 10]/[mean annual precipitation/1,000].

Total ecosystem C stocks 50 years into the future also decreased with regional climatic aridity, specifically, for every unit (°C mm^–1) increase in annual heat:moisture index (i.e., increase in aridity), projected C stocks decreased by about 3% across all treatments. However, harvesting intensity interacted with climatic aridity such that C losses due to increasing aridity were higher in forests where single trees were retained and lower in forests where large tree patches were retained. Note that since changing decay rates with temperature had no significant effect on projected C stocks across the harvesting treatments, it was not included in the best-fit model (Table 3). In addition, the use of the two growth multipliers did not significantly change the magnitude of the effects of harvesting intensity, aridity, and climate change on C stocks (total, above and belowground; Table 3).

TABLE 3

Table 3. Effect of harvesting intensity, climatic aridity (1981–2010 averages) and climate change (low-emissions RCP 2.6 and high-emissions RCP 8.5) on model projected Douglas-fir forest carbon stocks (2019–2069; log₁₀-transformed; no growth multiplier; n = 3 for all locations, except for Two-bit where n = 2 and Venables where n = 2 for the uncut control).

Projected aboveground C stocks of Douglas-fir forests decreased with increasing harvesting intensity, climatic aridity, and climate change scenarios (Figures 2C,D and Table 3). Clearcutting reduced future aboveground C stocks by half, whereas the retention of 60% of the stand area (i.e., large patch) only reduced aboveground C stocks by 33% compared to full retention. Similar to ecosystem C stocks, the negative effect of climatic aridity on aboveground C stocks was more pronounced in the high-intensity harvesting treatments (especially single tree retention) and less pronounced in the low-intensity treatments (patch retention). In addition, aboveground stocks were projected to slightly decrease (−7%; 95% confidence interval [−10.98, −1.81]) under the high-emission scenario (RCP 8.5) compared to the no-emission baseline (i.e., no climate change) across climatic aridity and harvesting treatments. This decline was primarily due to decay of snag and coarse woody debris pools. Under the low-emission scenario (RCP 2.6), however, aboveground C stocks were not projected to decrease significantly (−5%; 95% confidence interval [−9.13, 0.24]).

Belowground C stocks 50 years into the future were less sensitive to harvesting intensity than aboveground C stocks. In particular, clearcutting reduced future belowground C stocks by only 15% (compared to 50% aboveground) and low-intensity harvesting treatments did not have an impact on 2069 belowground stocks (Figures 2E,F and Table 3). However, projected belowground C stocks were as sensitive as aboveground stocks to climatic aridity (same beta coefficient). Similar to results on aboveground C stocks, climatic aridity had a higher impact on belowground C stocks in forests that experienced higher intensity treatments (especially single tree retention). Note that since temperature effects on decay had no significant effect on projected belowground C stocks across the treatments, it was not included in the final model (Table 3).

Discussion

Improved forest management may contribute to climate mitigation needed to hold warming to below 2°C. However, uncertainties persist about the effects of harvesting intensity on forest carbon sequestration, especially when considering interactions with regional climate and climate change. Our modeling study revealed that increasing harvesting intensity led to significant losses in total ecosystem C stocks 50 years post-harvest. In particular, forests that underwent clearcutting were projected to stock about 36% less C by 2069 than forests that were left untouched. In addition, C losses were greater in arid interior Douglas-fir forests than in humid, more productive forests, and growth multipliers and decay due to climate change scenarios had no significant effect on future total ecosystem C stocks.

Harvesting Intensity

Consistent with our first hypothesis that increasing harvesting intensity would result in declining projected C stocks, our study revealed a sharp, linear decrease in total ecosystem C stock from clearcut to large patch retention across all climate regions. We found that clearcut harvesting resulted in the largest reductions in total ecosystem C stocks compared to the lowest harvesting intensity treatments where remaining trees provided an initial source of aboveground C. Comparably, Paquin (2011) used CBM-CFS3 and showed that a 10% commercial thinning scenario resulted in boreal forest C projections that surpassed control scenarios (uncut/natural succession) 300 years after the beginning of the simulation. Similarly, over a 240-year simulation with multiple harvesting rotations, red spruce stands in Nova Scotia, Canada, experienced increased C storage when partially harvested and decreased C storage when clearcut (Taylor et al., 2008). These studies and the present study highlight that partially harvested forests are expected to maintain greater C stocks compared to clearcut forests.

In our study, total ecosystem C accumulation primarily occurred aboveground, which suggests that aboveground C stocks drive total ecosystem C recovery following harvesting. Similar trends were observed in Scots pine forests of Spain where heavily thinned stands contained significantly less (2.6% less) C (296 Mg C ha^–1) than moderately thinned stands (304 Mg C ha^–1), and differences were only significant in aboveground biomass (Ruiz-Peinado et al., 2016). In addition, over a 350-year period, unmanaged stands of red pine in Minnesota stopped accumulating C aboveground at year 100 whereas thinned stands continued to accumulate C (Powers et al., 2012). By contrast, Zhou et al. (2013) conducted a meta-analysis suggesting that, over time, there is no difference in aboveground C stocks among harvesting treatments of different intensities. Hence, it is possible that aboveground C stock would level off in our study if simulations were run further than 50 years.

Clearcutting was the only harvesting treatment that impacted projected belowground C stocks (15% reduction). In general, the effect of harvesting intensity on belowground C stocks is unclear. For example, results from a meta-analysis of 26 studies revealed that harvesting had little to no effect on soil C (Johnson and Curtis, 2001), whereas results from a different meta-analysis suggested that a decline in belowground C stocks of 11% was expected across all harvesting intensities (James and Harrison, 2016). In the forests represented in this study, 60% of forest floor C was lost on average 1-year post-harvest, whereas the mineral soil was unaffected (Simard et al., 2020). The greater loss in forest floor C in Simard et al. (2020) compared with previous meta-analyses may be explained by from greater mechanization of recent logging practices. Many studies stress the importance of further research on belowground C stocks as data below a depth of 60 cm in the mineral soil are scarce (Nave et al., 2010; Clarke et al., 2015). The dataset used to calibrate CBM-CFS3 in the present study considers C stocks within the mineral soil to a depth of 1 m, but in cases where data was not provided to a depth of 1 m, the last horizon observed was assumed to extend to a depth of 1 m (Shaw et al., 2005). The assumption that soil C below 1 m is inactive and unaffected by harvesting and climate change has been challenged; for example, deep soils have been shown to potentially to lose 17.7% of their C following harvest (James and Harrison, 2016).

Climate Change

We hypothesized that future C stocks would be negatively impacted by climate change scenarios. However, increasing temperatures that sped-up decay and/or growth multipliers had little to no effect on total ecosystem C projections. The only statistically significant effect was that the high-emission scenario (RCP 8.5) reduced aboveground C stocks by 7% compared to the no-emission baseline. This was due to more rapid decay of standing dead trees. Similar trends were observed by Garcia-Gonzalo et al. (2007) where results from an alternative growth and yield model (FINNFOR) revealed that total ecosystem C of Finland’s forests increased by only 1% over 100 years due to climate change. However, the use of other models such as LANDIS-II, showed that regions that are currently moisture-limited may become less productive, whereas cold-limited regions will likely become more productive (Dymond et al., 2016; Hof et al., 2017). Altogether, CBM-CFS3’s limitations in simulating climate change may be the cause for the lack of an effect of RCPs in our study since the model is unable to make predictions about changes in tree growth based on changes in temperature and precipitation regimes. To address this issue, we applied a growth multiplier to simulate increased or decreased productivity in tree growth (i.e., the ±17% growth multiplier), but this did not change the magnitude of the effect of climate change on C stocks.

Contrastingly, Metsaranta et al. (2011) found that total ecosystem C increased by 3.7% in an increased productivity scenario and decreased by 8.3% in a decreased productivity scenario. However, the growth multiplier in Metsaranta et al. (2011) was set to ±50% to capture a changing site index for a variety of species, whereas the growth multiplier in the present study was only set to ±17% to reflect changes in Douglas-fir site index. In addition, Metsaranta et al. (2011) ran simulations over an 80-year period compared to the 50-year period in the present study, further allowing C stocks to grow with the ±50% multiplier. Perhaps it is a combination of these two factors (growth multiplier and time-scale) that led to the lack of significance in climate change scenarios in our study. Based on our findings, short term changes in Douglas-fir growth due to climate change are not anticipated to have an impact on ecosystem carbon. However, longer-term studies might find more substantial impacts of these growth multipliers. Our results must be interpreted cautiously, however, because they do no account for widespread mortality that can occur due to amplifications of wildfire or insect disturbances with climate change.

Climatic Aridity

We found that historic aridity plays a significant role in projecting future above- and belowground C stocks for BC’s forests. Consistent with our hypothesis, projections suggest that more arid forest regions are expected to accumulate less C by year 50 compared to more humid forests. This result is expected because abundant precipitation is conducive to greater plant growth, and therefore greater C stocks within forest ecosystems. Note that across BC, winter precipitation is likely to increase compared to historical levels (Spittlehouse, 2008). Consequently, forest ecosystem C stocks may increase across BC as well.

More interestingly, regional climate interacted with harvesting intensity such that ecosystem C stocks in more humid sites were greater within seed tree retention forests compared to humid forests with higher retention harvesting types (small patch and large patch retention). Our result suggests that in more humid forests, the seed tree harvesting type allows forests to accumulate only 5% less C than uncut forests in year 50 of the simulation. This is likely the case because seed trees provide excellent opportunities for regeneration, as seen in black spruce forests of Canada’s boreal forests (Montoro Girona et al., 2018). This information informs potential strategies for maximizing C storage using forest management within BC. Humid climate zones are expected to decrease substantially (Wang et al., 2012), and if we implement seed tree systems into these zones, we might be better able to preserve C stocks. Alternatively, it might be favorable to opt for large patch retention harvesting systems within arid forests, which are expected to expand substantially through the 2080s (Wang et al., 2012).

Conclusion

Our modeling results suggest that variation between climate change scenarios plays a minimal role in affecting forest C stocks. For this reason, we recommend continued research that links field data and modeling experiments to develop a more comprehensive picture of how climate change will most likely impact our forests. Furthermore, we found that aboveground C stocks drive total ecosystem C and are heavily influenced by harvesting intensity. However, we note that there is a large gap in understanding of how harvesting and climate change impact belowground C stocks, which limits the use of carbon budget models. Regardless, in a time of rapid climatic change, it is vital that we protect forest C stocks by implementing low intensity harvesting systems until we have a more complete understanding of the impacts of climate change on BC’s forest C stocks.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AR, SS, CDy, and WR conceived the ideas and designed the methodology. AR and CDy carried out the modeling. CDe and BP analyzed the data. SS and WR collected the data. AR and CDe led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Funding

This research was funded by a Forest Enhancement Society of British Columbia Grant to WR, Forest Carbon Initiative Grant F20-01774 to SS, and Natural Sciences and Engineering Research Council Discovery Grant to SS.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to extend a special thanks to Stephen Kull and Werner Kurz at Natural Resource Canada for providing the CBM-CFS3 workshop and assistance modifying modeling inputs.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2022.934067/full#supplementary-material

References

Appiah Mensah, A., Holmström, E., Petersson, H., Nyström, K., Mason, E. G., and Nilsson, U. (2021). The millennium shift: investigating the relationship between environment and growth trends of Norway spruce and Scots pine in northern Europe. For. Ecol. Manag. 481:118727. doi: 10.1016/j.foreco.2020.118727