# Combining thinning and diverse plantings to adapt to climate-change-induced timber supply shortage in British Columbia

Valentine Lafond<sup>a,b</sup>, Adam D. Polinko<sup>a,c</sup>, Cosmin D. Man<sup>d</sup>, Caren C. Dymond<sup>e</sup>, Gregory Paradis<sup>a</sup>, and Verena C. Griess<sup>f</sup>

<sup>a</sup>Department of Forest Resources Management, Faculty of Forestry, University of British Columbia, Forest Sciences Centre, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada; <sup>b</sup>Université de Lorraine, AgroParisTech, INRAE, Silva, 54000 Nancy, France; <sup>c</sup>Department of Forestry, Mississippi State University, Box 9681, Mississippi State, MS 39762, USA; <sup>d</sup>Forsite Consultants Ltd., #330 42<sup>nd</sup> Street SW, Salmon Arm, BC V1E 4R1, Canada; <sup>e</sup>Forest Carbon and Climate Services Branch, Government of British Columbia, PO Box 9544, Stn Prov Govt, Victoria, BC V8W 1T7, Canada; <sup>f</sup>Institute of Terrestrial Ecosystems, Department of Environmental System Sciences, ETH Zurich, Zurich, Switzerland

**Corresponding author:** Verena C. Griess (email: verena.griess@usys.ethz.ch)

## Abstract

Forestry is an important component of Canada's economy with British Columbia (BC) contributing almost half to the national roundwood production. Yet, the country's timber supply and forest economy are threatened by climate change, with increased frequency and severity of natural disturbances and changes in forest productivity. Mountain pine beetle (Dendroctonus ponderosae) outbreaks are endemic in BC, but the latest climate change-driven outbreak has resulted in a cumulative loss of over half of all merchantable pine, leading to a mid-term timber supply shortage. In this study, we investigate the potential of commercial thinning and alternative planting regimes based on species diversification and assisted species migration to mitigate the anticipated decrease in timber supply in BC. We simulated the long-term effects of these management options in a case study area in interior BC, using a toolbox that combines management- and research-oriented data and models. We found that combining commercial thinning and species diversification has the best potential to mitigate future timber supply shortages in BC. We discuss the limits of this toolbox approach and identify research needs and recommendations for future studies aiming at modelling cumulative effects of management, climate change, and natural disturbances on timber supply.

**Key words:** forest management, timber supply, commercial thinning, species diversity, climate change adaptation, forest growth and yield model, forest estate model

## 1. Introduction

