Managing forest carbon and landscape capacities

Widespread impacts of a warming planet are fuelling climate change mitigation efforts world-wide. Decision makers are turning to forests, the largest terrestrial primary producer, as a nature-based contribution to mitigation efforts. Resource-based economies, however, have yet to include carbon (C) in their resource planning, slowing the implementation of these important measures for atmospheric greenhouse gas reduction. The realisation of forest mitigation potential depends greatly on our ability to integrate C-sequestration practices in our forest management applications. This requires robust C-estimates, an understanding of the natural potential for a specific landscape to sequester C, the current state of the landscape relative to this potential, and the evaluation of management practices as a tool to sequester forest C in the midst of all the other values forests offer humans. Discrepancies between models used in management decisions and C estimation are the first hurdle impeding the application of forest-based mitigation strategies. Here, we combine forest disturbance and management models with a well-established C model on an open-source simulation platform. We then use the modelling system to produce C estimates of the natural C-holding capacity (potential) and two management scenarios for a study area in BC, Canada. Our simulations provide an essential metric if forests are to be managed for C-sequestration: the natural landscape C-holding capacity. Our simulations also point to a decreasing trend in simulated C on the study area over time and to a bias of the current C-levels compared to the landscape C-holding capacity (477 vs 405.5 MtC). Our explanations for this bias may provide an avenue for improved current C-state estimates. We provide a framework and the information needed for the implementation of nature-based solutions using forests for climate change mitigation. This study is a step towards modelling systems that can unify scientifically based forest management and informed C-management.


