Duration of Climate Change Mitigation Benefits from Increasing Boreal Forest Harvest Age by 10 Years

Abstract

We presented a case study and associated method for stand scale assessment of the duration of the climate change mitigation benefit provided by increasing forest harvest age (i.e., the age a stand is harvested). We used stand yield curves and newly developed equations to estimate carbon stocks in various boreal forest ecosystem pools in Ontario, Canada. The proposed method was applied to forest identified as available for harvesting in management plans for three forest management units with a combined area of more than 1900 km². Our analysis indicated that a 10-year increase in harvest age did not provide a mitigation benefit (reduced carbon stock) in about half the available harvest area (45.5%, 61.9%, and 62.1% of the total available harvest area in the management units). Increasing the harvest age by 10 years resulted in a mitigation benefit lasting longer than 25 years for 15.1%, 16.0%, and 13.0% of the total available harvest area in the management units. The results suggest that increasing harvest age may have limited mitigation potential in Ontario’s managed boreal forests in the short-term but can reduce overall carbon stocks in the longer term.

1. Introduction

The continuing rise in atmospheric greenhouse gas (GHG) concentrations and the associated changes to Earth’s climate system necessitate reducing emissions and seeking opportunities to increase carbon removals from the atmosphere [1]. Forests contain vast amounts of carbon and are an integral part of the global carbon cycle [2]. The unique role of forestry stems from the fact that, unlike other sectors of the economy, the use of wood fiber can reduce GHG emissions from other activities, and tree growth removes and stores them from the atmosphere [3]. Consequently, managing forests to maintain or increase their carbon stocks while reducing their emissions is built into many climate change mitigation strategies [1,3].

Forest-based climate change mitigation actions can be broadly categorized as protection, management, and restoration [4]. Protection usually includes reduced deforestation (sometimes referred to as avoided conversion) [4,5]. The restoration of forest cover involves reestablishing forest cover in historically forested areas that were converted to another land use [6,7]. Improved forest management may refer to a suite of actions, such as increased utilization of harvested wood, faster regeneration of harvested sites, and reduced harvesting levels [8,9]. These three categories and their climate change mitigation potential were assessed in [10], in a comprehensive analysis of natural climate solutions for Canada.

Increasing forest harvest age (also referred to as extending rotation length) is one of the mitigation opportunities classified as improved forest management [1]. The literature on the effects of increasing harvest age on carbon stocks varies substantially in methodology. For example, some studies focused on stand scale effects [11,12,13], while others used a landscape scale approach [14,15]. Simulation horizons ranged from 35 [16] to 200 years [15]. Studies also differed in their level of detail of accounting for carbon storage in harvested wood products (HWP) by omitting either upstream [14] or landfill emissions [17].

Despite variations in methods, the results overall indicate that increasing harvest age (extending rotation length) increases combined carbon stocks in forests and HWP. However, the latter conclusion is somewhat compromised by the fact that most of the studies allow for reduced harvest volume as a result of extended rotation length. The only exception we found was [18], in which the authors simulated the effects of different levels of tree retention on forest carbon stocks, while requiring that wood production remained the same in all scenarios. When harvest volumes are permitted to decrease with extended rotation length, the results depend on other assumptions, primarily on the size of substitution benefits from using wood in place of non-woody materials. For example, if substitution effects are projected to be high, then increasing the rotation length may even reduce the combined carbon storage in forest and HWP (e.g., [13]). The opposite is true as well, as documented in [12], who found that in the absence of substitution effects, forest harvesting has higher climate change mitigation potential than a “no harvest” scenario only in the long term (>100 years). Meanwhile, substitution is unlikely to ever be fully realized because of the so-called rebound effect, i.e., the gap between the decreased use of resources expected from increased “eco-efficiency” and their actual utilization [19]. For HWP substitution effects, this means a gap between the expected use of non-woody materials (due to the increased use of HWP) and their actual utilization.

Studies of increased forest harvest age are affected by the temporal aspect of mitigation benefits (often somewhat loosely referred to as benefit permanence). While the long-term persistence of carbon dioxide (CO₂) in the atmosphere makes longer time horizons, during which benefits are sustained preferable, the risk of crossing irreversible climate tipping points makes shorter-term mitigation more relevant, even if it means merely shifting the burden of reducing CO₂ concentrations further in time [20]. In forests, the period of additional sequestration (or storage) is constrained by the stages of forest growth and associated carbon removal. Therefore, assessing the duration during which a mitigation benefit from a given forest-based solution is sustained is just as important as establishing whether the solution has a mitigation benefit.

In Ontario, Canada, harvesting is allowed in almost 28 million ha of managed forest [21,22]. Management plans for these forests allow for nearly 300,000 ha to be harvested annually, although the actual harvest level averaged about 44% of this available area during 2009–2019 [21]. Nevertheless, with an average annual harvest of 121,000 ha, opportunities exist for implementing climate change mitigation actions and, specifically, for changes to harvest age, should it indeed carry climate change benefits.

The study objective was to evaluate the change in forest carbon stocks that occur if harvesting is shifted to 10-year older age classes of trees, with the condition that the total volume harvested remained the same. We assumed constant harvest volumes despite changes in harvest age, to consider only forest carbon stocks and leave HWP outside the study scope. Our focus was on the duration of mitigation benefits after the harvest age was increased. For this analysis, we presented a methodology for constructing mitigation benefit profiles with increasing harvest age for individual forest stands. The methodology was then applied to a sample of forest management units in Ontario to demonstrate the opportunity for mitigation benefit duration with increasing harvest age.

2. Methods and Materials

The general approach consisted of the following stages: (1) developing equations to estimate carbon curves for individual stands (i.e., curves showing the amount of carbon in various forest ecosystem pools); (2) creating a method to construct mitigation benefit profiles of increasing harvest age for a stand with given yield and carbon curves; (3) producing composite yield curves that account for successional transitions; and (4) applying mitigation benefit profiles to a sample of Ontario’s forest management units to assess the area of opportunity for a range of time horizons during which a mitigation benefit is sustained.

2.1. Carbon Curve Equations

