ESTIMATING CARBON STOCKS AND FLUXES IN AN EXPERIMENTAL LOGGING TRIAL WITHIN BRITISH COLUMBIA’S INLAND TEMPERATE RAINFOREST by Maria Alejandra Chadid Hernandez B.Sc., University of the Andes, 2010 M.Sc. National University of Colombia, 2015 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN NATURAL RESOURCES AND ENVIRONMENTAL STUDIES UNIVERSITY OF NORTHERN BRITISH COLUMBIA May 2024 © Maria Chadid Hernandez, 2024 Abstract Temperate rainforests in British Columbia, Canada, are recognized for providing a suitable habitat for fauna and flora and storing a significant amount of carbon above and below ground. However, silviculture practices such as clearcuts have affected these forests, in some cases shifting their role from carbon sinks to carbon sources. To avoid some of the impacts of clearcutting, partial harvest practices have been implemented, which aim to leave some pre-harvest biological legacies behind to promote faster ecosystem recovery and provide a habitat for biodiversity. Carbon-related research in these systems has focused either on short-term studies or analyzing only a few pools, highlighting the need for long-term studies with a comprehensive analysis of carbon dynamics above and below ground. Here, I studied the long-term effect of partial harvest and clearcutting on forest carbon stocks by examining the carbon dynamics in multiple pools above (trees, downed woody material) and below ground (forest floor and roots) in the Northern Wetbelt, a longterm, replicated, research trial. I used historical databases and empirical data collected 19 years after harvest to estimate total forest carbon stocks and their change over time in mature and old-growth Interior Cedar-Hemlock forests. In addition, I identified the treatment that promotes the greatest carbon accumulation rates in the live trees given different canopy conditions. Within 19 years of harvest forest carbon stocks were much higher in partial harvest with high forest retention (449.49 ± 32.22 Mg C ha-1) and low retention (192.86 ± 14.63 Mg C ha-1) compared to clearcut conditions (136.29 ± 11.94 MgC ha-1). The carbon dynamics showed that pools can have opposite trends affecting differently the overall forest carbon stores. Specifically, while the carbon stocks in live trees increased, the coarse woody debris ii decreased. This interaction kept forest carbon stocks stable with no recovery toward the preharvest values within 19 years. Finally, the partial harvest treatment with higher forest retention had the canopy conditions that promoted the greatest carbon accumulation rates in live trees. Furthermore, the trees located at the edge of a managed opening showed a more rapid annual growth (0.0023 ± 0.0001 Mg C year-1) than the ones at the interior of the sampling plots (0.0017 ± 0.0001 Mg C year-1). This thesis contributes to forest management by providing long-term carbon estimates from different harvesting practices in mature and old-growth inland temperate rainforests to support decision-making toward more sustainable use of natural resources. iii Table of contents Abstract .................................................................................................................................... ii Table of contents .................................................................................................................... iv List of tables............................................................................................................................ vi List of figures ......................................................................................................................... vii Glossary ................................................................................................................................ viii Acknowledgements ................................................................................................................ ix 1. Chapter One. Introduction ............................................................................................. 1 1.1. 2. Chapter Two. Methods.................................................................................................... 5 2.1. 2.1.1. 2.2. 4. Study site ............................................................................................................................... 5 Experimental design ......................................................................................................... 6 Data collection ...................................................................................................................... 8 2.2.1. Live and dead standing trees............................................................................................. 8 2.2.2. Planted live trees in harvested areas ................................................................................. 9 2.2.3. Downed woody material ................................................................................................. 10 2.2.4. Forest floor sampling ...................................................................................................... 11 2.3. 3. Research questions ................................................................................................................ 2 Carbon assessments............................................................................................................. 12 2.3.1. Live and dead standing trees........................................................................................... 12 2.3.2. Live and dead roots......................................................................................................... 14 2.3.3. Downed woody material ................................................................................................. 15 2.3.4. Forest floor ..................................................................................................................... 16 2.3.5. Forest carbon accumulation rate ..................................................................................... 16 2.3.6. Carbon stock error propagation estimation ..................................................................... 17 2.3.7. Linear mixed-effect model to estimate differences in Carbon stock 2020 ..................... 17 2.3.8. Edge effect on carbon accumulation in live trees ........................................................... 17 Chapter Three. Results ................................................................................................. 19 3.1.1. Carbon stocks 19 years after harvesting.............................................................................. 19 3.1.2. Carbon dynamics through time ........................................................................................... 20 3.1.3. Canopy condition effects on individual live tree carbon accumulation rate ...................... 24 Chapter Four. Discussion.............................................................................................. 26 4.1. Carbon stocks 19 years after harvesting.............................................................................. 26 4.2. Carbon dynamics through time ........................................................................................... 27 iv 4.3. 5. Canopy condition effect on individual live tree carbon accumulation rate ......................... 30 Chapter Five ................................................................................................................... 31 5.1. Conclusions ......................................................................................................................... 31 5.2. Project limitations and future work ..................................................................................... 33 Literature cited...................................................................................................................... 34 Appendix A ...................................................................................................................................... 42 Equations to model heights for trees ............................................................................................ 42 Allometric equations to estimate biomass for trees ...................................................................... 43 CWD Carbon stock adjustments by the sampling effort .............................................................. 44 Appendix B: ..................................................................................................................................... 47 Areas per site and Treatments ...................................................................................................... 47 Appendix C....................................................................................................................................... 48 Carbon stocks 19 years after harvest ............................................................................................ 48 Carbon stocks in live trees 19 years after harvest......................................................................... 50 Carbon stocks pre- and post-harvesting........................................................................................ 51 Appendix D ...................................................................................................................................... 53 Modelling Carbon stocks in harvested areas ................................................................................ 53 Appendix E ....................................................................................................................................... 55 Canopy conditions’ effect on live tree mean carbon accumulation rate ....................................... 55 Mean carbon accumulation rate ~ Treatment ............................................................................... 55 Mean carbon accumulation rate ~ Canopy conditions within each treatment .............................. 56 Appendix F ....................................................................................................................................... 57 Linear mixed-effect model for carbon stocks in 2020 .................................................................. 57 v List of tables Table 1. Characteristics from each of the study sites. Source (Jull et al., 2001; Matsuzaki et al., 2013) ................................................................................................................................... 6 Table 2. The proportion of sampled transects per plot per treatment and site. The missing plots imply they were not found and surveyed in 2020. .................................................................. 44 Table 3. Areas (ha) for all the treatments at each study site. The proportion of harvested/forested areas was estimated based on these areas. ................................................ 47 Table 4.Mean and standard error (±95% confidence interval) for LT, LR, DR, DST, CWD, FWD, FF and total forest C stocks in four forest harvesting treatment units at three ICH study areas, Units are MgC ha-1. Data 19 years after harvest (2020). .............................................. 48 Table 5. Mean carbon (MgC ha-1) from live trees (LT) in 2020 disaggregated by the contribution from the forested and harvested areas at each Treatment/Site. .......................... 50 Table 6. Mean and standard error (±95% confidence interval) for five pools (LT, LR, DST, DR, CWD), and total forest C stocks in four forest harvesting treatment units at three ICH study areas. Units are MgC ha-1. Data available preharvest (1999), 1 year (2001) and 19 years after harvest (2020). ................................................................................................................ 51 Table 7. Projected carbon stocks (MgC ha-1) in planted trees in regeneration plots for 2020. Data derived from growth and yield curves using SORTIE 7.05 and regeneration survey from 2016. Carbon stocks have already been adjusted by the proportional increase and by the proportion of the harvested area. ............................................................................................ 53 Table 8. a) Mean, standard error Live tree carbon accumulation rate (MgC year-1) given three canopy conditions within each treatment; b) Live tree mean carbon accumulation rate per treatment. ................................................................................................................................ 55 vi List of figures Figure 1. Location of the study sites at the ICH Northern Wetbelt Trials ............................... 5 Figure 2. Experimental design in the study sites at the Northern Wetbelt trials in British Columbia: a) East Twin, b) Minnow, c) Lunate. All the PSPs are located within the forested areas of each treatment unit. All regeneration subplots are located within the harvested areas of each treatment unit................................................................................................................ 7 Figure 3. Left: Structure and dimensions of the PSP located in the forested area of each treatment. Right: Crew member collecting data from a live tagged tree. ................................. 9 Figure 4. Layout of the permanent sampling plot (PSP) from a partial harvest treatment. To its left, the gray square represents the associated planted regeneration plot in the harvested areas (gray background), with the four nested subplots (white circles). Figure adjusted from M. Jull et al., 2002. .............................................................................................................................. 10 Figure 5. Forest floor data collection. a) example of a subsample of forest floor (red square) with the mineral floor (blue square) which was ultimately excluded and not analyzed; b) FF’s sampling and subsampling design: yellow dots (every 2 m) correspond to each subsample collected with the core within each PSP; c) final labelled sample of forest floor containing five subsamples. d) Experimental design to collect forest floor samples in the different treatments. ................................................................................................................................................. 12 Figure 6. Contribution of different carbon pools to mean total ecosystem carbon (MgC ha-1) 19 years after harvest in different harvesting treatments in the Northern Wetbelt for A) across all three sites and B) Each study site. Treatments are clearcut (CC), Group retention (GR) Group selection (GS) and unharvest control (UN). Sites are East Twin, Lunate, Minnow. Colors on bars represent different carbon pools (MgC ha-1): LT (live trees), LR (live roots), DR (dead roots), DST (dead standing trees), CWD (coarse woody debris), FWD (fine woody debris), FF (forest floor). ........................................................................................................ 20 Figure 7. A) Total ecosystem carbon stocks over time (MgC ha-1) per treatment across all sites B). Carbon stocks change over time for CWD. C) Carbon stocks change over time for live trees. D) Forest floor carbon stocks 10 and 19 years after harvest. 10-year post-harvest data from Matsuzaki et al., 2013. All stocks are in MgC ha-1. ....................................................... 22 Figure 8. Mean and standard error forest carbon accumulation rate by treatment (Mg C ha-1 year-1) across all sites. CC (clearcut), GR, (group retention), GS (group selection), UN (unharvested). The rate was estimated based on the live trees (LT) dead standing trees (DST), live roots (LR), dead roots (DR) and coarse woody debris (CWD) from 2001 and 2020. ..... 23 Figure 9. A) Ratio of sampled trees by canopy condition in each treatment. The results were standardized to the area of the PSPs within each treatment across all sites. B) mean and standard error of live tree carbon accumulation rate (MgC tree-1 yr-1) grouped by treatment (GS, GR) and by canopy conditions: Edge (E) and Closed (C).............................................. 25 vii Glossary AIC Akaike Information Criterion BC British Columbia BEC Biogeoclimatic zone C Closed canopy CC Clearcut CDW Coarse woody debris DBH Diameter at breast height DR Dead roots DST Dead standing trees E Edge canopy FCrate Forest Carbon accumulation rate FF Forest floor FWD Fine woody debris GR Group retention GS Group selection H Height ICH Interior Cedar-Hemlock IWB Interior Wetbelt LR Live roots LT Live trees O Open canopy PSP Permanent sampling plot UN Unharvested viii Acknowledgements “Once you make a decision, the universe conspires to make it happen.” — Ralph Waldo Emerson. I would like to thank my family for their constant love and support. Even at the distance you provided me with everything I needed at all times. My friends (too many to name them) who always have stood by me giving me the strength to confront and overcome challenges. Special acknowledgement to Miguel Arias, who has been there since the beginning of this journey dealing with the rollercoaster of my emotions at every stage. Thanks for pushing me with love, patience and perseverance and I really hope I can call you Dr. sooner rather than later, it will pay off the sacrifice. Thanks to Thomas Beaurepaire for constantly helping me to work on a better version of myself. I would also like to give a special thanks to all my committee members Michelle Venter, Caren Dymond, Phil Burton and Oscar Venter; you have given me the support, sympathy, wisdom, experience and trust that I needed to succeed with this project, your contributions played a crucial role at every stage of the research. Thanks for helping me to build critical thinking and important knowledge in these complex and beautiful forests which I hope we do not lose. I want to thank Michael Jull and honor the memory of Bruce Rogers, who established the amazing trial of Northern Wetbelt; the data collected and your constant insights were highly valuable for this research. I want to thank for the patience and contributions to this project to Flora Krivak-Tetley with Carbon modelling, Peter Ott with all the statistical analysis, Doug Thompson with tree analysis and the CSL colleagues with their time during fieldwork and insights about the project. I want to acknowledge the support from UNBC, their facilities and their staff at every stage of the project. In addition to the Forest Carbon Initiative, the Ministry of Forests and UNBC for providing the funds. Finally, I want to acknowledge that the land on which I lived and worked is the traditional unceded territory of the Lheidli T'enneh. None of this would have happened if everyone and everything in the universe had not conspired for this. ix 1. Chapter One. Introduction Forests sequester significant amounts of carbon from the atmosphere and can contribute to natural climate solutions (Bartowitz et al., 2022; Hudiburg et al., 2019). During 2000-2007, forests around the globe absorbed 2.3 ± 0.5 PgC year–1 , and temperate forests contributed 0.8 ± 0.1 PgC year–1 to these C sinks (Pan et al., 2011). The Interior Wetbelt (IWB) is a major temperate rainforest ecosystem located on the western flanks of the Canadian Rocky Mountains (DellaSala et al., 2022). According to the biogeoclimatic ecosystem classification (BEC) system used in British Columbia (BC), the IWB is located in the wet and cool areas of the Interior Cedar-Hemlock (ICH) (Jull et al., 2001). The old growth stage is dominated by the prevalence of large trees of western redcedar (Thuja plicata) and western hemlock (Tsuga heterophylla), which store large amounts of carbon in above and below ground carbon pools (Matsuzaki et al., 2013; Mildrexler et al., 2020). In BC, the temperate rainforests are important ecosystems for multiple reasons: they contain ‘antique forests’ that are the outcome of patchy natural disturbances, having forests with uneven-aged stands (Coxson et al., 2020; DellaSala et al., 2011); they are biodiverse and sustain complex food web dynamics and important species such as the endangered deep-snow mountain caribou (Rangifer tarandus caribou) (COSEWIC, 2014; DellaSala et al., 2011); they are known for having large carbon stocks, with estimates from field plots ranging from 583 MgC ha-1 up to 1275 MgC ha-1, being similar to some of the world’s most carbon-dense temperate forests (DellaSala et al., 2022; Matsuzaki et al., 2013). However, according to the Red-listed Ecosystem Criteria, forests of the IWB have a critical ecosystem status (DellaSala et al., 2022). Historical and cumulative logging since 1970 in the IWB has negatively affected these forests (DellaSala et al., 2022). Clearcutting is the most common practice of harvesting (Beese et al., 2019), which removes all or nearly all of the forest cover and creates homogeneous stands with planted regeneration (British Columbia. Ministry of Forests. Forest Practices Branch, 2003). An alternative to clearcuts are partial harvest systems, which aim to conserve attributes from the original forest to promote a faster recovery while retaining habitat and resources for different taxa (Beese et al., 2019; Gustafsson et al., 2012). Much of the research on partial harvesting has focused on forest values such as biodiversity (Fedrowitz et al., 2014; Stevenson et al., 2006) and forest health (Thorpe & Thomas, 2007). These studies have generally found that partial harvest can support higher levels of species richness and abundance of forest species than clearcuts (Fedrowitz et al., 2014; Simard et al., 2020). Moreover, research has found that these systems store more aboveground carbon than clearcuts converted to plantations 1 year after-harvest (Simard et al., 2020). More specifically in the IWB, Matsuzaki et al., (2013), found that when accounting only for the carbon aboveground pools one-year after harvesting, clearcutting reduced carbon stocks by approximately 85%, relative to an unharvested stand, whereas partial harvest with high forest retention (>70%) only reduced the stocks by 12%. In these silviculture systems, the openings could emulate gaps and create forest edges, defined as the boundary in canopy composition (Saeed et al., 2019). These edges are exposed to a variety of microclimatic conditions and their effect may influence the forest’s capacity to store carbon depending on variables such as tree species, edge structure and management (Saeed et al., 2019). Higher carbon stocks have been reported in the proximity of an edge due to the interaction of increased light availability, and a higher basal area from young trees (Meeussen et al., 2021; Reinmann & Hutyra, 2017; Simard et al., 2020). Although logging is an important disturbance in the temperate rainforests and results in immediate lower forest carbon stocks (Matsuzaki et al., 2013), it is still unclear how partial harvesting affects the long-term forest carbon dynamics across multiple pools above and belowground. It is also unclear how these dynamics can be affected by other variables such as the edge effect that result from the managed openings after harvesting. The lack of long-term empirical data on ecosystem carbon stocks, fluxes, and mechanisms in managed forests limits our ability to manage for carbon. 1.1. Research questions To address the gaps in the literature, I responded to the following research questions: 1. What are the carbon stocks in aboveground (i.e. living and dead trees, coarse and fine wood debris) and belowground (live and dead roots and forest floor), before, the year after and 19 years after harvest? 2 a. How do these vary across treatments, i.e. Group Selection (GS; high retention), Group Retention (GR; low retention), Clearcut (CC), and Unharvested forest (UN)? What is the carbon distribution among the surveyed pools and how does it change from pre-harvest to post-harvest? Based on the one-year post-harvest data from a temperate rainforest in the IWB, Matsuzaki et al., (2013) reported significant reductions in the total carbon stocks from treatments with little or no tree retention (GR and CC) compared to those with high retention (GS) and unharvested (UN) stands. For the same study site nineteen years after harvest, I expected to find similar patterns between the same treatments. Furthermore, I expected no significant differences between GS and UN carbon stocks. More specifically, the total forest ecosystem carbon 19 years after harvest would follow this order of diminishing magnitude UN ≥ GS > GR > CC, implying that the higher the forest retention, the higher the carbon stocks. 2. What is the forest carbon accumulation rate between the year after harvest to 19 years after harvest in partial harvest systems and clearcuts compared to an unharvested forest? In a meta-analysis on the impacts of partial harvest on forest structure, Zhou et al., (2013) found a negative impact on the stocks in the aboveground biomass, which decreased with harvest intensity. Although the losses induced by partial harvest would be recovered in time, it is challenging to determine how long it would take to recover given the lack of long-term available data. By using the data from 2001 and 2020, I expected to find carbon recovery and gains (positive carbon accumulation) in the treatments with high retention (GS) and UN treatment while expecting losses (negative carbon accumulation) in the low forest retention (GR) and CC; in addition, I expect the lower rates of carbon losses in GR compared to clearcut. 3. What forest management in partial harvest treatments create the canopy configuration or conditions (edge and closed) that promote the highest carbon accumulation rate? a. Are canopy conditions in partial harvested treatments (GS and GR) different? b. Is the live tree carbon accumulation rate different between partial treatments (not considering canopy condition as a factor)? c. Is the live tree carbon accumulation rate for each canopy condition different within the partial harvesting treatments? 3 Based on Reinmann & Hutyra, (2017) and Meeussen et al., (2021) that have reported faster growth in trees near the edge than in the forest interior, I expected that the treatment that creates more edge, and has more trees at these edges (due to harvesting) should have a higher carbon accumulation rate (GR > GS). To address my research questions and gaps from the literature, I remeasured a long-term silvicultural trial 19 years post-harvest and used historical data to estimate multiple components of the forest carbon dynamics. Overall, the results will help to understand the impact of different harvesting practices on the ecosystem carbon dynamics over a ~20-year period. The long-term accounting of these stocks under different silvicultural practices will contribute to reliable GHG emission estimates from different harvesting practices on the IWB, and it will support decision-making towards more sustainable forest management as a naturebased solution for climate change. 4 2. Chapter Two. Methods 2.1. Study site The Northern Wetbelt Silvicultural Systems Project is a set of experimental harvesting trials that compare partial harvesting systems to clearcut and unharvested stands in the Interior Wetbelt, which is composed of Interior Cedar-Hemlock (ICH) and wetter portions of the Engelmann Spruce-Subalpine Fir biogeoclimatic zones. The experimental sites were established between 2000-2001 as part of the BC Ministry of Forests Silvicultural Systems Program (Wiensczyk, 2012). I collected data in the three Cedar-Hemlock sites: Lunate Creek, Minnow Creek and East Twin Creek (Figure 1). Figure 1. Location of the study sites at the ICH Northern Wetbelt Trials Inland temperate rainforest forests combine maritime and continental weather conditions (DellaSala et al., 2011) (Table 1). All sites are dominated by western redcedar (Thuja plicata) and western hemlock (Tsuga heterophylla). Other species present, but with lower abundance, 5 are hybrid white spruce (Picea engelmannii x glauca) and subalpine fir (Abies lasiocarpa) (Jull et al., 2001; Matsuzaki et al., 2013). Although podzolic soils predominate across all sites, luvisols are also present. Humus forms were mostly humimor /hemimor and mormoder (Jull et al., 2001; Matsuzaki et al., 2013). Table 1. Characteristics from each of the study sites. Source (Jull et al., 2001; Matsuzaki et al., 2013) Site East Twin Creek Minnow Creek Lunate Creek BEC subzone and variant ICHwk3 ICHwk3 ICHwk2 Elevation (m) 900-1050 900-1050 950-1200 Coordinates 53° 30’N 120° 53° 28’N 120° 53° 50’ N 121° 20’O 18’O 28’O Aspect NW SW N Area (ha) 26 39.2 72.4 Mean annual 2.3-3.2 2.3-3.2 1.9-2.8 1169 939 1120 Estimated stand age 300-350 300-350 450-500 Harvesting year 2000 2001 2001 temperature (˚C) Mean annual precipitation (mm) 2.1.1. Experimental design The Northern Wetbelt trial consists of a randomized block design with three to four silvicultural treatments applied at each site: 1) “clearcut” (CC), 0% post-harvest retention; 2) “group retention” (GR), 30% post-harvest retention; 3) “group selection” (GS), 70% postharvest retention; 4) unharvested (UN), no harvest control (Figure 2). Due to constraints in the treatment layout, the East Twin site is a randomized incomplete block with only three treatments (i.e. UN, GS, CC), All harvesting was conducted in the winter of 2000 and 2001 and planting was done in 2001-2002 (Jull et al., 2001; Matsuzaki et al., 2013). All harvest was conducted as complete removal in the sub-set of the treatment block being harvested, i.e. there was no thinning-type treatment. 6 Figure 2. Experimental design in the study sites at the Northern Wetbelt trials in British Columbia: a) East Twin, b) Minnow, c) Lunate. All the PSPs are located within the forested areas of each treatment unit. All regeneration subplots are located within the harvested areas of each treatment unit. 7 2.2. Data collection To estimate the carbon stocks of above and belowground carbon pools 19 years post-harvest I collected data on live and dead standing trees, lying dead wood, and forest floor. All field data collection was undertaken from June to September 2020. Where possible, I adopted the methods presented in the experiment’s establishment report (Jull et al., 2001). In addition, I used the historical datasets from pre-harvest (1999) and 1-year post-harvest (2001) collected by Jull et al., (2001) to estimate forest carbon stocks and examine the carbon dynamics at a pool level. The partial harvesting treatments (i.e. GS and GR) resulted in increased forest openings, creating forest edges and a variety of canopy conditions, which can be measured in different ways. One of these is with the canopy closure, which is “the proportion of the sky hemisphere obscured by vegetation when viewed from a single point” (Jennings, 1999). I conducted a simple visual assessment of the canopy closure, which requires no specialized equipment, and allowed me to estimate the canopy closure by comparing the area of the canopy with a standard scale (Jennings, 1999). Although this method has low accuracy, it was enough to create three categories for the canopy conditions of each tree. Based on the literature (Andrieu et al., 2018; Jennings, 1999; Land Owner Resources Centre, 2000; Lorber et al., 2018) and the treatment design in the stands, I defined three canopy conditions within each sampling plot: 1) Closed: refers to the remnant natural area within the PSP, with a higher density of trees, whose adjusted canopy cover >50%; 2) Edge: refers to the transition zone between the forest interior and the managed openings, where most of the times it overlapped with the boundaries (“edge”) of the permanent sampling plots, whose adjusted canopy cover ranged between 40-50%; 3) Open (O): refers to the trees located in managed openings or natural gaps, whose adjusted canopy cover <30%. The canopy conditions for each tree were visually established while conducting fieldwork using as reference the previous definitions and the stand characteristics (Figure 3). 2.2.1. Live and dead standing trees Rectangular Permanent Sampling Plots (PSPs) measuring 20m x 50m or 0.1 ha were systematically established in each treatment among all three sites in 2000 (Figure 2). There 8 were eight to ten plots established per treatment unit (Jull et al., 2001). Each PSP was divided into four sectors for ease of sampling (Figure 3). For each live tree (LT) or dead standing tree (DST) tagged in 2001 I recorded data on: 1) the tree species; 2) diameter at breast height (DBH) using a diameter tape (cm); 3) the tree condition class, which was assessed based on British Columbia’s wildlife tree classification (B.C. Ministry of Forests and Range & B.C. Ministry of Environment, 2010); 4) canopy conditions: Edge, Open, and Closed; 5) new trees (e.g. ingress or natural regeneration) with DBH > 4 cm were tagged, and the same data were recorded, adding the height (H) (meters), which was measured with an electronic hypsometer. Figure 3. Left: Structure and dimensions of the PSP located in the forested area of each treatment. Right: Crew member collecting data from a live tagged tree. 2.2.2. Planted live trees in harvested areas The harvested areas within the GS, GR and CC treatments were planted with a 50-50 mix of Western redcedar and Hybrid white spruce in 2001 or 2002. Planted regeneration plots were established in these areas (Figure 2c). These consisted of four nested subplots (plot radius = 3.99 m, plot area = 0.005 ha) (Figure 4), for which all trees were tagged and flagged (Jull et al., 2002). The original team led monitoring surveys in 2007 and 2016, collecting tree data within the regeneration plots, measuring DBH, H, and tree condition. 9 Figure 4. Layout of the permanent sampling plot (PSP) from a partial harvest treatment. To its left, the gray square represents the associated planted regeneration plot in the harvested areas (gray background), with the four nested subplots (white circles). Figure adjusted from M. Jull et al., 2002. 2.2.3. Downed woody material Downed woody material was grouped in two categories based on Stevenson et al., 2001: 1) Coarse woody debris (CWD), defined as dead pieces, with diameter less than 7.5 cm, a tilt angle less than 45° from vertical, intercepting the transect (Harmon & Sexton, 1996); and 2) Fine woody debris (FWD) pieces ranging between 1 cm and 7.4 cm diameter. Each plot had a set of three independently directed 24 m line transects with random azimuths. The original transects established to collect data from CWD in 2001 were used to remeasure the CWD pieces in 2020 (Jull et al., 2001), and I recorded the FWD pieces for the first time. FWD pieces were recorded on the same transects but for only half of the distance of the original transects (i.e. 12 m). For each CWD and FWD piece, the diameter and length was measured using calipers and a 30m measuring tape; I recorded data on the canopy conditions: Open (O), Edge (E), Closed (C); new pieces were recorded following the same protocols as Stevenson et al., (2001). The decay class for the CWD was assessed based on the Field Manual for Describing Terrestrial Ecosystems (Stevenson et al., 2001); for the FWD, I defined fewer categories based on the decomposition state: Sound, Intermediate, Decayed. Finally, for the FWD, I also measured size in three categories: Small (1 cm < diameter > 3.4 cm), Medium (3.5 cm < diameter > 5.4 cm), and Large (5.5 cm < diameter > 7.4 cm) 10 2.2.4. Forest floor sampling The forest floor (FF) or the O horizon of the soil profile was defined as the dead organic surface that lies on top of the mineral soil. It includes either recognizable or unidentified organic material (e.g. small branches and twigs, fallen needles, small roots, dead grasses) (Pearson et al., 2005; Serniak, 2017). The FF samples were collected using a fabricated stainless steel corer of 5 cm diameter and 30 cm length, which was mounted on a battery-powered drill. To reduce the influence of spatial variability, each FF sample consisted of five subsamples, which were pooled every two meters starting from the plot centre up and downhill (Figure 5b). Each core subsample was placed and sorted on a recipient trough (“bed”) to separate it from the mineral soil, which was excluded from the sample and not analyzed, and subsequently I recorded the FF’s thickness (cm) with a metric ruler (Figure 5a). Each sample was stored in a plastic bag with a unique ID (Figure 5c) until returned from the field, then was stored at the Enhanced Forestry Lab at UNBC, and then shipped to Victoria, BC, to the Analytical Laboratory to estimate its dry weight (in grams, g), and Carbon (g) and Nitrogen content (g). A total of 249 FF samples were collected. The number of samples varied according to the treatment unit: in UN and CC, I collected two samples per PSP. In the partial harvesting treatments (GS, GR) I collected three samples per PSP in the different canopy conditions (O, E, C) (Figure 5d). The sample processing included several key steps: Each sample was dried at a constant temperature of 70°C and then weighed. Subsequently, the sample underwent flail grinding to reduce its consistency to a particle size of less than 2mm. Additionally, a subsample was taken and ring-ground for carbon and nitrogen analysis. 11 Figure 5. Forest floor data collection. a) example of a subsample of forest floor (red square) with the mineral floor (blue square) which was ultimately excluded and not analyzed; b) FF’s sampling and subsampling design: yellow dots (every 2 m) correspond to each subsample collected with the core within each PSP; c) final labelled sample of forest floor containing five subsamples. d) Experimental design to collect forest floor samples in the different treatments. 2.3. Carbon assessments 2.3.1. Live and dead standing trees Carbon stocks were defined as the amount of carbon stored in both live and dead biomass above and below ground at one point in time (Hoover & Smith, 2023), whose estimations were in megagrams of carbon per hectare (MgC ha-1) . For live trees (LT) I estimated the biomass or the dry weight, using species-specific allometric equations (Ung et al., 2008) based on the height (H) and (DBH) (Appendix A, Allometric equations). Since heights for trees were not collected in 2020, I modelled the height for each tree, using species-specific DBH:H relationships derived from quadratic equations, using as a reference value the historical databases collected at the Northern Wetbelt trial (Jull et al., 2001) (Appendix A, modelling equations). The obtained biomass (kg) was calculated at a tree level and converted to megagrams. All trees within each PSP were summed to get to a plot level value, which was then converted to per hectare units by multiplying by an expansion factor derived from the PSP-specific dimensions, 12 and finally converted to carbon per hectare (MgC ha-1) by multiplying a factor of 0.5 (Pearson et al., 2005). For each treatment, I divided the harvested/unharvested portions into two strata: 1) retained areas, which corresponded to the forest that remained in the treatment unit after harvesting, and 2) harvested areas, which corresponded to harvested areas of each treatment unit (Figure 4, Appendix B, Table 3). To get the total carbon stock for the treatment within each stratum, carbon per hectare values were applied to these units and then divided by the total area of the treatment unit to get a stand-scale per hectare carbon value. Carbon stored in planted live trees within the harvested areas in the CC, GS and GR treatments was estimated using data collected in 2016 within the regeneration plots, whereas carbon for Dead Standing Trees (DST) was considered negligible based on field observation of the treatment units and exploratory database analysis. Carbon estimations for LT here followed the same protocols as the retained areas, estimating the biomass using species-specific allometric equations (Ung et al., 2008), adjusting by the plot area factor and multiplying by 0.5 (Pearson et al., 2005). To estimate 2020 carbon stores, the 2016 measurements were used to parameterize growth and yield curves within SORTIE 7.05 (Canham & Murphy, 2001) to estimate 2020 carbon stores (Appendix D). To estimate carbon in the Dead standing trees (DST), I used Equation 1, which incorporates decay-related variables that helps to reduce uncertainty when accounting for carbon in this pool (Russell et al., 2015): Equation 1. Carbon on Dead Standing Trees (DST). Adjusted from Russell et al., 2015. ࡰࡿࢀ࡯ = ࡰࡿࢀ࡮ࡵࡻ ∗ ࡰ࡯ࡾࡲࡿࡰࢀ࢑࢔ ∗ ࡯ − ࡯ࡻࡺ࡯ࡿࡰࢀ࢑࢔ Where, DSTC is dead standing carbon (MgC ha-1), DSTBIO is the dead tree biomass which was obtained from species-specific allometric equations (Ung et al., 2008), and DCRFDSTkn is the decay class reduction factor, which refers to the ratio of the density of a decomposed piece to that of a non-decomposed piece of a given species group (i.e. softwood versus hardwood) (Russell et al., 2015). The tree condition of all recorded DSTs was based on the British Columbia’s wildlife tree classification (B.C. Ministry of Forests and Range & B.C. Ministry of Environment, 2010) and used as a reference to obtain the values from Harmon et al., (2011). 13 C-CONCDSTkn values are the C concentrations for DST in a given species group (i.e. softwood versus hardwood), which change throughout the decomposition process (Russell et al., 2015), and were obtained from Harmon et al., (2013). DST estimates for all trees were summed for the plot, and then expanded based on the plot dimensions. 2.3.2. Live and dead roots Data from live and dead roots was not collected during fieldwork. However, to account for the carbon in the live roots the estimations were as a function of total live tree biomass on a per-unit-area basis (MgC ha-1). I used equations for conifer and deciduous species (Equation 2) from cold temperate and boreal forests (Li et al., 2003). Equation 2. Equations and error terms to estimate root biomass for softwood and hardwood (Li et al., 2003). ࡾ࡮࢙ = ૙. ૛૛૛ ࡸࢀ࢙ (0.222 ± 0.004) Estimation for softwood species ࡾ࡮ࢎ = ૚. ૞ૠ૟ ࡸࢀ૙.૟૚૞ (1.576 ± 0:267; 0.615 ± 0.035) Estimation for hardwood species ࢎ Where, ‫ܶܮ‬௦ and ‫ܶܮ‬௛ are total live tree biomass (MgC ha-1) of softwood and hardwood, respectively, which then was converted to carbon. The carbon remaining in the dead roots (MgC ha-1) was estimated using the Equation 2 for biomass and using a time dependent decay function (Equation 3) (Wang et al., 2012). This pool entails the combined contribution from the harvested and retained forested areas: 1) to estimate root biomass from the harvested trees in the CC, GS and GR treatments, I used the preharvest data which was further adjusted by the proportion of the harvested areas at each treatment; 2) in retained areas, root biomass was estimated from dead trees located in the postharvest sampling plots, with further adjustments by the proportion of the retained areas. Equation 3. Equation to estimate carbon remaining in decaying dead roots at different times since harvest (Wang et al., 2012). ࢅ = ࢄ ࢋ(ି૙.૚૙૝૝∗࢟ࢋࢇ࢘) , where Y, is the carbon remaining; X, is the initial quantity of carbon using 2001 as the reference year; e, is the base of natural logarithms (2.718); -0.1044, is the decay coefficient; and year is the number of years since harvest. 14 2.3.3. Downed woody material The CWD and FWD sampling followed the same methods as Jull et al., (2001). Volume per hectare was estimated from each single-line transect surveyed using the Equation 4 and Equation 5. Carbon (MgC ha-1) was estimated for each transect and was calculated using Equation 6, by multiplying the volume (m3*ha-1), the decay class reduction factor (DCRFCWDkn), the carbon concentration (C-CONCCWDkn) and the decay class densities (WDCWD) (Russell et al., 2015). Equation 4.Volume estimation (m3*ha-1) for round pieces (Marshall et al., 2000) ࡯ࢃࡰ ࢜࢕࢒࢛࢓ࢋ࢘࢕࢛࢔ࢊ (࢓૜ ି૚ ) ∗ ࢎࢇ ࣊૛ = ∗ ෍ ࢊ૛࢏ ૡ∗ࡸ ࢐ୀ૚ where d is the diameter of the round piece (cm); L is the length of the transect (24m). Equation 5. Volume estimation (m3*ha-1) for odd pieces (Marshall et al., 2000) ࡯ࢃࡰ ࢜࢕࢒࢛࢓ࢋ࢕ࢊࢊ (࢓૜ ∗ ࢎࢇି૚ ) = ෍ ࢐ୀ૚ ࢃ∗ࡴ ࡸ where W is the width (cm) of the odd-shape piece; H is the height (cm) and L is the length of the transect (24m). Equation 6. Carbon on downed woody debris (CWD), adjusted from Russell et al., 2015. ࡯ࢃࡰ࡯ = ‫ ࡯ ∗ ࢔࢑ࡰࢃ࡯ࡲࡾ࡯ࡰ ∗ ܔܗ܄‬− ࡯ࡻࡺ࡯࡯ࢃࡰ࢑࢔ ∗ ࢃࡰ࡯ࢃࡰ where CWDC is the carbon (MgC ha-1); Vol is the total volume (m3*ha-1) in a given transect; DCRFCWDkn is the decay class reduction factor of a given species group (obtained from Harmon et al., (2011); C- CONCCWDkn is the C concentration for CWD in a given species group and was obtained from Harmon et al., (2013); WDCWD is a species-specific decay class density and was obtained from Harmon et al., (2008). Each PSP had associated three CWD transects whose carbon was summed for a final value per plot; subsequently, based on the PSP values, I calculated the mean, SD, and SE from these plots by treatment and site. CWD carbon estimations were based on the 2001 databases from Jull et al., (2001) and the 2020 fieldwork. To estimate the pre-harvest (1999) values, I used the carbon stocks from the 15 unharvested treatment (control) at each site from 2001 and applied them to the other treatments. Comparing treatments with historical control data can provide valid estimates of outcomes in the absence of the harvesting treatments (Paulus et al., 2014). Since the transects are ~20 years old and had little to no maintenance, not all the transects and CWD pieces from the 2001 survey could be found in 2020. To compensate for the missing information and to account for differences in the sampling effort between the two years, I adjusted the carbon in each plot by dividing it by its respective proportion of sampled transects (Appendix A, Table 2). Carbon per plot was used to estimate the mean, SD, and SE by treatment and site. This procedure allowed more accurate stock estimation, considering variations in survey efforts. 2.3.4. Forest floor The carbon estimations (g C) from the FF laboratory results were transformed into MgC and the methods are consistent with Matsuzaki et al. (2013). To estimate the area sampled (m2), I accounted for the number of subsamples in each sample and the core cylinder (m2) area used to collect the FF. To extrapolate the stocks on a per-hectare basis, I used expansion factors indicating the area each sample represents. These were calculated as the area of one hectare (10000 m2) divided by the total area sampled. This value was subsequently multiplied by the carbon content obtained in the lab to get MgC ha-1, which was further summarized by treatment. 2.3.5. Forest carbon accumulation rate The forest carbon accumulation rate (FCrate) was defined as the average annual change in carbon stocks following harvest (Hoover & Smith, 2023). It was calculated from the remeasurement of the same pools in 2020 that were originally surveyed in 2001, expressed as Mg C ha -1 yr-1 (Equation 7). Equation 7. Estimation of the forest carbon accumulation rate ࡲ࡯࢘ࢇ࢚ࢋ = (ࡸࢀ࢚૛ + ࡰࢀ࢚૛ + ࡸࡾ࢚૛ + ࡰࡾ࢚૛ + ࡯ࢃࡰ࢚૛ ) − (ࡸࢀ࢚૚ + ࡰࢀ࢚૚ + ࡸࡾ࢚૚ + ࡰࡾ࢚૚ + ࡯ࢃࡰ࢚૚ ) ࢤ࢚ Where, LT and DT is the carbon stock in the live and dead trees, LR and DR is the carbon stock in live and dead roots respectively, and CWD is the carbon stock in the coarse woody debris. 16 t1 is year one post-harvest, t2 is year 19 post-harvest, and Δt is 18. All pools were estimated on a per hectare basis (MgC ha-1). 2.3.6. Carbon stock error propagation estimation Given the high variability from each pool, the combination of errors is not linear, therefore it cannot be summed in a straightforward calculation. To estimate the total error associated to each total carbon mean (MgC ha-1) per treatment, I used as input the standard error from each pool and applied the Monte Carlo method. This technique uses repetitive calculations to find the empirical distribution of the data by randomly changing the input (Anderson, 1976). 2.3.7. Linear mixed-effect model to estimate differences in Carbon stock 2020 I used a mixed-effect model to test for significant differences among the pools and the treatments 19 years after harvest, using the pool, treatment, and their interaction as fixed effects, while using the site and the treatment as random effects to account for other sources of variability. All the analyses were conducted in R (RStudio, 2020) and considered significant when α < 0.05. In addition, I used the estimate marginal means (emmeans), a post-hoc test derived from the best fitted model to identify whether a given interaction of treatment/pool was significant or not. To control for the normality of the residuals (heteroscedasticity), I included a function where the residual variance by pool was weighted (Appendix F). lme.allsites_lme3 <- lme(Carbon_mean ~ Treatment*Pool, random = ~1 | Site/Treatment, data=Carbon_dataset, weights = varIdent(form = ~ 1 | Pool) ) 2.3.8. Edge effect on carbon accumulation in live trees To assess the edge effect on carbon accumulation in live trees, I accounted for the canopy conditions of the edge and closed, from each tree within the PSP in the partial harvesting treatments (i.e. GS and GR). I removed the trees in Edge from natural gaps (as natural disturbance), to assess only the edge effect from gaps due to harvesting. In addition, I removed from the analysis the trees from the open canopy conditions since these were mostly planted trees located in the managed harvesting opening gaps To calculate the individual tree carbon accumulation rate (MgC yr-1), I used the Equation 8 and the estimated carbon (MgC) per tree from the PSPs in 2001 and 2020 from both partial 17 harvest treatments. These values were subsequently used to estimate the mean carbon accumulation rate per treatment and per canopy condition within each treatment. Equation 8. Estimation of individual live trees' carbon accumulation rate ࡸࢀ࢘ࢇ࢚ࢋ(࢔) = (ࡸࢀ࢚૛࢔ − ࡸࢀ࢚૚࢔ ) ࢤ࢚ where LTt1n is the carbon stock (MgC) one year post harvest in each individual tree (n); LTt2n is the carbon stock 19 years post- harvest from the same tree (n); Δt is 18 years. The results were standardized to the area of the PSPs within each treatment across all sites. To test whether the canopy conditions within each treatment and the treatments themselves affected the live tree carbon accumulation differentially, I used the Kruskal-Wallis test, a nonparametric analog to the one-way ANOVA to compare the groups of interest (Hui, 2018). Furthermore, I applied a Dunn test, a post-hoc test for multiple pairwise comparisons (Dinno, 2015) and compared 1) the mean rates between the treatments without accounting for the canopy conditions, and 2) the mean rates given in two canopy conditions within each treatment. 18 3. Chapter Three. Results 3.1.1. Carbon stocks 19 years after harvesting As expected, 19 years after harvest, UN and GS stored the most total ecosystem carbon (475.3 ± 19.7 and 449.5 ± 32.2 MgC ha-1, respectively) (Figure 6a, Table 4), with no statistical differences in all measured pools (Appendix F). The stocks were followed by GR and CC (192.9 ± 14.6 MgC ha-1 and 36.3 ± 11.9, respectively) (Appendix C, Table 4), where partial harvest and UN showed significant differences in the pools LT, LR, DST and DR compared to CC (Appendix F). UN had higher stocks than GS across Minnow and East Twin, but not in the oldest site, Lunate, where 19 years after harvesting the GS stored more carbon than UN (571.45 ± 62.17 MgC ha1 and 514.61 ± 31.13 MgC ha-1, respectively) (Figure 6b), which was likely due to higher pre- harvest carbon stores within the areas subjected to GS. Across all three sites, live trees contributed an average of 59% of the total ecosystem carbon in unharvested forests. As expected, as the levels of retention decreased so did the LT contribution, where it accounted for 63% in GS, and 28% in GR (Figure 6a; Appendix C, Table 4). In treatments with little or no forest retention (GR and CC), CWD made the largest contribution to the total carbon stores (42% and 45%, respectively) (Figure 6b). Forest floor stores remained fairly consistent across all treatments 19 years after harvest, with a larger contribution to total ecosystem carbon in CC with 45% (Appendix C, Table 4). 19 Figure 6. Contribution of different carbon pools to mean total ecosystem carbon (MgC ha-1) 19 years after harvest in different harvesting treatments in the Northern Wetbelt for A) across all three sites and B) Each study site. Treatments are clearcut (CC), Group retention (GR) Group selection (GS) and unharvest control (UN). Sites are East Twin, Lunate, Minnow. Colors on bars represent different carbon pools (MgC ha -1): LT (live trees), LR (live roots), DR (dead roots), DST (dead standing trees), CWD (coarse woody debris), FWD (fine woody debris), FF (forest floor). 3.1.2. Carbon dynamics through time Contrary to my prediction, none of the three harvested treatments showed an increase in total ecosystem carbon from one year to 19 years after harvest (Table 6, Figure 7a). This might be partially explained by the opposing dynamics across the most important pools, meaning the gain in the live tree pool did not exceed the loss in others, notably the CWD. One year after harvest, CWD increased proportionally to the amount harvested in each TU, with CC storing the highest stocks in this pool. However, in 2020 this pool showed an important decrease especially in CC with over a 50 percent reduction, indicating a large source of ecosystem carbon losses (Table 6, Figure 7b). On the other hand, live tree C stocks in the harvested treatments experienced a decrease from harvesting, but by 2020 all treatments showed an increase in this pool. The increase in these stocks is mostly attributed to the growth from the mature live trees in the retained areas, while 20 planted trees in the regeneration areas contributed very little of these stocks (Appendix C, Table 5, c). Forest floor carbon stores remained fairly consistent across all treatments between year 10 (first year pool was measured) (Matsuzaki et al., 2013) and year 19 after harvest. CC had the largest proportional increase of 0.4 (from 37.5 to 52.5 MgC ha-1). GS and GR also experienced slightly lower proportional increases (0.35 and 0.25, respectively), while the unharvested control showed almost no change (Figure 7Figure 7d). 21 Figure 7. A) Total ecosystem carbon stocks over time (MgC ha-1) per treatment across all sites B). Carbon stocks change over time for CWD. C) Carbon stocks change over time for live trees. D) Forest floor carbon stocks 10 and 19 years after harvest. 10-year post-harvest data from Matsuzaki et al., 2013. All stocks are in MgC ha-1. 22 Converting total carbon stock from different survey years (2001 and 2020) to carbon accumulation rate, I found that losses in CC (-7.1 ± 1.75 Mg C ha -1 yr-1) were at least five times the rate compared to GR (-1.3 ± 1.21 Mg C ha -1 yr-1), GS (-1.3 ± 2.5 Mg C ha -1 yr-1) and UN (-0.4 ± 1.6 Mg C ha -1 yr-1) (Figure 8). This negative rate was driven by a loss of carbon in CWD, DR, and DST that was not compensated by the carbon gain in LT and LR (Figure 7b), reinforcing my findings about the opposite effects from some pools in the total ecosystem carbon, leading to no recovery in their stocks towards their pre-harvest values within a 19 year window. Figure 8. Mean and standard error forest carbon accumulation rate by treatment (Mg C ha-1 year-1) across all sites. CC (clearcut), GR, (group retention), GS (group selection), UN (unharvested). The rate was estimated based on the live trees (LT) dead standing trees (DST), live roots (LR), dead roots (DR) and coarse woody debris (CWD) from 2001 and 2020. 23 3.1.3. Canopy condition effects on individual live tree carbon accumulation rate I found both partial harvest treatments created all canopy conditions but in different proportions; mature live trees were mostly in closed canopy but in different ratios: in UN it was 99.5%, GS 62%, and GR 61% (Figure 9a). As expected, GR had the highest proportion of trees in the edge (38%), followed by GS (36%) and UN (0.4%), however GR also had fewer individual trees contributing to a lower FCrate (Appendix E, Table 8a). Contrary to expectations, I found the high retention treatment (GS) created the gap conditions that promoted the greatest carbon accumulation rate in LT, which as a result, contributed to higher total carbon stocks in this treatment. Among both of these treatments, I found significant differences, where GR had a lower mean rate (0.0007 ± 0.00003 MgC yr-1) than GS (0.002 ± 0.00007 MgC yr-1), and GS had a higher rate than the unharvested stands (0.0017 ± 0.00005 MgC ha-1 Yr-1), Kruskal-Wallis = Chi square =182.19, df = 2, p < 0.01) (Appendix E, Table 8b; Figure 9b). I also found significant differences in the rate between the canopy conditions within the treatments, where trees located at the edge had higher accumulation rates than trees in closed canopy (Dunn post-hoc test, p < 0.001) (Appendix E, Table 8a; Figure 9b). Altogether, my results show that, as predicted, harvesting gaps had a positive edge effect on tree growth, with trees at the edge growing faster than the ones in the interior, leading to enhanced FC rates in both treatments. However, the rates are also affected by the overall amount of retained trees whereby GS benefits the most by showing the highest rate at the stand level. 24 Figure 9. A) Ratio of sampled trees by canopy condition in each treatment. The results were standardized to the area of the PSPs within each treatment across all sites. B) mean and standard error of live tree carbon 25 accumulation rate (MgC tree-1 yr-1) grouped by treatment (GS, GR) and by canopy conditions: Edge (E) and Closed (C). 4. Chapter Four. Discussion 4.1. Carbon stocks 19 years after harvesting My research aims to provide insights into the long-term forest carbon dynamics under different silviculture practices, up to 19 years after harvest in the IWB. To do this, I measured and examined carbon stocks from different pools and estimated the carbon accumulation rates, including individual tree responses to canopy conditions and analyze historical data. My results support that the carbon stored in these temperate rainforests can fluctuate depending on the forest management practice. In unharvested carbon forests, I found that total carbon stocks are 426.5 ± 19.6 MgC ha-1. After 20 years, the total C stocks in high retention (GS) were reduced by 18% of pre-harvest levels, the low retention by 59% and CC by 82%. Perhaps most surprisingly, I found net negative carbon flux in all harvest treatments, even in light harvest, 19 years after logging occurred. Although the forest carbon stocks from the partial harvest treatments have not shown a recovery, they have been quite stable with significantly less emissions in high retention. In the same forests stands, Matsuzaki et al. (2013) found that UN and GS stored more carbon than GR and CC immediately after harvest. Their results also highlighted the importance of these temperate rainforests as carbon stores and that harvesting led to a significant reduction in LT stocks by removing large trees, the most dynamic pool in carbon sequestration. In addition, recent research from other field plot calculations in the managed stands from the IWB region reported an average carbon estimate from the live above and belowground biomass of 433.81 Mg C ha-1 (DellaSala et al., 2022). Other ICH stands in the interior of BC reported total ecosystem carbon estimates from 260.8 MgC ha-1 to 564.7 MgC ha-1, with live trees as the pool with the largest percentage of storage (Roach et al., 2021). My findings are comparable to those reported elsewhere, with differences possibly due to the methods used to estimate carbon, forest stand characteristics (stand age, site productivity), uncertainty associated with differences in the collection and management of data, etc. My results highlight that silvicultural practices including partial harvest have a negative effect on the forest carbon stocks in mature to old-growth forest stands in rainforests from the interior 26 Northern Wetbelt. This reduction is proportional to the amount of forest harvested and within ~20 yr after harvest, it shows no recovery towards their pre-harvest stores. 4.2. Carbon dynamics through time One significant aspect of my study that adds to previous research is the contribution of different carbon pools through time and their interaction. I found there were opposing trends of increasing stocks in the live trees and the decreasing stocks in the CWD, with the net result of no recovery of ecosystem carbon stores compared to the pre-harvest condition. Understanding the dynamics and the interactions between the pools is important to gain insights into the ability of mature rainforests to recover carbon stocks after alternative logging practices. In unharvested forests, I found that 59% of the total forest carbon stock is stored in the LT. In addition, across all sites and harvesting treatments, this pool showed an increase in their stocks compared to the one year after harvest, where my results showed that high retention of mature live trees is key to maintaining and increasing the carbon stocks up to 20 years of harvest. Because of a large tree’s allometry, even smaller increases in diameter correspond to significant increments in the carbon stocks. Mildrexler et al., (2020) reported that while a tree with 25cm DBH could store 90-121 kg of aboveground carbon, a tree with 50 cm DBH could hold between 541-583 kg of AGC. The patterns found in the LT in the Wetbelt trials are aligned with other research across the province. For instance, Simard et al. (2020) reported that live trees accounted for 50% of the total carbon stocks in different forests across British Columbia, which experienced reduced stocks in proportion to harvesting intensity. Specifically, in the forests of the IWB, DellaSala et al. (2022) found that about half of the carbon is sequestered in LT, and the rest is grouped in the dead pools and soils, highlighting the importance of leaving CWD in forests after harvest for climate mitigation. In the Northern Wetbelt trials, using one-year after harvest data, Matsuzaki et al., (2013) reported a similar pattern, with live trees accounting for > 50% of the ecosystem carbon whose stocks were reducing proportionally to the amount of forest harvested. CWD on the other hand, showed an opposite trend to the LT with decreasing stocks over time. Due to harvesting, the partial harvest treatments and CC increased the carbon stocks in this pool. However, ~20yr later this pool has decreased the stocks particularly in CC losing 63% 27 of the stores. Based on the one-year after harvest data for the Northern Wetbelt trial, Matsuzaki et al., (2013) although reported lower stocks compared to the ones from this research, they found the same pattern with CC having the greatest stocks, whereby differences between the stocks are likely due to the carbon estimation method. In addition, Farnell et al., (2020) found similar patterns in the CWD from an ITR, with clearcuts having a larger proportion of CWD pieces one-year after harvesting compared to unharvested stands. However, 27 years after harvest, the pieces’ decay class varied with the retention level, with CC having a smaller proportion of less decayed pieces than UN stands. Environmental factors such as differences in the temperature and humidity between treatments could also explain the more rapid decomposition and therefore the carbon depletion in CC than in the other treatments, where research has demonstrated a positive effect of temperature on decomposition rates. For instance, Magnússon et al., (2016) found decreasing carbon stocks in CWD under increasing temperatures in subalpine forests. Moreover, the CWD pool is not only relevant for its carbon storing capacity, but also for other values such as providing suitable habitat for wildlife. Farnell et al., (2020) reported that inland temperate rainforest stands with high levels of retention (>70%) had similar favorable habitat conditions as unharvested areas for the American marten (Martes americana), a CWDdependent mammal. Furthermore, Magnússon et al., (2016) found that CWD may also enhance the microbial community’s abundance and diversity, which helps decompose these pieces and can contribute to the soil C stocks. Therefore, if enhancing multiple values in forest management is a goal, partial harvesting can achieve this through the maintenance of high CWD stocks. Although the carbon stored below ground (in the forest floor) tends to be slow to respond to disturbances, it had important contributions to the total forest carbon and was also affected by harvesting practices. In the case of the Northern Wetbelt trial, there is no one-year after-harvest data from this pool to estimate the stocks. Instead, I compared my FF results to the values 10 years after harvest from Matsuzaki et al., (2013), where my stocks were slightly higher, suggesting a carbon input from other pools (e.g. CWD, snags, roots, etc.) and a valuable contribution to the total carbon, particularly in CC with 45%. Nave et al., (2010) found that FF is more susceptible to carbon loss due to logging than mineral soils, where FF decreased by 30 28 ± 6% following timber harvesting; in addition, this loss was smaller in coniferous/mixed stands (-20%) as opposed to hardwoods (-36%). Contrary to my prediction, I found negative FCrate, for all treatments, implying overall carbon losses. The rate varied among the treatments and depended on the forest stocks from 2001 and 2020. Based on the contribution and the dynamic from each individual pool, these results highlight the importance of estimating carbon flux and stocks using as many pools as possible. Using only a few of them (e.g. live trees) could lead to biased or inaccurate results. When comparing the FCrate among the treatments, I found GR and GS had smaller net carbon losses than CC but higher than UN. This may be partially explained by tree retention per se, where live trees are crucial in contributing to the total forest carbon stocks. In addition, since partial harvest results in multi-cohort stands from natural and planted regeneration, the carbon accumulation in LT is likely enhanced by the growth of younger trees, highlighting the potential influence of the forest age on carbon accumulation rates (Hoover & Smith, 2023). Hoover & Smith (2023) found higher rates of carbon accumulation in younger forests compared to older stands, where the rates ranged from 3 tC ha -1 yr-1 to 0.03 tC ha -1 yr-1, whereas in the oldest stands the highest accumulation rate reached was 1.33 tC ha -1 yr-1. In addition, Gray et al., (2016) reported in the Pacific Northwest National Forests, USA, that carbon accumulation rate was affected by the site productivity, the plant community type, and the forest age, reaching the greatest values in forests with larger proportions of younger stands and declining as the forests mature while reaching their stocks’ maximum capacity. Despite having an increase in the LT and have no apparent disturbances from forest management, the unharvested treatment also showed net carbon losses. The negative flux was driven by decomposition of the dead carbon pools, which includes DST, CWD, and DR, where it is likely these stands had higher than normal amounts of carbon stocks at the time the experiment was established, possibly due to a widespread disturbance in the forests, such as windthrow, within the preceding decades, and then experiencing carbon decrease from CWD, DST and DR pools. Another potential factor contributing to faster decay in these pools is global warming, where Magnússon et al., (2016) reported positive correlations between increasing temperatures and wood decay. A final interacting factor in the carbon stocks trend is the site 29 and stand-level differences across all study sites, where Gray et al., (2016) reported differences in carbon given the forest age. Altogether, the unharvested treatment had negative FCrate but it had the largest carbon stocks. These findings are aligned on the one hand, with Luyssaert et al. (2008) who highlighted the critical role of old-growth forests as carbon sinks; on the other hand, with Gray et al., (2016) and D’Amato et al., (2011), who suggested that despite large trees are important carbon stocks, they play a minor role in additional C accumulation, which can be relevant when managing forests for a particular objective/value. 4.3. Canopy condition effect on individual live tree carbon accumulation rate The high tree retention (GS) creates the best conditions that benefit carbon sequestration, especially near the edge. These results, although not intuitive, might be explained by the interaction of multiple variables such as the amount of forest retained, presence of large trees, stem density, basal area, forest age, and size of the gaps, among others (Meeussen et al., 2021; Saeed et al., 2019; Ziter et al., 2014). In the Northern Wetbelt trial, higher retention of large trees contributes to a more significant carbon accumulation rate, which could be potentially enhanced by the contribution of the ones located near the edge, whose rates would be higher due to access of more resources such as light availability (Reinmann & Hutyra, 2017; Ziter et al., 2014). My results show that higher landscape fragmentation and lower forest retention can weaken the forest carbon stores to a point where the enhanced growth from the LT at the edges cannot compensate for the loss. These gaps can also alter microclimatic variables (e.g. temperature, humidity), increasing the trees’ stress levels, particularly during drought and hot seasons (Reinmann & Hutyra, 2017). Gap impacts forest in other ways as well, for instance Coates & Burton (1997) found that different gap sizes come with a particular combination of light regimes and other microclimate attributes which leads to vegetation responses and affects seed germination and seedling mortality. Meeussen et al. (2021) found that higher aboveground carbon stocks near the forest edges than the interior across Europe were attributed to increased stem density. In the temperate broadleaf 30 forests in the USA, Reinmann & Hutyra, (2017) consistently found more forest growth at the edge than the interior, however, the magnitude of that difference was subject to annual variations in climate, with heat stress as an important explanatory variable. Projections predict that by the end of the century, the forest edge growth could be reduced by one-third due to heat stress. Therefore, in the case of these rainforests having multicohort stands, high forest retention would lead to a higher carbon accumulation rate in live trees, particularly at the edge from gaps due to harvesting. Nevertheless, given the site’s environmental conditions, a ~20-year window is not long enough timeframe for the trees in the partial harvest treatments to recover their preharvest stocks. 5. Chapter Five 5.1. Conclusions In the mature to old-growth temperate rainforests from the Northern Wetbelt region, the stands have developed from complex interactions of natural disturbances through centuries, allowing them to build up very large carbon stocks across the pools, highlighting one of the many values of these rare and important ecosystems. This research is one of the few to report field data on the longer-term effects of timber harvesting on the forest carbon stocks in the above and belowground pools for the ITR in the Wetbelt region. Quantifying carbon stocks using empirical data demonstrates how harvesting practices affect pools differentially and highlights the crucial role of the temperate rainforests in storing large amounts of carbon, while using these ecosystems as a climate-based solution when accounting for carbon stocks contributions to Canada’s National Determined Contribution (NDC) under the Paris Climate Agreement. As expected, all aboveground pools from partial harvest treatments stored more carbon stocks than clearcut systems. However, within ~20-years, although all treatments kept the carbon stocks relatively stable, they did not recover their pre-harvest carbon stocks. The opposite effects of the increasing carbon stocks in LT and decreasing in the CWD, along with the net losses from the FCrate (including in the UN), highlight the complexity of these ecosystems and 31 provide insights about the susceptibility and recovery time to disturbances such as harvesting in these stands. When managing these forests for economic purposes, logging practices that aim for high forest retention (GS), although they show some potential for forest management, my results show that it will take more than 20 years to recover the carbon stocks to pre-harvest values. In addition to these recovery times due to harvesting, it is crucial to account for other variables such as structural stage, environmental conditions, and natural disturbances (e.g. wildfires and insects) (Burton & Boulanger, 2018; DellaSala et al., 2022), whose interaction could lead to prolonged recovery times towards the pre-harvest carbon stocks. Live trees had positive carbon accumulation rates despite the net losses in forest carbon accumulation rates. When accounting for the C accumulation in the partial harvest treatments, the amount of retained trees was more important than the canopy conditions created by the edge effect. These results underscore the importance of retaining large, mature trees for a longterm carbon reservoir, while playing crucial roles in sustaining complex food webs and biodiversity. These values should be considered when managing temperate rainforests with similar characteristics. Finally, the long-term research trials throughout the province of BC are critical study sites since they provide baseline data that can produce high-value information to support decisionmaking toward more sustainable forest management. Results from the Interior Wetbelt research trials suggest that given the importance of these inland temperate rainforests for multiple values and its critical status, a precautionary approach for management should be taken, aiming to protect the complexity and diversity inherent in these ecosystems. Actions at the landscape scale could include the retention of Old Growth Management Areas (OGMAs), riparian areas, visual corridors and rare ecosystems (Stevenson et al., 2011). 32 5.2. Project limitations and future work One of the limitations was time for data collection and my inability to include more pools such as mineral soil. The mineral pools are important contributors in temperate forests; however they are the most stable pool in response to harvesting (Nave et al., 2010). Collecting data from more pools and data availability for pre- and post-harvest for all pools could have provided more insights into the long-term effects of silviculture practices on the carbon dynamics at the IWB. Due to time limitations, I did not survey planted trees located in the openings of the partial harvest and CC treatments. Instead, I used a database with the last survey in 2016 from Jull et al. (2001). Collecting data from the planted trees during the same survey year rather than modelling carbon based on previous surveys could provide more accurate outcomes. In addition, collecting data of heights in the mature trees in the PSP instead of modelling them would potentially improve the accuracy of the biomass and carbon content. Fieldwork lasted 3 ½ months with two permanent members and occasional volunteers. It is likely that ‘more hands-on-ground’ could have made a difference regarding the amount of data collected. For instance, many of the CWD transects were not re-located in 2020, though much effort was afforded to this. However, this option came with its challenges and limitations, mostly due to strict COVID-19 policies. Overall, this project answered relevant research questions, filled gaps, and provided information to support decision-making. Among the challenges during the discussion and comparison with similar research was accounting for the initial ecosystem condition and history of these ITR, which are unique stands and have been historically unmanaged. Many of the study sites from other research have been influenced by past forest management or large disturbances that have depleted the legacy of carbon pools, affecting the stocks, going forward. Therefore, it is important to acknowledge the history of the stand, to make more reliable comparisons. Finally, adding more analysis will allow a better understanding of the carbon dynamics: for instance, accounting for the wood harvesting products allows tracking of carbon outside the ecosystem and could be further linked to economic analysis to support decision-making. Other analyses accounting for the canopy condition as an experimental parameter might explore the growth response as a function of the initial tree diameter and the species. 33 Literature cited Anderson, G. . (1976). Error propagation by the Monte Carlo method in geochemical calculations. Geochimica et Cosmochimica Acta, 40(12), 1533–1538. https://doi.org/10.1016/0016-7037(76)90092-2 Andrieu, E., Cabanettes, A., Alignier, A., Van Halder, I., Alard, D., Archaux, F., Barbaro, L., Bouget, C., Bailey, S., Corcket, E., Deconchat, M., Vigan, M., Villemey, A., & Ouin, A. (2018). Edge contrast does not modulate edge effect on plants and pollinators. Basic and Applied Ecology, 27(November), 83–95. https://doi.org/10.1016/j.baae.2017.11.003 B.C. Ministry of Forests and Range, & B.C. Ministry of Environment. (2010). Field Manual for Describing Terrestrial Ecosystems. 2nd Edition. In Environment (second edi). Bartowitz, K. J., Walsh, E. S., Stenzel, J. E., Kolden, C. A., & Hudiburg, T. W. (2022). Forest Carbon Emission Sources Are Not Equal: Putting Fire, Harvest, and Fossil Fuel Emissions in Context. Frontiers in Forests and Global Change, 5(May), 1–11. https://doi.org/10.3389/ffgc.2022.867112 Beese, W. J., Deal, J., Dunsworth, B. G., Mitchell, S. J., & Philpott, T. J. (2019). Two decades of variable retention in British Columbia: a review of its implementation and effectiveness for biodiversity conservation. Ecological Processes, 8(1), 33. https://doi.org/10.1186/s13717-019-0181-9 British Columbia. Ministry of Forests. Forest Practices Branch. (2003). Silvicultural Systems Handbook for British Columbia (Issue March). http://www.for.gov.bc.ca/HFP/Pubssilvsystems.htm%0AFurther Burton, P. J., & Boulanger, Y. (2018). Characterizing combined fire and insect outbreak disturbance regimes in British Columbia, Canada. Landscape Ecology, 33(11), 1997– 2011. https://doi.org/10.1007/s10980-018-0710-4 Canham, C. D., & Murphy, L. E. (2001). SORTIE-ND (7.05). Institute of Ecosystem Studies. http://www.sortie-nd.org/software/index.html Coates, K. D., & Burton, P. J. (1997). A gap-based approach for development of silvicultural 34 systems to address ecosystem management objectives. Forest Ecology and Management, 99(3), 337–354. https://doi.org/10.1016/S0378-1127(97)00113-8 COSEWIC. (2014). Caribou (Rangifer tarandus) COSEWIC assessment and status report 2014: chapter 1. https://www.canada.ca/en/environment-climate- change/services/species-risk-public-registry/cosewic-assessments-statusreports/caribou-2014/chapter-1.html Coxson, D., Goward, T., & Werner, J. R. (2020). The Inland Temperate Rainforest and Interior Wetbelt Biomes of Western North America. In Encyclopedia of the World’s Biomes (Issue January, pp. 88–102). Elsevier. https://doi.org/10.1016/B978-0-12-409548-9.12055-X D’Amato, A. W., Bradford, J. B., Fraver, S., & Palik, B. J. (2011). Forest management for mitigation and adaptation to climate change: Insights from long-term silviculture experiments. Forest Ecology and Management, 262(5), 803–816. https://doi.org/10.1016/j.foreco.2011.05.014 DellaSala, D. A., Alaback, P., Spribille, T., von Wehrden, H., & Nauman, R. S. (2011). Just what are Temperate and boreal rainforest? In D. A. DellaSala (Ed.), Temperate and Boreal Rainforests of the world: ecology and conservation (pp. 1–41). DellaSala, D., Keith, H., Sheehan, T., Strittholt, J., Mackey, B., Connolly, M., Werner, J. R., & Fredeen, A. L. (2022). Estimating carbon stocks and stock changes in Interior Wetbelt forests of British Columbia, Canada. Ecosphere, 13(4), 1–18. https://doi.org/10.1002/ecs2.4020 Dinno, A. (2015). Nonparametric Pairwise Multiple Comparisons in Independent Groups using Dunn’s Test. The Stata Journal: Promoting Communications on Statistics and Stata, 15(1), 292–300. https://doi.org/10.1177/1536867X1501500117 ESRI. (2011). ArcGIS Desktop: Release 10.5.1 (10.5.1). Environmental Systems Research Institute. :purple_heart: 1. https://support.esri.com/en/products/desktop/arcgis- desktop/arcmap/10-5 Farnell, I., Elkin, C., Lilles, E., Roberts, A. M., & Venter, M. (2020). The effects of variable retention forestry on coarse woody debris dynamics and concomitant impacts on 35 American marten habitat after 27 years. Canadian Journal of Forest Research, 50(9), 925–935. https://doi.org/10.1139/cjfr-2019-0417 Fedrowitz, K., Koricheva, J., Baker, S. C., Lindenmayer, D. B., Palik, B., Rosenvald, R., Beese, W., Franklin, J. F., Kouki, J., Macdonald, E., Messier, C., Sverdrup‐Thygeson, A., & Gustafsson, L. (2014). REVIEW: Can retention forestry help conserve biodiversity? A meta‐analysis. Journal of Applied Ecology, 51(6), 1669–1679. https://doi.org/10.1111/1365-2664.12289 Gray, A. N., Whittier, T. R., & Harmon, M. E. (2016). Carbon stocks and accumulation rates in Pacific Northwest forests: role of stand age, plant community, and productivity. Ecosphere, 7(1), 1–19. https://doi.org/10.1002/ecs2.1224 Gustafsson, L., Baker, S. C., Bauhus, J., Beese, W. J., Brodie, A., Kouki, J., Lindenmayer, D. B., Lõhmus, A., Pastur, G. M., Messier, C., Neyland, M., Palik, B., Sverdrup-Thygeson, A., Volney, W. J. A., Wayne, A., & Franklin, J. F. (2012). Retention Forestry to Maintain Multifunctional Forests: A World Perspective. BioScience, 62(7), 633–645. https://doi.org/10.1525/bio.2012.62.7.6 Harmon, M. E., Fasth, B., Woodall, C. W., & Sexton, J. (2013). Carbon concentration of standing and downed woody detritus: Effects of tree taxa, decay class, position, and tissue type. Forest Ecology and Management, 291, 259–267. https://doi.org/10.1016/j.foreco.2012.11.046 Harmon, M. E., & Sexton, J. (1996). Guidelines for Measurements of Woody Detritus in Forest Ecosystems. Publication No. 20. (Vol. 20, p. 73). Harmon, M. E., Woodall, C. W., Fasth, B., & Sexton, J. (2008). Woody detritus density and density reduction factors for tree species in the United States: a synthesis. In USDA General Technical Report NRS-29. USDA Forest Service Northern Research Station. https://doi.org/10.2737/NRS-GTR-29 Harmon, M. E., Woodall, C. W., Fasth, B., Sexton, J., & Yatkov, M. (2011). Differences between standing and downed dead tree wood density reduction factors: A comparison across decay classes and tree species. In U.S. Department of Agriculture, Forest Service, 36 Northern Research Station, Research Paper NRS-15 (Vol. 15). https://www.fs.usda.gov/treesearch/pubs/38699 Higgins, J., Thomas, J., Chandler, J., Cumpston, M., Li, T., Page, M., & Welch, V. (2023). Selecting studies and collecting data. In J. Higgins, J. Thomas, J. Chandler, M. Cumpston, T. Li, M. Page, & V. Welch (Eds.), Cochrane Handbook for Systematic Reviews of Interventions version 6.4. Available from www.training.cochrane.org/handbook. Hoover, C. M., & Smith, J. E. (2023). Aboveground live tree carbon stock and change in forests of conterminous United States: influence of stand age. Carbon Balance and Management, 18(1), 7. https://doi.org/10.1186/s13021-023-00227-z Hudiburg, T. W., Law, B. E., Moomaw, W. R., Harmon, M. E., & Stenzel, J. E. (2019). Meeting GHG reduction targets requires accounting for all forest sector emissions. Environmental Research Letters, 14(9), 095005. https://doi.org/10.1088/1748- 9326/ab28bb Hui, E. G. M. (2018). Learn R for applied statistics: With data visualizations, regressions, and statistics. In Learn R for Applied Statistics: With Data Visualizations, Regressions, and Statistics. Apress Media LLC. https://doi.org/10.1007/978-1-4842-4200-1 Jennings, S. (1999). Assessing forest canopies and understorey illumination: canopy closure, canopy cover and other measures. Forestry, 72(1), 59–74. https://doi.org/10.1093/forestry/72.1.59 Jull, M., Stevenson, S., & Rogers, B. (2002). Northern Wetbelt Silvicultural Systems Project ESTABLISHMENT REPORT. 2nd Approximation (Issue April 2002). Jull, M., Stevenson, S., Rogers, B., Sandborn, P., & Eastham, A. (2001). Establishment report (First Approximation): Summary of initial Pre-harvest conditions and harvest treatments. Land Owner Resources Centre. (2000). Conserving the Forest Interior: A Threatened Wildlife Habitat. In Extension Notes Ontario. http://www.lrconline.com/Extension_Notes_English/pdf/forInterior.pdf Li, Z., Kurz, W. A., Apps, M. J., & Beukema, S. J. (2003). Belowground biomass dynamics in the Carbon Budget Model of the Canadian Forest Sector: recent improvements and 37 implications for the estimation of NPP and NEP. Canadian Journal of Forest Research, 33(1), 126–136. https://doi.org/10.1139/x02-165 Lorber, J., Thomas-Van Gundy, M. A., & Croy, S. (2018). Characterizing effects of prescribed fire on forest canopy cover in the George Washington and Jefferson National Forests. In Research Paper NRS-31 (Issue June). https://doi.org/10.2737/NRS-RP-31 Luyssaert, S., Schulze, E. D., Börner, A., Knohl, A., Hessenmöller, D., Law, B. E., Ciais, P., & Grace, J. (2008). Old-growth forests as global carbon sinks. Nature, 455(7210), 213– 215. https://doi.org/10.1038/nature07276 Magnússon, R. Í., Tietema, A., Cornelissen, J. H. C., Hefting, M. M., & Kalbitz, K. (2016). Tamm Review: Sequestration of carbon from coarse woody debris in forest soils. Forest Ecology and Management, 377, 1–15. https://doi.org/10.1016/j.foreco.2016.06.033 Marshall, P. L., Davis, G., & LeMay, V. . (2000). Using Line Intersect Sampling for Coarse Woody Debris. In For. Serv. Tech. Rep. TR-003 (Vol. 128, Issue 3). https://www.webpages.uidaho.edu/for373new/pdfs/for373/lineintersectsampling_tr003. pdf Matsuzaki, E., Sanborn, P., Fredeen, A. L., Shaw, C. H., & Hawkins, C. (2013). Carbon stocks in managed and unmanaged old-growth western redcedar and western hemlock stands of Canada’s inland temperate rainforests. Forest Ecology and Management, 297, 108–119. https://doi.org/10.1016/j.foreco.2012.11.042 Meeussen, C., Govaert, S., Vanneste, T., Haesen, S., Van Meerbeek, K., Bollmann, K., Brunet, J., Calders, K., Cousins, S. A. O., Diekmann, M., Graae, B. J., Iacopetti, G., Lenoir, J., Orczewska, A., Ponette, Q., Plue, J., Selvi, F., Spicher, F., Sørensen, M. V., … De Frenne, P. (2021). Drivers of carbon stocks in forest edges across Europe. Science of The Total Environment, 759, 143497. https://doi.org/10.1016/j.scitotenv.2020.143497 Mildrexler, D. J., Berner, L. T., Law, B. E., Birdsey, R. A., & Moomaw, W. R. (2020). Large Trees Dominate Carbon Storage in Forests East of the Cascade Crest in the United States Pacific Northwest. Frontiers in Forests and Global Change, 3(November), 1–15. https://doi.org/10.3389/ffgc.2020.594274 38 Nave, L. E., Vance, E. D., Swanston, C. W., & Curtis, P. S. (2010). Harvest impacts on soil carbon storage in temperate forests. Forest Ecology and Management, 259(5), 857–866. https://doi.org/10.1016/j.foreco.2009.12.009 Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., Phillips, O. L., Shvidenko, A., Lewis, S. L., Canadell, J. G., Ciais, P., Jackson, R. B., Pacala, S. W., McGuire, A. D., Piao, S., Rautiainen, A., Sitch, S., & Hayes, D. (2011). A Large and Persistent Carbon Sink in the World’s Forests. Science, 333(6045), 988–993. https://doi.org/10.1126/science.1201609 Paulus, J. K., Dahabreh, I. J., Balk, E. M., Avendano, E. E., Lau, J., & Ip, S. (2014). Opportunities and challenges in using studies without a control group in comparative effectiveness reviews. Research Synthesis Methods, 5(2), 152–161. https://doi.org/10.1002/jrsm.1101 Pearson, T., Walker, S., & Brown, S. (2005). Sourcebook for Land use, Land-use change and forestry projects. In Winrock International and the BioCarbon Fund of the World Bank 57 (2005) (Vol. 21, Issue 3). http://wbcarbonfinance.org/docs/Background_LULUCF_Sourcebook_compressed.pdf Reinmann, A. B., & Hutyra, L. R. (2017). Edge effects enhance carbon uptake and its vulnerability to climate change in temperate broadleaf forests. Proceedings of the National Academy of Sciences, 114(1), 107–112. https://doi.org/10.1073/pnas.1612369114 Roach, W. J., Simard, S. W., Defrenne, C. E., Pickles, B. J., Lavkulich, L. M., & Ryan, T. L. (2021). Tree Diversity, Site Index, and Carbon Storage Decrease With Aridity in DouglasFir Forests in Western Canada. Frontiers in Forests and Global Change, 4(June). https://doi.org/10.3389/ffgc.2021.682076 RStudio, T. (2020). RStudio: Integrated Development Environment for R (4.0.3). RStudio, PBC. http://www.rstudio.com/ Russell, M. B., Fraver, S., Aakala, T., Gove, J. H., Woodall, C. W., D’Amato, A. W., & Ducey, M. J. (2015). Quantifying carbon stores and decomposition in dead wood: A review. 39 Forest Ecology and Management, 350, 107–128. https://doi.org/10.1016/j.foreco.2015.04.033 Saeed, S., Yujun, S., Beckline, M., Chen, L., Zhang, B., Ahmad, A., Mannan, A., Khan, A., & Iqbal, A. (2019). Forest edge effect on biomass carbon along altitudinal gradients in Chinese Fir ( Cunninghamia lanceolata ): A study from Southeastern China. Carbon Management, 10(1), 11–22. https://doi.org/10.1080/17583004.2018.1537517 Simard, S. W., Roach, W. J., Defrenne, C. E., Pickles, B. J., Snyder, E. N., Robinson, A., & Lavkulich, L. M. (2020). Harvest Intensity Effects on Carbon Stocks and Biodiversity Are Dependent on Regional Climate in Douglas-Fir Forests of British Columbia. Frontiers in Forests and Global Change, 3. https://doi.org/10.3389/ffgc.2020.00088 Stevenson, S. K., Armleder, H. M., Arsenault, André, Coxon, D., DeLong, C. S., & Jull, M. (2011). Changing communities, changing values, changing uses. In British Columbia’s Inland Rainforest Ecology, Conservation, and Management (p. 456). UBC Press. Stevenson, S. K., Jull, M. J., & Rogers, B. J. (2006). Abundance and attributes of wildlife trees and coarse woody debris at three silvicultural systems study areas in the Interior CedarHemlock Zone, British Columbia. Forest Ecology and Management, 233(1), 176–191. https://doi.org/10.1016/j.foreco.2006.06.023 Stevenson, S., Rogers, B., & Jull, M. (2001). Field Data Collection Protocol Manual For internal use only : Northern Interior Wetbelt Silviculture Systems Project 2001. Thorpe, H. C., & Thomas, S. C. (2007). Partial harvesting in the Canadian boreal: Success will depend on stand dynamic responses. The Forestry Chronicle, 83(3), 319–325. https://doi.org/10.5558/tfc83319-3 Ung, C.-H., Bernier, P., & Guo, X.-J. (2008). Canadian national biomass equations: new parameter estimates that include British Columbia data. Canadian Journal of Forest Research, 38(5), 1123–1132. https://doi.org/10.1139/X07-224 Wang, G. G., Van Lear, D. H., Hu, H., & Kapeluck, P. R. (2012). Accounting Carbon Storage in Decaying Root Systems of Harvested Forests. AMBIO, 41(3), 284–291. https://doi.org/10.1007/s13280-011-0161-5 40 Wiensczyk, A. M. (2012). Status of British Columbia’s Longterm Silvicultural Systems Research Trials. http://www.forrex.org/sites/default/files/publications/full_issues/fr-1201.pdf Zhou, D., Zhao, S. Q., Liu, S., & Oeding, J. (2013). A meta-analysis on the impacts of partial cutting on forest structure and carbon storage. Biogeosciences, 10(6), 3691–3703. https://doi.org/10.5194/bg-10-3691-2013 Ziter, C., Bennett, E. M., & Gonzalez, A. (2014). Temperate forest fragments maintain aboveground carbon stocks out to the forest edge despite changes in community composition. Oecologia, 176(3), 893–902. https://doi.org/10.1007/s00442-014-3061-0 41 Appendix A Equations to model heights for trees I used quadratic equations to predict the height (H) of each tree whose DBH was measured in the summer of 2020. The equations are species-specific and are based on the historical data taken in 2001 provided by Jull et al., 2001. For a total of N trees, I modeled DBH vs heights and chose the model with highest R2 and lowest Akaike Information Criterion (AIC) value. Common name Scientific name Species code Subalpine fir Abies lasiocarpa Bl Western red cedar Thuja plicata Cw Western hemlock Tsuga heterophylla Hw Douglas maple Acer glabrum Mr Engelmann spruce Picea engelmannii Se Spruce hybrid Picea engelmannii x Picea glauca Sx Unknown U • Bl H = -0.00868(DBH)2 + 1.14156DBH -1.83397 • Cw H = -0.002377(DBH)2 + 0.610526(DBH) + 2.062073 • Hw H = -0.006642(DBH)2 + 0.960497(DBH) -2.017930 • Mr H = 0.04261(DBH) + 4.80009 • Se H = -0.004862(DBH)2 + 0.921360(DBH) -0.282847 • Sx H = -0.007334(DBH)2 + 1.072208(DBH) -1.063502 42 Allometric equations to estimate biomass for trees Biomass and carbon estimates from live and dead trees were calculated on a per-unit-area basis (MgC ha-1). Biomass for each individual tree was estimated from historical databases when the research trial was established (Jull et al., 2001), in combination with species-specific allometric equations (Ung et al., 2008) using height (H) and diameter at breast height (DBH). • Bl Biomass = (0.022*DBH^1.6469*H^1.1714)+(0.0061*DBH^1.8603*H^0.7693)+(0.0265*DBH^3. 6747*H^-1.5958)+(0.0509*DBH^2.9906*H^-1.2271) • Cw Biomass = (0.0188*DBH^1.3376*H^1.5293)+(0.0002*DBH^2.4369*H^1.1315)+(0.0611*DBH^1 .9208)+(0.1097*DBH^1.553) • Hw Biomass = (0.0113*DBH^1.9332*H^1.1125)+(0.0019*DBH^2.3356*H^0.6371)+(0.0609*DBH^2 .0021)+(0.2656*DBH^2.0107*H^-0.7963) • Mr Biomass = (0.0353*DBH^2.0249*H^0.7048)+(0.009*DBH^1.8677*H^0.7144)+(0.0448*DBH^2. 6855*H^-0.5911)+(0.0869*DBH^1.8541*H^-0.5491) • Se Biomass = (0.0133*DBH^1.3303*H^1.6877)+(0.0086*DBH^1.6216*H^0.8192)+(0.0428*DBH^2 .7965*H^-0.7328)+(0.0854*DBH^2.4388*H^-0.7630) • Sx Biomass = (0.0276*DBH^1.6868*H^1.0953)+(0.0101*DBH^1.8486*H^0.5525)+(0.0313*DBH^2 .9974*H^-1.0383)+(0.1379*DBH^2.3981*H^-1.0418) • U Biomass = (0.0283*DBH^1.8298*H^0.9546)+(0.012*DBH^1.6378*H^0.7746)+(0.0338*DBH^2. 6624*H^-0.5743)+(0.1699*DBH^2.3289*H^-1.1316) 43 CWD Carbon stock adjustments by the sampling effort To account for differences in the sampling effort between the survey years (2001, 2020), the carbon stocks in each plot in 2020 were adjusted by the proportion of the number of sampled transects per plot (Table 2), which is the total number of transects surveyed in 2020 divided by the total surveyed in 2001 (Equation 9). Equation 9. The proportion of sampled transects per PSP. ܲܵܲᇱ ‫= ݈݀݁݌݉ܽݏ ݊݋݅ݐݎ݋݌݋ݎܲ ݏ‬ ܰ‫݋‬. ܶ‫ ݏݐܿ݁ݏ݊ܽݎ‬2020 ܰ‫݋‬. ܶ‫ ݏݐܿ݁ݏ݊ܽݎ‬2001 Table 2. The proportion of sampled transects per plot per treatment and site. The missing plots imply they were not found and surveyed in 2020. Site Treatment Plot Proportion sampled 2020 CC East Twin GS UN 1 0.67 2 0.67 3 0.67 4 0.67 5 0.33 6 1.00 7 0.67 1 0.33 2 0.33 3 0.67 4 0.67 5 1.00 6 0.33 1 0.67 2 1.00 3 1.00 4 1.00 44 CC GR Lunate GS UN CC Minnow GR GS 5 1.00 6 1.00 2 0.67 7 0.33 1 1.00 5 0.67 6 0.33 7 0.67 8 0.67 2 1.00 4 1.00 8 1.00 3 1.00 4 1.00 5 1.00 6 0.67 7 0.67 1 0.67 2 0.33 3 0.67 4 0.67 5 0.67 2 0.33 3 0.33 4 0.33 1 1.00 2 0.67 3 1.00 4 1.00 5 0.67 45 UN 6 1.00 7 0.50 8 0.50 1 0.67 2 0.67 3 1.00 4 0.67 5 0.67 6 0.67 46 Appendix B: Areas per site and Treatments Summary of the forested, harvested portions and total areas per site and per treatment. These areas were estimated using ArcMap 10.5.1 (ESRI, 2011) from the base maps provided by Jull et al., (2001) (Figure 2). These areas were used as reference for estimating the proportion of harvested and forested areas at each site/treatment, which were used to adjust the carbon stocks at specific pools. Table 3. Areas (ha) for all the treatments at each study site. The proportion of harvested/forested areas was estimated based on these areas. Site All sites East Twin Lunate Minnow Treatment CC GR GS UN CC GS UN CC GR GS UN CC GR GS UN Total Harvested Area (ha) area (ha) 30.8 28 41.7 42 7.5 8.7 11 16.3 19 22 21 7 9 11 10 0 4.8 33.3 42 7.5 2.1 0 16.3 15.5 3.9 0 7 7.7 2.5 0 Retained area (ha) 30.8 23.2 8.4 0 0 6.6 11 0 3.5 18.2 21 0 1.3 8.6 10 Proportion Proportion harvested forested 1 0.8 0.2 0 1 0.2 0 1 0.8 0.2 0 1 0.9 0.2 0 0 0.2 0.8 1 0 0.8 1 0 0.2 0.8 1 0 0.1 0.8 1 47 Appendix C Carbon stocks 19 years after harvest Table 4.Mean and standard error (±95% confidence interval) for LT, LR, DR, DST, CWD, FWD, FF and total forest C stocks in four forest harvesting treatment units at three ICH study areas, Units are MgC ha-1. Data 19 years after harvest (2020). Site All Sites Mean Carbon pool (MgC ha-1) Unharvested Group Selection Group Retention Clearcut Live trees -LT 280.68 ± 31.64 [x] 281.85 ± 63.06 [x] 54.77 ± 7.31 [y] 14.50 [y] c Live roots -LR 62.32 ± 7.04 [x] 62.58 ± 14 [x] 12.16 ± 1.62 [x] 3.22 [y] c Dead roots -DR 2.21 ± 0.42 [x] 3.45 ± 0.43 [x] 5.92 ± 0.27 [y] 9.72 ± 0.6 [z] Dead Standing trees -DST 15.14 ± 6.86 [x] 11.87 ± 5.05 [x] 2.16 ± 0.7 [x] b Coarse Woody Debris -CWD 66.17 ± 23.74 [x] 41.64 ± 14.91 [x] 81.89 ± 30 [x] 54.26 ± 21.08 [x] Fine woody Debris -FWD 1.62 ± 0.6 [x] 1.28 ± 0.38 [x] 1.04 ± 0.93 [x] 1.26 ± 0.46 [x] Forest floor -FF 47.14 ± 4.71 [x] 46.82 ± 4.94 [x] 34.91 ± 3.92 [x] 52.34 ± 13.34 [x] Total Carbon 475.29 ± 19.74 449.49 ± 32.22 192.86 ± 14.63 136.29 ± 11.94 Live trees -LT 318.15 ± 66.45 393.74 ± 142.87 66.35 ± 11.28 4.53 c Live roots -LR 70.66 ± 14.78 87.41 ± 31.72 14.73 ± 2.5 1c Dead roots -DR 2.23 ± 0.95 3.88 ± 0.69 6.46 ± 0.31 10.45 ± 0.98 Dead Standing trees -DST 9.97 ± 8.36 12.09 ± 9.26 2.86 ± 0.79 b Coarse Woody Debris -CWD 54.74 ± 26.27 28.76 ± 8.5 67.85 ± 44.02 62.94 ± 133.95 Fine woody Debris -FWD 2.4 ± 2.38 1.41 ± 0.53 1.09 ± 1.63 0.86 ± 2.18 Forest floor -FF 56.46 ± 7.81 44.16 ± 5.89 36.38 ± 6.1 61.66 ± 34.2 Total Carbon 514.61 ± 31.13 571.45 ± 62.17 195.71 ± 19.43 141.44 ± 58.35 Live trees -LT 284.39 ± 37.22 243.83 ± 42.83 a 18.39 c Live roots -LR 63.13 ± 8.26 54.13 ± 9.51 a 4.08 c Dead roots -DR 1.41 ± 0.62 3.77 ± 0.97 a 9.08 ± 0.76 East Dead Standing trees -DST 11.02 ± 6.84 7.07 ± 6.66 a b Twin Coarse Woody Debris -CWD 40.77 ± 29.39 37.08 ± 41.24 a 50.97 ± 27.45 Fine woody Debris -FWD 1.3 ± 0.43 1.04 ± 0.82 a 1.36 ± 0.64 Forest floor -FF 37.61 ± 4.64 45.26 ± 6.93 a 34.45 ± 12.23 Total Carbon 439.64 ± 20.08 392.18 ± 24.88 a 118.33 ± 12.31 Live trees -LT 245.98 ± 35.91 196.52 ± 33.34 51.54 ± 8.91 17.05 c Live roots -LR 54.61 ± 7.97 43.65 ± 7.37 11.45 ± 1.98 3.78 c Minnow Dead roots -DR 2.99 ± 0.56 2.62 ± 0.44 5.23 ± 0.36 9.44 ± 1.56 24.43 ± 19.63 16.02 ± 12.57 1.39 ± 1.06 b 49.9 ± 22.4 105.3 ± 37.81 55.37 ± 50.15 Lunate Dead Standing trees -DST Coarse Woody Debris -CWD 101.11 ± 60.59 48 Fine woody Debris -FWD 1.29 ± 0.68 1.41 ± 0.72 0.95 ± 2.36 1.29 ± 1.38 Forest floor -FF 50.55 ± 12.66 53.92 ± 18.43 33.44 ± 6.26 62.45 ± 19.26 Total Carbon 480.96 ± 29.93 364.03 ± 19.68 209.3 ± 16.68 149.4 ± 19.62 Note: Means across “All sites” sharing the same letter (x,y,z) are not statistically significant using 95% confidence intervals generated from a linear mixed-effect (see Appendix F) a. GR treatment is not present at East Twin b. Assumed as negligible (see Appendix D) c. No error reported (see Appendix D) 49 Carbon stocks in live trees 19 years after harvest Table 5. Mean carbon (MgC ha-1) from live trees (LT) in 2020 disaggregated by the contribution from the forested and harvested areas at each Treatment/Site. Site Treatment GR GS All sites CC UN GS East Twin CC UN GR GS Lunate CC UN GR GS Minnow CC UN Strata Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration Retained Regeneration C mean (MgC ha-1) 39.78 14.99 279.91 1.94 0.00 14.50 280.68 0 242.39 1.44 0 18.39 284.39 0.00 56.96 9.39 392.17 1.57 0 4.53 318.15 0.00 30.25 21.29 193.77 2.75 0.00 17.05 245.98 0.00 Total C mean (MgC ha-1) 54.77 281.85 14.50 280.68 243.83 18.39 284.39 66.34 393.74 4.53 318.15 51.55 196.51 17.05 245.98 50 Carbon stocks pre- and post-harvesting Table 6. Mean and standard error (±95% confidence interval) for five pools (LT, LR, DST, DR, CWD), and total forest C stocks in four forest harvesting treatment units at three ICH study areas. Units are MgC ha-1. Data available preharvest (1999), 1 year (2001) and 19 years after harvest (2020). Site Mean Carbon pool (MgC ha-1) Live trees LT Live roots LR Dead roots DR All Sites Dead Standing trees -DST Unharvested Group Selection Group Retention Clearcut 1999 266.24 ± 41.07 315.03 ± 51.92 225.39 ± 22.33 308.27 ± 40.57 2001 262.04 ± 32.13 253.36 ± 55.24 34.44 ± 6.88 0±0 2020 280.68 ± 31.64 281.85 ± 63.06 54.77 ± 7.31 14.50 c 1999 59.1 ± 9.12 69.94 ± 11.53 50.04 ± 4.96 68.26 ± 9.04 2001 58.19 ± 7.14 56.46 ± 12.33 7.65 ± 1.53 0±0 2020 62.32 ± 7.04 62.58 ± 14 12.16 ± 1.62 3.22 c 1999 7.33 ± 1.23 3.71 ± 0.88 3.48 ± 0.8 2.36 ± 0.52 2001 2.12 ± 1.02 13.02 ± 1.14 37.05 ± 1.69 63.62 ± 3.92 2020 2.21 ± 0.42 3.45 ± 0.43 5.92 ± 0.27 9.72 ± 0.6 1999 30.15 ± 10.96 15.56 ± 7.68 14.37 ± 7.32 9.61 ± 4.64 2001 24.24 ± 12.15 12.49 ± 11.06 2.41 ± 1.41 0±0 2020 15.14 ± 6.86 11.87 ± 5.05 2.16 ± 0.7 b d 86.53 ± 26.01 86.53 ± 26.01 86.53 ± 26.01 86.53 ± 26.01 2001 86.53 ± 26.01 89.76 ± 31.42 98.38 ± 33.47 145.33 ± 61.84 2020 66.17 ± 23.74 41.64 ± 14.91 81.89 ± 30 54.26 ± 21.08 1999 449.35 ± 24.14 490.77 ± 28.67 379.81 ± 16.66 475.04 ± 23.51 2001 433.12 ± 21.02 424.87 ± 31.72 179.92 ± 16.15 208.95 ± 29.8 2020 426.52 ± 19.61 401.39 ± 32.13 156.91 ± 14.51 81.7 ± 10.09 1999 305.28 ± 94.5 424.15 ± 79.42 242.84 ± 28.31 332.73 ± 75.38 2001 319.94 ± 75.81 349.09 ± 160.7 43.94 ± 9.91 0±0 2020 318.15 ± 66.45 393.74 ± 142.87 66.35 ± 11.28 4.53 c 1999 67.77 ± 20.98 94.16 ± 17.63 53.91 ± 6.28 73.87 ± 16.74 2001 71.07 ± 16.87 77.5 ± 35.68 9.75 ± 2.2 0±0 2020 70.66 ± 14.78 87.41 ± 31.72 14.73 ± 2.5 1c 1999 9.62 ± 2.42 5.72 ± 1.73 3.49 ± 1.10 1.86 ± 0.5 2001 0.73 ± 0.32 16.86 ± 1.54 39.97 ± 1.95 68.41 ± 6.4 2020 2.23 ± 0.95 3.88 ± 0.69 6.46 ± 0.31 10.45 ± 0.98 Dead Standing trees -DST 1999 44.8 ± 27.54 26.65 ± 19.1 16.09 ± 12.12 7.89 ± 6.06 2001 8.43 ± 6.31 8.08 ± 9.33 1.86 ± 1.25 0±0 2020 9.97 ± 8.36 12.09 ± 9.26 2.86 ± 0.79 b Coarse Woody Debris -CWD 1999 d 100 ± 65.36 100 ± 65.36 100 ± 65.36 100 ± 65.36 2001 100 ± 65.36 110.51 ± 81.15 77.07 ± 37.55 215.27 ± 157.21 2020 54.74 ± 26.27 28.76 ± 8.5 67.85 ± 44.02 62.94 ± 133.95 1999 527.47 ± 49.1 650.68 ± 44.33 416.32 ± 29.86 516.35 ± 42.27 2001 500.16 ± 42.42 562.03 ± 77.73 172.59 ± 16.58 283.68 ± 66.79 Coarse Woody Debris -CWD Total Carbon Live trees LT Live roots LR Lunate Year Dead roots DR Total Carbon 1999 51 2020 455.75 ± 30.94 525.88 ± 62.12 158.24 ± 19.25 78.92 ± 56.65 1999 266.15 ± 46.37 311.16 ± 78.63 a 280.68 ± 59.59 2001 260.91 ± 38.29 245.54 ± 41.24 a 0±0 2020 284.39 ± 37.22 243.83 ± 42.83 a 18.39 c 1999 59.09 ± 10.29 69.08 ± 17.46 a 62.31 ± 13.23 2001 57.92 ± 8.5 54.51 ± 9.15 a 0±0 2020 63.13 ± 8.26 54.13 ± 9.51 a 4.08 c 1999 4.69 ± 0.95 2.34 ± 1.03 a 3.13 ± 1.03 2001 1.31 ± 0.45 11.23 ± 1.08 a 59.45 ± 4.98 2020 1.41 ± 0.62 3.77 ± 0.97 a 9.08 ± 0.76 Dead Standing trees -DST 1999 21.85 ± 11.35 10.88 ± 12.38 a 15.43 ± 11.72 2001 18.54 ± 18.73 2.07 ± 1.46 a 0±0 2020 11.02 ± 6.84 7.07 ± 6.66 a b Coarse Woody Debris -CWD 1999 d 60.11 ± 41.02 60.11 ± 41.02 a 60.11 ± 41.02 2001 60.11 ± 41.02 70.6 ± 74.8 a 103.72 ± 41.47 2020 40.77 ± 29.39 37.08 ± 41.24 a 50.97 ± 27.45 1999 411.89 ± 24.83 453.56 ± 35.51 a 421.85 ± 30.01 2001 398.79 ± 24.35 383.96 ± 34.17 a 163.17 ± 17.66 2020 400.72 ± 19.98 345.89 ± 24.71 a 82.52 ± 11.24 1999 220.77 ± 88.04 208.83 ± 33.39 202.13 ± 35.56 307.76 ± 140.01 2001 205.27 ± 28.77 159.49 ± 31.94 23.36 ± 8.5 0±0 2020 245.98 ± 35.91 196.52 ± 33.34 51.54 ± 8.91 17.05 c 1999 49.01 ± 19.54 46.36 ± 7.41 44.87 ± 7.89 67.62 ± 31.41 2001 45.57 ± 6.39 35.42 ± 7.07 5.19 ± 1.89 0±0 2020 54.61 ± 7.97 43.65 ± 7.37 11.45 ± 1.98 3.78 c 1999 4.99 ± 1.23 1.71 ± 0.56 2.67 ± 0.89 0.90 ± 0.55 2001 4.33 ± 2.99 10.88 ± 1.87 33.16 ± 2.26 61.82 ± 10.21 2020 2.99 ± 0.56 2.62 ± 0.44 5.23 ± 0.36 9.44 ± 1.56 1999 21.35 ± 14.96 7.97 ± 6.12 12.08 ± 11.13 4.21 ± 7.13 2001 45.73 ± 30.29 26.72 ± 33.86 2.55 ± 2.39 0±0 2020 24.43 ± 19.63 16.02 ± 12.57 1.39 ± 1.06 b d 97.24 ± 43.72 97.24 ± 43.72 97.24 ± 43.72 97.24 ± 43.72 2001 97.24 ± 43.72 83.37 ± 28.53 119.68 ± 61.44 91.67 ± 49.77 2020 101.11 ± 60.59 49.9 ± 22.4 105.3 ± 37.81 55.37 ± 50.15 1999 393.37 ± 39.44 362.12 ± 22.48 358.99 ± 22.57 477.74 ± 54.47 2001 398.15 ± 24.85 315.88 ± 23.36 183.94 ± 26.36 153.49 ± 20.63 Live trees LT Live roots LR Dead roots DR East Twin Total Carbon Live trees LT Live roots LR Dead roots DR Minnow Dead Standing trees -DST Coarse Woody Debris -CWD Total Carbon 1999 2020 429.12 ± 29.44 308.7 ± 18.07 174.9 ± 16.46 85.65 ± 18.13 a. GR treatment is not present at East Twin b. Assumed as negligible (see Appendix D) c. No error reported (see Appendix D) d. Values of unharvested plots area assumed to be the same in 1999 and 2001 because CWD was not sampled in 1999. Also, I extrapolated the 1999 unharvested CWD as proxy values for the GS, GR, and CC. Although the stocks may had been different if sampled, this proxy of comparing treatments with historical control data (Paulus 52 et al., 2014) allowed me to use the available stocks to gain deeper insights of the carbon dynamics through time, rather than removing the entire pool due to the absence of pre-harvest data. Appendix D Modelling Carbon stocks in harvested areas I used SORTIE 7.05 (Canham & Murphy, 2001) to estimate projected carbon stocks in the harvested areas (i.e. planted regeneration plots) for the GS, GR and CC treatments at all sites and integrate these values into the total carbon stored in 2020 in the Wetbelt treatments. To calibrate the initial model parameters, I used the site characteristics, species composition and stem densities from the original regeneration survey in 2016. In addition, reference curves with different planting density stands were created (i.e. 1200 and 1800 stems per hectare). Each site-treatment was run independently. These were 5-year runs accounting for 0-4 timesteps. The output is the mean modelled carbon (MgC ha-1) from stems>4 cm DBH for 2020 (Table 7). To avoid overestimation in the carbon stocks, these values were adjusted multiplying by: 1) a “proportional increase” (derived from the proportion of modelled C measured in the empirical data); 2) proportion of the harvested area at each site-treatment (Table 3; see methods). Table 7. Projected carbon stocks (MgC ha-1) in planted trees in regeneration plots for 2020. Data derived from growth and yield curves using SORTIE 7.05 and regeneration survey from 2016. Carbon stocks have already been adjusted by the proportional increase and by the proportion of the harvested area. Site/TU East Twin Lunate Minnow Year GS CC GS GR CC GS GR CC 2020 1.44 18.39 1.57 9.39 4.53 2.75 21.29 17.05 Given the parameters to run the model, only the mean was available for this carbon stock. In the partial harvest treatments, these values were added to the 2020 empirical values for the live tree pool from the forested areas, whose error was estimated and accounted for in the total carbon stocks. However, it is not the case in the CC, whereby it was not possible to report an associated error to this pool (Table 6). 53 Since the carbon in the live root estimation relies on the live trees, same error situation applies to the live root pool for 2020 (Table 6). To estimate the carbon in the live tree and live root pools in 2020 for all sites and compare only by treatment (TU), I used the combined mean formula (Equation 10) and the final estimations from the modelled data by site/TU (Higgins et al., 2023). No associated error is reported due to the nature of the data source. Equation 10. Formulae for combining groups (Higgins et al., 2023). ‫= ݊ܽ݁݉ ܾ݀݁݊݅݉݋ܥ‬ ܰଵ ‫ܯ‬ଵ + ܰଶ ‫ܯ‬ଶ + ܰଷ ‫ܯ‬ଷ ܰଵ + ܰଶ + ܰଷ ܹℎ݁‫݁ݎ‬, N is the sample size and M is the carbon mean (MgC ha-1) from each site (East Twin, Lunate and Minnow). Since the number of dead standing trees (DST) from the regeneration survey in 2016 was insignificant (23 of 1324 trees), I did not model carbon projections for this pool in 2020 and considered it as negligible (Table 6). 54 Appendix E Canopy conditions’ effect on live tree mean carbon accumulation rate Table 8. a) Mean, standard error Live tree carbon accumulation rate (MgC year-1) given three canopy conditions within each treatment; b) Live tree mean carbon accumulation rate per treatment. Treatment a) Canopy Mean Rate Standard Standard conditions (MgC year-1) deviation error n GR Closed 0.0005 0.0008 0.0000 489 GR Edge 0.0008 0.0011 0.0001 303 GS Closed 0.0017 0.0021 0.0001 574 GS Edge 0.0023 0.0022 0.0001 332 UN Closed 0.0017 0.0020 0.0001 1348 UN Edge 0.0028 0.0022 0.0009 6 b) Mean Rate Standard Standard (MgC year-1) deviation error GR 0.0007 0.0009 0.00003 802 GS 0.0020 0.0021 0.00007 930 UN 0.0017 0.0020 0.00005 1355 Treatment n To explore whether the partial harvest treatments had a significant influence on the carbon accumulation in the live trees, I performed a Kruskal Wallis and Dunn’s tests (post-hoc) between 1) Live trees mean carbon accumulation rate ~ treatment, and 2) Live trees mean carbon accumulation rate ~ canopy conditions (O, E, C). Mean carbon accumulation rate ~ Treatment • Kruskal-Wallis test Kruskal-Wallis chi-squared = 182.19, df = 2, p-value < 2.2e-16 • Dunn’s tests (post-hoc) # A tibble: 3 × 9 .y. group1 group2 n1 n2 statistic * p p.adj p.adj.signif 55 1 Rate GR 2 Rate GR GS UN 792 906 12.9 5.32e-38 1.60e-37 **** 792 1354 10.7 6.15e-27 1.84e-26 **** 3 Rate GS UN 906 1354 -3.41 6.61e- 4 1.98e- 3 ** Mean carbon accumulation rate ~ Canopy conditions within each treatment • Kruskal-Wallis test Kruskal-Wallis chi-squared = 8.88, df = 1, p-value = 0.002883 • Dunn’s tests (post-hoc) $UN # A tibble: 1 × 9 .y. group1 group2 n1 n2 statistic * 1 Rate C E 1348 6 p p.adj p.adj.signif 1.49 0.137 0.137 ns $GS # A tibble: 1 × 9 .y. group1 group2 n1 n2 statistic * 1 Rate C E 574 332 p p.adj p.adj.signif 4.43 0.00000952 0.00000952 **** $GR # A tibble: 1 × 9 .y. group1 group2 n1 n2 statistic * 1 Rate C E 489 303 p p.adj p.adj.signif 3.85 0.000117 0.000117 *** 56 Appendix F Linear mixed-effect model for carbon stocks in 2020 Because of the differential contribution of each pool to the total ecosystem carbon and the fact that each pool had its own residual variance, it made more sense to use a mixed-effect model across all sites to test for significant differences among the pools and the treatments 19 years after harvest. I fitted a linear mixed-effect model (estimated using REML and nlminb optimizer) to predict mean of Carbon with Treatment and Pool: lme.allsites_lme3 <- lme(Carbon_mean ~ Treatment*Pool, random = ~1 | Site/Treatment, data=Carbon_dataset, weights = varIdent(form = ~ 1 | Pool) ) The model included Site and treatment as random effects (formula: ~1 | Site/Treatment). Because of the large variation in carbon among the pools, this weighted model is the best fit since the “weights” component gives a unique residual variance to each pool. The model's explanatory power related to the fixed effects alone (marginal R 2) is 1.00. The model's intercept corresponding to the treatment and pool is at 58.55 (t(42) = 14.58, p = 0.000). Within this model: Linear mixed-effects model fit by REML Data: Carbon_dataset AIC BIC logLik 477.2214 547.2187 -201.6107 Random effects: Formula: ~1 | Site (Intercept) StdDev: 1.004952e-05 Formula: ~1 | Treatment %in% Site (Intercept) Residual 57 StdDev: 7.558637e-05 20.68942 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Pool Parameter estimates: CWD DR DT FF FWD LR LT 1.00000000 0.04540317 0.27517868 3.01996093 0.01187068 0.81155678 2.97016120 Fixed effects: Carbon_mean ~ Treatment * Pool Value Std.Error DF t-value p-value (Intercept) 58.54727 4.01451 42 14.583911 0.0000 Treatment1 -33.11375 6.69085 5 -4.949109 0.0043 Treatment2 -14.49872 7.68720 5 -1.886085 0.1180 Treatment3 30.51674 6.69085 5 4.560965 0.0061 Pool1 3.23413 6.69182 42 0.483296 0.6314 Pool2 -53.50383 4.02186 42 -13.303244 0.0000 Pool3 -47.77522 4.27631 42 -11.172061 0.0000 Pool4 71.99154 16.65949 42 4.321354 0.0001 Pool5 -57.33969 4.01501 42 -14.281319 0.0000 Pool6 -10.48721 5.91568 42 -1.772782 0.0835 Treatment1:Pool1 27.76342 11.15304 42 2.489315 0.0168 Treatment2:Pool1 39.29271 12.81386 42 3.066422 0.0038 Treatment3:Pool1 -53.71664 11.15304 42 -4.816324 0.0000 Treatment1:Pool2 37.41608 6.70311 42 5.581902 0.0000 Treatment2:Pool2 14.70097 7.70128 42 1.908899 0.0631 Treatment3:Pool2 -32.18893 6.70311 42 -4.802091 0.0000 Treatment1:Pool3 22.34170 7.12719 42 3.134715 0.0031 Treatment2:Pool3 16.35548 8.18851 42 1.997369 0.0523 Treatment3:Pool3 -25.96924 7.12719 42 -3.643688 0.0007 Treatment1:Pool4 11.09056 27.76581 42 0.399432 0.6916 Treatment2:Pool4 -11.30928 31.90049 42 -0.354517 0.7247 58 Treatment3:Pool4 45.70662 27.76581 42 1.646148 0.1072 Treatment1:Pool5 32.95834 6.69169 42 4.925265 0.0000 Treatment2:Pool5 14.31379 7.68817 42 1.861795 0.0696 Treatment3:Pool5 -30.43633 6.69169 42 -4.548377 0.0000 Treatment1:Pool6 -14.45965 9.85946 42 -1.466576 0.1499 Treatment2:Pool6 18.07140 11.32766 42 1.595334 0.1181 Treatment3:Pool6 -0.77789 9.85946 42 -0.078898 0.9375 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.59291411 -0.63394429 -0.02183757 0.49339911 2.03294586 Number of Observations: 77 Number of Groups: Site Treatment %in% Site 3 11 anova(lme.allsites_lme3) numDF denDF F-value p-value (Intercept) 1 42 442.0063 <.0001 Treatment 3 5 0.5443 0.673 Pool 6 Treatment:Pool 42 81.7568 <.0001 18 42 12.0422 <.0001 $emmeans (Post-hoc) Pool = CWD: Treatment emmean SE df lower.CL upper.CL CC 56.431 11.945 2 5.036 107.83 GR 86.575 14.630 2 23.629 149.52 GS 38.581 11.945 2 -12.814 89.98 UN 65.538 11.945 2 14.142 116.93 59 Pool = DR: Treatment emmean SE df lower.CL upper.CL CC 9.346 0.542 2 7.012 11.68 GR 5.246 0.664 2 2.388 8.10 GS 3.371 0.542 2 1.038 5.70 UN 2.211 0.542 2 -0.122 4.54 Pool = DT: Treatment emmean SE df lower.CL upper.CL CC 0.000 3.287 2 -14.143 14.14 GR 12.629 4.026 2 -4.693 29.95 GS 15.320 3.287 2 1.177 29.46 UN 15.140 3.287 2 0.997 29.28 Pool = FF: Treatment emmean SE df lower.CL upper.CL CC 108.516 36.074 2 -46.696 263.73 GR 104.731 44.181 2 -85.364 294.83 GS 206.762 36.074 2 51.550 361.97 UN 102.147 36.074 2 -53.065 257.36 Pool = FWD: Treatment emmean SE df lower.CL upper.CL CC 1.052 0.142 2 0.442 1.66 GR 1.023 0.174 2 0.275 1.77 GS 1.288 0.142 2 0.678 1.90 UN 1.467 0.142 2 0.857 2.08 Pool = LR: Treatment emmean SE df lower.CL upper.CL 60 CC 0.487 9.694 2 -41.224 42.20 GR 51.633 11.873 2 GS 77.799 9.694 2 36.089 119.51 UN 62.322 9.694 2 20.612 104.03 0.548 102.72 Pool = LT: Treatment emmean SE df lower.CL upper.CL CC 2.203 35.479 2 -150.449 154.86 GR 46.504 43.452 2 -140.457 233.46 GS 280.327 35.479 2 127.674 432.98 UN 280.676 35.479 2 128.024 433.33 Degrees-of-freedom method: containment Confidence level used: 0.95 $contrasts Pool = CWD: contrast estimate SE df t.ratio p.value CC - GR -30.1443 18.887 5 -1.596 0.4563 CC - GS 17.8496 16.893 5 1.057 0.7276 CC - UN -9.1066 16.893 5 -0.539 0.9456 GR - GS 47.9939 18.887 5 2.541 0.1667 GR - UN 21.0378 18.887 5 1.114 0.6978 GS - UN -26.9561 16.893 5 -1.596 0.4565 Pool = DR: contrast estimate SE df t.ratio p.value CC - GR 4.1001 0.858 5 4.781 0.0183 CC - GS 5.9745 0.767 5 7.790 0.0021 CC - UN 7.1347 0.767 5 9.302 0.0009 GR - GS 1.8744 0.858 5 2.186 0.2458 61 GR - UN 3.0347 0.858 5 3.539 0.0581 GS - UN 1.1602 0.767 5 1.513 0.4947 Pool = DT: contrast estimate SE df t.ratio p.value CC - GR -12.6288 5.197 5 -2.430 0.1883 CC - GS -15.3195 4.649 5 -3.296 0.0745 CC - UN -15.1398 4.649 5 -3.257 0.0776 GR - GS -2.6907 5.197 5 -0.518 0.9512 GR - UN -2.5110 5.197 5 -0.483 0.9596 GS - UN 0.1797 4.649 5 0.039 1.0000 Pool = FF: contrast estimate CC - GR SE df t.ratio p.value 3.7848 57.037 5 0.066 0.9999 CC - GS -98.2465 51.016 5 -1.926 0.3252 CC - UN 6.3690 51.016 5 0.125 0.9992 GR - GS -102.0314 57.037 5 -1.789 0.3754 GR - UN 2.5842 57.037 5 0.045 1.0000 GS - UN 104.6155 51.016 5 2.051 0.2845 Pool = FWD: contrast estimate SE df t.ratio p.value CC - GR 0.0295 0.224 5 0.132 0.9991 CC - GS -0.2358 0.201 5 -1.176 0.6652 CC - UN -0.4153 0.201 5 -2.071 0.2783 GR - GS -0.2653 0.224 5 -1.183 0.6612 GR - UN -0.4448 0.224 5 -1.984 0.3056 GS - UN -0.1795 0.201 5 -0.895 0.8084 62 Pool = LR: contrast estimate SE df t.ratio p.value CC - GR -51.1461 15.328 5 -3.337 0.0714 CC - GS -77.3123 13.710 5 -5.639 0.0091 CC - UN -61.8353 13.710 5 -4.510 0.0231 GR - GS -26.1662 15.328 5 -1.707 0.4083 GR - UN -10.6892 15.328 5 -0.697 0.8940 GS - UN 15.4770 13.710 5 1.129 0.6899 Pool = LT: contrast estimate SE df t.ratio p.value CC - GR -44.3004 56.097 5 -0.790 0.8565 CC - GS -278.1234 50.174 5 -5.543 0.0098 CC - UN -278.4731 50.174 5 -5.550 0.0098 GR - GS -233.8230 56.097 5 -4.168 0.0316 GR - UN -234.1727 56.097 5 -4.174 0.0314 GS - UN -0.3497 50.174 5 -0.007 1.0000 Degrees-of-freedom method: containment P value adjustment: tukey method for comparing a family of 4 estimates 63