Introduction
Forests are the largest terrestrial primary producer, absorbing the largest proportion of atmospheric CO 2 (Friedlingstein et al 2020). They have considerable potential to help mitigate the impacts of changing climatic conditions (Griscom et al 2017, Cook-Patton et al 2020. To realise this potential however, C needs to be an integral part of forest management decisions. Impediments to including C in forest management decisions vary per jurisdiction, but have a common trait: discrepancies between models used in management decisions and C estimation. The sheer complexity of the forest C cycle makes it challenging for measuring, monitoring, and projecting changes to forest C (Boisvenue and White 2019, Forests in focus 2021). Further, the geographic extent and spatiotemporal variability in the system prohibits, in most cases, a complete census of forest C (FAO 2019), which implies that sampling and modelling are therefore essential to forest management and planning. Projecting productivity via modelling is the basis for sustainable forest management and modelling of species dynamics supports silvicultural interventions. Modelling efforts have also resulted in forest C and its variations being estimated with a diversity of models (Xia et al 2017, Walker et al 2022. Despite these efforts, gaps remain in our theoretical and empirical understanding of C allocation (Merganičová et al 2019).
The call for integration of forest-based solutions in emission reduction policy is increasing (Anderegg et al 2020, Matthews et al 2020, Lamb et al 2021 with policies being put in place for integration of forest-based intervention in the mitigation efforts (Environment and Climate Change Canada 2020). Their application, where management actions aimed at increased C storage are implemented, are still restricted to academic publications (Liu et al 2016, Pilli et al 2021, Sleeter et al 2022, Zhao et al 2022 5 . Even in jurisdictions where forest C is quantified, most do not include forest C when planning greenhouse gas reductions (Lamb et al 2021), nor is C included in forest management decisions 6 . This is partly because our ability to manipulate the C content of forests rests within the realm of forest management, and models supporting those decisions differ in structure and spatial extent from forest C models (e.g. Fortin andLavoie 2022 vs Walker et al 2022). Nonetheless, forest-based interventions need C information and decision makers face these modelling discrepancies amongst other challenges (Grassi et al 2017, Markusson 2022. Sequestration potentials have been estimated using C-models (Griscom et al 2017, Smyth et al 2020, Drever et al 2021. In Canada for example, forest and specifically improved forest management is estimated to offer the 4th largest mitigation opportunities for 2030 among natural climate solutions (Drever et al 2021). Even with a well-established forest C model, as is the case for Canada (Environment and Climate Change Canada 2020), forest C management is, as of yet, not incorporated into resource planning. Effective forest-based solutions require (a) robust estimates of current C, (b) an assessment of the land's potential and unrealised potential for C-sequestration, and (c) integration of C estimates and potential into land-management decision schema (Walker et al 2022).
Landscape characteristics partly determine the amount of C that can be maintained, i.e. its potential. This amount varies through time, and is bounded by environmental and site conditions (Boisvenue and Running 2006, Running 2012, Walker et al 2022. 5 No C in resource values: www2.gov.bc.ca/gov/content/industry/ forestry/managing-our-forest-resources/forest-stewardship-plans. 6 Example: www2.gov.bc.ca/gov/content/industry/forestry/ managing-our-forest-resources/timber-supply-review-andallowable-annual-cut.
The potential or C-holding capacity is the mass of C stored in a landscape under prevailing historical environmental conditions and natural disturbance regimes (Keith et al 2010, Liang et al 2017 excluding anthropogenic interventions. A disturbance regime is a historic pattern (e.g. frequency and extent) of natural processes such as fire, insects, wind, and mass movement that affect the ecosystems and landscapes in a particular area. The C-holding capacity is a theoretical value that reflects the long-term stabilization of the landscape C under the dominant natural disturbance. If the dominant disturbance is fire, this Cpotential is linked to the historic fire return interval (FRI). Current potential to hold C is a legacy of past environmental conditions. It will change with changing environmental conditions as those drive the disturbance regimes. Most forest-based economies have legislated sustainable forestry practices 7 , but planning horizons are short, varying from five to ten years 8 . Given the short forest management time horizon, the notion of current capacity is necessary for informed C-management decisions. Longer-term planning will need to be informed by projections and probabilities of changing disturbance regimes under changing environmental conditions, a field of research not yet ripe for integration in shorter-term forest management decisions. Many C-centric studies compare various management approaches but omit the notion of current C-holding capacity (Pilli et al 2016, Krug 2019. Other studies assess the potential of the landscape to hold C but do not link it to the decisionmaking process involved in forest resource management (Keith et al 2010, Wei et al 2014, Domke et al 2020. Given current political and environmental motivations to maximize our climate change mitigation capacities, quantitative assessments of C in the landscape, of management interventions C-effects, and of landscape C-holding capabilities are in urgent need. Here, we develop a modelling system, combining forest disturbance and management simulations with a well-established C model for a large study area in western Canada. We then apply this system to (a) estimate C, and (b) calculate the current Cholding potential of a landscape. We then compare current potential to current conditions, and two scenarios of resource removal and discuss our findings. We follow modelling best practices of McIntire et al (2022) resulting in our approach being transparent, repeatable and reusable. We provide an avenue for implementing forest-based greenhouse gas reduction interventions, and transparent tools to evaluate and eventually improve forest C management.

Study area
We selected a study area in British Columbia, Canada, a forest-based economy 9 where a reputable approach is used for forest C estimation and reporting (Kurz et al 2009, Stinson et al 2011, Metsaranta et al 2017 10 . Our study area is located at the northeastern corner of the province of British Columbia (figure 1). It encompasses five forest management units in the provincial forest management system with 31 Mha of 53.4 Mha region under provincial forest management agreements. The topography of the units forms a gradient of increasing relief from east to west. Under provincial legislation (Forest Act 1996 11 ), these lands are actively managed for the sustainable production of resources including wood products, with a combined 19 Mm 3 yr −1 of annual allowable harvest 12 .

Context
Even in places like Canada, where the Federal Government reports annually on greenhouse gas emissions under the United Nations Framework Convention on Climate Change 13 for all managed forests, C modelling is not linked to forest management modelling. In this case, C modelling uses similar growth information as forest management models, however, incompatibilities in the software and information processing exist, and a mismatch in spatial resolution between the modelling systems results in C information being currently external to the modelling used in forest resource management. Growth modelling has emerged from the history of resource management for wood supply. Modelling of forest C was developed for the international greenhouse gas reporting, with Canada submitting its first report 2006 10 . In our example, the Chief Forester of BC determines annual allowable cut (AAC) level 14 , informed by growth and yield modelling which represents forest dynamics. In the case of our study area, these models are the basis for the C-estimates produced by the Canadian Forest Service, who is responsible for reporting on the 9 www.nrcan.gc.ca/our-natural-resources/forests/state-canadasforests-report/forest-industry-contribute/16517. 10 www.nrcan.gc.ca/climate-change-adapting-impacts-andreducing-emissions/climate-change-impacts-forests/forestcarbon/13085. 11 www.bclaws.gov.bc.ca/civix/document/id/complete/statreg/ 96157_00 12 www2.gov.bc.ca/gov/content/industry/forestry/managing-ourforest-resources/timber-supply-review-and-allowable-annualcut/allowable-annual-cut-timber-supply-areas. 13 unfccc.int/process-and-meetings/transparency-and-reporting/ greenhouse-gas-data/ghg-data-unfccc/ghg-data-from-unfccc. 14 Example: Fort Nelson Timber Supply Area news.gov.bc.ca/ releases/2019FLNR0187-001442. greenhouse gases (GHG) balance. For GHG reporting, growth models inputs and outputs are processed differently and applied at a coarser scale than that needed on forest management decisions (Kurz et al 2009). Here we start with the model currently used for reporting, the Carbon Budget Model of the Canadian Forest Sector (CBM-Kurz et al 2009), adapting it so it can provide C-estimates compatible with forest management decisions. We then link the Cmodelling with a fire model to enable the calculation of the C-holding capacity of the landscape, since fires are the dominant natural disturbances in our study area . Finally, we add a forest extraction model to assess scenarios of resource extraction and their effects on C.
Our modelling efforts strive towards an open reusable approach, as described in McIntire et al (2022). All three models were re-implemented or linked via SpaDES (Spatial Discrete Event Simulation-McIntire et al 2019), an open source platform that enables model connections and expansion to a spatially explicit format. This platform offers a solid framework to build interoperable and reusable models (e.g. Micheletti et al 2021) in the R language, a free software environment 15 . The source code used for this study is publicly available on GitHub 16 and follows FAIR data standards (Findable, Accessible, Interoperable and Reusable data; Stall et al 2019).

Forest carbon modelling
An executable for CBM is freely available 17 , but runs on Microsoft Access, which in many cases limits a user's ability to run, modify and understand the model and data processing. We therefore re-implemented CBM in R, on the SpaDES platform (Supplementary information Appendix 1). The implementation of CBM in R and SpaDES enabled modifications to the spatial resolution of simulations, making it more flexible and able to match a spatial resolution of forest management models. It also renders the model and data processing transparent, enabling us to address the incompatibilities in the software and information processing between the C model and the forest management models. CBM has parameters for all managed forests of Canada (Stinson et al 2011), which we use for our study area and make available for all of Canada in our modelling system. CBM derives live aboveground biomass estimates, which drive its growth and post-disturbance dynamics, using field-measured yield information from the forest industry (m 3 ha −1 over time) translated into biomass using statistical equations (Boudewyn et al 2007). Other C pools in CBM are estimated via process 15 www.r-project.org/. 16 github.com/cboisvenue/spadesCBM_RIA.git. 17 www.nrcan.gc.ca/climate-change/impacts-adaptations/climatechange-impacts-forests/carbon-accounting/carbon-budgetmodel/13107. representation during an initialisation phase. The initialisation phase aims at stabilizing C in dead organic matter pools. Parameters and processes for CBM are described in detail in Kurz et al (2009). CBM tracks C in 21 pools for each unit (stand or pixel), plus three pools for emissions (CO 2 , CH 4 , CO), and one for C removed from the system via harvested wood products. The flexibility and transparency enabled by R and SpaDES (McIntire et al 2019, 2022) permitted a direct link to the management and fire simulator resulting in the management and C models being driven by the same growth models, sharing input data and data processing.

Resource extraction modelling
Two management scenarios were simulated using a SpaDES adaptation of the ws3 model, a Python-based open source software package (Paradis and LeBel 2018). Ws3 reproduces industry-standard forest estate modelling functionality, used by foresters to determine AAC. We used ws3 to generate plausible patterns of harvesting and regeneration for the managed portions of each timber supply area (TSA). The spades_ws3 module (Supplementary information Appendix 2) wraps ws3 in R scripts to yield a SpaDES-compatible module that we used in our study to combine with CBM.

Fire modelling
We used a SpaDES implementation of the fire model described in Armstrong and Cumming (2003) (Supplementary Information Appendix 3). This model treats fire as a three-stage stochastic process of ignition, escape, and spread adapted to the specific landscape characteristics. To eliminate discrepancies between the fire simulation and the C-estimation, these processes had to be parameterized to simulate the FRI used in the CBM initialization phase. Disturbance-return intervals used for CBM initialization differ for each terrestrial ecozone and vary between 125 and 175 years in our study area (see Kurz et al 2009).

Data
We maintained FAIR data standards (Stall et al 2019) as much as possible for our input data and used publicly available parameters for our simulations. Our data preparation modules 18 , which we developed for each scenario, provide links to all publicly available data.
We used the default Canadian parameters for CBM (Kurz et al 2009, Stinson et al 2011 that are specific to our study area. British Columbia has an open access policy and we used their Vegetation Resource Inventory to populate our forest stand ages for each of the scenarios and across models. All forest inventory stands in CBM are represented by the growth curve of the dominant species, reflecting the predominantly stand-replacing impacts of boreal disturbance regimes. Stands can contain hardwood and softwood components, each associated with individual growth information. The growth curves used for forest management are not publicly available. We simplified the forest management growth curves (described in Supplementary Information Appendix 2) so we could provide all inputs to our system, making this modelling system and simulations transparent. We provide a link in our data preparation module to the simplified curves14. Disturbance information was specific to each scenario and is described in the following section.

Simulations scenarios 2.4.1. Landscape C-holding capacity
Disturbed stands can be expected to be emitters of atmosphere greenhouse gases either because of combustion or because of decay of organic material or both. Under a disturbance regime, a landscape can maintain a relatively stable amount of C (Keith et al 2010, Liang et al 2017. To estimate this C-holding capacity for our study area, we simulated fire at return intervals that matched the initializing phase used in CBM, with no other disturbances, starting from the British Columbia Vegetation Resource Inventory ages and dominant species in 2020, for 520 years. This accounted for the stochasticity in fire ignition, escape and spread, as well as the stabilization of the C on the landscape. The location of these fire disturbances were the product of our fire model (see section 2.3) and account for the characteristics of our study area. The representation of the C pool impacts of the disturbances in each scenario are those described in CBM and specific to our location in Canada (see Kurz et al 2009, Stinson et al 2011.

Present day
We simulated a representation of present-day conditions for our study area by simulating the same landscape, using the same modelling framework but for the time-period 1985-2015. By simulating a presentday scenario, we have a representation of current conditions through the lens of our modelling system which we use as a comparison point. Growth information for this scenario matched that used in other scenarios. Landscape conditions were derived from the BC Vegetation Resource Inventory with stand replacing disturbances locations from National Forest Information System 19 (Hermosilla et al 2016). As per all C-implications of disturbances in our scenarios, 19 opendata.nfis.org/mapserver/nfis-change_eng.html. disturbances followed their specific C transactions prescribed in CBM for this location in Canada (see Kurz et al 2009, Stinson et al 2011.

Harvesting scenarios
We simulated two resource removal scenarios one higher ('base') and one lower ('less') from the year 2020-2099. This time horizon represents multiple planning cycles in forest management 8 . The base scenario simulates harvesting the maximum sustainable even-flow volume at every time step. As a comparison point, and to corroborate the general pattern of increased resource extraction corresponding to a decrease in C storage from the scientific literature (Harmon and Marks 2002, Bradford 2011, Sharma et al 2013, Simard et al 2020, the less scenario simulates harvesting half of the volume from the base scenario. Both scenarios extract resources from the landscape via whole-stand disturbance over the horizon period as estimated by spades_ws3. Spades_ws3 schedules these events along with an expected level of wildfires, supplying the location of each stand-replacing disturbance to the modelling system. Under the less scenario, we expected that more C would be maintained on the landscape.

Results
Our simulations estimate a C-holding capacity of 405.5 MtC for our study area (table 1, figure 3 panel (a)). While all scenarios showed a decrease in total carbon from the start to the end of simulations (Table 1), they all started with landscapelevel C above the C-holding capacity (table 1). The distribution of C in the landscape differs between scenarios with the present-day and harvest scenarios having fewer pixels with low tC/ha but more with mid-range tC/ha than the C-holding capacity scenario (figure 2). The present-day scenario ended simulations with landscape-level C higher than C-holding capacity (477.2 MtC), while the two harvest scenarios ended with C lower than C-holding capacity. Harvesting less retained more carbon on the landscape (∆ 33.7 MtC) and although estimates from both harvesting scenarios kept less C on the landscape than the estimated C-holding capacity, the lesser-level of resource removal was close to C-holding capacity.
The average influx of C per ha under the C-holding capacity scenario stabilises around the 480th year of simulation at ∼0.60 MtC ha yr −1 (table 2). The less harvest scenario C uptake has the same pattern as the present-day simulations, with a decreasing level of C-uptake through the simulation horizon. The base harvest case maintains a higher net C-uptake. Current levels of resource removals as represented in our simulations, are lower than the representation of maximum sustainable removal of wood products (base). When only looking at the

Discussion
Our modelling system linked C, disturbance and management models to produce C estimates, C-holding capacity, and two levels of resource removal in our study area (figure 3). C-holding capacity simulations display a negative exponential age class distribution (figure 4(a); Wagner 1978), indicating that our goal of stabilizing landscape-C was achieved. How much C a landscape can naturally hold is an essential metric if forests are to be managed for C-sequestration (Griscom et al 2017, Drever et al 2021, Walker et al 2022. Our simulations show a decreasing trend in Cretention, with all initialized as well as the presentday (2015) landscapes having higher C than the Cholding capacity. Given the conceptual similarities between the C-holding capacity simulations and the initialization process (i.e. both repeat a grow-burn cycle), their results should be similar. Fire suppression or management are not accounted for in CBM initialization process further indicating that both should yield similar results. They differ, however, in that the CBM-initialization process burns stands at regular intervals equal to the FRI and our C-holding simulations burn at stochastic intervals with a mean return interval of FRI. This results in an upwards biased initial C estimate because of the mathematical reality of asymmetric C accumulation in growth curves on either side of the FRI age. In other words, if a stochastic fire burns 80 years younger than FRI, the stand will have not accumulated X amount of C, while the amount of C gained due to 80 years of extra growth beyond FRI age is a × X, where a < 1. This effect seems to be insensitive to the particular distribution of FRI used as it will occur as long as the growth is slowing at FRI age (i.e. negative exponential or normally distributed around a mean of FRI; see Supplementary Information Appendix 4). This lower C storage (greater loss) due to more burns younger than FRI age has recently been noticed in nearby northern Canadian forests (Baltzer et al 2021). In our C-holding capacity simulations, burns will stochastically occur during the pre-FRI age, resulting in a C accumulation that is cut short compared to the FRI age. Because the growth curves prescribe slowing of growth past FRI age, the stochastically longer-than-FRI-age inter-burn periods do not make up for the many instances of pre-FRI burns (see Supplementary Information Appendix 4). We also see a narrower age distribution (figures 4(a) and (c)), which is likely either the result of these forests being outside of their natural range due to forest management that uses a fixed rotation age or inaccurate age estimates especially of the older trees (Cumming et al 2000). This also may contribute to why initialized landscapes have higher C than the C-holding capacity. The initialization step has been necessary because we have lacked data on the highly variable soil and dead C-pools. Varying the burn interval in the CBM-initialization process has been shown to be a source of variation  (Metsaranta et al 2017) and initialization processes in other contexts have been shown to be the source of systematic errors (Oberpriller et al 2021). With the increase in availability of soils and belowground C pools data (Sothe et al 2022), initialization processes may eventually become unnecessary. The comparison of the landscape-state to its C-holding capacity will not only require more exploration of this initialization process, but will probably necessitate data on Ccontent of dead organic layers.
In our simulations, both harvesting scenarios have the same initial landscape-C values. Even if these were wrong, the comparison of resulting landscape-C would still be valid. Our results show that the higher extraction in our base-scenario results in less landscape-C than the lower extraction level (∆ 33.6 MtC) matching other published trends (Harmon and Marks 2002, Bradford 2011, Sharma et al 2013, Simard et al 2020. Some mitigation scenarios suggest storing C outside the system, such as in durable harvested wood products, and keeping the system at a higher influx (higher C-influx in the base scenario table 2, Smyth et al 2020). Assessing this scenario would require adding a wood product and life cycle analysis model to our modelling system, an opportunity presented by our modular and flexible system. C-storage outside the system may justify a higher harvest level through C-sequestration in wood products, but may have other undesired consequences such as lowering landscape resilience (Ibáñez et al 2019, Albrich et al 2020 and soil C impacts (Prescott et al 2021). Further, higher influx into the system does not necessarily equal C-residency (Körner et al 2005).
Both harvesting scenarios result in landscape-C below C-holding capacity, although the less harvest scenario is closer to the C-holding capacity (∆ <5 MtC over the 53.4 Mha). Assuming that our landscape C-holding capacity approximates the actual holding capacity and that our harvesting simulations represent the correct trends, this simulation study suggests that, if only landscape-C is considered, extracting a lower amount of resource may match the landscape-C to its C-holding capacity. However, we have not yet determined the consequences of maintaining a landscape above or below its natural Cholding capacity (for example, for an increased Cstorage elsewhere). Our current landscape (2015) has higher C than the C-holding capacity. This could suggest that the landscape is holding more C than what it could under the assumed-constant natural disturbance regime, potentially due to its historical anthropogenic interventions. It could also be stabilizing after too-high C-values from the initialization phase (see above). Given the uncertainty about the initialization process, to really determine if a landscape can hold more than its C-holding capacity and the consequences of maintaining a landscape in that state, more research is required. Our system is wellpositioned to advance this research.
Our simulations indicate that the stable-state of the forests may be better captured by improved fire modelling. A stable-state, however, may be an irrelevant concept under changing environmental conditions which would require forecasts of evolving disturbance regimes in our landscapes (Walter et al 2019). Further, most forested regions are not pristine areas devoid of human influence, but are places where people have played a role in shaping the trajectory of forest composition, and thus C storage, over time (Lewis et al 2015), be it only via fire suppression (Hurteau et al 2019).

Summary and concluding remarks
We successfully adapted CBM to match the spatial scale of the management model, combined CBM with a fire model that uses the characteristics of the landscape to simulate the current disturbance regime until the landscape-level C stabilized, and standardized the software and information processing across all three models. This modelling system attempts to follow the reproducible, transparent and reusable approach as per McIntire et al (2022) rendering it readily adaptable to other study areas and available to iterative updating (Dietze et al 2018). Given that forest management and C models are jurisdiction-specific, our example is most easily reusable in BC, can be adapted to the rest of Canada, but could be adapted to jurisdictions with only the jurisdiction-specific information modified, rather than the whole workflow. At a minimum, it provides a transparent blueprint to linking models for C-inclusion in management decisions.
The realisation of forest mitigation potential depends greatly on our ability to integrate Csequestration practices in our forest management applications. This requires robust C-estimates, an understanding of the potential for a specific landscape to sequester C, the current landscape-state relative to this potential, and the evaluation of management practices as C-sequestration tools in the midst of all the other commodities forests offer humans. Modelling efforts must strive to disentangle forest responses. Our study presents a solid basis for such an effort. Improved representation of our systems' current state will require improved data sets, the ability to iteratively update and improve simulations (as per McIntire et al 2022), which may even enable projections of the potential states of forests (Trouillier et al 2020) while using models that can support applied decisions. Our study is a step towards modelling systems that can support scientifically based decisionmaking and informed C-management.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// github.com/cboisvenue/spadesCBM_RIA.git.

Funding
This research was undertaken as a project in the Science-based improvements to the National Forest Carbon Monitoring, Accounting and Reporting System, in the work package Forest carbon tools and assessments for mitigation and reporting, in the Forest Climate Change program of the Canadian Forest Service.

Statement of authorship
C B conceived of the idea, completed the joint simulations, and wrote the manuscript; E M and C B completed the analyses and jointly developed the discussion; E M and A C are lead developers of the SpaDES toolkit, tested and helped adapt models to the modelling system; G P developed the harvesting module (ws3), participated in the conception of the research, provided harvest scenario simulations and growth representation; I E implemented and provided the C-holding capacity fire simulations and provided GIS support; all authors contributed substantially to revisions.