Canada encompasses approximately 362 million ha of forest, with Canada's forest product-based industry contributing CAD $34.8 billion to the nation's gross domestic product in 2021 (NRCan 2022). About 40% of the industrial roundwood harvested in Canada comes from British Columbia (BC), with three quarters of this volume coming from interior BC (Government of British Columbia 2021; NRCan 2022). The mountain pine beetle (MPB; *Dendroctonus ponderosae*) outbreaks of 1975–1985 and 2002–2012 caused major mortality in the lodgepole pine (*Pinus contorta*) dominated forests of interior BC (Alfaro et al. 2015). MPB is endemic in BC but it is estimated that the most recent and climate-change driven outbreak affected 18.3 million ha of forests (i.e., >20% of BC's total land area) and killed more than half of BC's merchantable lodgepole pine growing stock (>750 million m³) (Walton 2013). This is already having profound negative socioeconomic and environmental effects that are expected to last for decades (Dhar et al. 2016). The allowable annual cut is pro-

jected to fall by about 32% in the next 50 years, resulting in a cumulative present value loss of CAD $57.37 billion (or 1.34%) in GDP and a CAD $90 decline in welfare (Corbett et al. 2016). Moreover, the cumulative effects of MPB outbreaks and wildfires have turned BC forests into a carbon source for the last 2 decades (Giles-Hansen and Wei 2022; Kurz et al. 2008).

Climate change is expected to increase the intensity and severity of abiotic and biotic disturbances (Seidl et al. 2017), alter forest productivity (Boisvenue and Running 2006), and change the distribution of species (Searle and Chen 2017). In parts of BC, forest fires are expected to increase in severity while some forest insects are expected to expand their range northward and upward (Haughian et al. 2012), therefore increasing the risk of forest dieback and economic losses. In addition, the geographic distribution of ecological niches of BC's key ecosystems have already shifted by over 20% since the 1970s (Wang et al. 2012) and are predicted to further expand their potential habitat upward and northward. In some cases, the anticipated rate of habitat expansion is expected
to exceed the natural migration capacity of some tree popu-
lations, presenting a risk of species maladaptation or forest
dieback that would result in uneconomic production, which
could however be potentially mitigated with assisted migra-
tion (Aitken et al. 2008).

Several forest management options for climate change
adaptation involve changing species composition, which may
require selecting species better adapted to major distur-
bances or future climatic conditions (either inside the local
species pool or outside—via assisted migration), switching
from mono-specific to mixed-species forests, and more gen-
erally increasing species diversity (Park et al. 2014; Hof et al.
2017; Messier et al. 2019). The rationale behind using species
diversification (i.e., increasing the number of tree species in a
stand) as a climate change adaptation option relies on the as-
sumption that, because different species usually exhibit vari-
able susceptibility to a given threat (e.g., altered climatic con-
ditions or disturbance regimes), diverse forests are likely to
be less impacted (i.e., more resistant) and will recover faster
(i.e., more resilient) than less diverse forests (Bauhus et al.
2017). Assisted migration can take various forms, including
moving populations (genotypes) within the current species
range, expanding species range (short- to medium-range mi-
gration) or moving species far from their current range (long-
distance migration, introduction of exotic species) (Williams
and Kasten Dumroese 2013). In BC, assisted migration has
been implemented for several years through northward and
upward seed transfers aiming at compensating for past cli-
mate change (Ying and Yanchuk 2006), and seed transfer
guidelines have been further updated to anticipate future cli-
mate change (O'Neill et al. 2017). It is increasingly acknowl-
edged that new planting regimes based on species diversi-
fication, assisted migration or a combination of both have
the potential to increase forest resilience to climate change
and insect outbreaks in North America (Williams and Kasten
Dumroese 2013), and more specifically in BC (Dymond et al.
2014, 2015; Devisscher et al. 2021).

In addition to changes in species composition, several silvi-
cultural practices have been recommended as climate change
adaptation options to increase the resistance and resilience
of stands to natural disturbances, including thinning (Park
et al. 2014; Moreau et al. 2022). Thinning is as an intermedi-
ate treatment where some trees are removed before a final
harvest, which can be used to obtain revenue earlier than
during the final harvest (commercial thinning), to improve
the growth of residual trees, increase wood quality, or ob-
tain certain species or sizes needed for specialty products
(Ashton and Kelty 2018). Thinning has also been advocated
as a possible management option to decrease forest vulnera-
bility to disturbances such as drought, wildfire, insects, and
pathogens (Park et al. 2014; Moreau et al. 2022). Currently,
thinning is rarely implemented in BC public forests (which
cover 95% of the forested area in the Province) due to multi-
ple factors: limited road access, high management costs, and
an unfavorable tenure system that limits silvicultural invest-
ments (Zhang and Pearse 1996; Mitchell et al. 2017). There
is, however, a renewed interest in commercial thinning as
an option to mitigate the mid-term timber supply shortage
induced by the unbalanced age class distribution inherited

from past MPB outbreaks in Western Canada, as it allows to
access timber volume in immature stands or in areas man-
aged for non-timber objectives, and shortens the time to mer-
chantability (Griess et al. 2019; Lieffers et al. 2023).

Assisted migration, species diversification, and thinning
have thus been advocated as options to potentially re-
duce climate change-related risks to timber supply in BC
(Spittlehouse and Stewart 2003; Devisscher et al. 2021), but
their combined effectiveness has not yet been investigated.
A likely reason is the lack of adequate tools that can simu-
late the combined impacts of climate change, natural distur-
bances, and silviculture on timber supply, while still provid-
ing reliable yield estimates that can be compared to the ex-
isting timber supply analyses. One option to overcome this
problem is to combine multiple tools, each specialized for
one component of the problem, into what has been coined as
meta-modelling or toolkit framework (Sturtevant et al. 2007).

We applied a toolkit framework to a case study area in
north-western BC and aimed to assess the potential benefits
of: (i) commercial thinning as a way to mitigate the timber
supply impact of severe insect outbreaks; (ii) novel planting
regimes based on assisted migration and species diversifica-
tion as a way to mitigate future climate change impacts; and
(iii) combined management strategies of commercial thin-
ning and novel planting regimes.

2. Materials and methods

2.1. Study area

The Copper Pine Creek study area is located in the south-
west of Bulkley timber supply area (TSA) in north-western BC
(Figs. 1a and 1b). It covers approximately 96 000 ha, with a
broad elevation range (from 460 to 2500 m) and is transi-
tional between coastal and interior climates. Average annual
temperatures range between 0.5 and 3.1 °C, and average an-
nual precipitation ranges from 500 to 1600 mm (historical
average values for the period 1961–1990; see Hof et al. 2017).
The Copper Pine Creek area encompasses three main bio-
geoclimatic ecosystem classification (BEC; Pojar et al. 1987)
zones and five variants which were used to determine cli-
mate change modifiers of growth and planting regimes (Fig.
1c). At the lowest elevation, we find the Sub-Boreal Spruce
zone (SBS), which is composed of mixed forests dominated
by hybrid spruce (Picea engelmannii x Picea glauca), lodgepole
pine, and subalpine fir (Abies lasiocarpa), with Douglas-fir be-
coming increasingly more prevalent on dryer sites. At inter-
mediate elevations, the Interior Cedar Hemlock (ICH) zone
is dominated by subalpine fir, western hemlock (Tsuga hetero-
phylla), lodgepole pine, hybrid spruce, with Douglas-fir and
western larch (Larix occidentalis) on dryer and mesic variants.
At the highest elevations, the Engelmann Spruce—Subalpine
Fir zone (ESSF) is dominated by subalpine fir, lodgepole pine,
and Engelmann spruce forests. Above 1500 m elevation, the
forest transitions to parkland and alpine zones. More than
66% of the forested area is covered by the ESSF zone, 25% by
the SBS zone, and 10% by the ICH zone. Conifers are highly
dominant in the landscape, with mainly subalpine fir leading
stands (69%), followed by spruce (11%), pine (8%) and hem-
**Fig. 1.** Location of the Copper Pine Creek study area (a, in British Columbia; b, within Bulkley timber supply area) and description of the biogeoclimatic zone variants (c) and land base (d). ESSF = Engelmann spruce—subalpine fir zone; mc = moist cold variant; wv = wet very cold variant; ICH = interior cedar hemlock zone; mc1 = moist cold variant 1; SBS = sub-boreal spruce zone; dk = dry cool variant; mc2 = SBS moist cold variant 2; non-TSA/NF = forested area not included in the Bulkley TSA (e.g., community forests) or non-forested area; THLB = timber harvesting land base; non-THLB = area part of the TSA but not of the THLB. Maps were created using ArcMap 10.05 and assembled from geographical data publicly available on BC Data Catalogue (https://catalogue.data.gov.bc.ca/). All data sources are also available in the dataset published with this study (Lafond et al. 2024).

<figure><img src="image_1.png" /></figure>

lock (5%) leading stands (Fig. 2a). Broadleaves cover 6% of the forested area. The age class distribution (Fig. 2b) is the result of past natural disturbances regimes (e.g., wildfires and MPB outbreaks), combined with a relatively low historic management intensity.

Approximately 66% of the Copper Pine Creek area is non-forested land (e.g., alpine zone, wetlands, water bodies) or non-public forests (e.g., private land, community forests), while the remaining 34% constitutes the productive forested land base where management activities could take place (Fig. 1d; see Table A2). Approximately 18% of the study area was harvest-eligible land, a subset of the productive forested land base after subtracting reserves, physically inoperable areas,

non-commercial stands, future wildlife tree retention areas,
and future right-of-way allowances.

2.2. Toolkit modelling framework

We used a toolkit framework (Sturtevant et al. 2007) com-
bining the best data and stand- or landscape-scale models
available in BC to simulate each component of our forest
management and timber supply analysis problem (Fig. 3). The
framework combines different data sources as well as out-
puts from stand-level growth and yield models and from a
climate-sensitive model into a forest estate model to sched-
ule management actions and optimize timber supply at the
landscape level, as described below (see also Appendix B).
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Area distribution by leading species (a, top) and by age class (b, bottom) in the timber harvest land base (black) and in the rest of the case study area (non-timber harvest land base, grey).</figcaption></figure>

### 2.2.1. Input data and stand stratification

Geographical data included administrative boundaries, forest inventory data, and areas with non-timber management objectives (Fig. 3; Appendix A). These data were overlaid to define the land base (Fig. 1d; see also Table A2) and to link the resulting polygons (about 22 K) to yield curves and management objectives (see below; see also Appendix B).

Similar polygons in terms of existing and future stand characteristics (i.e., management area, site productivity, MPB impacts) were stratified into 198 unique analysis units (AUs; Appendix A), each AU being linked to one yield curve (see 2.2.2; see also Appendix D). Each forest polygon was therefore linked to one current AU that represented the current forest condition, as well as to multiple future AUs to represent all possible future outcomes considered in this study. Finally, each polygon was linked to all applicable management objectives (see 2.2.3; see also Appendix C).

### 2.2.2. Stand level models

Following a classic approach to timber supply analyses in BC and Canada (Penner 2020; Metsaranta et al. 2024), forest growth and yield models were used to derive yield curves (i.e., volume as a function of stand age) that were used as inputs for the forest estate model Patchworks™ (Fig. 3; see Appendix D). In this study, forest growth and yield curves (GYCs) were developed using reference growth and yield models supported by the BC Government (Penner 2020; Metsaranta et al. 2024): variable density yield projection (VDYP; Government of British Columbia 1993), tree and stand simulator (TASS III; Government of British Columbia 2018b) and table interpolation program for stand yield (TIPSY; Government of British Columbia 2018a).

Reference GYC under baseline conditions were generated for each AU (see Appendix D) using VDYP and TIPSY (for natural and planted stands, respectively), which are the two models routinely used to develop yield curves for timber supply analyses in BC (Penner 2020; Metsaranta et al. 2024).
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Schematic description of the toolkit modelling framework displaying information transfer between models (black ovals) via a pipeline approach. GYC = growth and yield curves.</figcaption></figure>

Since RDYP cannot be used to simulate stand response to commercial thinning treatments (Penner 2020), we used TASS III to simulate post-thinning responses for each AU. Following Griess et al. (2019), yield curve modifiers were computed as the age-dependent ratios between yields predicted for unthinned and thinned stands by TASS III for each AU (Fig. D1). These curve modifiers were then used to pro-rate post-CT outcomes of natural and managed yield curves developed with the previous two growth and yield models Fig. 3).

At the time of analysis, none of the three growth and yield models available for BC were able to simulate the impacts of climate change (Penner 2020; Metsaranta et al. 2024). Yield curves under climate change for baseline composition (i.e., species listed in stocking standards) and alternative planting regimes (see Table 2) were developed using curve modifiers derived from the analysis of yield predicted by LANDIS-II in the case study are Hof et al. 2017 (see Appendix A). We analyzed changes in biomass production predicted by Hof et al. (2017) under climate change to derive climate change curve modifiers, which were ratios between predicted yields under baseline and climate change conditions (see Appendix D). Following Dymond et al. (2020), we selected the most severe climate change scenario (HADGEM A1B, +7 °C in temperature and -200 mm in precipitation over 100 years), which was the only scenario that projected negative climate change impacts in Hof et al. (2017) for the whole landscape (results varied within the landscape, with both negative impacts due to

maladaptation and positive impacts due to improved growing conditions) and thus the only one that called for climate change adaptation measures. Although this was the most severe scenario used in Hof et al. (2017), it falls well within the range of current model projections for the area.

We selected two alternative planting regimes designed for climate change adaptation from Hof et al. (2017) : (i) Douglas-fir enrichment, where Douglas-fir is added to the species mixture defined by the stocking standards in the SBSmc2 and ICHmc1 BEC variants as a conservative approach to assisted migration and diversification; and (ii) diversification, where two to four species are added to the species mixture depending on the BEC variant (Table 1). Note that for each planting scenario, species were planted in equal proportion on each 1 ha site on the landscape during the initiation phase and stand dynamics were simulated by the ForCS model component of Landis II, where stand productivity was modulated by interspecific competition.

For each AU, we derived climate change curve modifiers by calculating time-dependent ratios between yields predicted by Landis II under baseline climate and climate change for the baseline planting regime (i.e., stocking standards). Similarly, we calculated curve modifiers for alternative planting regimes under climate change as the ratios between yields predicted for each of the alternate planting regimes and yields predicted for species composition corresponding to stocking standards. Then, the climate change and planting regime curve modifiers (i.e., ratios) were multiplied with the
<table><caption>**Table 1.** Planting regime per biogeoclimatic ecosystem classification (BEC) zone (modified after Hof et al. (2017)) under different scenarios.</caption><thead><tr><th>BEC zone →</th><th>ESSF (mc&wv)</th><th>SBSmc2</th><th>ICHmc1</th><th>SBSdk</th></tr></thead><tbody><tr><td>Planting regime ↓</td><td>Coldest and cold fir-spruce</td><td>Moderate pine</td><td>Moderate fir-Hemlock</td><td>Warmest spruce-pine-aspen</td></tr><tr><td>Baseline (stocking standards)</td><td>Hybrid spruce, subalpine fir</td><td>Hybrid spruce, lodgepole pine</td><td>Hybrid spruce, western hemlock, subalpine fir</td><td>Hybrid spruce, lodgepole pine, Douglas-fir*, Western larch*</td></tr><tr><td>Douglas-fir enrichment</td><td>Stocking standards</td><td>Stocking standards + Douglas-fir</td><td>Stocking standards + Douglas-fir</td><td>Stocking standards</td></tr><tr><td>Diversification</td><td>Stocking standards + Douglas-fir, lodgepole pine, whitebark pine</td><td>Stocking standards + Douglas-fir, subalpine fir, western hemlock, black spruce</td><td>Stocking standards + Douglas-fir, lodgepole pine, western larch, western redcedar</td><td>Stocking standards + subalpine fir, black spruce</td></tr></tbody></table>

Baseline (stocking standards), Douglas-fir enrichment and diversification. Asterisks (*) indicate species that are already considered as a climate change adaptation option under current stocking standards (see Appendix A). ESSF = Engelmann Spruce—Subalpine Fir Zone (includes ESSFmc moist cold variant and ESSFw wet very cold variant); ICHmc1 = Interior Cedar Hemlock Zone, moist cold variant 1; SBS = Sub-Boreal Spruce Zone; SBSmc2 = Sub-Boreal Spruce, moist cold variant 2; SBSdk = Sub-Boreal Spruce, dry cool variant.

reference yield curves (from VDYP or TIPSY) to derived modified yield curves under climate change. Depending on the status of each stand, the switch to modified yield curves was conducted either gradually or fully: (i) gradually to yield curves of existing stands and first generations of future stands (i.e., linear interpolation between 0—no effect in 2017 and 1—full impact of climate change reached by year 100 of the 250-year planning horizon) and (ii) fully to yield curves of second + generations of future stands.

### 2.2.3. Forest estate model

Landscape scale forest development and management over time was simulated using the forest estate model Patchworks™ (Spatial Planning Systems 2020), which allows to simulate forest growth and forestry operations, find optimal management solutions, and predict timber supply under a given set of management objectives, rules, and constraints (Fig. 3; see Appendix B). Patchworks™ is a spatially-explicit forest management optimization model that is widely used for strategic and tactical forest management and timber supply planning while considering operational constraints in Canada. It uses a heuristic simulated annealing algorithm, which allows the consideration of a range of complex objectives while producing close-to-operational reality and close-to-optimal solutions (Rouillard and Moore 2008).

For this study, we modified the forest estate model developed by Griess et al. (2019) for Bulkley TSA. We restricted the geographic extent of the model and input data to match the Copper Pine Creek study area, added new yield curves for the different scenarios (derived from the application of the curve modifiers to the base case curves, as explained above), and modified some modeling assumptions (e.g., minimum harvest age definition; see Appendix B). The full dataset and the forest model developed for the case study area can be found on Borealis (Lafond et al. 2024). We simulated forest management operations (see 2.3.3) and analyzed timber flow over a 250-year planning horizon, which represents two rotations and thus allows the forest to reach equilibrium under the prescribed management. For each simulation, we set the following

constraints: (i) all non-timber objectives in managed areas are met (see Appendix C); (ii) harvest levels do not change by more than ± 10% per decade (each decade compared to previous) throughout the 250 year planning horizon (but exceptions up to ± 20% were allowed for the first 2 decades if it helped reducing the mid-term timber supply shortage); and (iii) timber harvesting land base (THLB) growing stock (i.e., standing volume) does not decline in the last 100 years of the 250 years planning horizon. We did not simulate natural disturbances in the managed forested area, to make sure differences in simulation outputs only reflect differences in management scenarios (see below). Note that this is also consistent with timber supply analyses, where natural disturbances are often not simulated on the THLB explicitly, but their impact are considered through removing buffer areas from the THLB.

### 2.3. Forest management

Baseline forest management in the case study area followed the classic silvicultural system implemented in most of BC's public land: clearcutting with reserves, followed by artificial regeneration using a restricted list of locally authorized species (i.e., stocking standards; Table 1), and no intermediate treatments (i.e., no thinning). Additionally, in this study we investigated alternative management scenarios enabling the implementation of commercial thinning and the use of alternative planting regimes for climate change adaptation (Table 2).

Forest management rules for the baseline scenario (see 2.3.4) simulated in this study were based on management practices described in the latest timber supply review (TSR) for the Bulkley TSA, with some exceptions (see Appendix B). Notably, minimum harvest ages (MHAs) were set at the age of maximum growth (i.e., age of culmination of mean annual increment, CMAI) for each stand to maximize harvest rate.

CT was added as an optional treatment, which could occur 20–10 years before MHA was reached and allow the removal of 40% of the merchantable volume (with a minimum of 40 m³/ha), but would make thinned stands unavailable for final harvest for 40 years (see Appendix B). The decision to
<table><caption>**Table 2.** Scenario description.</caption><thead><tr><th>Scenario</th><th>Climate</th><th>Management</th><th>Commercial thinning</th><th>Planting regime</th></tr></thead><tbody><tr><td>Baseline</td><td>Baseline</td><td>Baseline</td><td>No</td><td>Stocking standards</td></tr><tr><td>CC</td><td>Climate change</td><td>Baseline</td><td>No</td><td>Stocking standards</td></tr><tr><td>CT</td><td>Baseline</td><td>Thinning</td><td>Yes</td><td>Stocking standards</td></tr><tr><td>CC-CT</td><td>Climate change</td><td>Thinning</td><td>Yes</td><td>Stocking standards</td></tr><tr><td>CC-DF</td><td>Climate change</td><td>DF</td><td>No</td><td>Douglas-fir enrichment</td></tr><tr><td>CC-DIV</td><td>Climate change</td><td>Diversification</td><td>No</td><td>Diversification</td></tr><tr><td>CC-DF-CT</td><td>Climate change</td><td>DF &amp; Thinning</td><td>Yes</td><td>Douglas-fir enrichment</td></tr><tr><td>CC-DIV-CT</td><td>Climate change</td><td>Diversification &amp; Thinning</td><td>Yes</td><td>Diversification</td></tr></tbody></table>

CT = commercial thinning, CC = climate change, DF = Douglas-fir, DIV = diversification.

thin was left to the optimization model, with the objective of maximizing harvest rate under the constraints described above.

After a final cut, all stands were planted (Table 1). Under the baseline planting regime, species composition for each BEC zone followed the list of authorized species for each area (i.e., stocking standards). Note that Table 1 describes the generic species composition per BEC variant, while species composition was specific to each AU (see Appendix A, Table A3) and could slightly vary within a BEC variant. Note also that natural regeneration (natural ingress) is also accounted for in yield curves of regenerated stands. Two alternative planting regimes were also considered as climate change adaptation options (Table 1): Douglas-fir enrichment, and diversification. Douglas-fir appeared to be a better candidate for climate change adaptation in the case study area than other species (Hof et al. 2017) because the study area is located at the northern edge of the natural range of Douglas-fir; therefore, it is likely better adapted to future climatic conditions. Furthermore, Douglas-fir has a high economic value and is not a MPB host. The diversification planting regime was based on existing plans for diversification benchmarks, which were developed through a stakeholder consultation, with the objective of maintaining or increasing species diversity and promoting resilient landscapes (Mah et al. 2014). Alternative planting regimes were modeled in this study by replacing the baseline management yield curves by yield curves modified after Hof et al. (2017). See Section 2.2.2 for implementation details.

We compared a total of eight scenarios that resulted from variations in climatic conditions and forest management (Table 2). To assess the impact of climate change, we simulated future timber supply with baseline management under baseline climatic conditions and under climate change. To assess the effect of commercial thinning on the future timber flow, we additionally simulated scenarios with thinning under both climatic conditions. To assess the climate change adaptation impact of alternative planting regimes, we added two regeneration scenarios: Douglas-fir enrichment and diversification. Finally, we combined commercial thinning and alternative planting regimes in two scenarios.

# 3. Results

## 3.1. Baseline management: timber supply shortage and climate change impact

Under baseline management and baseline climatic conditions, simulations projected a reduction in harvested volume by 20%–30% for the next 100 years (Fig. 4a Baseline), which is caused by an unbalanced age distribution that partially resulted from past MPB outbreaks (1975–1985 and 2002–2012; Fig. 2b), as well as by constraints on harvesting age (MHA). Harvest levels were projected to decrease gradually during the first 50 years of the 250-year planning horizon (from 46 413 m³/year at year 0 to 32 760 m³/year at year 50, i.e., a 29.4% reduction), and remained at lowered—and almost constant—levels until year 80, before increasing to 40 365 m³/year on average for the remainder of the simulation time frame (resulting in harvest levels ~13% lower than initial harvest levels). Meanwhile, growing stock in the THLB (Fig. 4b Baseline) was projected to decrease from ~3.7 million m³ at year 0 to a minimum of ~2.9 million m³ by year 50 (–22%), before increasing again to 3.2 million m³ by year 100 (14% lower than initial growing stock). Growing stock did not decline in the final 100 years of the simulation, as a result of the sustainability criterion in the model design.

Under baseline management with climate change, harvest levels were projected to be lower than under baseline climate for the entire planning horizon (Fig. 4a CC). The difference increased over time, from –11.3% in the short-term (first 3 decades) to –16.9% in the mid-term (years 40–100) and to –24% in the long-term (after year 100). While harvest levels under climate change showed a similar mid-term decrease compared to harvest levels under baseline climate (from 40 764 m³/year at year 0 to 27 770 m³/year at year 50, i.e., a 31.9% reduction over the first 50 years), they were not projected to recover in the long-term, leading to a long-term 24.3% timber supply shortage (30 665 m³/year in average after year 100). In the meantime, THLB standing volume (Fig. 4b CC) decreased over time, from 3.7 million m³ at the beginning of the simulation to 2.8 million m³ by year 50 (–24%), and below 2.7 million m³ after year 90–100 (–28%).
**Fig. 4.** Harvest levels (a, top) and growing stock (b, bottom) under the different climate and management scenarios: baseline climate and management (baseline, black); baseline climate and commercial thinning (CT, dashed black); baseline management under climate change (CC, red); commercial thinning under climate change (CC-CT, dashed red); Douglas-fir enrichment under climate change (CC-DF, green); Douglas-fir enrichment and commercial thinning under climate change (CC-DF-CT, dashed green); diversification under climate change (CC-DIV, blue); diversification and commercial thinning under climate change (CC-DIV-CT, dashed blue). Note that all alternative management scenario under climate change must be compared to the CC scenario, not the baseline. See Table 2 for detailed scenario description.

<figure><img src="image_4.png" /></figure>

## 3.2. Commercial thinning for timber supply shortage mitigation

Under baseline climatic conditions (Fig. 4a CT), CT allowed higher harvest levels throughout the 250-years planning horizon as compared to baseline management (by + 5.6%), with a stronger effect in the short- (+8.3% for years 0–30) and mid-term (+8.5% in average for years 40–100; +10.5% for years 50–80) than in the long-term (+3.7% in average after year 100). Yet, because of the unbalanced age structure, harvest levels still showed a 27.8% decrease over time (reduction from 50 113 m³/year at year 0 to 36 200 m³/year for years 40–80). Nevertheless, this reduction was slightly lower than under the baseline management scenario (–29.4%, from 46 413 m³/year at year 0 to

32 760 m³/year for years 40–80). Moreover, these minimum harvest levels were 22% lower than initial levels under baseline management (46 413 m³/year at t0), i.e., if not considering the initial increase in harvest level under the thinning scenario. Still, because of the landbase history, THLB growing stock (Fig. 4b CT) was projected to follow a similar trend than under baseline management, decreasing from 3.7 million m³ to a minimum of 2.8 million m³ for years 30–60 (–25%), before increasing again to 3.1 million m³ by year 100 (17% lower than initial growing stock). Thinning thus allowed the model to partially mitigate the reduction in harvest levels in the mid-term without further reducing the growing stock in the landbase (only 0.1 million m³ lower than without thinning; Fig. 4b CT).
Under severe climate change (Fig. 4a CC-CT), the benefits of commercial thinning scenario were more pronounced than under baseline climate (by + 10.4% on average throughout the planning horizon when compared to baseline management under climate change; Fig. 4a CC). This increase was especially pronounced in the short term (+15.4% to + 24.6% for years 0–30). In the mid term, thinning increased harvest levels by 6.5% only (as compared to + 8.5% under baseline climate; Fig 4a CT), before rising to + 9.7% in the long-term (after year 100). Comparing scenarios CT and CC-CT directly (Fig. 4a) shows that thinning, by allowing a strong increase in harvest levels in the short-term, mitigated the short-term timber supply shortage attributed to climate change, but did not mitigate further reductions (–19.4% in the mid term, –19.2% in the long term). In the meantime, THLB growing stock (Fig. 4b) decreased over time, from 3.7 million m³ at year 0 to a minimum of 2.4 million m³ by year 50 (–34%), before increasing again to 2.5 million m³ after year 70 (–32%). Thinning reduced growing stocks by 0.2–0.3 million m³ when compared to growing stock trend under baseline management and climate change (Fig. 4b CC and CC-CT).

## 3.3. Alternative planting regimes for climate change adaptation

Alternative planting regimes with Douglas-fir enrichment led to higher harvest levels than with current stockings standards under climate change (Fig. 4a CC-DF and CC-DF-CT vs. to CC and CC_CT), with an average of + 21.3% over the whole planning horizon (independently of the thinning treatment). Under baseline harvesting management (Fig. 4a CC-DF), harvest levels increased by 5.0–6.1% in the short-term (years 0–30), 9.8–15.8% in the mid-term (years 40–100), and 27.3% in the long-term (after year 110) compared to baseline management under climate change (Fig. 4a CC).

If combined with thinning (Fig. 4a CC-DF-CT), harvest levels followed a similar trend, with a relatively higher increase in the mid term: they increased by 1.2–2.7% in the short term (years 0–30), 18.6% in average in the mid-term (years 40–100), and 26.3% in the long term (after year 100). Douglas-fir enrichment also allowed mitigating the reduction in THLB growing stock (Fig. 4b CC-DF vs. CC), which lowered from 3.7 million m³ at year 0 to 3.0 million m³ at year 50 (–19%) in the mid term before increasing again to 3.4 million m³ in the long-term (–7% lower than initial stocking), while it kept decreasing (to 2.8 million m³ and 2.7 million m³, respectively) with stocking standards. Here also, allowing commercial thinning treatments led to lower THLB growing stocks (Fig. 4b CC-DF-CT). Yet, the combination with Douglas-fir enrichment allowed the growing stock to partly recover in the long term (as compared to with commercial thinning only; Fig. 4b CC-CT), as it increased again from 2.8 million m³ in year 50 (–24%) to 3.2 million m³ (–14%) by year 100, therefore reaching similar levels than under baseline management and climate scenario (Fig. 4b Baseline).

Alternative planting regimes with diversification led to even higher harvest levels than Douglas-fir enrichment, with an average increase of 43.3% or 44.8% over the whole planning horizon (Fig. 4a CC-DIV and CC-DIV-CT, respectively) when

compared to scenarios with stocking standards under climate change (Fig. 4a CC). Without commercial thinning (Fig. 4a CC-DIV), harvest levels were 6.2–15.5% higher than with stockings standards in the short-term (0–30 years), 32.7% higher in the mid-term (40–100 years), and 57.0% in the long-term (after year 100). If combined with commercial thinning (Fig. 4a CC-DIV-CT), harvest levels followed a similar trend, except for a slight reduction for the first decade (–2.3%). They were on average higher by 2.3% in the short term (years 0–30), but then gradually over time, from 18% at year 40 to 52.4% at year 90 (35.8% on average for the mid-term period), and stabilized at about 54% in the long-term (after year 100).

Even though diversification did not prevent the mid-term decrease in THLB growing stock, it allowed to better mitigate the decrease than with standard stockings or Douglas-fir enrichment, and further allowed to fully recover growing stock in the long-term, reaching or even exceeding initial levels (Fig. 4b). Under baseline management scenario (Fig. 4b CC-DIV), THLB growing stock lowered from 3.7 million m³ to 3.3 million m³ by year 50 (–12%), before increasing to 4.0 million m³ by year 100, thus exceeding initial stocking levels (+9%). Under commercial thinning scenarios (Fig. 4b CC-DIV-CT), THLB growing stocks were constantly lower, and decreased from 3.7 million m³ to 3.1 million m³ by year 50 (–16%), yet they still allowed to fully recover initial stocking levels (3.7 million m³) by year 100.

# 4. Discussion

This study examined the potential benefits of CT and alternative planting regimes as climate change adaptation options in a case study area in BC. We demonstrated that both strategies can contribute to the mitigation of timber supply shortages associated with past MPB outbreaks and future climate change impacts. Combining both options had the highest harvest rates over the planning horizon under a severe climate change scenario due to the additive nature of both interventions.

## 4.1. Mitigating the midterm timber supply shortage with commercial thinning

By allowing timber extraction from immature stands, commercial thinning increased harvest levels over the entire 250-year planning horizon and mitigated the projected reduction in timber supply in the Copper Pine Creek area in interior BC. This is consistent with the simulation study conducted by Griess et al. (2019) in the same area but at a larger scale (for the entire Bulkley TSA), where commercial thinning reduced both the amplitude and duration of the timber supply shortage caused by past MPB outbreaks. Our study additionally projected the effect of a severe climate change scenario and showed its negative impacts on timber supply, which worsened the expected shortage. We also showed that, not only did commercial thinning maintain its timber supply shortage mitigation potential under climate change, but also it had an even stronger positive effect on harvest levels in the short and long term. This is likely an effect of the forest estate model simulating thinning early in the planning hori-
zon to accommodate future climate change associated tim-
ber supply shortage, as well as accessing volume in stands
affected by model constraints (e.g., visual quality and water-
shed protection). Here, thinning provides an opportunity to
access volume early to mitigate both the MPB-induced mid-
term shortage and future volume losses induced by climate
change.

Post-thinning growth response at the tree and stand level is complex to quantify and model. Thinning response is dependent on species composition, the type and intensity of thinning, stand age, and density (Ashton and Kelty 2018; Bose et al. 2018). Because thinning was implemented late in the rotation in this study, it was not expected to trigger significant growth responses at the tree level (Bose et al. 2018), which simplified the modelling of stand growth after thinning and limited simulation biases induced by modelling assumptions. Some studies assume a post-thinning “boost” that increases both tree- and stand-level growth (Franklin et al. 2009; Nabuurs et al. 2000). Other studies assume no impact on volume (i.e., yield curve slopes remain the same, growing stock is simply reduced according to the thinning intensity). Those assumptions can highly impact simulation results and estimated benefits from thinning, but we still lack evidence to support one assumption over the other. In the current study, we likely underestimated the benefits of thinning for climate change adaptation, as the model used thinning as a way to limit future volume losses by accessing volume earlier, but did not simulate the effect of thinning on resistance and resilience to biotic or abiotic disturbances. While there is strong evidence of the benefits of thinning to reduce forest vulnerability to these natural disturbances, thinning effects may vary over time and between contexts, calling for more empirical-based studies (Moreau et al. 2022). Simulating these effects in our modeling framework would require developing yield curves that reflect the increase in growing space combined with stochastic disturbances. Additional research efforts are necessary to predict the probability of and individual tree response to disturbance. More generally, wide experiments are needed to better understand—and therefore better simulate—the combined effects of silviculture (e.g., thinning), climate change, and natural disturbances on forest dynamics (Moreau et al. 2022; Seidl et al. 2017).

## 4.2. Mitigating negative climate change impacts with assisted migration and species diversification

This study also showed the potential benefits of alternative planting regimes focusing on (i) assisted migration and (ii) species diversification. Both options allowed limiting negative climate change effects, allowing to maintain or even increase forest productivity and timber supply under climate change.

Assisted migration, here implemented by adding Douglas-
fir to the set of species used for regeneration in two for-
est types or adding two or more species in the species di-
versification scenarios, appeared as an opportunity for cli-
mate change adaptation in our case study area located in
north-western BC. Although Douglas-fir is already considered

as a candidate species for climate change adaptation in the
warmest forest types (BEC zone) included in this study, we
extended its introduction to colder areas and showed the po-
tential benefits of this proactive climate change adaptation
strategy. In a previous study in this area, Hof et al. (2017)
showed that introducing ponderosa pine was only relevant
in the warmest situations, and that introducing Douglas-fir
or lodgepole pine exhibited more consistent positive effects.
However, given the increased risk of MPB and *Dothistroma* nee-
dle blight outbreaks under climate change in BC (Haughian
et al. 2012), lodgepole pine enrichment does not appear as
a viable climate change adaptation strategy. Instead, actively
replacing the currently dominant MPB-host lodgepole pine
by a mix of species including Douglas-fir was advocated as
a way to increase forest resilience to MPB outbreaks in two
different areas in interior BC (Dymond et al. 2015, 2014).

More generally, the implementation of assisted migration
is often advocated together with favoring mixed species
stands, which are more resistant and/or resilient than mono-
specific forests to disturbances (Park et al. 2014; Messier et
al. 2019). In this study, the species diversification scenario,
which relied on the addition of two to four species to the
current species composition and included the assisted mi-
gration of several other species, yielded the best outcomes
in terms of future timber supply under climate change.
Yet, we did not evaluate the increase in resistance and re-
silience at any scale in this study, and the positive effect
of species diversification on timber supply under climate
change is the combined result of the LANDIS-II yield curves
simulating (i) increased climatic fitness of the species mix-
ture in this scenario and (ii) increased productivity in mixed-
stands in general (i.e., transgressive overyielding), as shown
in Hof et al. (2017). Biodiversity—productivity relationships
have been demonstrated in forests across multiple scales and
biomes (Paquette and Messier 2011; Forrester and Bauhus
2016; Liang et al. 2016). Though empirical and model-based
evidence of overyielding in mixed-stands is clear, the con-
text under which they occur are not, and the processes in-
volved are complex (Forrester and Bauhus 2016). Interactions
between species are driven by neighborhood effects (Fichtner
et al. 2018) that include positive, neutral, and negative in-
teractions. The sum of these interactions is what forms the
basis for the generalizations of stand level to global scale
biodiversity-productivity relationships. Gradients of produc-
tivity across multiple scales as well as the spatial arrange-
ment, stand density, structural diversity, and composition of
species (diversity and evenness) within stands likely influ-
ences the relationship and may even cause a negative rela-
tionship between diversity and productivity. The mixing ef-
fects derived from LANDIS-II assume that stands are fully
stocked with an even initial distribution of species (which
will later evolve as the most competitive species dominate
and develop more biomass), which may overestimate the
net effect of species-mixing in our study compared to in-
situ planted forests in BC. Future studies of timber supply
should aim to more closely examine the interactions between
species mixtures, disturbances, and climate change impacts.
## 4.3. Study limitations and future directions

Like with any modelling study, the reliability and predictive capacity of our toolkit modelling framework is conditioned by several simplifications, assumptions, and methodological limits. Our simulations were for example conditioned by the age constraint implemented for clearcutting (MHA), which is consistent with current forest management practices but might have negatively influenced the projected mid-term timber supply shortage, and thus reinforced the benefits associated with commercial thinning in our study. Future studies might consider removing the MHA constraint to directly compare the effect of thinning versus early clearcuts, as well as relax other constraints (e.g., even-flow) and allow surge harvest to allow post-disturbance salvage logging (Lieffers et al. 2023). However, this will likely have cascading impacts on other ecosystem services like recreation, carbon sequestration, or biodiversity conservation, which would need to be assessed.

Other simplifications and modelling assumptions need to be considered when interpreting the results of this study. The impact of future natural disturbances (MPB outbreaks, wildfires) on the timber supply was not simulated, despite evidence that cumulative effects from increased disturbances regimes have turned BC forests into a carbon source over the last decades (Kurz et al. 2008). In addition, our approach of using outputs from Hof et al. (2017) to modify reference yield curves involved several simplifications, including the use of a single forest model and a single climatic scenario, which ignored any uncertainty associated with modelling assumptions of ecological processes or climate change scenarios, which are known to strongly influence model outputs (Huber et al. 2021). Moreover, we assumed that changes in biomass predicted by Landis II for fully stocked stands can be directly translated into equivalent changes in volume yields predicted by VDYP and TIPSY for natural and planted stands (respectively), thus ignoring interactions of treatment or climate change effects with stand density and competition, as well as differences between modelled processes and outputs. Finally, our modelling framework relies on a series of models feeding into one another, each with their own modeling assumptions, which increases the risk of cumulating error in model projections. Though debatable, this methodological approach was used as a way to overcome the lack of climate sensitive growth and yield models that are both applicable in a large range of ecological conditions and validated by the provincial government for use in timber supply analyses in BC (Penner 2020). This issue is not specific to BC, and several barriers to the development of climate sensitive growth and yield models for Canadian forests and their transition from research to operational tools have been identified, including the complexity of the governmental approval of growth and yield models, which limits both model development and the use of ensemble modelling approaches in timber supply analyses (Metsaranta et al. 2024). There is yet hope that the ongoing national-scale and multi-partner initiative led by the Canadian Forest Service will succeed at fostering the development of climate sensitive growth and yield models, and Metsaranta et al. (2024) made several recommendations to

support this initiative, including the use of open data and open-source code principles. We further recommend that future timber supply studies should not only consider the uncertainty associated with different climate change scenarios, disturbance regimes, and modelling assumptions of ecological processes, but also include an uncertainty analysis to assess the amplitude and relative importance of these different sources of uncertainty (e.g., Huber et al. 2021). In summary, there is a critical need for more data to empirically evaluate model assumptions as well as identify and quantify projection biases (Courbaud et al. 2015).

Finally, though our study showed potential opportunities associated with thinning, diversification, and assisted migration to mitigate timber supply shortages as well as to limit climate change impacts, one of the largest assumptions in our study is that the implementation of these strategies begin at time zero, thanks to the abovementioned “crystal ball” effect. Though our analysis is a snapshot of a potential reality, forests are inherently complex adaptive systems that may shift developmental trajectories in response to both small- or large-scale disturbances (Messier et al. 2015), which are by nature stochastic. In addition, it is more likely that managers will not take a proactive but a reactive approach to management and lose years or decades in opportunities to mitigate these events because of the multiple barriers to the implementation of adaptation measures, in particular commercial thinning. In Western Canada, the public opinion seems to support climate change adaptation (Peterson St-Laurent et al. 2018), but the perceived acceptability of different regeneration strategies seems to decrease with the degree of associated human intervention, from mixed species stands to the introduction of exotic species (Moshofsky et al. 2019). Future studies might examine a delay or consequence of a delayed reaction, the sensitivity of yield projections on timber supply and climate change adaptation as well as proactive management approaches (e.g., fuzzy sets, Ducey and Larson 1999).

## Acknowledgements

The authors would like to thank Kathleen Coupland and Christina Howard for their help with this project and their inputs on an earlier version of the manuscript.

## Article information

### History dates

Received: 20 September 2023

Accepted: 24 May 2024

### Copyright

© 2024 The Author(s). This work is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
## Data availability

Landscape level forest data, growth and yield data, metadata, and models for modelling timber supply are available for download from Borealis, the Canadian Dataverse Repository (https://doi.org/10.5683/SP3/L2Y27G).

## Author information

### Author contributions

Conceptualization: ADP, CDM, CCD, VCG

Data curation: CDM

Formal analysis: VL, ADP, CDM

Funding acquisition: GP, VCG

Methodology: CDM

Project administration: VCG

Resources: CCD

Supervision: VCG

Visualization: VL, CDM

Writing – original draft: VL, ADP

Writing – review & editing: VL, ADP, CDM, CCD, GP, VCG

### Competing interests

The authors declare there are no competing interests.

### Funding information

This work was generously supported by the Department of Natural Resources Forest Innovation Program—Canadian Wood Fibre Center at Natural Resources Canada (CWFC1820-034).

## References

Aitken, S.N., Yeaman, S., Holliday, J.A., Wang, T., and Curtis-McLane, S. 2008. Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol. Appl. **1**: 95–111. doi:10.1111/j.1752-4571.2007.00013.x. PMID: 25567494.

Alfaro, R.I., van Akker, L., and Hawkes, B. 2015. Characteristics of forest legacies following two mountain pine beetle outbreaks in British Columbia, Canada. Can. J. For. Res. **45**(10): 1387–1396. doi:10.1139/cjfr-2015-0042.

Ashton, M.S., and Kelty, M.J. 2018. The practice of silviculture: applied forest ecology. John Wiley & Sons.

Bauhus, J., Forrester, D.I., Gardiner, B., Jactel, H., Vallejo, R., and Pretzsch, H. 2017). Ecological stability of mixed-species forests. In Mixed-species forests: Ecology and management. pp. 337–382. doi:10.1007/978-3-662-54553-9_7

Boisvenue, C., and Running, S.W. 2006. Impacts of climate change on natural forest productivity–evidence since the middle of the 20th century. Global Change Biol. **12**(5): 862–882. doi:10.1111/j.1365-2486.2006.01134.x.

Bose, A.K., Weiskittel, A., Kuehne, C., Wagner, R.G., Turnblom, E., and Burkhart, H.E. 2018. Does commercial thinning improve stand-level growth of the three most commercially important softwood forest types in North America? Forest Ecol. Manag. **409**: 683–693. doi:10.1016/j.foreco.2017.12.008.

Corbett, L.J., Withey, P., Lantz, V.A., and Ochuodho, T.O. 2016. The economic impact of the mountain pine beetle infestation in British Columbia: provincial estimates from a CGE analysis. Forestry, **89**(1): 100–105. doi:10.1093/forestry/cpv042.

Courbaud, B., Lafond, V., Lagarrigues, G., Vieilledent, G., Cordonnier, T., Jabot, F., and De Coligny, F. 2015. Applying ecological model evaluation: lessons learned with the forest dynamics model Samsara2. Ecol. Modell. **314**: 1–14. doi:10.1016/j.ecolmodel.2015.06.039.

Devisscher, T., Spies, J., and Griess, V.C. 2021. Time for change: learning from community forests to enhance the resilience of multi-value forestry in British Columbia, Canada. Land Use Policy, **103**: 105317. doi:10.1016/j.landusepol.2021.105317.

Dhar, A., Parrott, L., and Hawkins, C.D. 2016. Aftermath of mountain pine beetle outbreak in British Columbia: stand dynamics, management response and ecosystem resilience. Forests, **7**(8): 171. doi:10.3390/f7080171.

Ducey, M.J., and Larson, B.C. 1999. A fuzzy set approach to the problem of sustainability. Forest Ecol. Manag. **115**(1): 29–40. doi:10.1016/S0378-1127(98)00433-2.

Dymond, C.C., Tedder, S., Spittlehouse, D.L., Raymer, B., Hopkins, K., McCallion, K., and Sandland, J. 2014. Diversifying managed forests to increase resilience. Can. J. For. Res. **44**(10): 1196–1205. doi:10.1139/cjfr-2014-0146.

Dymond, C.C., Spittlehouse, D.L., Tedder, S., Hopkins, K., McCallion, K., and Sandland, J. 2015. Applying resilience concepts in forest management: a retrospective simulation approach. Forests, **6**(12): 4421–4438. doi:10.3390/f6124377.

Dymond, C.C., Beukema, S., Nitschke, C.R., Coates, K.D., and Scheller, R.M. 2016. Carbon sequestration in managed temperate coniferous forests under climate change. Biogeosciences, **13**(6): 1933–1947. doi:10.5194/bg-13-1933-2016.

Dymond, C.C., Giles-Hansen, K., and Asante, P. 2020. The forest mitigation-adaptation nexus: economic benefits of novel planting regimes. Forest Policy Econ. **113**: 102124. doi:10.1016/j.forpol.2020.102124.

Fichtner, A., Härdtle, W., Bruelheide, H., Kunz, M., Li, Y., and Von Oheimb, G. 2018. Neighbourhood interactions drive overyielding in mixed-species tree communities. Nat. Commun. **9**(1): 1144. doi:10.1038/s41467-018-03529-w. PMID: 29559628.

Forrester, D.I., and Bauhus, J. 2016. A review of processes behind diversity—Productivity relationships in forests. Curr. Forestry Rep. **2**: 45–61. doi:10.1007/s40725-016-0031-2.

Franklin, O., Aoki, K., and Seidl, R. 2009. A generic model of thinning and stand density effects on forest growth, mortality and net increment. Ann. Forest Sci. **66**(8): 815. doi:10.1051/forest/2009073.

Giles-Hansen, K., and Wei, X. 2022. Cumulative disturbance converts regional forests into a substantial carbon source. Environ. Res. Lett. **17**(4): 044049. Environmental Research Letters, **17**(4), p.044049. doi:10.1088/1748-9326/ac5e69.

Government of British Columbia. 2018a. TIPSY, version 4.4 [online]. Forestry Analysis and Inventory Branch. BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/field-forms-and-software/software-download [accessed 09 July 2023].

Government of British Columbia. 2018b. Tree & Stand Simulator (TASS) [online]. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/tree-and-stand-simulator-tass [accessed 09 July 2023].

Government of British Columbia. 1993. Variable Density Yield Projection (VDYP) [online]. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/variable-density-yield-projection-vdypl [accessed 09 July 2023].

Government of British Columbia. 2021. 2121 B.C. Forest Sector—Statistics Summary. Available from https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/forest-industry-economics/economic-state/2021 bc forest sector_-_statistics_summary.pdf [accessed 09 July 2023].

Griess, V.C., Man, C.D., Polinko, A.D., and Spies, J. 2019. Mitigating midterm timber supply shortage using commercial thinning operations. A case study from British Columbia, Canada. Forest Ecol. Manag. **443**: 1–18. doi:10.1016/j.foreco.2019.04.003.

Haughian, S.R., Burton, P.J., Taylor, S.W., and Curry, C. 2012. Expected effects of climate change on forest disturbance regimes in British Columbia. J. Ecosyst. Manag. **13**(1). doi:10.22230/jem.2012v13n1a152.

Hof, A.R., Dymond, C.C., and Mladenoff, D.J. 2017. Climate change mitigation through adaptation: the effectiveness of forest diversification by novel tree planting regimes. Ecosphere, **8**: e01981. doi:10.1002/ecs2.1981.
<figure><img src="image_5.png" /></figure>

Huber, N., Bugmann, H., Cailleret, M., Bircher, N., and Lafond, V. 2021. Stand-scale climate change impacts on forests over large areas: transient responses and projection uncertainties. Ecol. Appl. **31**(4): e02313. doi:10.1002/eap.2313. PMID: 33630399.

Kurz, W.A., Stinson, G., Rampley, G.J., Dymond, C.C., and Neilson, E.T. 2008. Risk of natural disturbances makes future contribution of Canada's forests to the global carbon cycle highly uncertain. Proc. Natl. Acad. Sci. USA, **105**(5): 1551-1555. doi:10.1073/pnas.0708133105.

Lafond, V., Polinko, A., Man, C., Dymond, C.C., Paradis, G., and Griess, V.C. 2024. Replication data for: combining thinning and diverse plantings to adapt to climate-change induced timber supply shortage in British Columbia. Borealis. doi:10.5683/SP3/L2Y27G.

Liang, J., Crowther, T.W., Picard, N., Wiser, S., Zhou, M., Alberti, G., et al. 2016. Positive biodiversity-productivity relationship predominant in global forests. Can. J. For. Res. **53**(6309): aaf8957.

Lieffers, V.J., Pinno, B.D., Johnson, A., Hossain, K.L., and Gooding, T. 2023. Intensive management increases flexibility in managing wood supply. Can. J. For. Res. **53**(7): 544-555. doi:10.1139/cjfr-2022-0317.

Mah, S., and Astridge, K. 2014. Landscape-level ecological tree species benchmarks pilot project: first approximation benchmarks in five British Columbia timber Supply areas. Technical Report-Ministry of Forests, Lands and Natural Resource Operations, British Columbia(082).

Messier, C., Bauhus, J., Doyon, F., Maure, F., Sousa-Silva, R., Nolet, P., et al., 2019. The functional complex network approach to foster forest resilience to global changes. Forest Ecosyst. **6**(1): 1-16. doi:10.1186/s40663-019-0166-2.

Messier, C., Puettmann, K., Chazdon, R., Andersson, K.P., Angers, V.A., Brotons, L., et al., 2015. From management to stewardship: viewing forests as complex adaptive systems in an uncertain world. Conserv. Lett. **8**(5): 368-377. doi:10.1111/conl.12156.

Metsaranta, J.M., Fortin, M., White, J.C., Sattler, D., Kurz, W.A., Penner, M., et al., 2024. Climate sensitive growth and yield models in Canadian forestry: challenges and opportunities. Forest. Chronicle, **100**(1): 1-19. doi:10.5558/tfc2024-005.

Mitchell, S., Larson, B., Griess, V., Simard, S., Delong, D., Nelson, H., et al., 2017. UBC rethink silviculture discussion paper. University of British Columbia, Vancouver 16. doi:10.14288/1.0394382. http://hdl.ha ndle.net/2429/76032.

Moreau, G., Chagnon, C., Achim, A., Caspersen, J., D'Orangeville, L., Sánchez-Pinillos, M., and Thiffault, N., 2022. Opportunities and limitations of thinning to increase resistance and resilience of trees and forests to global change. Forestry, **95**(5): 595-615.

Moshofsky, M., Gilani, H.R., and Kozak, R.A., 2019. Adapting forest ecosystems to climate change by identifying the range of acceptable human interventions in western Canada. Can. J. For. Res. **49**(5): 553-564. doi:10.1139/cjfr-2018-0076.

Nabuurs, G.J., Schelhaas, M.J., and Pussinen, A., 2000. Validation of the European forest Information Scenario Model (EFISCEN) and a projection of Finnish forests. Silva Fennica, **34**(2): 167-179. doi:10.14214/sf.638.

NRCan. 2022. The State of Canada's Forests. Annual Report 2022. Natural Resources Canada, Canadian Forest Service, Ottawa.

O'Neill, G., Wang, T., Ukrainetz, N., Charleson, L., McAuley, L., Yanchuk, A., and Zedel, S., 2017. A proposed climate-based seed transfer system for British Columbia (No. 099). Government of British Columbia, Victoria, BC.

Park, A., Puettmann, K., Wilson, E., Messier, C., Kames, S., and Dhar, A., 2014. Can boreal and temperate forest management be adapted to the uncertainties of 21st century climate change? Crit. Rev. Plant Sci. **33**(4): 251-285. doi:10.1080/07352689.2014.858956.

Paquette, A., and Messier, C. 2011. The effect of biodiversity on tree productivity: from temperate to boreal forests. Global Ecol. Biogeogr. **20**(1): 170-180. doi:10.1111/j.1466-8238.2010.00592.x.

Penner, M. 2020. Growth and yield Modelling Systems for BC: A review. *Forest Analysis Ltd.: Huntsville, ON, Canada*. Technical report prepared for the Government of British Columbia. Ministry of Forests, Lands, Natural Resource Operations and Rural Development, Forest Analysis and Inventory Branch. Available from https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-in-dustry/forestry/stewardship/forest-analysis-inventory/growth-yield/modelreview_finalreport.pdf [assessed 09 July 2023].

Peterson St-Laurent, G., Hagerman, S., Kozak, R., and Hoberg, G. 2018. Public perceptions about climate change mitigation in British Columbia's forest sector. PLoS One, **13**(4): e0195999. doi:10.1371/journal.pone.0195999. PMID: 29684041.

Pojar, J., Klinka, K., and Meidinger, D.V. 1987. Biogeoclimatic ecosystem classification in British Columbia. Forest Ecol. Manag. **22**(1-2): 119-154. doi:10.1016/0378-1127(87)90100-9.

Rouillard, D., and Moore, T. 2008. Patching together the future of forest modelling: implementing a spatial model in the 2009 Romeo Malette Forest Management Plan. Forest. Chronicle, **84**(5): 718-730. doi:10.5558/tfc84718-5.

Scheller, R.M., Domingo, J.B., Sturtevant, B.R., Williams, J.S., Rudy, A., Gustafson, E.J., and Mladenoff, D.J. 2007. Design, development, and application of LANDIS-II, a spatial landscape simulation model with flexible temporal and spatial resolution. Ecol. Modell. **201**(3-4): 409-419. doi:10.1016/j.ecolmodel.2006.10.009.

Searle, E.B., and Chen, H.Y. 2017. Persistent and pervasive compositional shifts of western boreal forest plots in Canada. Global Change Biol. **23**(2): 857-866. doi:10.1111/gcb.13420.

Seidl, R., Thom, D., Kautz, M., Martin-Benito, D., Peltoniemi, M., Vacchiano, G., et al. 2017. Forest disturbances under climate change. Nat. Clim. Change, **7**(6): 395-402. doi:10.1038/nclimate3303.

Spatial Planning Systems. 2020. Spatial Planning Systems [online]. Available from http://www.spatial.ca [accessed 09 July 2023].

Spittlehouse, D., and Stewart, R.B. 2003. Adaptation to climate change in forest management. BC J. Ecosyst. Manag. **4**: 1-11.

Sturtevant, B.R., Fall, A., Kneeshaw, D.D., Simon, N.P., Papaik, M.J., Berninger, K., et al. 2007. A toolkit modeling approach for sustainable forest management planning: achieving balance between science and local needs. Ecol. Soc. **12**(2). doi:10.5751/ES-02102-120207.

Walton, A. 2013. Provincial-Level Projection of the Current Mountain Pine Beetle Outbreak: Update of the Infestation Projection Based on the Provincial Aerial Overview Surveys of Forest Health Conducted from 1999 through 2012 and the BCMPB Model (Year 10). BC Ministry of Forests, Lands and Natural Resource Operations: Victoria, BC, Canada, 2013. [online] Available from http://www.for.gov.bc.ca/ftp/hre/external/ publish/web/bc mpb/year10/BCMPB.v10.BeetleProjectionUpdate.pdf [accessed on 09 July 2023].

Wang, T., Campbell, E.M., O'Neill, G.A., and Aitken, S.N. 2012. Projecting future distributions of ecosystem climate niches: uncertainties and management applications. Forest Ecol. Manag. **279**: 128-140. doi:10.1016/j.foreco.2012.05.034.

Williams, M.I., and Dumroese, R.K. 2013. Preparing for climate change: forestry and assisted migration. J. Forestry, **111**(4): 287-297. doi:10.5849/jof.13-016.

Ying, C.C., and Yanchuk, A.D. 2006. The development of British Columbia's tree seed transfer guidelines: purpose, concept, methodology, and implementation. Forest Ecol. Manag. **227**: 1-13. doi:10.1016/j.foreco.2006.02.028.

Zhang, D., and Pearse, P.H. 1996. Differences in silvicultural investment under various types of forest tenure in British Columbia. Forest Science **42**(4): 442-449. doi:10.1093/forestscience/42.4.442.

## Appendix A. Forest data and stand stratification into analysis units

### Geographical data

Geographical data included all administrative boundaries, forest inventory data, and areas with non-timber management objectives listed in Table A1. These data were overlaid and consolidated into a resultant file free of topological errors (21, 903 unique spatially-defined polygons), which served as one of the key inputs to the forest estate model. The resultant file was used to define the land base (i.e., non-forest, productive land base, and harvest-eligible land base; see Table A2) and to link each polygon to yield curves and management objectives (see Appendix B and Appendix D). All
<figure><img src="image_6.png" /><figcaption>**Table A1.** Data source names and descriptions.</figcaption></figure>

<table><thead><tr><th>Feature class name</th><th>Description</th></tr></thead><tbody><tr><td>Pine_creek_bdry</td><td>Copper Pine Creek boundary—study area boundary</td></tr><tr><td>Own</td><td>Ownership, used to identify lands excluded from the analysis (private, community forest, non-public lands)</td></tr><tr><td>IR</td><td>Indian reserves, excluded from the analysis</td></tr><tr><td>SRMP_BulkleyValley</td><td>Bulkley valley strategic resource management plan boundary</td></tr><tr><td>LUP_Consolidated</td><td>Consolidated land use plans boundaries used for non-timber objectives—landscape corridors and core ecosystems</td></tr><tr><td>LU</td><td>Landscape units, administrative units used for non-timber objectives—landscape-level biodiversity, landscape corridors, core ecosystems</td></tr><tr><td>Watershed_FishSensitive</td><td>Fisheries sensitive watersheds, used to set-up non-timber objectives for fisheries sensitive watersheds</td></tr><tr><td>Rec_Poly</td><td>Recreation polygons, reserves deferred from harvesting</td></tr><tr><td>TSM</td><td>Terrain stability mapping, includes predictive ecosystem mapping (PEM), used to identify unstable terrain, red and blue listed species, and define analysis units (link to yield curves).</td></tr><tr><td>BEC</td><td>Biogeoclimatic zones, used to set-up landscape-level biodiversity objectives and implement wildfire disturbances.</td></tr><tr><td>VRI</td><td>Vegetation resource inventory, used to define the land base (forestd and non-forested), low sites and problem forest types, and analysis units (link to yield curves)</td></tr><tr><td>blocks</td><td>Consolidated cutblocks of previous harvest not included in VRI. Used to complement VRI for recent harvest events.</td></tr><tr><td>FWA_water</td><td>Fresh water atlas—lakes, wetlands, rivers, stream network. Used to complement VRI to identify non-forested land base and prepare riparian buffer reserves.</td></tr><tr><td>VLI</td><td>Visual landscape inventory, used to set-up visual quality objectives</td></tr><tr><td>Rec_Lines_buff</td><td>Buffered recreation lines, reserves deferred from harvesting</td></tr><tr><td>Roads</td><td>Buffered roads, used to complement VRI to identify non-forested land base.</td></tr><tr><td>riparian_buffer</td><td>Consolidated buffered riparian (lakes, wetlands, rivers, stream network), deferred from harvesting.</td></tr><tr><td>Resultant_v1_Cpine</td><td>Resultant for copper pine creek that consolidated all above datasets and was used to define the land base, analysis units (link to yield curves), post-disturbance outcomes, and management objectives. The productive forested land base was transferred to patchworks.</td></tr></tbody></table>

data sources were publicly accessed from Data BC and compiled in August, 2017.

## Land base definitions

Table A2 identifies elements of the landscape that were removed from the total area to define the productive forested land base (PFLB), from the PFLB to define the timber harvesting land base (THLB), and from the THLB to define the long-term THLB. Layers of spatial data described in Table B1 were used to identify non-crown (i.e., non-public) lands (e.g., community forests, Indian reserves), non-forested lands (e.g., roads, waterbodies), protected areas (e.g., parks, reserves, wildlife and riparian habitat), unstable terrain, and low productivity sites.

## Stand stratification into analysis units

Analysis units (AUs) represent areas covered by similar stands. The definition of similarity between stands depends on the timber characteristics and non-timber objectives, the forest manager, or the analyst. Typically, some stratification criteria are used to aggregate similar stands into AUs to simplify the analysis and ease the interpretation of results. In this analysis, the criteria to stratify AUs included: management era, BEC variant, and site series (Table A3). The management area included four categories: (1) existing natural stands—any stand with no logging history or any stand 40 years and older

(AU < 100), (2) existing managed (i.e., artificially regenerated) stands—stands younger than 40 years and with logging history (AU 200 series), (3) existing managed stands impacted by Dothistroma needle blight (AU 300 series), and (4) future managed stands (AU 100 series).

The BEC variant and site series were sourced from the predictive ecosystem mapping, a type of expert system applied to large-scale classification and mapping of ecosystems across BC but with higher resolution than forest cover inventory.
<table><caption>**Table A2. Land base definition.**</caption><thead><tr><th>Factor</th><th>Total area (ha)</th><th>Effective area (ha)</th><th>% of total area</th><th>% of PFLB</th></tr></thead><tbody><tr><td>Total area</td><td>96 244</td><td></td><td>100.0%</td><td></td></tr><tr><td>Less:</td><td></td><td></td><td></td><td></td></tr><tr><td>Community forests</td><td>32 537</td><td>32 537</td><td>33.8%</td><td></td></tr><tr><td>Woodlots</td><td>709</td><td>709</td><td>0.7%</td><td></td></tr><tr><td>Private</td><td>10 049</td><td>10 049</td><td>10.4%</td><td></td></tr><tr><td>Alpine</td><td>9498</td><td>5160</td><td>9.9%</td><td></td></tr><tr><td>Water bodies</td><td>5238</td><td>2388</td><td>5.4%</td><td></td></tr><tr><td>Non-vegetated</td><td>1578</td><td>840</td><td>1.6%</td><td></td></tr><tr><td>Vegetated not treed</td><td>8383</td><td>3820</td><td>8.7%</td><td></td></tr><tr><td>Non-FMLB</td><td>28 680</td><td>6906</td><td>29.8%</td><td></td></tr><tr><td>Wetlands</td><td>1376</td><td>186</td><td>1.4%</td><td></td></tr><tr><td>Existing Roads</td><td>3113</td><td>857</td><td>3.2%</td><td></td></tr><tr><td><strong>Productive forested land base (PFLB)</strong></td><td></td><td><strong>32 792</strong></td><td><strong>34.1%</strong></td><td><strong>100.0%</strong></td></tr><tr><td>Less:</td><td></td><td></td><td></td><td></td></tr><tr><td>Crown (public) reserves</td><td>111</td><td>111</td><td>0.1%</td><td>0.3%</td></tr><tr><td>Special management zone 1</td><td>74</td><td>74</td><td>0.1%</td><td>0.2%</td></tr><tr><td>Unstable terrain</td><td>730</td><td>730</td><td>0.8%</td><td>2.2%</td></tr><tr><td>Low sites (spruce)</td><td>107</td><td>106</td><td>0.1%</td><td>0.3%</td></tr><tr><td>Low sites (hemlock)</td><td>571</td><td>542</td><td>0.6%</td><td>1.7%</td></tr><tr><td>Low sites (fir)</td><td>8378</td><td>8071</td><td>8.4%</td><td>24.6%</td></tr><tr><td>Problem forest types</td><td>2120</td><td>2034</td><td>2.1%</td><td>6.2%</td></tr><tr><td>Recreation reserves</td><td>1285</td><td>843</td><td>0.9%</td><td>2.6%</td></tr><tr><td>Mountain goat</td><td>1316</td><td>375</td><td>0.4%</td><td>1.1%</td></tr><tr><td>Red and blue listed</td><td>12</td><td>3</td><td>0.0%</td><td>0.0%</td></tr><tr><td>Riparian</td><td>2913</td><td>1631</td><td>1.7%</td><td>5.0%</td></tr><tr><td><strong>Timber harvesting land base (THLB)</strong></td><td></td><td><strong>18 270</strong></td><td><strong>19.0%</strong></td><td><strong>55.7%</strong></td></tr><tr><td>Less:</td><td></td><td></td><td></td><td></td></tr><tr><td>Wildlife tree betention (1%-34%)</td><td>599</td><td>0.6%</td><td>1.8%</td><td></td></tr><tr><td>Future roads (3.9%)</td><td></td><td>713</td><td>0.7%</td><td>2.2%</td></tr><tr><td><strong>Long-term THLB</strong></td><td></td><td><strong>16 959</strong></td><td><strong>17.6%</strong></td><td><strong>51.7%</strong></td></tr></tbody></table>

<table><caption>**Table A3. Analysis unit definitions.**</caption><thead><tr><th>AU</th><th>Management era</th><th>BEC</th><th>Site series</th><th>Moisture class</th><th>Species composition*</th><th>MHA</th><th>Future AU</th><th>THLB (ha)</th><th>Non-THLB (ha)</th></tr></thead><tbody><tr><td>6</td><td>Natural</td><td>ESSFmc</td><td>02, 03</td><td>Dry</td><td>Bl50 Sx25 Pl22 At3</td><td>140</td><td>106</td><td>53</td><td>10</td></tr><tr><td>7</td><td>Natural</td><td>ESSFmc</td><td>01, 04</td><td>Fresh</td><td>Bl69 Sx22 Pl8 At1</td><td>150</td><td>107</td><td>2, 177</td><td>1, 542</td></tr><tr><td>8</td><td>Natural</td><td>ESSFmc</td><td>06, 07</td><td>Moist</td><td>Bl73 Sx25 Pl2</td><td>140</td><td>108</td><td>184</td><td>113</td></tr><tr><td>9</td><td>Natural</td><td>ESSFmc</td><td>08, 09, 10</td><td>Wet</td><td>Bl74 Sx23 Pl3</td><td>150</td><td>109</td><td>33</td><td>87</td></tr><tr><td>14</td><td>Natural</td><td>ESSFwv</td><td>02, 03</td><td>Dry</td><td>Bl56 Sx23 Pl14 Hw6 At1</td><td>150</td><td>114</td><td>1, 093</td><td>588</td></tr><tr><td>15</td><td>Natural</td><td>ESSFwv</td><td>01, 04</td><td>Fresh</td><td>Bl72 Hw12 Hw12 Pl3 Ba1</td><td>150</td><td>115</td><td>4, 259</td><td>5, 837</td></tr><tr><td>16</td><td>Natural</td><td>ESSFwv</td><td>05, 06</td><td>Moist</td><td>Bl72 Sx15 Hw9 Ba2 Ba2</td><td>140</td><td>116</td><td>2, 487</td><td>1, 776</td></tr><tr><td>17</td><td>Natural</td><td>ESSFwv</td><td>07, 09</td><td>Wet</td><td>Bl68 Sx20 Hw4 Hw4 Hm3 At1</td><td>140</td><td>117</td><td>146</td><td>304</td></tr><tr><td>18</td><td>Natural</td><td>ICHmc1</td><td>02</td><td>Dry</td><td>Pl34 Bl27 Hw17 Sx14 At6 Ep2</td><td>90</td><td>118</td><td>5</td><td>47</td></tr><tr><td>19</td><td>Natural</td><td>ICHmc1</td><td>01, 03, 04, 05</td><td>Fresh</td><td>Bl33 Hw32 Sx22 Pl8 At4 Ba1</td><td>130</td><td>119</td><td>1, 445</td><td>990</td></tr><tr><td>20</td><td>Natural</td><td>ICHmc1</td><td>06</td><td>Wet</td><td>Bl56 Sx16 Hw15 Pl13</td><td>130</td><td>120</td><td>19</td><td>44</td></tr><tr><td>25</td><td>Natural</td><td>MHmmp</td><td>01, 03, 05</td><td>Fresh</td><td>Bl71 Hw10 Ba7 Hm6 Sx5 Pl1</td><td>150</td><td>125</td><td>8</td><td>24</td></tr><tr><td>28</td><td>Natural</td><td>SBSdk</td><td>02, 03</td><td>Dry</td><td>Sx54 Pl25 At20 Bl1</td><td>110</td><td>128</td><td>8</td><td>5</td></tr><tr><td>29</td><td>Natural</td><td>SBSdk</td><td>01, 05, 06</td><td>Fresh</td><td>Sx48 Pl26 At17 Bl5 Ep2 Ep2</td><td>100</td><td>129</td><td>693</td><td>1, 634</td></tr><tr><td>30</td><td>Natural</td><td>SBSdk</td><td>07, 08</td><td>Moist</td><td>Sx54 At24 Pl16 Bl5 Ep1</td><td>110</td><td>130</td><td>25</td><td>195</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Table A3. (concluded).</figcaption></figure>

<table><thead><tr><th>AU</th><th>Management era</th><th>BEC</th><th>Site series</th><th>Moisture class</th><th>Species composition*</th><th>MHA</th><th>Future AU</th><th>THLB (ha)</th><th>Non-THLB (ha)</th></tr></thead><tbody><tr><td>33</td><td>Natural</td><td>SBSmc2</td><td>01, 03, 05, 06, 09</td><td>Fresh</td><td>Sx37 Bl36 Pl20 At4 Hw2 Ba1</td><td>120</td><td>133</td><td>2, 190</td><td>908</td></tr><tr><td>34</td><td>Natural</td><td>SBSmc2</td><td>07, 10</td><td>Moist</td><td>Sx41 Bl31 Pl25 At3</td><td>130</td><td>134</td><td>264</td><td>136</td></tr><tr><td>35</td><td>Natural</td><td>SBSmc2</td><td>12</td><td>Wet</td><td>Sx53 Bl24 Pl22 At1</td><td>130</td><td>135</td><td>26</td><td>2</td></tr><tr><td>206</td><td>Managed</td><td>ESSFmc</td><td>02, 03</td><td>Dry</td><td>Pl80 Bl20</td><td>100</td><td>206</td><td>67</td><td>1</td></tr><tr><td>207</td><td>Managed</td><td>ESSFmc</td><td>01, 04</td><td>Fresh</td><td>Sw40 Bl40 Pl20</td><td>110</td><td>207</td><td>493</td><td>14</td></tr><tr><td>208</td><td>Managed</td><td>ESSFmc</td><td>06, 07</td><td>Moist</td><td>Sw45 Bl35 Pl20</td><td>120</td><td>208</td><td>42</td><td>5</td></tr><tr><td>214</td><td>Managed</td><td>ESSFwv</td><td>02, 03</td><td>Dry</td><td>Pl80 Bl20</td><td>130</td><td>214</td><td>378</td><td>40</td></tr><tr><td>215</td><td>Managed</td><td>ESSFwv</td><td>01, 04</td><td>Fresh</td><td>Sw50 Bl30 Pl20</td><td>120</td><td>215</td><td>276</td><td>55</td></tr><tr><td>216</td><td>Managed</td><td>ESSFwv</td><td>05, 06</td><td>Moist</td><td>Sw50 Bl30 Pl20</td><td>130</td><td>216</td><td>199</td><td>17</td></tr><tr><td>217</td><td>Managed</td><td>ESSFwv</td><td>07, 09</td><td>Wet</td><td>Bl60 Sw40</td><td>180</td><td>217</td><td>21</td><td>0</td></tr><tr><td>229</td><td>Managed</td><td>SBSdk</td><td>01, 05, 06</td><td>Fresh</td><td>Pl80 Sw20</td><td>70</td><td>229</td><td>258</td><td>3</td></tr><tr><td>230</td><td>Managed</td><td>SBSdk</td><td>07, 08</td><td>Moist</td><td>Sw80 Pl20</td><td>80</td><td>230</td><td>10</td><td>4</td></tr><tr><td>233</td><td>Managed</td><td>SBSmc2</td><td>01, 03, 05, 06, 09</td><td>Fresh</td><td>Pl50 Sw35 Bl15</td><td>90</td><td>233</td><td>936</td><td>59</td></tr><tr><td>234</td><td>Managed</td><td>SBSmc2</td><td>07, 10</td><td>Moist</td><td>Sw50 Bl30 Pl20</td><td>90</td><td>234</td><td>81</td><td>15</td></tr><tr><td>318</td><td>Dothistroma delay 10 years</td><td>ICHmc1</td><td>02</td><td>Dry</td><td>Pl80 Hw20</td><td>80</td><td>218</td><td>23</td><td>3</td></tr><tr><td>319</td><td>Dothistroma delay 10 years</td><td>ICHmc1</td><td>01, 03, 04, 05</td><td>Fresh</td><td>Sw50 Hw20 Bl15 Pl15</td><td>90</td><td>219</td><td>371</td><td>59</td></tr><tr><td>320</td><td>Dothistroma delay 10 years</td><td>ICHmc1</td><td>06</td><td>Wet</td><td>Sw55 Bl25 Hw20</td><td>100</td><td>220</td><td>1</td><td>1</td></tr><tr><td>106</td><td>Future</td><td>ESSFmc</td><td>02, 03</td><td>Dry</td><td>Pl80 Bl20</td><td>100</td><td>106</td><td>0</td><td>0</td></tr><tr><td>107</td><td>Future</td><td>ESSFmc</td><td>01, 04</td><td>Fresh</td><td>Sw40 Bl40 Pl20</td><td>110</td><td>107</td><td>0</td><td>0</td></tr><tr><td>108</td><td>Future</td><td>ESSFmc</td><td>06, 07</td><td>Moist</td><td>Sw45 Bl35 Pl20</td><td>120</td><td>108</td><td>0</td><td>0</td></tr><tr><td>109</td><td>Future</td><td>ESSFmc</td><td>08, 09, 10</td><td>Wet</td><td>Sw55 Bl45</td><td>140</td><td>109</td><td>0</td><td>0</td></tr><tr><td>114</td><td>Future</td><td>ESSFwv</td><td>02, 03</td><td>Dry</td><td>Pl80 Bl20</td><td>130</td><td>114</td><td>0</td><td>0</td></tr><tr><td>115</td><td>Future</td><td>ESSFwv</td><td>01, 04</td><td>Fresh</td><td>Sw50 Bl30 Pl20</td><td>120</td><td>115</td><td>0</td><td>0</td></tr><tr><td>116</td><td>Future</td><td>ESSFwv</td><td>05, 06</td><td>Moist</td><td>Sw50 Bl30 Pl20</td><td>130</td><td>116</td><td>0</td><td>0</td></tr><tr><td>117</td><td>Future</td><td>ESSFwv</td><td>07, 09</td><td>Wet</td><td>Bl60 Sw40</td><td>180</td><td>117</td><td>0</td><td>0</td></tr><tr><td>118</td><td>Future</td><td>ICHmc1</td><td>02</td><td>Dry</td><td>Pl80 Hw20</td><td>70</td><td>118</td><td>0</td><td>0</td></tr><tr><td>119</td><td>Future</td><td>ICHmc1</td><td>01, 03, 04, 05</td><td>Fresh</td><td>Sw50 Hw20 Bl15 Pl15</td><td>80</td><td>119</td><td>0</td><td>0</td></tr><tr><td>120</td><td>Future</td><td>ICHmc1</td><td>06</td><td>Wet</td><td>Sw55 Bl25 Hw20</td><td>90</td><td>120</td><td>0</td><td>0</td></tr><tr><td>125</td><td>Future</td><td>MHmmp</td><td>01, 03, 05</td><td>Fresh</td><td>Ba70 Hm15 Bl15</td><td>160</td><td>125</td><td>0</td><td>0</td></tr><tr><td>128</td><td>Future</td><td>SBSdk</td><td>02, 03</td><td>Dry</td><td>Pl80 Sw20</td><td>90</td><td>128</td><td>0</td><td>0</td></tr><tr><td>129</td><td>Future</td><td>SBSdk</td><td>01, 05, 06</td><td>Fresh</td><td>Pl80 Sw20</td><td>70</td><td>129</td><td>0</td><td>0</td></tr><tr><td>130</td><td>Future</td><td>SBSdk</td><td>07, 08</td><td>Moist</td><td>Sw80 Pl20</td><td>80</td><td>130</td><td>0</td><td>0</td></tr><tr><td>133</td><td>Future</td><td>SBSmc2</td><td>01, 03, 05, 06, 09</td><td>Fresh</td><td>Pl50 Sw35 Bl15</td><td>90</td><td>133</td><td>0</td><td>0</td></tr><tr><td>134</td><td>Future</td><td>SBSmc2</td><td>07, 10</td><td>Moist</td><td>Sw50 Bl30 Pl20</td><td>90</td><td>134</td><td>0</td><td>0</td></tr><tr><td>135</td><td>Future</td><td>SBSmc2</td><td>12</td><td>Wet</td><td>Sw100</td><td>140</td><td>135</td><td>0</td><td>0</td></tr><tr><td>9000</td><td>other non THLB</td><td>ALL</td><td>ALL</td><td>ALL</td><td>ALL</td><td>0</td><td>9000</td><td>0</td><td>3</td></tr></tbody><tfoot><tr><td colspan="10">Total</td><td>18, 270</td><td>14, 521</td></tr></tfoot></table>

Thus, for finer scale estimates of site productivity, four moisture classes, relative to BEC variant and site series (a classification system of moisture gradient), were used in this analysis.

Each AU was attributed a specific yield curve, derived from reference growth and yield models for BC (VDYP and TIPSY for natural and planted stands, respectively) and optionally modelled to reflect the impacts of thinning and climate change

(see Appendix D). Each forest polygon was therefore linked to one current AU that represented the current forest condition, as well as to multiple future AUs to represent all possible future outcomes considered in this study (each with its specific yield curve).

Finally, each polygon was linked to all applicable management objectives, i.e., policies related to timber flow and con-

<footnote>* Species composition is represented as species code (letters, see Table C2) followed by percentage composition (i.e., Pl50 means 50% of Lodgepole pine). AU = analysis unit; BEC = biogeoclimatic zones (see Fig. 1-c). MHA = minimum harvest age; THLB = timber harvesting land base.</footnote>
<figure><img src="image_8.png" /></figure>

servation measures to maintain landscape-level biodiversity, fish and wildlife habitat, and visual quality (see 2.2.3; see also **Appendix C**). Given that the management objectives designated areas for specific management objectives overlapped, it was possible for one polygon to contribute to multiple management objectives (Figure C1).

# Appendix B. Forest estate modelling with patchworks

## Forest estate model description

Patchworks™ is a spatially-explicit heuristic (simulated annealing algorithm) forest estate model that addresses a range of complex objectives to converge to solutions that are close-to-optimal and close-to-operational reality. Unlike classic linear programming optimization forest estate models, patchworks will always produce solutions from a relatively complex matrix of aspatial and spatial objectives. Patchworks can be set-up both, as an area- and unit-restricted model, linked to a road network to solve transportation problems, and with maximum formulation flexibility. Each spatially-explicit polygon can be programmed to follow one a-priori established development path (e.g., only clearcutting), or multiple development paths (e.g., possibly thinning followed by clearcutting, possibly only clearcutting followed by thinning of future stands and then clearcutting), or a combination of the two. Patchworks is widely used by professional foresters to support management decisions and prepare forest management plans.

## Input data and simulation setup

Like any of the spatially-explicit forest estate models, patchworks requires a spatial dataset that links each polygon to yield curves (existing and future) and management objectives (timber and non-timber), the actual yield curves and definition of management objectives, and definition of silvicultural systems and other disturbance events (i.e., when and where they can occur and post-disturbance outcomes). In this study, the spatial dataset was the productive forested land base subset of the resultant file obtained when combining different geographical data (administrative boundaries, forest inventory data, and areas with non-timber management objectives). Similar forest polygons were stratified into 198 unique analysis units, and each analysis unit was linked to one yield curve per scenario, which were derived from reference models used in BC (see Section 2.2). Disturbance events simulated for the harvest-eligible productive land base were: (1) clearcutting followed by artificial regeneration using a restricted list of locally authorized species (stocking standards) and no intermediate treatments (i.e., no commercial thinning); (2) commercial thinning followed by clearcutting. See below, paragraph on forest management.

A set of management objectives were defined for each polygon. Timber objectives were programmed to produce the maximum sustainable harvest rate possible throughout the planning horizon. Thus, given a surplus of the older stands, it was possible to allocate a relatively high harvest rate at the start of the planning horizon, then decline, by down to 10% per decade, to a lower mid-term harvest rate, before increasing, by up to 10% per decade, to a long-term even-flow harvest rate. In addition, the growing stock on the harvest-eligible land base was maintained in a non-declining state over the last 100 years of the planning horizon. Non-timber objectives either led to areas being set as reserves (i.e., removed from the
<figure><img src="image_9.png" /></figure>

land base eligible to harvest) or were implemented as harvest
constraints (see Appendix C).

For all scenarios, the planning horizon was set to 250 years,
25 periods of 10-year length each. The following approach
was implemented to find a close-to-optimal solution: (1) set
a relatively high harvest target (i.e., approximately 1.5 times
the theoretical long-run sustained yield subject to each stand
being harvested at the age of maximum growth) and let the
model run for 1 million iterations (allowed the heuristic al-
gorithm to freely explore the harvest options to avoid local
optimums at later simulation stages); (2) activated all objec-
tives (except non-declining growing stock) and let the model
reach a relatively stable maximum harvest rate (typically un-
evenly distributed over the planning horizon); (3) activated
the merchantable growing stock objective—non-declining in
the last 100 years of the 250-year planning horizon as a sur-
rogate for harvest sustainability (i.e., harvest equals growth);
and (4) manually fit the best harvest flow (maximum ini-
tial with a lower mid-term and a non-declining long-term,
±10%/decade change) and let the model run until the over-
all objective function improvement between 500 000 consec-
utive iterations was < 0.000000001%.

The solving time ranged between 24 and 48 h (100 + mil-
lion iterations) on a computer with Windows 7 (64-bit), Intel
Core i7 (3.4 GHz), and 32 GB RAM.

Forest management

All forest management operations simulated in this study
follow the dominant silvicultural system in BC, which is
clearcutting with reserves. Forest management rules for the
baseline scenario simulated in this study were based on man-
agement practices described in the latest TSR document for
the Bulkley TSA (FLNRORD 2012), with some exceptions. The
main difference relates to the rules used to determine the
MHA : in the TSR, the MHA for clearcuts was determined
as the age at which a minimum merchantable volume of
150 m³/ha was reached (merchantability levels were set to
12.5 cm diameter at breast height (DBH) for pine and 17.5 cm
DBH for all other species); in this study, the MHA was set at
the age of CMAI, i.e., the age of maximum growth, with the
objective of maximizing harvest rate. In addition, and con-
trary to the previous version of the forest estate model (Griess
et al. 2019), the CMAI was determined annually in this study
(i.e., CMAI was determined from an annual yield curve using
linear interpolation, as opposed to decadal yield curve). This
change was meant to provide a more accurate estimate of the
age of CMAI, and a more accurate representation of the op-
erability windows (harvesting and thinning). Moreover, the
age of CMAI and operability windows were recalculated ev-
ery time a yield curve was changed (i.e., harvest, thinning,
climate change). This dynamic change of operability windows
allowed for a better comparison between scenarios.

Following a clearcutting event, the regeneration assump-
tions described in Table 1 were applied, with various species
mixture depending on the BEC zone and scenario: under
baseline management, they followed stocking standards in
terms of species composition and density (Government of
British Columbia, n.d.); under the Douglas-fir enrichment

and diversification scenarios, the species composition was al-
tered to add one or several species (respectively) to stocking
standards.

An optional commercial thinning treatment was added to some scenarios. The window of opportunity for commercial thinning was defined as the age range 20–10 years prior to the age at which MHA was reached. The CT intensity was set to a removal of 40% of merchantable volume accounting for a minimum of 40 m³/ha to ensure stand stability as well as minimum economic requirements (Ministry of Forests 1999a, 1999b). The growth of residual stand was adjusted with the corresponding post-thinning factor curves. A thinned stand was excluded from harvesting for 40 years following commercial thinning, to allow for stands to respond to the treatment.

**Appendix C. Non-timber objectives**

Non-timber objectives can be met by either setting aside
some area, which is thus excluded from the managed area
eligible for harvest (THLB), or implementing specific man-
agement practices in managed areas, here in the form of
harvest constraints built in the forest estate model. Har-
vest constraints are defined by minimum or maximum val-
ues associated with a stand metric (e.g., age, area) or an in-
dex (e.g., equivalent clearcut area). The following non-timber
objectives were considered in the case study area (Fig. C1)
: landscape-level biodiversity objectives, fisheries sensitive
watershedsual quality objectives (VQO), and an ecosystem
network composed of landscape corridors (LCOR) and core
ecosystems (COEC).

Landscape-level biodiversity

The landscape level biodiversity objectives were imple-
mented to maintain certain levels of young, mature, and
old seral stands within each unique combination of land-
scape unit, biogeoclimatic ecosystem classification (BEC) sub-
zone, and natural disturbance type (Table C1, adapted from
latest timber supply review of Bulkley timber supply area),
following the British Columbia Forest Practices Code (1995).
Landscape units were introduced since 1995 as key manage-
ment units to ensure consistent biodiversity and manage-
ment across BC: 5, 000 to 100, 000 ha in size including pri-
vate and public lands, defined by ecological boundaries, wa-
tersheds or groups of watersheds rather than administrative
boundaries (Fenger 1996). In each landscape unit, biodiver-
sity conservation targets are defined for each BEC zone and
natural disturbance type, following recommendations from
the in the Biodiversity Guidebook (British Columbia Forest
Service and British Columbia Environment 1995). There are
five natural disturbance types (NDT) in BC, defined by the
frequency of disturbances and their impact on forest stand
structure: the first three types describe ecosystems with rare
(NDT1), infrequent (NDT2), or frequent (NDT3) stand initi-
ating disturbances; NDT4 describes ecosystems with stand-
maintaining wildfires; and NDT5 describes tundra and alpine
areas.

In this study, landscape-level biodiversity objectives were
modelled as legally established seral stage distribution tar-
gets applicable by landscape unit (LU), BEC subzone, and
**Fig. C1.** Area distribution by key non-timber objectives. Note that each polygon can contribute to several objectives. VQO = visual quality objectives, LCOR = landscape corridors, COEC = core ecosystems, FSW = fisheries sensitive watersheds.

<figure><img src="image_10.png" /><figcaption>**Table C1.** Landscape-level biodiversity objectives</figcaption></figure>

<table><thead><tr><th>Landscape Unit</th><th>BEC subzone</th><th>NDT</th><th>Old Seral</th><th>Mature + Old Seral</th><th>Young Seral</th><th>PFLB (ha)</th></tr></thead><tbody><tr><td>Bulkley</td><td>ESSFmc</td><td>NDT2</td><td>9%>250 years</td><td>14%>120 years</td><td>None</td><td>204</td></tr><tr><td>Bulkley</td><td>ESSFwv</td><td>NDT1</td><td>28%>250 years</td><td>54%>120 years</td><td>17%<40 years</td><td>526</td></tr><tr><td>Bulkley</td><td>ICH mc 1</td><td>NDT2</td><td>13%>250 years</td><td>46%>100 years</td><td>27%<40 years</td><td>1,088</td></tr><tr><td>Bulkley</td><td>SBS dk</td><td>NDT3</td><td>10%>140 years</td><td>None</td><td>None</td><td>1,491</td></tr><tr><td>Bulkley</td><td>SBS mc 2</td><td>NDT3</td><td>10%>140 years</td><td>None</td><td>None</td><td>94</td></tr><tr><td>Copper</td><td>ESSFmc</td><td>NDT2</td><td>9%>250 years</td><td>28%>120 years</td><td>36%<40 years</td><td>2,379</td></tr><tr><td>Copper</td><td>ESSFwv</td><td>NDT1</td><td>19%>250 years</td><td>36%>120 years</td><td>22%<40 years</td><td>15,301</td></tr><tr><td>Copper</td><td>ICH mc 1</td><td>NDT2</td><td>13%>250 years</td><td>46%>100 years</td><td>27%<40 years</td><td>1,989</td></tr><tr><td>Copper</td><td>MH mm 2</td><td>NDT1</td><td>19%>250 years</td><td>36%>120 years</td><td>22%<40 years</td><td>149</td></tr><tr><td>Copper</td><td>SBS mc 2</td><td>NDT3</td><td>11%>140 years</td><td>23%>100 years</td><td>54%<40 years</td><td>1,681</td></tr><tr><td>Kitseguecla</td><td>ESSFwv</td><td>NDT1</td><td>19%>250 years</td><td>36%>120 years</td><td>22%<40 years</td><td>523</td></tr><tr><td>Kitseguecla</td><td>ICH mc 1</td><td>NDT2</td><td>13%>250 years</td><td>46%>100 years</td><td>27%<40 years</td><td>5</td></tr><tr><td>Telkwa</td><td>ESSFmc</td><td>NDT2</td><td>9%>250 years</td><td>28%>120 years</td><td>36%<40 years</td><td>2,079</td></tr><tr><td>Telkwa</td><td>ESSFwv</td><td>NDT1</td><td>19%>250 years</td><td>36%>120 years</td><td>22%<40 years</td><td>172</td></tr><tr><td>Telkwa</td><td>SBS dk</td><td>NDT3</td><td>11%>140 years</td><td>23%>100 years</td><td>54%<40 years</td><td>1,287</td></tr><tr><td>Telkwa</td><td>SBS mc 2</td><td>NDT3</td><td>11%>140 years</td><td>23%>100 years</td><td>54%<40 years</td><td>2,752</td></tr><tr><td>Trout Creek</td><td>ESSFwv</td><td>NDT1</td><td>19%>250 years</td><td>36%>120 years</td><td>22%<40 years</td><td>247</td></tr><tr><td>Trout Creek</td><td>ICH mc 1</td><td>NDT2</td><td>9%>250 years</td><td>31%>100 years</td><td>36%<40 years</td><td>80</td></tr><tr><td>Trout Creek</td><td>SBS dk</td><td>NDT3</td><td>11%>140 years</td><td>23%>100 years</td><td>54%<40 years</td><td>89</td></tr></tbody></table>

NDT (Table C1). For the study area, there were 19 reporting units, each with three sets of targets corresponding to the seral requirement. These targets are defined as the percentage of area that should be covered by “old,” “mature” and “young” stands. Here, “old” is defined as > 250 years in all subzones except SBSdk/mc2, and as > 140 years in the SB-Sdk/mc2; “Mature” is defined as > 120 years in the MHmm2 and ESSFmc/mk/wv, as > 100 years in the ICHmc1/mc2 and SBSdk/mc2, and as > 80 years in the CWHws2; “Young” is defined as <= 40 years in all subzones.

## Fisheries sensitive watersheds

Fisheries sensitive watersheds are meant to recognize the integral between upland conditions and their influence on maintaining aquatic conditions necessary to sustain healthy

fish populations. A fishery sensitive watershed is legally defined under BC forest legislation as a watershed with significant fisheries values and sensitivity. The disturbance level within fisheries sensitive watersheds was controlled via equivalent clearcut area index, an indicator inverse to the hydrological recovery of stands used to measure the clearcut area impact on watershed hydrological features (Table C2).

One fisheries sensitive watersheds (Cumming Creek) overlapped with the study area that required a maximum equivalent clearcut area of 30%. The curves of equivalent clearcut area were developed for each future managed analysis unit during the yield development based on height/hydrological recovery relationship (Table C2).
<table><caption>**Table C2.** Relationship between stand height, hydrological recovery, and equivalent clearcut area.</caption><thead><tr><th>Height (m)</th><th>Hydrological Recovery (%)</th><th>Equivalent clearcut area (%)</th></tr></thead><tbody><tr><td>&lt;=1</td><td>0</td><td>100</td></tr><tr><td>&gt;1 and &lt;=2</td><td>1</td><td>99</td></tr><tr><td>&gt;2 and &lt;=3</td><td>8</td><td>92</td></tr><tr><td>&gt;3 and &lt;=4</td><td>18</td><td>82</td></tr><tr><td>&gt;4 and &lt;=5</td><td>31</td><td>69</td></tr><tr><td>&gt;5 and &lt;=6</td><td>46</td><td>54</td></tr><tr><td>&gt;6 and &lt;=7</td><td>61</td><td>39</td></tr><tr><td>&gt;7 and &lt;=8</td><td>76</td><td>24</td></tr><tr><td>&gt;8 and &lt;=9</td><td>91</td><td>9</td></tr><tr><td>&gt;9</td><td>100</td><td>0</td></tr></tbody></table>

## Visual quality objectives

The purpose of visual quality objectives (VQO) is to ensure that the scenic qualities of a forested hillside continue to meet the expectations of the public and the tourism industry while providing opportunities for harvesting to the forest industry. These scenic areas are typically steep forested hillsides which are important to the tourism industry and public social values; they are spatially defined in the visual landscape inventory dataset. Management of the area does not exclude timber harvesting but requires harvesting practices to be carried out in a manner whereby the designated objective for the area continues to be met. There are five designated objectives, each with different level of alteration allowed: (1) preservation where no alteration is allowed, (2) retention, where up to 1.5% alteration is allowed, (3) partial retention up to 7%, (4) modification up to 18%, and (5) maximum modification up to 35%. The alteration is defined as the amount of area allowed under a green-up height within each spatially-defined VQO polygon.

Areas with visual quality objectives (VQO) are treated differently depending on the objectives. Where the objective is set to preservation (P), the polygons are excluded from the THLB and set as visual reserves. For other VQO objectives, constraints are set on the maximum percentage of area that can be visually altered in each VQO polygon (i.e., % of area cleared where regeneration has not achieved the visually effective green-up state defined by minimum targets in terms of stocking, heights and coverage).

In this study, only VQO polygons with objectives of retention, partial retention, and modification overlapped with the study area. For these polygons, the alteration allowed was thus capped at 1.5%, 7%, and 18%, respectively. Because slope data was not available, topography could not be considered in the definition of green-up heights and VQO targets, calling for a simplified approach that was likely less constraining than the analysis done in the timber supply review for the area. The visually effective green-up height was set to 6 m for all polygon, and the alteration was thus tracked in each VQO polygon as the forested area under this threshold height of 6 m. The average age where the 6 m green-up height was reached was determined from all future stands within each VQO polygon

(as the area-weighted average). The average ages ranged between 10 and 30 years (i.e., it took 10 to 30 years, depending on site productivity, to reach 6 m in height). For the study area, there were 39 reporting units, each with a maximum alteration corresponding to the designated objective (1.5% to 18%). Note that part of the areas with VQO were removed from the area eligible for harvest (i.e., were part of the NTHLB in Fig C1) because they overlapped with other non-timber objectives that required to set the area aside as reserves.

## Ecosystem network: core ecosystems and landscape corridors

The design of forest ecosystem networks is a long-term landscape-scale management approach applied in BC to meet conservation objectives by preserving sensitive areas and improving landscape connectivity. Ecosystem networks are composed of core ecosystems and landscape corridors: core ecosystems are areas identified spatially as important for the preservation of rare forest structures (e.g., old seral or mature forests), rare ecosystems (including non-forested areas), wildlife habitat, water quality (e.g., headwaters or lakes) or cultural heritage; while landscape corridors ensure habitat connectivity within the landscape.

In this study, the landscape corridors and core ecosystems were spatially-explicit management zones defined legally by the Bulkley Valley Sustainable Resource Management Plan (Government of British Columbia 2005). These areas were designed to provide habitat connectivity within the landscape, reduce habitat fragmentation, and protect regionally significant, sensitive, and vulnerable wildlife and plant species communities. In the case of landscape corridors (LCOR in Fig. C1), a minimum of 70% of the forested area had to be older than 80 years for each unique combination of the landscape corridor ID and landscape unit. In the case of the core ecosystems (COEC in Fig. C1), a maximum of 5% of the forested area could be under 50 years for each unique combination of core ecosystem ID and landscape unit. For the study area, there were 12 reporting units for core ecosystems and 34 for landscape corridors. Operationally, there will be no harvesting in areas identified as rare ecosystems, cultural heritage resources, and other important forest structures.

## Appendix D. Growth and yield predictions

### Pre-existing natural stands—Variable density yield projection (VDYP)

VDYP (Variable Density Yield Projection; Government of British Columbia 1993) is an empirical stand level model that projects the average yield of naturally regenerated and unmanaged stands, using inventory data as input (incl. species composition, age, stand height, basal area, trees/ha and BEC zone).

In this study, VDYP (batch v7.30a, build: 299) was used to extract the yield curves for each vegetation resource inventory (VRI) polygon that overlapped with the study area. Then,
<figure><img src="image_11.png" /></figure>

for each AU corresponding to an existing natural stand (AU series < 100 in Table A3), an area-weighted yield curve (with live and dead components) was consolidated from all VDYP yields overlapping with the AU in question.

For stands impacted by mountain pine beetle (and identified in the VRI, with information on the year of attack and percentage of mortality), an area-weighted dead yield value was calculated for each AU and maintained for 15 years following the MPB attack year (also determined as an area-weighted average year for each AU). Past this 15 years “shelf life” post MPB attack, the dead volume component was removed and replaced by a regeneration layer similar to the original yield, delayed by 10 years, and adjusted with the estimated MPB impact (i.e., dead percentage). Finally, live and dead yield components were included in the forest estate model and tracked separately for management and reporting purposes.

## Planted stands—Table interpolation program for stand yields (TIPSY)

TIPSY (Table Interpolation Program for Stand Yield; Government of British Columbia 2018a), the most widely used model in timber supply analyses in BC, is a metamodel of TASS (see below). TIPSY predicts growth and yield using a database of stand yield tables generated by TASS for single-species, even-aged stands. For mixed species stands, yield tables are computed as the weighted average of pure-species yield tables.

In this study, the TIPSY model (batch v4.3.2) was used to derive yields for managed stands (i.e., existing and future stands artificially regenerated stands after harvest events, AU 100, 200, and 300 series in Table A3), using regeneration assumptions detailed in Table D1 (reproduced from latest timber supply review for Bulkley timber supply area (FLNRORD 2012)). Endemic losses were accounted via operational adjustment factors (OAF). There were two OAFs: OAF 1 affected the magnitude of the yield curve, was constant across all ages and was associated with uneven stocking or gaps, whereas the impact of OAF 2 accelerated with age and represented the impact of decay, waste and breakage. The values for OAF 1 were different while OAF 2 was set to 0.95 for all AUs. Yield curves for the *Dothiostroma* stands (AU 300 series in Table A3) were identical to the yield curves for the existing managed stands (AU 200 series), except that they were delayed by 10 years to account for regeneration delays caused by Dothistroma needle blight (i.e., yield at age 50 in AU 200 series was the yield at age 60 in 300 series—*Dothiostroma* stands).

## Thinned stands—tree and stand simulator (TASS III)

TASS (Tree And Stand Simulator; Government of British Columbia 2018b) is a biologically based, spatially explicit, individual tree, and crown-driven model used by the BC government to simulate managed forests in BC. The public version of TASS III (v. 4.0.93 beta) was used to develop post-thinning curve modifiers. For each AU, two sets of yield curves were modelled: one unaltered set that estimated yield curves in normal conditions (i.e., unthinned) and another set that ap-

plied a commercial thinning removing 40% of merchantable volume 20 years prior to the minimum harvest age (MHAs, see Table A3).

The ratio between the two curves resulted in nonlinear factor curves that were applied to each AU for commercially thinned stands in the forest estate model (Fig. D1). By multiplying their reference yield curve (derived either from VDYP or TIPSY) with the corresponding growth modifier curve, commercially thinned stands had their residual yield curves adjusted to reflect post-thinning responses modelled by TASS III.

Note that the thinning response was not assessed for every year of the 10-year opportunity window (20 to 10 years prior to MHA), but only once for a thinning date set at 20 years prior to MHA. We expected this simplification to have little impact since the 10-year opportunity window is quite short, and the 10-year time step of the forest estate model limits possible misalignment with post-thinned yield curves.

The limitation of TASS III was that it only included two applicable species to the study area: white spruce and lodgepole pine. Consequently, the actual species composition of each AU (Table D2) was mapped to the two available species in TASS III based on similarities between species (Table D2), following the approach by Griess et al. (2019).

## Yield curves under climate change and different planting regimes (Landis II)

The Forest Carbon Succession extension to LANDIS-II (Scheller et al. 2007; Dymond et al. 2016) is a widely used forest landscape modelling framework capable of simulating ecological processes over time at landscape level, including climate change effects. Previous studies applied LANDIS-II in the Copper-Pine Creek valley (Dymond et al. 2020, 2016; Hof et al. 2017), and their outcomes were used as a basis for modelling the impact of climate change and alternative species composition in this study.

Following the same approach as for simulating post-thinning response, curve modifiers were derived from LANDIS-II to develop yield curves under climate change for baseline composition (i.e., even mixture of all species listed in stocking standards) and alternative planting regimes (i.e., Douglas-fir enrichment and species diversification), based on changes in biomass projected by LANDIS-II in the case study area (Hof et al. 2017).

Curve modifiers for climate change and alternative planting regimes were defined as the age-dependent ratio between the yield under baseline climatic and planting scenario (base scenario in Fig. D2) and the yield under each of the alternative scenarios predicted by LANDIS-II (climate change, Douglas-fir enrichment and diversification). These ratios were then multiplied to the reference yield curves used in this study (derived from VDYP and TIPSY, see above). The figure below shows the resulting modified yield curves under baseline climate, climate change and different planting regimes for five example AUs located in different ecological conditions (BEC variants; Fig. D2).

Climate change impacts on yield curves were however modelled differently for existing stands (i.e., already present
<table><caption>**Table D1. TIPSY regeneration assumptions.**</caption><thead><tr><th>AU</th><th>BEC</th><th>RM</th><th>RP</th><th>Density</th><th>RD</th><th>OAF1</th><th>Species Composition*</th><th>SI Spp1</th><th>GW Spp1</th><th>SI Spp2</th><th>GW Spp2</th><th>SI Spp3</th><th>GW Spp3</th><th>SI Spp4</th><th>GW Spp4</th></tr></thead><tbody><tr><td>106</td><td>ESSF</td><td>P</td><td>0.1</td><td>1,000</td><td>1</td><td>0.85</td><td>Pl80 Bl20</td><td>13.5</td><td>7</td><td>13.3</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>106</td><td>ESSF</td><td>N</td><td>0.9</td><td>8,900</td><td>1</td><td>0.85</td><td>Pl80 Bl20</td><td>13.5</td><td></td><td>13.3</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>107</td><td>ESSF</td><td>P</td><td>1</td><td>1,200</td><td>1</td><td>0.85</td><td>Sw40 Bl40 Pl20</td><td>15.3</td><td>8</td><td>13.7</td><td></td><td>16</td><td>7</td><td></td><td></td></tr><tr><td>108</td><td>ESSF</td><td>P</td><td>1</td><td>1,200</td><td>1</td><td>0.85</td><td>Sw45 Bl35 Pl20</td><td>11.7</td><td>8</td><td>14.9</td><td></td><td>18.3</td><td>7</td><td></td><td></td></tr><tr><td>109</td><td>ESSF</td><td>P</td><td>1</td><td>1,000</td><td>1</td><td>0.95</td><td>Sw55 Bl45</td><td>12.6</td><td>8</td><td>10.8</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>114</td><td>ESSF</td><td>P</td><td>0.5</td><td>1,200</td><td>1</td><td>0.85</td><td>Pl80 Bl20</td><td>10</td><td>7</td><td>10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>114</td><td>ESSF</td><td>N</td><td>0.5</td><td>2,500</td><td>1</td><td>0.85</td><td>Pl80 Bl20</td><td>10</td><td></td><td>10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>115</td><td>ESSF</td><td>P</td><td>0.5</td><td>1,200</td><td>1</td><td>0.95</td><td>Sw50 Bl30 Pl20</td><td>14.8</td><td>8</td><td>11.4</td><td></td><td></td><td>7</td><td></td><td></td></tr><tr><td>115</td><td>ESSF</td><td>N</td><td>0.5</td><td>2,200</td><td>1</td><td>0.95</td><td>Sw50 Bl30 Pl20</td><td>14.8</td><td></td><td>11.4</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>116</td><td>ESSF</td><td>P</td><td>1</td><td>1,200</td><td>1</td><td>0.85</td><td>Sw50 Bl30 Pl20</td><td>12.4</td><td>8</td><td>12</td><td></td><td></td><td>7</td><td></td><td></td></tr><tr><td>117</td><td>ESSF</td><td>P</td><td>0.3</td><td>1,000</td><td>1</td><td>0.85</td><td>Bl60 Sw40</td><td>10</td><td></td><td>10</td><td>8</td><td></td><td></td><td></td><td></td></tr><tr><td>117</td><td>ESSF</td><td>N</td><td>0.7</td><td>3,500</td><td>1</td><td>0.85</td><td>Bl60 Sw40</td><td>10</td><td></td><td>10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>118</td><td>ICH</td><td>P</td><td>0.2</td><td>1,000</td><td>1</td><td>0.85</td><td>Pl80 Hw20</td><td>19.6</td><td>7</td><td>17.1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>118</td><td>ICH</td><td>N</td><td>0.8</td><td>5,300</td><td>1</td><td>0.85</td><td>Pl80 Hw20</td><td>19.6</td><td></td><td>17.1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>119</td><td>ICH</td><td>P</td><td>0.5</td><td>1,200</td><td>1</td><td>0.85</td><td>Sw50 Hw20 Bl15 Pl15</td><td>23.3</td><td>8</td><td>18.5</td><td></td><td>18.18</td><td></td><td>22.4</td><td>7</td></tr><tr><td>119</td><td>ICH</td><td>N</td><td>0.5</td><td>2,200</td><td>1</td><td>0.85</td><td>Sw50 Hw20 Bl15 Pl15</td><td>23.3</td><td></td><td>18.5</td><td></td><td>18.8</td><td></td><td>22.4</td><td></td></tr><tr><td>120</td><td>ICH</td><td>P</td><td>1</td><td>1,000</td><td>1</td><td>0.85</td><td>Sw55 Bl25 Hw20</td><td>21</td><td>8</td><td>15</td><td></td><td>15</td><td></td><td></td><td></td></tr><tr><td>125</td><td>MH</td><td>P</td><td>0.15</td><td>900</td><td>1</td><td>0.85</td><td>Ba70 Hm15 Bl15</td><td>12.6</td><td></td><td>12.6</td><td></td><td>12.6</td><td></td><td></td><td></td></tr><tr><td>125</td><td>MH</td><td>N</td><td>0.85</td><td>6,800</td><td>1</td><td>0.85</td><td>Ba70 Hm15 Bl15</td><td>12.6</td><td></td><td>12.6</td><td></td><td>12.6</td><td></td><td></td><td></td></tr><tr><td>128</td><td>SBS</td><td>P</td><td>0.3</td><td>1,200</td><td>1</td><td>0.85</td><td>Pl80 Sw20</td><td>16.1</td><td>7</td><td>14.9</td><td>8</td><td></td><td></td><td></td><td></td></tr><tr><td>128</td><td>SBS</td><td>N</td><td>0.7</td><td>3,600</td><td>1</td><td>0.85</td><td>Pl80 Sw20</td><td>16.1</td><td></td><td>14.9</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>129</td><td>SBS</td><td>P</td><td>0.5</td><td>1,200</td><td>1</td><td>0.95</td><td>Pl80 Sw20</td><td>20.3</td><td>7</td><td>18.8</td><td>8</td><td></td><td></td><td></td><td></td></tr><tr><td>129</td><td>SBS</td><td>N</td><td>0.5</td><td>2,400</td><td>1</td><td>0.95</td><td>Pl80 Sw20</td><td>20.3</td><td></td><td>18.8</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>130</td><td>SBS</td><td>P</td><td>1</td><td>1,000</td><td>1</td><td>0.85</td><td>Sw80 Pl20</td><td>20</td><td>8</td><td>21.7</td><td>7</td><td></td><td></td><td></td><td></td></tr><tr><td>133</td><td>SBS</td><td>P</td><td>0.5</td><td>1,200</td><td>1</td><td>0.95</td><td>Pl50 Sw35 Bl15</td><td>19.3</td><td>7</td><td>19.3</td><td>8</td><td>15.6</td><td></td><td></td><td></td></tr><tr><td>133</td><td>SBS</td><td>N</td><td>0.5</td><td>2,200</td><td>1</td><td>0.95</td><td>Pl50 Sw35 Bl15</td><td>19.3</td><td></td><td>19.3</td><td></td><td>15.6</td><td></td><td></td><td></td></tr><tr><td>134</td><td>SBS</td><td>P</td><td>1</td><td>1,000</td><td>1</td><td>0.95</td><td>Sw50 Bl30 Pl20</td><td>18.6</td><td>8</td><td>19</td><td></td><td>18.8</td><td>7</td><td></td><td></td></tr><tr><td>135</td><td>SBS</td><td>P</td><td>1</td><td>400</td><td>1</td><td>0.85</td><td>Sw100</td><td>12</td><td>8</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

<footnote>*Species composition is represented as species code (see meaning of letters thereafter) followed by percentage composition (i.e., Pl40 means 40% of Lodgepole pine). At: trembling aspen, Ba: amabilis fir, Bl: subalpine fir, Cw: western red cedar, Ep: paper birch, Hm: mountain hemlock, Hw: western hemlock, Pl: lodgepole pine, Sw: white spruce, Sx: spruce hybrid, Yc: yellow cedar. AU = Analysis Unit; BEC zones = Biogeographical zones (see Fig. 1-c); RM = Regeneration Method (N = Natural, P = Planted); RP = Regeneration Percentage (% of natural or planted); Density = planting density (t/ha); RD = Regeneration Delay (years); OAF1 = Operational Adjustment Factor 1; SI Sp1 = Site Index for species 1 (ordered by decreasing species proportion), *idem* for SI Sp2, SI Sp3 and SI Sp4; GW Sp1 = Genetic Worth (or gain) for species 1 (ordered by decreasing species proportion), *idem* for GW Sp2, GW Sp3 and GW Sp4.</footnote>
<figure><img src="image_12.png" /><figcaption>**Table D2.** Species composition conversion in TASS III.</figcaption></figure>

<table><thead><tr><td>Spp</td><td>Prorated Spp in TASS 3</td><td>Species name</td><td>Genus species authority</td></tr></thead><tbody><tr><td>At</td><td>Pl</td><td>trembling aspen</td><td>Populus tremuloides Michx.</td></tr><tr><td>Ba</td><td>Sx</td><td>amabilis fir</td><td>Abies amabilis (Douglas ex Loudon) Douglas ex Forbes</td></tr><tr><td>Bl</td><td>Sx</td><td>subalpine fir</td><td>Abies lasiocarpa (Hook.) Nutt.</td></tr><tr><td>Cw</td><td>Sx</td><td>western red cedar</td><td>Thuja plicata Donn ex D. Don</td></tr><tr><td>Ep</td><td>Pl</td><td>paper birch</td><td>Betula papyrifera Marshall var. papyrifera</td></tr><tr><td>Hm</td><td>Pl</td><td>mountain hemlock</td><td>Tsuga mertinsiana (Bong.) Carrière</td></tr><tr><td>Hw</td><td>Pl</td><td>western hemlock</td><td>Tsuga heterophylla (Raf.) Sarg.</td></tr><tr><td>Pl</td><td>Pl</td><td>lodgepole pine</td><td>Pinus contorta Douglas ex Loudon</td></tr><tr><td>Sw</td><td>Sx</td><td>white spruce</td><td>Picea glauca (Moench) Voss</td></tr><tr><td>Sx</td><td>Sx</td><td>spruce hybrid</td><td>Picea glauca x Picea engelmannii Parry ex Engelm.</td></tr><tr><td>Yc</td><td>Sx</td><td>yellow cedar</td><td>Callitropsis nootkatensis (D. Don) Oerst. Ex D.P. Little</td></tr></tbody></table>

at the beginning of the planning horizon), the first generation of future stands of (i.e., stands regenerated after harvest of existing stands), and further generations of future stands (i.e., stands regenerated after harvest of first generation of future stands).

For existing stands, climate change impacts applied gradually, with no impact at the beginning of the simulation and full impact by year 100 of the simulation, which is consistent with the climate change projections (HADGEM A1B predicts + 7 °C in temperature and -200 mm in precipitation over 100 years) and impacts predicted by Landis II. Growth modifier curves derived from Landis II projections were therefore multiplied by a gradual curve (linear interpolation between 0 at age 0 and 1 at age 100) to simulate gradual change in climatic conditions. This curve was further shifted to match the current age of each stand (within each AU) at the beginning of simulation. Thus, each stand within each AU had its original natural yield changed, with the change starting with the current age relative to year 2017 (i.e., year 2017 was year zero of the planning horizon), and then incrementally changing as stand aged for the next 100 years. After the first 100 years of simulations, climate change impacts remained constant (equal to full impact simulated at year 100 of the planning horizon).

For the first generation of future stands, climate change impacts were also applied gradually (but without the age shift), with no impact at age 0 and full impact by age 100. Note that, due to technical limitations, we did not adjust the timing of the gradual application of climate change impacts to the actual date of stand initiation to reflect current climatic conditions (i.e., stand regenerated at year 50 of the planning horizon experienced no impact at age 0). This implied lags of varying importance in the simulation of climate change impacts for some of the first generation stands. Note that most existing stands are transitioned to future stands in the first 50 years of the planning horizon.

For the further generations of future stands, climate change impacts are simulated in full from age 0 (i.e., full values of growth modifiers are used to modulate reference yield curves), as these stands are regenerated after year 100 in the planning horizon. This was done by transitioning these stands to new AUs, with climate change yield curves.
<figure><img src="image_13.png" /><figcaption>**Fig. D2.** Example growth and yield curves with baseline, climate change, Douglas-fir enrichment and diversification modifiers for AUs located in five BEC zones in the case study area.</figcaption></figure>

**Additional references**

British Columbia Forest Service and British Columbia Environment. 1995. Biodiversity Guidebook. Queens Printer, British Columbia. Available from https://www.for.gov.bc.ca/hfd/library/documents/bib19715.pdf [accessed 09 July 2023].

Fenger, M. 1996. Implementing biodiversity conservation through the British Columbia Forest Practices Code. Forest Ecol. Manag. 85: 67-77. doi:10.1016/S0378-1127(96)03751-6.

FLNRORD. 2012. Bulkley timber supply area timber supply review: updated data package following completion of the timber supply analysis. Ministry of Forests, Lands, and Natural Resource Operations. Available from https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/stewardship/forest-analysis-inventory/tsr-annual-allowable-cut/bulkley_tsa_data_package.pdf [accessed 09 December 2023].

Government of British Columbia. n.d. Stocking Standards [online]. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/silviculture/stocking-standards [accessed 09 December 2023].

Government of British Columbia. 2018a. TIPSY, version 4.4 [online]. Forestry Analysis and Inventory Branch. BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/field-forms-and-software/software-download [accessed 09 July 2023].

Government of British Columbia. 2018b. Tree & Stand Simulator (TASS) [online]. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/tree-and-stand-simulator-tass [accessed 09 July 2023].

Government of British Columbia. 2005. Bulkley Valley Sustainable Resource Management Plan. Available from https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/natural-resource-use/land-water-use/crown-land/land-use-plans-and-objectives/skeena-region/bulkleyvalley-srmp/bulkley_srmp.pdf [accessed 09 December 2023].

Government of British Columbia. 1993. Variable Density Yield Projection (VDYP) [online]. Available from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory/growth-and-yield-modelling/variable-density-yield-projection-vdy [accessed 09 July 2023].

Ministry of Forests. 1999a. Guidelines for developing stand density management regimes. Ministry of Forests, British Columbia, Canada.

Ministry of Forests. 1999b. Guidelines for commercial thinning. Ministry of Forests, British Columbia, Canada.