Forest stand level carbon curves describing the amount of carbon in a forest ecosystem at a given stand age were developed to emulate the algorithms and parameters of FORCARB-ON2, the Ontario adaptation of the forest carbon budget model FORCARB2 [23,24]. FORCARB2 was developed by the USDA Forest Service and used to generate regional and national-scale estimates of carbon stocks in U.S. forests and HWP. The Ontario version of the model followed the same structure as the original model, but with an added module to simulate the effects of fire disturbance and a substantially more detailed description of carbon stocks and fluxes associated with HWP [23]. To estimate carbon stocks in Ontario’s forests, FORCARB-ON2 includes an updated set of parameter values developed for the application of FORCARB2 in the northeastern U.S. [24]. FORCARB-ON2 has been used in numerous applications to estimate current carbon stocks and to predict future carbon stocks in forests and HWP at forest management unit and provincial scales for various forest management scenarios [25,26,27].

FORCARB-ON2 provides estimates of forest ecosystem carbon stocks in six pools: live trees (above- and belowground), standing dead trees (above- and belowground), down dead wood (DDW, which includes logs and branches ≥ 76 mm diameter as well as stumps), understory vegetation, forest floor (dead organic matter above the mineral soil horizon, including branches and logs < 76 mm diameter, litter, and humus), and soil [25]. As in the original FORCARB2 model, carbon stock estimates in various pools were derived using the relationships with net stand merchantable volume (for carbon in live and standing dead trees) with live tree biomass (for DDW), with stand age (for forest floor and understory vegetation), and with forest type (for soil). These relationships were based on regression equations developed using the USDA Forest Service’s Forest Inventory and Analysis database [28] and were designed to describe the average carbon stocks in forest stands of given species composition and age. While this approach worked well in large scale applications, it is not suitable for estimating carbon stocks in an individual stand. In particular, the estimation of carbon stocks in dead organic matter pools (standing dead trees, DDW, and forest floor) requires more direct links among the pools over time. Therefore, for this study, we developed a new set of equations (hereafter referred to as dynamic) predicting carbon stocks in dead organic matter pools, as follows.

The approach to developing dynamic carbon stock equations consisted of (1) selecting a set of growth and yield curves predicting the amount of net merchantable volume in a stand at a given age; (2) for each curve, calculating carbon stocks in all pools using the original FORCARB-ON2 equations; and (3) fitting a new dynamic equation to the carbon estimates obtained. Since growth and yield curves in Ontario are commonly estimated in 10-year increments, parameter values were assumed to describe changes in carbon stocks over 10-year periods. Some parameters of the dynamic equations were pre-set to values reported in the literature, while others were estimated by fitting the equation to estimates calculated from the original equations; in the latter case, estimated parameter values were averaged over all selected growth and yield curves.

In the first step, carbon stocks in total live tree biomass, B, were estimated using equations from [24,25] (see Table S1 [29,30,31,32] in the Supplementary Materials). The second step was to develop an equation describing the dynamics of carbon stocks in standing dead trees. These dynamics were modeled using the equation:

$$StD\left(t\right)=s_1·B\left(t-10\right)+s_2·\left(B\left(t\right)-B\left(t-10\right)|if >0\right)+s_3·StD\left(t-10\right)$$

(1)

where StD denotes carbon stocks in standing dead trees, t is stand age, and s₁, s₂, and s₃ are parameters. The first term describes an increase in standing dead trees’ stock due to the base mortality of trees, while the second corresponds to an additional transfer of carbon from live to standing dead tree pools due to self thinning; consequently, the second term is proportional to the 10-year difference in total live tree biomass carbon if that difference is positive, and to zero otherwise. The third term corresponds to the fraction of carbon retained by standing dead trees over the preceding 10 years (i.e., standing dead trees that have not fallen). Decay of standing dead trees (as well as DDW and forest floor pools) is commonly described using a negative exponential equation (although other approaches are documented in the literature, e.g., see [33]:

$$StD\left(t\right)=StD\left(0\right)·\exp \left(-k_{StD}·t\right)$$

(2)

where StD(0) is a tree’s carbon when it died and k_StD is the annual decay rate (year⁻¹); then, parameter s₃ in Equation (1) is s₃ = exp(−10∙k_StD).

Table 1. Estimated decay rates of standing dead trees (k_StD) in North American boreal forests.

Common Name	Scientific Name	DBH (cm)	kStD (1/Year)	Location	Source
Softwoods
Balsam fir, red pine	Abies balsamea L. Mill., Picea rubens Sarg.	2.5–42	0.0866 a	Maine, USA	[34]
Firs, hemlocks, Norway and black spruce b		40	0.0495 a	British Columbia, Canada	[35]
Lodgepole pine	Pinus contorta Dougl. Ex Loud.	<25	0.1733 c	Oregon, USA	[36]
		>25	0.0866 c
		20	0.0815 a	Oregon, USA	[37]
		40	0.0660 a
		3–23	0.0693 a	Washington, USA	[38]
		23–41	0.0462 a
		NA	0.0210	Colorado, USA	[33]
		40	0.0433 a	British Columbia, Canada	[35]
		7.5–30	0.0760 d	Montana, USA	[39]
Red pine	Pinus resinosa Ait.	NA	0.0147 e	Ontario, Canada	[40]
White pine	Pinus strobus L.	NA	0.0408 e	Ontario, Canada	[40]
Deciduous
Trembling aspen, cottonwood, and paper birch		40	0.0248 a	British Columbia, Canada	[35]
Hardwoods (mainly oak)		12–24	0.1779	Arkansas, USA	[41]
		>25	0.0528
Hardwoods f		>10	0.1173 g	Wisconsin, USA	[42]
Red maple, white birch	Acer rubrum L., Betula papyrifera Marshall	15	0.0990	Maine, USA	[34]

^a Estimated from half-life. ^b Except subalpine fir. ^c Estimated from the percentage of trees standing 8 years after death. ^d Estimated from the percentage of trees standing 15 years after death. ^e Decay rate estimated from weighted average half-life for 5 decay classes; half life for each decay class calculated from the 5-year probability of survival. ^f Includes black ash, green/white ash, American basswood, paper birch, balsam fir, eastern hemlock, ironwood, red maple, sugar maple. ^g Estimated from the percentage of trees standing after 4.5 years since death.

includes estimates of the decay rate of standing dead trees from the studies involving tree species present or similar to those present in Ontario.

The median (chosen over the mean to avoid the undue influence of possible outliers) decay rates listed in Table 1 are 0.0677 and 0.0990 for softwoods and hardwoods, respectively, which translates into 10-year mass retention rates of 0.5084 and 0.3715. Therefore, parameter s₃ in Equation (1) was assumed equal to 0.50 and 0.37 for softwoods and hardwoods, respectively, while parameters s₁ and s₂ were estimated by means of linear regression using values of B(t) and StD(t) calculated using the original FORCARB-ON2 equations.

In the third step, we modeled changes in down dead wood carbon stocks using the equation:

$$DDW\left(t\right)=d_1·B\left(t-10\right)+d_2·\left(1-s_3\right)·StD\left(t-10\right)+d_3·DDW\left(t-10\right)$$

(3)

where DDW(t) is carbon in down dead wood at stand age t, while d₁, d₂, and d₃ are parameters. The first term on the right side of (3) corresponds to carbon transfer from the live tree biomass pool and reflects contributions to down dead wood from fallen large branches of live trees. Following the approach used in FORCARB2, we assumed DDW contributions from large branches are a constant proportion of live-tree biomass [24] and set parameter d₁ to 0.01. The second term corresponds to the fraction of carbon lost by standing dead trees over the preceding 10 years. Parameter d₃ in Equation (3) is analogous to s₃ in (1), in that it shows the amount of carbon retained in the down dead wood pool over the preceding 10 years. Similar to standing dead trees, down dead wood decay is most often modeled using a negative exponential function.

Table 2. Estimated decay rates of down dead wood (k_DDW) in North American boreal forests. Multiple values for a species reported by a single source reflect either different methods of estimation or multiple sites included in the study.

Common Name	Scientific Name	kDDW (1/Year)	Location	Source
Softwoods
Black spruce	Picea mariana (Mill.) BSP	0.0210	Quebec, Canada	[43]
		0.0600	Manitoba, Canada	[44]
		0.0500
Eastern hemlock	Tsuga canadensis (L.) Carr	0.0210	Wisconsin, USA; Michigan, USA	[45]
Engelmann spruce	Picea engelmannii Parry ex Engelm.	0.0054	Alberta, Canada	[46]
		0.0025
Jack pine	Pinus banksiana Lamb.	0.0420	Minnesota, USA	[47]
Lodgepole pine	Pinus contorta Dougl. ex Loud.	0.0507	Alberta, Canada	[48]
		0.0171	Alberta, Canada	[46]
		0.0299
		0.0153
		0.0045
		0.0035
		0.0720	Alberta, Canada	[49]
		0.0210	British Columbia, Canada	[50]
		0.0180
		0.0120	Wyoming, USA	[51]
Red pine	Pinus resinosa Ait.	0.0550	Minnesota, USA	[47]
Red spruce	Picea rubens Sarg.	0.0330	New Hampshire, USA	[52]
Sitka spruce	Picea sitchensis (Bong.) Carr.	0.0110	Washington, USA	[53]
Subalpine fir	Abies lasiocarpa (Hook.) Nutt.	0.0286	Alberta, Canada	[48]
		0.0520	Alberta, Canada	[49]
Western hemlock	Tsuga heterophylla (Raf.) Sarg.	0.0124	Washington, USA	[54]
White spruce	Picea glauca (Moench) Voss	0.0271	Alberta, Canada	[48]
		0.0710	Minnesota, USA	[47]
		0.0240	Alberta, Canada	[49]
Hardwoods
American basswood	Tilia americana L.	0.0750	Wisconsin, USA	[55]
Bigtooth aspen	Populus grandidentata Michx.	0.0900	Michigan, USA	[56]
Sugar maple	Acer saccharum Marsh.	0.0750	Wisconsin, USA	[55]
Trembling aspen	Populus tremuloides Michx.	0.0800	Minnesota, USA	[47]
		0.0800
White ash	Fraxinus americana L.	0.0495	Minnesota, USA	[57]
Hardwoods a		0.0365 b	Ontario, Canada	[40]
		0.0330 c

^a Mostly sugar maple, with important components of American beech, yellow birch, and eastern hemlock. ^b Unmanaged stands; decay rate calculated from half-life time. ^c Managed stands; decay rate calculated from half-life time.

shows the rates of negative exponential decay as reported in the literature for relevant tree species.

Using the same approach as for standing dead trees, we estimated the median decay rate from Table 2 and translated it to a decadal retention rate to set the parameter d₃ in Equation (3) to 0.75 and 0.50 for softwoods and hardwoods, respectively. Finally, parameter d₂ was estimated using linear regression by fitting Equation (3) to the values of B(t), StD(t), and DDW(t). The second term could have been simplified to contain only one parameter as a multiplier for StD(t − 10), but we used the explicit form to show that this source of carbon transfer to the down dead wood pool is the fraction of the standing dead tree pool lost during the preceding decade. Note that B(t) and DDW(t) were calculated using the original FORCARB-ON2 equations, while StD(t) was estimated from B(t) using Equation (1) and its parameters, estimated in the first step of the current methods.

In the fourth step, we estimated parameters of the following equation, which is used to simulate changes in forest floor carbon:

$$\begin{array}{ll}FF\left(t\right)& =f_1·B\left(t-10\right)+f_2·\left(1-d_2\right)·\left(1-s_3\right)·StD\left(t-10\right)\\ & +f_3·\left(1-d_3\right)·DDW\left(t-10\right)+f_4·FF\left(t-10\right)\end{array}$$

(4)

where FF(t) is carbon in the forest floor at stand age t, while f₁, f₂, f₃, and f₄ are parameters. Terms on the right side of Equation (4) correspond to carbon transfers from live tree biomass, standing dead trees, down dead wood, and forest floor, respectively. Parameter f₄ in Equation (4) is similar to s₃ in (1) and d₃ in (3) and shows the amount of carbon retained in the forest floor over the preceding 10 years, modeled using a negative exponential function. To reduce the number of parameters to be estimated by fitting Equation (4) to data for the forest floor, we reviewed the literature on decay rates of forest floor components, summarized in

Table 3. Estimated decay rates of forest floor carbon (k_FF) in North American boreal forests (where “Years” indicates study period).

Common Name	Scientific Name	kFF (1/Year)	Years	Location	Source
Leaves/needles
Softwoods
Black spruce	Picea mariana (Mill.) BSP	0.1703	6	Canada a	[58]
Douglas-fir	Pseudotsuga menziesii (Mirb.) Franco	0.1294	6	Canada a	[58]
		0.1833	5	British Columbia, Canada	[59]
Jack pine	Pinus banksiana Lamb.	0.1446	6	Canada a	[58]
Lodgepole pine	Pinus contorta Dougl. ex Loud.	0.2811	4	British Columbia, Canada b	[60]
		0.2100	5	British Columbia, Canada	[59]
		0.1316	3	Alberta, Canada c	[61]
		0.1500	Comp d	Wyoming, USA	[62]
Red pine	Pinus resinosa Ait.	0.2550	6	Massachusetts, USA	[63]
Tamarack	Larix laricina (Du Roi) K. Koch	0.1223	6	Canada a	[58]
Western red cedar	Thuja plicata Donn ex D. Don	0.1058	6	Canada a	[58]
White spruce	Picea glauca (Moench) Voss	0.2582	5	British Columbia, Canada	[59]
		0.1828	3	Alberta, Canada	[61]
Pine	Pinus spp. e	0.0740	10	Alaska, USA	[64]
		0.0905	10	Alaska, USA	[64]
		0.0438	10	Colorado, USA	[64]
		0.2037	10	New Hampshire, USA	[64]
		0.1309	10	Washington, USA	[64]
		0.0285	10	Wisconsin, USA	[64]
Hardwoods
Beech	Fagus grandifolia Ehrh.	0.1122	6	Canada a	[58]
Black oak	Quercus velutina Lam.	0.3100	6	Massachusetts, USA	[63]
Red maple	Acer rubrum L.	0.2250	6	Massachusetts, USA	[63]
Sugar maple	Acer saccharum Marsh.	0.1175	10	Alaska, USA	[64]
		0.1184	10	Alaska, USA	[64]
		0.0652	10	Colorado, USA	[64]
		0.1464	10	New Hampshire, USA	[64]
		0.1636	10	Washington, USA	[64]
		0.0674	10	Wisconsin, USA	[64]
Trembling aspen	Populus tremuloides Michx.	0.1446	6	Canada a	[58]
		0.2733	3	British Columbia, Canada b	[60]
		0.3290	5	British Columbia, Canada	[59]
White birch	Betula papyrifera Marsh.	0.1798	6	Canada a	[58]
		0.2100	5	British Columbia, Canada	[59]
Yellow birch	Betula alleghaniensis Britt.	0.2700	6	Massachusetts, USA	[63]
Samples from F and H horizons, surface-placed		0.1563	3	British Columbia, Canada	[65]
Wood stakes, surface-placed
Ramin	Gonystylus bancanus (Miq.) Kurz	0.0408	10	USA f	[66]
Trembling aspen	Populus tremuloides Michx.	0.0209	10	Montana, USA g	[67]
Western hemlock	Tsuga heterophylla (Raf.) Sarg.	0.0604	6	Canada h	[68]
Wood stakes, buried
Trembling aspen	Populus tremuloides Michx.	0.0432	10	Montana, USA	[67]
Forest floor material
Samples from F horizon, buried		0.0424	4	British Columbia, Canada	[60]
		0.0946	4	British Columbia, Canada	[60]
		0.0959	4	British Columbia, Canada	[60]

^a Average of decomposition rates observed at 18 upland forest sites across Canada. ^b Average decomposition rate across 10 sites located in boreal white and black spruce, sub-boreal spruce, interior Douglas-fir, and interior cedar–hemlock subzones in BC, Canada. ^c Average decomposition rate over three sites. ^d Observations on foliage in various stages of decomposition were used to construct an 8-year-long chrono-sequence. ^e Pine species varied across study site and included eastern white pine (Pinus strobus L.), red pine (Pinus resinosa Ait.), and slash pine (Pinus elliotii Engelm.). ^f Average decomposition rate over 5 sites located in boreal, temperate conifer, and temperate deciduous forests. ^g Average decomposition rate over two sites. ^h Average of decomposition rates observed at 16 upland forest sites across Canada (two sites from [58] with a shorter observation period not included).

.

If decay rates reported in the literature did not use a negative exponential model, then the rates in Table 3 were calculated from the mass of forest floor material remaining at the end of the observation period. As is evident from Table 3, the number of studies on decay rates differs among forest floor components. Foliar decay is studied relatively more often than woody components of the litter (L layer of forest floor); similarly, fewer studies were conducted on the decay rate of F and H horizons than of forest floor litter. To estimate a composite decay rate for forest floor carbon, we combined rates for individual components using their fractions in the respective layer, and then combined components based on their fraction in the forest floor.

Studies focused on assessing the above fractions are relatively infrequent and reported estimates vary considerably. For example, Puhlick et al. [69] studied forest floor composition in central Maine, USA, in uneven-aged conifer-dominated stands (that included balsam fir, red spruce, eastern hemlock, northern white-cedar, and eastern white pine mixed with maples, birches, and aspens) and estimated fractions of non-woody and woody materials in the litter layer (O_i horizon) at 94% and 6%, respectively. However, Mukhortova and Bezkorovainaya [70] estimated the fraction of needles and leaves at only about 20% of the litter layer of larch-dominated stands in central Siberia. Similar variation is evident in estimates of the fraction of the litter layer in the overall forest floor mass. Prescott et al. [71] studied litter accumulation in three stands (dominated by lodgepole pine, white spruce, and fir-spruce, respectively) in the Kananaskis Valley in the Rocky Mountains of southwestern Alberta (Canada) and concluded that “only 1–3% of the total litter mass was in the L layer.” Similar estimates were reported in [72] in a study of litter accumulation and decomposition on 16 sites in predominantly hardwood forests (containing a minor component of coniferous species) representing four climatic regions in Maine, USA, with the average fraction of litter in the overall forest floor mass ranging from 3.3–7%. Slightly higher litter layer fractions relative to overall forest floor mass were observed in [73] (11–18% in two different larch ecosystems in central Siberia) and [74] (14–27% in beech forests in southern France). The results of [69] were at the highest end of the range, with the litter layer constituting 39% of the total forest floor mass.

Given variations observed for litter mass, we used median values for decay rates of foliage and woody components in the litter layer (from Table 3) and applied these rates to permutations of the fraction of foliage in the litter layer and the fraction of the litter layer in the overall mass of forest floor, using the lowest and highest estimates of these fractions from the literature cited above. The resulting mass retention rate over 10 years ranged from 33.53% to 64.65%, with a median value of 47.90%. The latter rate was used as parameter f₄ when fitting Equation (4). We calculated the values of parameters f₂ and f₃ using the ratio between mass loss in standing dead trees and down dead wood that can be attributed to fragmentation and mineralization. Lambert et al. [75] studied the decay of balsam fir boles in New Hampshire, USA, and attributed 37% of mass loss to mineralization and 63% to fragmentation. We estimated the fraction of losses due to mineralization by fitting a negative exponential function to the 10th year data from Oregon, USA [54], on volume and density losses to arrive at 24.3% and 28.9% for Douglas fir and western hemlock snags, respectively. Based on these estimates, we assumed that parameter f₃ equaled 0.7 (i.e., 70% loss due to fragmentation), while parameter f₂ was calculated from the estimate of d₂ using the same fragmentation loss rate of 70% for standing dead trees (i.e., total fragmentation losses minus the fraction transferred to down dead wood). For comparison, a similar fractional loss due to fragmentation was used in [76] (38.4% calculated for 10-year losses due to decay and physical transfer of carbon in snag stems).

To estimate the remaining parameters in Equations (1), (3), and (4) we used yield curves developed for so-called standard forest units (standard FU), defined as an aggregation of forest stands with similar species composition that develop in a similar manner and are managed under the same silvicultural system [77]. Curves were selected to represent managed provincial forests across Ontario, Canada (Figure 1). We followed the approach used in FORCARB-ON2 to estimate carbon stocks in live vegetation and dead organic matter pools using volume densities for the following species groups: P—jack, red, and white pine; SF—black and white spruce and fir; AB—aspen-birch; and MB—maple-beech. To calculate carbon stocks for each curve, each of the 16 major Ontario tree species represented in the standard FUs was assigned to one of the species groups (P, SF, AB, or MB); in addition to the above listed tree species, eastern white cedar and hemlock and “other conifers” (species category in FU) were assigned to the SF group, balsam poplar to AB, and red oak and yellow birch to MB. The species group volume density (m³ ha⁻¹) was estimated as the sum of net merchantable volumes for individual species in the group divided by the fraction of the species group volume in the total stand volume. The species group volume density was converted to live tree carbon using FORCARB-ON2 equations, and the resulting carbon stocks were summed to provide stand totals weighted by the fraction of a species group’s volume in total stand volume.

We applied the same approach to estimate parameters in Equations (1), (3), and (4) for the four species groups (P, SF, AB, or MB). For P, SF, and AB species groups, we used 50 yield curves, while for MB it was 31 yield curves. Fewer yield curves were available for MB because of the more limited geographic distribution of these species in Ontario’s managed forests. Finally, to assess the accuracy of Equations (1), (3), and (4), we calculated relative standard errors for 115 yield curves for Ontario’s standard FUs that were not used to estimate parameters. For a given stand age, we calculated the relative difference between the estimates produced by FORCARB-ON2 and newly developed dynamic equations, divided by those of FORCARB-ON2; the square root of the mean relative difference averaged over the number of ages covered by the yield curve provided the relative standard error for the given yield curve.

2.2. Mitigation Benefit Profiles with Increasing Harvest Age

To assess the effects of harvest age on forest carbon, our approach was as follows. Assume that we need to harvest total merchantable volume M from a sufficiently large forest stand. If harvested at age A (referred to as base age), the area S_A needed to produce M equals:

$$S_A=\frac{M}{V\left(A\right)}$$

(5)

where V(A) is the net merchantable volume of a stand at age A per unit of area. If harvest age was increased by dA to A + dA, then the area S_A+dA needed to produce M would equal:

$$S_{A+dA}=\frac{M}{V\left(A+dA\right)}$$

(6)

Consequently, S_A and S_A+dA are related as:

$$S_{A+dA}=S_A·\frac{V\left(A\right)}{V\left(A+dA\right)}$$

(7)

Let us assume that the stand area is no less than S_A + S_A+dA; for simplification, assume that the stand is of post-fire origin (to avoid the need for additional assumptions about initializing dead organic matter pools that, in a post-harvest stand, depend on carbon stocks in the preceding stand). Let us also assume that the harvested area is immediately regenerated using the same tree species and density as were in the pre-harvest stand. If volume M was harvested at base age A, then total forest carbon stocks t years after harvest would equal:

$$CTOT\left(A+t\right)=S_{A+dA}·C\left(A+t|0\right)+S_A·C(t|A)$$

(8)

where CTOT(A + t) is the total forest ecosystem (including all pools) carbon stock in area S_A + S_A+dA, C(t) is the total forest ecosystem carbon stock at stand age t per unit of stand area, and the variable after the vertical divider specifies the age of the stand before the most recent harvest that is needed to initialize carbon stocks in dead organic matter pools. If volume M was harvested at age A + dA, then total forest carbon stocks at the same point in time as in (8) would equal:

$$CTOT\left(A+t\right)=S_A·C\left(A+t|0\right)+S_{A+dA}·C(t-dA|A+dA)$$

(9)

In both (8) and (9), the first and second terms on the right side of the equation correspond to unharvested and harvested areas, respectively; note that the variable after the vertical divider for unharvested areas (S_A+dA in (8) and S_A in (9)) is set to zero because the stand is of post-fire origin. Thus, the scenario (harvesting at age A or A + dA) having higher total carbon stocks (CTOT) at time A + t is determined by comparing the quantities produced by Equations (8) and (9). To simplify this comparison, we can substitute S_A+dA from (7) into (8) and (9) and then divide both by S_A (division by a positive number does not affect which quantity is larger). Then, the comparison would be between the quantities:

$$\frac{V\left(A\right)}{V\left(A+dA\right)}·C\left(A+t|0\right)+C\left(t|A\right)\langle ?\rangle C\left(A+t|0\right)+\frac{V\left(A\right)}{V\left(A+dA\right)}·C(t-dA|A+dA)$$

(10)

If the left side of (10) is larger than the right, then harvesting area S_A at base age A results in higher total forest carbon stocks at time A + t, and vice versa. Here, increasing harvest age provides a climate change mitigation benefit at time A + t if total carbon stocks in the scenario with increased harvest age are higher at time A + t than those in the base harvest age scenario, i.e., right side of (10) is greater than the left side; the length of time t, during which the latter inequality is true, constitutes the duration of climate change mitigation benefit. Conversely, increasing harvest age carries no climate change mitigation benefit at time A + t if total carbon stocks in the scenario with increased harvest age are lower at time A + t than those in the base harvest age scenario (i.e., right side of (10) is less than the left side). Note that determining which scenario results in higher stocks is unaffected by the assumption that the stand is large enough to contain S_A + S_A+dA. The same comparison as in (10) is arrived at by means of simple transformations if we assume that total stand area was either S_A or S_A+dA (whichever is larger) and then estimate total stocks in that area for both scenarios. To restate the approach, the determination of the mitigation benefit with a dA delay in harvest depends on the stand carbon with vs. without dA delay in the “base” forest area (i.e., the harvested area at the base age A), as well as in the area that has to be harvested if base harvest age A was delayed by dA.

2.3. Composite Yield Curves and Areas of Mitigation Opportunity

To investigate whether increasing harvest age increases carbon stocks (i.e., provides a climate change mitigation benefit), we used the most recent management plans for three forest management units (FMU). Forestry activities in each FMU are governed by an approved forest management plan (FMP) that is updated on a ten-year cycle. FMPs and reports on completed forest operations are submitted annually and are available for public viewing (https://nrip.mnr.gov.on.ca/s/fmp-online?language=en_US accessed 4 August 2022). The three selected FMUs are in the boreal zone of Ontario’s managed forest (Figure 1) and cover 1.9 million ha of forested area (1,112,000, 378,000, and 425,000 ha for Hearst Forest, Lakehead Forest, and Pic River Forest, respectively). The latest FMPs for these FMUs were approved in 2019 for Hearst Forest, 2020 for Lakehead Forest, and 2013 for Pic River Forest. Each plan includes yield curves for FUs present in the FMU and area available for annual harvesting by FU; yield curves allow net merchantable volume to be estimated for each parcel of land available for harvesting and a carbon curve to be derived for each FU.

Yield curves included in FMPs provide volume estimates for stand ages 5–255 years. The fate of older stands (either approaching 255 years or surviving past this age) is defined by FMU-specific forest successional transitions; for each FU, one or more transitions are possible, and these are specified as the fraction of area that, at a given age, succeed from one FU into another, along with the age of the post-transition FU. To assess carbon stocks in forests past the age of 255 years, we constructed composite yield curves by combining yield curve for the original FU (pre-transition curve) with yield curves for all post-transition FUs, using fractions of the succeeding area as weights; if in turn a fraction of the area occupied by a post-transition FU could undergo another transition before the original stand reached age 325, these secondary transitions were accounted for as well. In addition to allowing carbon stocks to be estimated past the age of 255, composite curves ensured that carbon stocks in older stands were not unduly affected by a commonly observed artifact that can occur in yield curves used in FMPs, in which a steep decline in net merchantable volume occurs past the age of maximum yield. An example of successional transition rules and construction of a composite yield curve is provided in the Supplementary Materials.

We calculated FU carbon stocks for stand ages 5–325 years using the above-described composite yield curves and then constructed a “mitigation benefit profile” to examine the effect of increasing harvest age on carbon stocks. Mitigation benefit profiles compared the two sides of inequality (10) for each curve and assessed the length of time t, during which increasing a given base harvest age A by a given increment dA would result in higher total forest carbon stocks. For the base age ranges we used, the minimum and maximum ages of stands available for harvesting as specified in the FMPs were, respectively, 75–225 for Hearst Forest, 55–185 for Lakehead Forest, and 75–205 for Pic River Forest. To reduce the number of possible combinations, we used one 10-year increment as dA. Finally, we combined these carbon profiles with the FMP-defined areas available for harvesting to estimate the areas of mitigation opportunity, i.e., areas in which increasing base harvest age by 10 years would increase total forest carbon stocks.

3. Results

Estimated parameter values for Equations (1), (3), and (4) are presented in Table S1. The verification procedure showed good agreement between estimates from the equations developed in this study and those from FORCARB-ON2. The relative standard error for tested curves ranged from 2.48–14.75%; the average relative standard error for 115 curves was 6.82%. An example of estimates produced by newly developed versus FORCARB-ON2 equations is shown in Figure 2 for a jack pine-dominated stand growing on a productive site in northwestern Ontario (estimates of soil carbon stocks not shown). As seen from Figure 2, the main difference between estimates is caused by the forest floor at older ages. Estimates produced by the new equations are driven by stocks in other pools contributing to the forest floor (live trees, standing dead trees, down dead wood), while the original estimates were based entirely on stand age, resulting in the same slowly increasing curve regardless of the actual stocks in other carbon pools. Consequently, if the standard relative error between the new and old estimates was calculated for each yield curve to age 200 years (instead of 255), then the average over 115 tested curves would drop to 4.95%. According to Watkins [78], forests over 180 years old account for just 0.4% of the total area occupied by managed forests in Ontario. Thus, for most of Ontario’s managed forests, the potential error introduced by switching from the original FORCARB-ON2 equations to those developed in this study is less than 5%.

Composite yield curves were constructed for all 36 forest units present in the three FMUs: 13 each for Hearst Forest and Lakehead Forest and 10 for Pic River Forest; for the description and species composition of the forest units, see Table S2 in Supplementary Materials. The number of successional transitions per forest unit ranged from one to eight. Figure 3 provides an example of a composite curve for a jack pine-black spruce forest unit (with admixture of trembling aspen) for Pic River Forest. As seen in Figure 3, the net merchantable volume of a stand declines sharply after reaching maximum at age 125 (dashed green line). This decline is because an individual yield curve reflects the volume of the original stand, ignoring subsequent successional transitions. To avoid these artificially large changes in net merchantable volume, the composite curve accounts for all possible transitions and is estimated as the sum of area-weighted post-transition yield curves. Note that, in Figure 3, yield curves for post-transitional FUs are weighted by the fraction of area that undergoes a given transition.

Mitigation benefit profiles for increasing base harvest age by 10 years were constructed for each composite yield curve. An example of a carbon profile for a jack pine-black spruce FU is shown in Figure 4; the mitigation profiles for all 36 analyzed FUs is included in Supplementary Materials (see Figure S1). For each possible harvest age, the green horizontal bar corresponds to the duration of time during which carbon stocks would be higher if harvesting was delayed 10 years; the brown bar shows the duration of time for which the opposite is true. For example, the length of the green bar for a base harvest age of 75 years corresponds to 33.6 years, meaning that if harvesting was delayed from 75 years until stand age 85, carbon stocks would be higher until the stand reaches an age of 108.6 (= 75 + 33.6), after which the pattern reverses and the stand has less carbon. In other words, the mitigation benefit of increasing the harvest age from 75 to 85 for this forest unit lasts 33.6 years, beyond which the forest contains less carbon than if harvesting were not delayed. As seen from Figure 4, the duration of the mitigation benefit due to increasing harvest age by 10 years decreases with increasing base harvest age. At base harvest age 115, not only is mitigation benefit nil, but increasing harvest age by 10 years reduces forest carbon stocks. This pattern persists until harvest age 155, at which time delaying the base harvest age by 10 years would again result in a mitigation benefit, lasting from 15.5 years (if harvesting was delayed from stand age 155 to 165 years) to 29.4 years (if harvesting was delayed from stand age 185 to 195 years).

We used mitigation benefit profiles for individual forest units to assess the length of mitigation benefit obtained by increasing harvest age for stand areas of different age class available for harvesting in each FMU. For example, the total forest area available annually for harvesting in Hearst Forest is 7771 ha. In this area, delaying harvesting by 10 years on 3544 ha (45.5%) would reduce forest carbon stocks, i.e., in these stands, delaying harvesting would cause a net increase in atmospheric carbon. The breakdown of the remaining 54.5% of the annual available harvest area by the duration of mitigation benefit is shown in Figure 5a; in 21.1% of the available harvest area, increasing the base harvest age by 10 years would provide a 15- to 20-year long mitigation benefit; in 18.2% of the area, the benefit would last between 20 and 25 years, etc. Estimates for the other two FMUs are 4298 ha and 3941 ha (total harvest area available annually), with 2660 ha (61.9%) and 2445 ha (62.1%) of the harvest area providing no mitigation benefit (i.e., forest carbon stocks are reduced) from increasing the base harvest age by 10 years in the Lakehead Forest and Pic River Forest, respectively. The breakdown of harvest area by duration of the mitigation benefit due to increasing base harvest age by 10 years for the Lakehead and Pic River Forests is shown in Figure 5b,c, respectively.

4. Discussion

4.1. Duration of Mitigation Benefits

The approach presented in this study allows for the construction of a temporal profile of the mitigation benefits for a given yield curve by assessing, for each possible harvest age, the duration of the benefit. This approach was demonstrated using a harvest age increase of 10 years. The effects of increasing harvest age by more than 10 years can be analyzed using the same approach. However, 10 years was a practical choice given that both yield curves and forest age structure are provided in 10-year increments in the study jurisdiction’s forest management plans. These results indicate that, depending on the age-carbon profile of individual stands, increasing stand age produced either an increase or a decrease in stand carbon over time.

Total forest carbon stocks in both scenarios (with and without harvesting delay) is equal to the sum of stocks in the unharvested area and the area regenerating after harvesting (Equations (8) and (9)). The dynamics of this sum depends on the rate of carbon accumulation in the regenerating area and on the base harvest age. As the age of unharvested area increases, its total carbon stocks start declining (Figure 2). Once the break-even point between carbon stocks in scenarios with and without base harvest age increase is reached (i.e., the two sides of Equation (10) are equal), climate change mitigation benefit turns into carbon loss and vice versa. Thus, the duration of the mitigation benefit for a given stand depends on the regeneration rate and base harvest age.

The pattern shown in Figure 4 is common to almost all FUs (see Figure S1) and can be explained qualitatively based on the position on the yield curve at the harvest age (Figure 3). If harvesting takes place when the stand volume is rising, increasing the harvest age allows the wood volume to be obtained by harvesting a smaller area, resulting in a mitigation benefit. The picture reverses as the base harvest age approaches the age of maximum stand volume (the first maximum shown in Figure 3), hence the duration of the mitigation benefit becomes shorter as the base harvest age increases, eventually turning into carbon losses (i.e., at a certain point, not only there is no mitigation benefit, but increasing harvest age creates a carbon source compared to the base case). The decline in stand volume at older stand ages is caused by tree mortality; once the understory trees grow and fill the gaps in the main canopy, the volume starts increasing again (the right maximum in Figure 3). Consequently, as the base harvest age approaches the low point between the two volume maxima, delaying harvesting may again result in a mitigation benefit. The pattern in Figure 4 is not entirely determined by the relative volume density at the base and increased harvest ages and the respective harvest areas. The accurate detection of mitigation benefits (their presence and duration) is determined by the sign of inequality (10); other factors affecting it include the carbon stocks in other forest ecosystem pools, steepness of the yield curve prior to and post the point of maximum stand volume, etc. However, there is a general correlation between the presence and duration of a mitigation benefit and the position of the base harvest on the yield curve that helps explain the pattern shown in Figure 4.

Carbon accumulation in the regenerating area may be accelerated with enhanced regeneration; for example, Drever et al. [10] assumed 10% enhancement of forest yield in harvested areas, achieved by protecting advanced regeneration during harvest, planting seedlings well-adapted to future growing conditions, and moderately fertilizing. However, such measures would constitute a climate change mitigation action that is separate from increasing harvest age. To avoid confounding the effects of two actions, we assumed that the regeneration rate in the post-harvest area was identical to that of the original stand.

Our application of mitigation benefit profiles to three Ontario FMUs showed limited opportunity for this type of climate change mitigation. Results were consistent across the FMUs evaluated, with only about half or less of the available harvest area producing any mitigation benefit from increasing base harvest age by 10 years. More importantly, for most of this area, the benefit of delaying harvesting is short-lived; the area with mitigation benefit lasting longer than 25 years amounted to only 15.1, 16.0, and 13.0% of the total available harvest area for Hearst Forest, Lakehead Forest, and Pic River Forest, respectively.

4.2. Computational Approaches

Using composite curves rather than those provided in forest management plans was necessary for two reasons. First, it allowed volume to be accounted for from trees that emerge during stand succession and corrected an anomalous abrupt volume reduction that may occur because yield curves track the net merchantable volume of the original stand only. Second, it extended yield curves past age 255 (the maximum age specified in FMPs), which allowed us to use a 100-year-long horizon for assessing mitigation benefits for harvest ages as old as 205 years. A similar approach to constructing composite curves was used in [79], in the study of management scenarios in parts of northeastern Ontario with long natural fire cycles. Note that stand ages at later stages of composite curves (in post-transitional stands) refer not to the age of trees, but to time since the last stand-replacing disturbance.

We developed equations estimating carbon stocks in all forest ecosystem pools. These equations reflect transfers of matter between various pools, with the values of more than half the parameters based on field studies of decay rates of dead organic matter pools. The new equations also provided an explicit connection among some of the pools; for example, forest floor carbon stocks now depend on the stocks in other pools, as opposed to the previous formulation in which values for forest floor stocks were generated using an average age-driven regression equation. Since the original equations used in FORCARB2 were fitted to data from the extensive USDA Forest Service’s Forest Inventory and Analysis database [28], testing our values against estimates produced by the original equations is validation by proxy, in which models are tested using previously established relationships between original data and surrogate variables [80].

4.3. Direct and Indirect Effects of Harvest Age Adjustment

In addition to the relevance of these results to finding climate solutions that provide long-term mitigation benefits, this study is relevant to Canada’s commitment to achieve net-zero GHG emissions by 2050 [81]. A net-zero target for 2050 puts more emphasis on actions with a benefit horizon that goes beyond the middle of this century; for example, Drever et al. [10] explicitly stated that this goal prompted them to estimate the mitigation potential of natural climate solutions to 2050. Our analysis also showed the length of time during which a cumulative mitigation benefit accrues—the fact that a cumulative benefit comes to zero indicates that the net flux for the mitigation action becomes positive (i.e., emissions outweigh removals) even earlier. Thus, a consideration of fluxes would further reduce the area of forest where increasing harvest age could result in a mitigation benefit in 2050.

As previously discussed, most studies allow harvest volume to vary while studying the effects of extending rotation length, with the notable exception of [18], whose authors stated, “we are not aware of any study that keeps the wood production constant, as we did in this study” (p. 43). In addition, aside from the previously mentioned issues with substitution effects and carbon stocks in HWP in use and landfills, allowing harvest volume to change ignores the leakage of harvesting to other jurisdictions, in which case the local effects of delayed harvesting are likely overstated, or the possibility of increased emissions due to substitution with non-wood materials that have higher production emissions. Leakage refers to phenomena in which reduced emissions in one jurisdiction are offset by increases elsewhere [82]. According to some estimates, up to 80% of reduced harvesting may be replaced by increased harvesting elsewhere [83]; the average amount of leakage estimated in [84] based on their meta-analysis of 46 studies was 39.6%. Thus, studies ignoring leakage may predict a local mitigation potential from increasing rotation length, but they do not indicate a global effect. Accounting for potential leakage or substitution with non-wood materials would undoubtedly reduce climate change mitigation benefits assessed in such studies.

Assuming harvest volumes remain constant regardless of changes in harvest age alleviated us of the need to consider carbon stocks in HWP and associated issues (such as substitution effects and leakage). However, it made it difficult to compare our results with the other published literature on the effects of harvest age. By and large, other studies indicate that increasing harvest age provides a mitigation benefit through an increase in forest carbon stocks. Our results suggest that constraining harvest volume production at a constant level reduces such projections to a 50-50 split between areas with and without mitigation benefits from increased harvest age.

Implementing increased harvest age may disrupt wood supply in some jurisdictions, as stands scheduled for harvesting are “placed on hold” until they reach the desired delay in harvest age. Such a situation is less likely to arise in Ontario, where less than 50% of the available area is being harvested, making it likely that older stands are available for harvesting immediately. However, other factors may influence the decision to increase harvest age, even in areas where it would lead to a mitigation benefit. First, increasing harvest age may lead to changes in forest age structure beyond thresholds allowed by the rules of sustainable management. Second, it may increase costs by forcing harvesting to more remote stands. Both factors were beyond the scope of our study. We also did not attempt to account for possible changes in forest growth or decomposition caused by changing climate; however, this omission is not expected to have a substantial effect on our results due to the short life span of mitigation benefits (i.e., in more than 84% of available harvest area in the three FMUs considered, the mitigation benefit lasts less than 25 years) and the slow growth of forests in Ontario.

The focus of this study was the duration of mitigation benefit rather than the benefit size. Next, we plan to address the potential scale of mitigation benefits by considering the effects of the continuous implementation of harvest age changes at the landscape scale. Although additional research is needed to more fully assess the effects of increasing harvest age, our study findings suggest that it may have limited climate change mitigation potential in Ontario’s managed boreal forests.