Evaluating the terrestrial carbon dioxide removal potential of improved forest management and accelerated forest conversion in Norway

As a carbon dioxide removal measure, the Norwegian government is currently considering a policy of large‐scale planting of spruce (Picea abies (L) H. Karst) on lands in various states of natural transition to a forest dominated by deciduous broadleaved tree species. Given the aspiration to bring emissions on balance with removals in the latter half of the 21st century in effort to limit the global mean temperature rise to “well below” 2°C, the effectiveness of such a policy is unclear given relatively low spruce growth rates in the region. Further convoluting the picture is the magnitude and relevance of surface albedo changes linked to such projects, which typically counteract the benefits of an enhanced forest CO2 sink in high‐latitude regions. Here, we carry out a rigorous empirically based assessment of the terrestrial carbon dioxide removal (tCDR) potential of large‐scale spruce planting in Norway, taking into account transient developments in both terrestrial carbon sinks and surface albedo over the 21st century and beyond. We find that surface albedo changes would likely play a negligible role in counteracting tCDR, yet given low forest growth rates in the region, notable tCDR benefits from such projects would not be realized until the second half of the 21st century, with maximum benefits occurring even later around 2150. We estimate Norway's total accumulated tCDR potential at 2100 and 2150 (including surface albedo changes) to be 447 (±240) and 852 (±295) Mt CO2‐eq. at mean net present values of US$ 12 (±3) and US$ 13 (±2) per ton CDR, respectively. For perspective, the accumulated tCDR potential at 2100 represents around 8 years of Norway's total current annual production‐based (i.e., territorial) CO2‐eq. emissions.


| INTRODUC TI ON
Balancing CO 2 emissions and removals by the latter half of the 21st century aligns with ambitious efforts ascribed under the Paris Agreement to limit global mean temperature increases to "well below 2°C" from pre-industrial levels. Aggressive mitigation of total emissions by 2050 coupled with substantial carbon dioxide removal is required to increase the chance of limiting warming to 1.5°C and 2°C by 2100 (Roe et al., 2019, Rogelj et al., 2018, and land sector measures alone could sustainably contribute around 30% (~15 Gt CO 2 -eq./year) of the required mitigation needed by mid-century to limit the increase to 1.5°C with zero overshoot (Roe et al., 2019). Terrestrial carbon dioxide removal (tCDR) is considered one of the safest and lowest cost large-scale climate change mitigation measures (Shepherd, 2009), and within the land sector, the mitigation potential of tCDR alone amounts to ~10 Gt CO 2 -eq./ year by 2050 ("roadmap" potential of Roe et al., 2019). Among these, afforestation and reforestation-or the planting of forests on nonforested or recently deforested areas-are a tCDR strategy with a combined estimated annual CO 2 removal potential between 1.5 and 10.1 Gt CO 2 -eq./year (Smith et al., 2020), with up to 4.9 Gt CO 2 -eq./ year obtainable by 2050 at a cost of under 200 US$/t CO 2 (Doelman et al., 2019).
Other tCDR strategies involving forests include improved forest management (IFM; Griscom & Cortez, 2013;Putz et al., 2008), which could give an additional tCDR potential of between 0.4 and 2.1 Gt CO 2 -eq./year without adverse effects resulting from changes to land use (Smith et al., 2020). IFM includes activities such as reduced impact logging, logged to protected forest, extended rotation age/cutting cycle, and low-productive to highproductive forest (LtHP)-with the latter including (but not limited to) a change in tree species (Verified Carbon Standard, 2019). More recent bottom-up estimates put the combined mitigation potential from both aff-/reforestation and IFM at around 5.75 Gt CO 2 -eq./ year ("roadmap" potential of Roe et al., 2019)-which is around 12% of the total mitigation needed by 2050 to limit warming to 1.5°C with no overshoot (66% confidence).
As a relatively low-cost tCDR measure, the Norwegian government is considering a policy of large-scale planting of Norway spruce (Picea abies (L) H. Karst) on mainland areas where spruce is expected to grow better than the naturally occurring species (Norwegian Environment Agency, 2019). Climate modeling studies suggest, however, that tCDR benefits of coniferous forestation projects in high-latitude regions like Norway are largely outweighed by decreases to the surface albedo (Arora & Montenegro, 2011, Bala et al., 2007, Betts, 2000, Devaraju et al., 2015, as dark coniferous canopies can "mask" surfaces with higher albedos-particularly in late winter and spring when snow is present and incident solar radiation is increasing in intensity. Conclusions from earlier modeling studies have spurred debate about the merits of coniferous forestry projects in boreal regions (e.g., Bastin et al., 2019;Veldman et al., 2019) like Norway, which has delayed their implementation and prompted new research.
New research, however, should be mindful of the caveats or limitations of global modeling studies. A first and obvious limitation is their coarse horizontal spatial resolution which is needed to reduce computational expense. Both forest productivity and surface albedo are sensitive to environmental conditions that may vary widely within a typical climate model grid cell, which often exceeds 100 km × 100 km. This is especially true for topographically complex regions like Norway where near-surface air temperature, precipitation, and incident solar radiation vary significantly at kilometer scales (Erlandsen, Tallaksen, & Kristiansen, 2019;Lussana, Saloranta, et al., 2018;Lussana, Tveito, & Uboldi, 2018b;Olseth & Skartveit, 1986). A second major limitation of many climate modeling studies, particularly those global in scope, is that it is difficult to assimilate local knowledge and information needed to identify areas considered viable candidates for forestry projects (Montenegro et al., 2009). For instance, using information about historical deforestation as a guide for constraining the geographic location and extent of simulated reforestation, Pongratz, Reick, Raddatz, Caldeira, and Claussen (2011) demonstrated that the counteracting albedo change forcing was notably reduced when compared to the findings of earlier climate modeling studies of idealized deforestation (e.g., Bala et al., 2007, Bathiany, Claussen, Brovkin, Raddatz, & Gayler, 2010, Betts, 2000, Claussen, Brovkin, & Ganopolski, 2001, Sitch et al., 2005. While the spatial extent of historical deforestation can inform about the potential extent of reforested area (as in Pongratz et al., 2011), local insight into the spatial patterns of recent and current land usage can provide additional context needed to help further constrain the identification of candidate reforestation area according to local decision-making criteria. A third and final limitation of global modeling studies is that important parameters governing CO 2 exchange and surface albedo are calibrated to perform optimally at global scales, meaning that spurious predictions may arise in specific regions. Global climate models still suffer from biased predictions of surface albedo in high-latitude forests (Li et al., 2016;Wang et al., 2016), and the parameterizations of soil and vegetation processes remain a persistent source of uncertainty in global climate model estimates of terrestrial carbon cycling (Booth et al., 2012;Todd-Brown et al., 2013;Trugman, Medvigy, Mankin, & Anderegg, 2018).
Here, we carry out a comprehensive and empirically rooted analysis to quantify the tCDR potential of large-scale spruce planting in Norway on areas which have recently undergone-or are currently undergoing-the transition to a forest dominated by deciduous broadleaved tree species (Betula spp.) of lower productivity. For areas which have already undergone the forest transition, spruce planting projects fall under the IFM category "LtHP." These areas are limited to post-abandoned secondary forests-or those regenerating largely through a natural process after abandonment of alternative land uses such as agriculture or pasture (FAO, 2003). For areas currently undergoing a process of transition but which do not yet meet the forest definition (FAO), we consider spruce planting projects as a form of accelerated forest conversion (AFC) rather than aff-/reforestation which, to our knowledge, is neither a formal IFM activity nor an activity typically considered in mitigation studies involving forestry.
We make use of Norway's National Forest Inventory (NFI) system to both identify the candidate IFM + AFC areas in Norway and initialize a set of empirical models employed to simulate transient developments in ecosystem carbon pools and surface albedo over time on those areas. The CO 2 -stock equivalence of the albedo change-driven radiative forcing is added to the simulated change in atmospheric CO 2 stock to arrive at net-accumulated tCDR estimates for 2050, 2100, and 2150 after scaling to the national level according to the total AFC + IFM area the NFI plots represent. We start by describing our methods in Section 2 followed by a presentation and analysis of results in Section 3, concluding with a discussion in Section 4 surrounding our study's merits and limitations.

| MATERIAL S AND ME THODS
Our overall methodological framework for predicting future forest carbon and surface albedo states at the NFI plots is depicted schematically in Figure 1.
The framework employs several empirical models or functions developed and calibrated within Norway, many of which using data from the Norwegian NFI. This includes those employed to model forest growth, forest mortality, biomass/volume expansion, litter production and turnover, and surface albedo-which are described in more detail in Sections 2.2-2.6.

| NFI subset and spatial representativeness
To estimate the total candidate area for AFC + IFM projects in Norway, we use the network of permanent plots operated by the Norwegian NFI. The NFI system in Norway can be characterized as a single-phase, permanent, systematic, and stratified survey employing an interpenetrating panel design where one-fifth of the sample plots evenly distributed across the country (the so-called "panel") using a Latin square design are measured each year. NFI sample plots are placed at intersections of grid lines to ensure a systematic plot distribution. NFI sampling densities vary between three strata, where a 3 × 3 km (Easting × Northing) grid is used in the lowlands including Finnmark county, a 3 × 9 km grid is used in the mountains not located in Finnmark, and a 9 × 9 km grid is used in the mountainous area of Finnmark county. The area of a stratum-or A h -is estimated by multiplying the proportion of points on the grid belonging to the stratum h with Norway's land area. The representation factor-also known as the design weight or the inverse of the sampling probability-determines how much area of Norway one sample plot represents. The representation factor of a sample plot is given by A h /n h , where n h is the number of sample plots on the grid that is specific to the stratum. That is, when a land use category is identified by a given set of attributes, the total area covered by this land use category can be determined by summing the representation factor for all plots identified with those attributes. This is a robust and well-accepted system for deriving national level land resource estimates, and as such, is relied on by many nations to report greenhouse gas emissions (GHGs) of the land use land use change and forestry (LULUCF) sector (Grassi et al., 2017;Petrescu et al., 2020) Among non-forested plots in the NFI database, we consider a plot as candidate for AFC provided the edaphic and climatic conditions are supportive of sustaining forest of medium to high productivity (i.e., site index). NFI plots classified as "Cropland," "Grassland," or "Settlements" are only considered candidates for AFC provided they are no longer managed or utilized as such (i.e., they have very recently been abandoned). As the majority of non-forested plots in the NFI are not visited by field personnel, their candidacy for AFC with respect to these two aforementioned criteria is assessed by assimilating information from two map-based products. The first provides information about the potential site productivity for the planting of Norway spruce (Bjørdal & Bjørkelo, 2006), while the second provides information about the status, continuity, and change history for agricultural landscapes within Norway (Stokstad, Fjellstad, Eiter, & Dramstad, 2017)-or in other words-whether a non-forested NFI plot is at an early state of natural forest transition.
Altogether, 1,122 plots meet the above IFM + AFC candidate area criteria (Figure 2), corresponding to a total representative area ( ∑ h A h ∕n h ) of just under 1 million hectares (9,619 km 2 ; Table 1). The majority of the representative NFI plots (~80%) are located in the coastal areas of western and northern Norway ( Figure 2).
Plots are classified as one of six UNFCCC LULUCF land use categories used by Norway in national GHG inventory reporting (Table 1). While most plots belong to the category "Forest," notable area also exists for categories "Cropland" and "Other." The latter F I G U R E 2 Spatial distribution of the included NFI plots, climate statistical summary (present day, or 2008-2017 mean), and topographic statistical summary stratified by land use land use change and forestry category. Presentday climate and topographic statistics are means of 1 km 2 grid cells. CRO, cropland; DJFMA SC, mean normalized difference snow index snow cover percentage for months December-January-February-March-April; FOR, forest; GRA, grassland; OTH, other; P ann , annual precipitation sum; SET, settlements; SW ↓,ann , annual mean incident radiation at the tilted surface; T ann , annual mean 2 m air temperature; WET, wetland constitutes-for this subset-areas recorded by the NFI as "other wooded land" following the FAO/OECD definition, and other open areas potentially capable of sustaining forest growth such as the coastal heathlands, which cover substantial areas along the west coast. A few plots in our dataset are categorized as "Grassland"-or grass-dominated areas intensively grazed by domestic animals which may or may not have sufficient trees to meet the FAO definition of forest (i.e., minimum 10% crown cover). Other LULUCF categories ("Wetland" and "Settlements") constitute only a minor percentage of the plots identified as potential candidates for AFC/IFM projects.
Regarding attribution, we consider tCDR as IFM for spruce planting simulations carried out on NFI plots classified as "Forest," and AFC for spruce planting simulations carried out on plots classified as non-forest at the start of the simulation (described further in Section 2.2). It is important to reiterate that candidate areas for AFC projects only include lands no longer under grazing or other active management-as for IFM-but with the difference being that these areas do not yet meet the definition of a forest (i.e., they have only very recently been abandoned). Thus, when simulating carbon and albedo trajectories in our baseline or no AFC/IFM scenario (described in Section 2.2), we do not suppress natural forest expansion by broadleaved deciduous pioneers onto these areas.

| Transient developments in future forest states
Simulating the long-term development in forest state on our subset of NFI plots is carried out using SiTree and with no management interventions or treatments (including final harvests). The future forest state simulated by SiTree essentially provides the foundation for our estimates of carbon bound in living biomass and soils (Section 2.4), as well as surface albedo (Section 2.5), which is depicted schematically in Figure 1.
SiTree uses the most recently published growth and mortality functions for Norway (Bollandsås, Buongiorno, & Gobakken, 2008;Bollandsås & Naesset, 2009). In SP, an initial stocking density of 2,000 spruce and 50 birch per hectare is applied at all IFM and AFC plots (Rindal et al., 2013). In BAU, both stand-level (i.e., age, site index) and tree-level (i.e., species, diameter at breast height [DBH], height) observations at the NFI plots are used to initialize the growth models employed by SiTree for IFM ( Figure 1). For ingrowth-or trees that grow beyond 50 mm at each 5 year interval-a 1-nearest neighbor imputation is carried out using the last three Norwegian NFIs (i.e., the last two remeasurement periods) as a reference database. For example, to estimate ingrowth for a certain plot (target plot) during the simulation, we look for a similar plot in similar conditions (i.e., density, site productivity, main species) in the reference database; once the similar plot (i.e., reference plot) is found, we assign its ingrowth to the target plot. In other words, the same type (species), number, and size of trees are assigned to the target plot. The effect of climate change on forest development is included through a climate-sensitive site index model calibrated for SiTree outputs a tree list with diameter, height, and species when alive and at the time of death, as well as a list of plot characteristics that includes stand age at every simulated time step. Volume and biomass are calculated from tree diameter and height using species-specific individual tree volume equations (Braastad, 1966;Brantseg, 1967;Vestjordet, 1967). Biomass is estimated using species-specific allometric equations (Marklund, 1988;Smith, Granhus, & Astrup, 2016) with tree diameter and tree height as independent variables. Biomass transfers from tree and ground vegetation to the soil pool are described in Section 2.4.

| Transient climate changes
Ensemble estimates of future near-surface air temperature, precipitation, and surface snow pack are provided through the Norwegian Centre for Climate Services (https://klima servi cesen ter.no/) at 1 km spatial resolution and daily time steps (Hanssen-Bauer et al., 2017).
The atmospheric variables originate from combined statistical downscaling and bias correction of nine regional climate model simulations from the EURO-CORDEX archive (numbers 1-6 and 8-10 in table 1 of Wong, Haddeland, Lawrence, & Beldring, 2016). We limit our analysis to simulations under the representative concentration pathway (RCP) 4.5 for two reasons. First, the downscaled and biascorrected EURO-CORDEX data for Norway are only available for RCP4.5 and RCP8.5, and the current global emission trajectory more closely resembles RCP4.5 than RCP8.5, the latter of which being more of a worst case than a business-as-usual scenario (Hausfather & Peters, 2020). Second, RCP4.5 is the only RCP compatible with TA B L E 1 Number of NFI plots by land use classification and total representative area in Norway. The category "Other" comprises "other wooded land" and shrub/heathland areas potentially capable of sustaining forest growth  Thomson et al., 2011). Coherent estimates of future surface snow pack (as water equivalent) are derived from the hydrological model HBV (Beldring, Engeland, Roald, Saelthun, & Voksø, 2003) which simulates snow accumulation and depletion employing a degree-day approach with land cover-specific melt factors. Daily means are used to compute monthly means for 5 year timeslices centered over the simulation time step (e.g., monthly means for 2018-2022 centered over 2020 as the first time step, means of 2023-2027 centered over 2025 as the second time step, etc.). Because our climate data time series ceases at 2100, we use the 10 year monthly means for the period 2091-2100 for each succeeding time step remaining in the simulation period (see Section S1 for a summary of climate statistics at select time steps).
The climate dataset described above is augmented here with coherent, high-resolution projections of downwelling solar radiation incident on the tilted surface (SW sfc ↓ ). First, the mean downwelling solar radiation flux incident on the horizontal surface (SW sfc ↓ ) is derived by downscaling monthly surface downwelling shortwave radiation available from the EURO-CORDEX-ensemble following the quantile mapping methods described in Wong et al. (2016) Table S1) and constrained by empirical mass loss, the parameters determining humus dynamics are primarily constrained by soil C stock data. Application to Norwegian conditions is further described in  and .

| Estimating carbon pools at stand level
Yasso07 inputs include litter (including large residues) as described in the Norwegian GHG inventory (Norwegian Environment Agency et al., 2020), mean annual temperature, annual precipitation, and annual temperature amplitude. Litter includes foliage, branch, and root biomass shed from living trees (Table S2) and natural tree mortality as estimated from the SiTree biomass development and plot characteristics (Table S3; Antón-Fernández & Astrup, 2012), ground vegetation biomass (Table S4; , and turnover (Peltoniemi, Mäkipää, Liski, & Tamminen, 2004; Table S5) where belowground biomass for shrubs, herbs, and grasses in the IFM and AFC simulations (after transition to forest) is assumed to be twice the aboveground biomass. Assumed dimensions for litter are 0 cm (non-woody), 2 cm (fine-woody), and 10 cm (coarse-woody). The litter input estimates and model application are done for each NFI plot individually. Each plot is initialized by a spin-up to equilibrium in five pools (e.g., the three labile + insoluble + humus) using the average litter input 1990-2016 grouped by site index and dominant tree species (Table S6). Forested plots are then subject to a pre-simulation or "backcast" of 30 years (6 × 5 year time steps) as described in  and . This backcast is based on land use prior to NFI plot establishment and national standards for harvest volumes and ages for the relevant time period.
The SiTree-Yasso07 framework was originally developed for forested areas. In the current study, however, non-forested plots have been included as described in Section 2.   (Table S7) and assuming litter quality as for forest herbs and grasses (Table S1). where C is carbon in terrestrial pool i and 44/12 is the ratio of CO 2 to C molar mass. For assessing transient CO 2 radiative forcing (RF CO 2 ) over time (Section 2.5), terrestrial carbon stocks (as CO 2 ) are first converted to fluxes by differentiating Equation (2) with respect to time.

| Estimating surface albedo at stand level
We employ the empirical models of Bright and Astrup (2019)-developed in Norway-which are based on satellite optical remote sensing and which predict the monthly surface albedo as a function of land cover type, forest structure (if forest), monthly near surface air temperature, and monthly snow cover. These models provide estimates of the monthly mean "black-sky" and "white-sky" surface albedos (entire shortwave broadband) at local solar noon. Because the actual surface albedo is a linear combination of the two that is weighted by the fraction of diffuse solar radiation, we estimate the diffuse fraction using an empirical function (Erbs, Klein, & Duffie, 1982) that relates the monthly clearness index (the ratio of terrestrial-toextraterrestrial solar radiation) to diffuse solar radiation at the surface. Because the black-sky albedo depends on solar geometry and varies over the diurnal cycle (being higher in the mornings and evenings than at midday), the true daily mean of the black-sky albedo will typically be higher than the black-sky albedo registered by the satellite product at the local solar noon (i.e., when zenith angles are lowest). However, Wang et al. (2015) recently compared monthly mean observations of the true daily mean and the local solar noon albedos for a variety of land cover types and found differences of ~1% irrespective of the month or land cover type. Because the primary objective here is to quantify the albedo difference between two land cover types (i.e., spruce forest and other forest or non-forest), this ~1% deviation is irrelevant to consider.
The albedo models of Bright and Astrup (2019) require inputs of snow cover percentage based on the normalized difference snow index (NDSI) as provided by the version 6 MOD10 Snow Cover product (Hall & Riggs, 2016). Since these inputs are not available for our forward-looking simulations, we develop an empirical function through ordinary least squares regressions that give NDSI snow cover percentage as a function of land cover type and monthly mean SWE (see Section S3 for details).
For forests, the models of Bright and Astrup (2019)  are applied, with results weighted by the relative volume share of each tree species simulated at each NFI plot ( Figure S4). For the AFC simulations, the forest models are also applied for two reasons: (a) model parameters for forests with zero volume are similar to the non-forest model parameters; and (b) the initial non-forest land use type transitions to a forest within the first two decades after the start of the simulation. This decision has little effect on the simulated surface albedo differences between SP and BAU (see Figure S7).

| Radiative forcing and metrics for surface albedo change
Estimates of the local annual mean instantaneous shortwave RF (in W/m 2 ) at the top-of-the-atmosphere following a surface albedo change (or difference) are based on a simple radiative transfer parameterization described in Bright and O'Halloran (2019): where Δ m,t is a surface albedo change (or difference) in month m and year t and T m,t is the all-sky monthly mean clearness index (unitless) computed with SW sfc ↓ (i.e., T = SW sfc ↓ ∕SW toa ↓ ) in month m and year t. After normalization to the Earth's total surface area, we convert the annual RF Δ to annual CO 2 -equivalent emissions (or removals) following the time-dependent emissions equivalent (TDEE) method described in Bright et al. (2016) that makes use of a CO 2 impulseresponse function (described below) to account for the time-dependent abundance of CO 2 in the atmosphere following emissions. See Section S4 for additional details describing the TDEE measure and (2) ΔCO 2 (TH) = − 44 12 a comparison to two alternative methods for converting RF Δ into "CO 2 -equivalence." Annual RF Δ is also compared directly to the annual RF CO 2 that results from changes to CO 2 emissions or removals by terrestrial carbon sinks. RF CO 2 is estimated from changes to atmospheric CO 2 abundance computed as a convolution integral between emissions (or removals) and a CO 2 impulse-response function: where t is the time dimension, tʹ is the integration variable, e(tʹ) is the CO 2 emission (or removal) rate (in kg)-or difference in atmospheric CO 2 flux between the SP and BAU scenarios, y CO 2 is the multi-model mean CO 2 impulse-response function described (Joos et al., 2013; Myhre & Shindell, 2013) for a CO 2 background concentration of 389 ppmv, and k CO 2 is the radiative efficiency per kg CO 2 emitted upon the same background concentration (i.e., 1.76 × 10 -15 W/m 2 kg -1 ) which is assumed constant and time-invariant given the small size of the perturbation (i.e., "e").
We include the albedo change effect in our definition of tCDR: where ΔCO 2 is the difference in atmospheric CO 2 stock between SP and BAU (i.e., Equation 2) at a given TH, and ∑ TDEE is the atmospheric CO 2 stock equivalence of the surface albedo difference at the same TH (i.e., the solution to Equation S3 summed up to TH).

| Net present value of tCDR
We estimate net present value (NPV) for each scenario and TH after revenue and expense estimates are first carried out for each plot (n = 1,122) and time step: where R are revenues, E are expenses, t is the time step, and r is a discount rate of 4%. Expenses include those associated with stand establishment and the operating costs associated with hypothetical future harvesting at the same three policy THs (i.e., 2050, 2100, 2150) for which tCDR estimates are quantified. Establishment costs are obtained at the county level and include the costs of clearing, site preparation, planting, and supplementary planting within the first 5 years, and one or more rounds of tending until the stock has exceeded 5 cm in DBH. Operating costs are estimated using a cost-calculation model described in Fønhus (2018) and include the costs of harvesting, forwarding, and road transport to industry. Harvesting costs are estimated on the basis of the simulated tree size, simulated tree density, and terrain conditions, while forwarding costs depend on the simulated volume per hectare, forwarding distance, and terrain characteristics.
Revenues for each NFI plot and time step are estimated following the method of Blingsmo and Veidahl (1992), where the value of each tree is based on saw and pulpwood prices at the municipality level (mean for 2013-2017) and the simulated tree sizes (i.e., DBH and height outputs from SiTree). In cases where there are no price statistics from municipalities, county-level price statistics are used instead. Revenues and expenses are compiled in Norwegian crowns (NOK) and converted to US dollar ($) using an exchange rate of 8.5 NOK/US$.
The difference in NPV at a given TH per metric ton tCDR is then computed, treating the sign of tCDR as positive with respect to the terrestrial biosphere (i.e., where net CO 2 -eq. atmospheric removals are positive).

| RE SULTS
Only the difference in outcomes between the SP and BAU scenarios is presented henceforth. The reader is referred to Section S5 for absolute results for each scenario, and for the disaggregation of net terrestrial CO 2 emissions/removals by pool (i.e., LBC and SOC). Given that the total representative area for "Grassland," "Settlements," and "Wetland" constitutes only 8.9% of the total AFC area potential (and 1.2% of IFM + AFC), results for these land use categories are merged and presented henceforth as "GSW."

| Instantaneous RF and atmospheric CO 2 -eq. flux differences at stand level
Means and standard deviations for RF CO 2 and RF Δ per hectare and land use category are presented in Figure 3. Local RF Δ is plotted as the right-hand y-axis for the reader's convenience. Starting with IFM, a large and sustained positive RF CO 2 of ~3 × 10 -14 W/m 2 is seen for the first ~20 years (Figure 3a, blue curve) which is attributed primarily to a large reduction in the LBC pool at the start of SP (the initial stock change for step t 1 − t 0 is treated as a CO 2 emission pulse at t 1 ; see Figures  For the AFC cases, RF Δ and RF CO 2 dynamics and magnitudes are similar to the IFM case over the medium to longer terms (Figure 3bd). Over the short term (2020-2060), however, the magnitude of the positive RF CO 2 and negative RF Δ is lower compared to the IFM case given smaller differences in initial carbon stocks ( Figure S5) and surface albedos ( Figure S7).
Converting RF Δ to TDEE (Equation S3) and plotting it together with the simulated CO 2 flux difference between scenarios (SP minus BAU) yields Figure 4.
Converting deciduous-dominant forest to spruce forest (Figure 4a) results in a mean CO 2 emission of 35 (±25) t ha −1 year −1 at 2025-of which 28 t ha −1 year −1 stems from the reduction of the LBC pool in SP, 4 t ha −1 year −1 as an emission from the SOC pool in SP, and 3 t ha −1 year −1 from foregone CO 2 removals in BAU ( Figure S6). At the same time step, however, the clear-felling of existing forest greatly increases the annual mean surface albedo by around 4% (Figure S7) resulting in a CO 2 -equivalent removal effect of around −16 (±11) t ha −1 year −1 . It is important to note here that while TDEE is an instantaneous measure, the sign and magnitude of the CO 2 equivalent effect in time do not necessarily correspond to that of RF Δ (and Δ ). As the planted spruce forest matures, the albedo difference switches from positive to negative over the medium term due to enhanced surface masking by forest canopies, resulting in a CO 2 -eq. emission effect of up to 21 (±10) t ha −1 year −1 at 2050. Over the longer term, the mean forest CO 2 sink is enhanced (Figure 4a, blue curve) leading to sustained atmospheric CO 2 removals, while surface albedo differences decrease F I G U R E 3 (a)-(d) Mean RF Δ and RF CO 2 from differences in mean atmospheric CO 2 -eq. fluxes ("ΔCO 2 sink") and surface albedo ("Δ ") between spruce planting (SP) and business-as-usual (BAU) scenarios per land use land use change and forestry category and hectare (±1 SD; "1σ"). "Δ " is the time-dependent CO 2 emission equivalence of the surface albedo difference between SP and BAU; "ΔCO 2 sink" is the difference in the atmospheric CO 2 flux between SP and BAU. For "∆NET," 1σ (gray shaded interval) is computed assuming independence between Δ and ΔCO 2 sink. Vertical dashed green lines indicate the three time horizons (TH) of interest. Note that local RF Δ is plotted on the right-hand y-axes. AFC, accelerated forest conversion; GSW, grasslands, settlements, wetlands; IFM, improved forest management ( Figure S7), resulting in a mean net CO 2 -eq. removal that is dominated by the terrestrial CO 2 sink (Figure 4a, black curve). Beyond 2100, however, forest mortality and respiration in soils increase in SP relative to BAU which, for some plots, results in net CO 2 -eq. emissions (Figure 4a, upper bound of 1σ gray shaded area).
For AFC, the temporal dynamics of emissions and removals from Δ and ΔCO 2 -sink on all AFC land use categories (Figure 4b-d) resemble those of the IFM case, although emissions and removals are lower in the short term as differences in forest age and thus ecosystem productivity between SP (i.e., planted spruce forest) and BAU (i.e., naturally regenerated deciduous-dominant forest) are more negligible.

| Transient development in atmospheric CO 2 -eq. stock at stand level
Figure 5 presents mean results at stand level for the six land use categories expressed as an atmospheric CO 2 -eq. stock difference F I G U R E 4 (a)-(d) Differences in mean atmospheric CO 2 -eq. fluxes between spruce planting (SP) and business-as-usual (BAU) scenarios per land use land use change and forestry category and hectare (±1SD; "1σ"). For "∆NET," 1σ (gray shaded interval) is computed assuming independence between Δ and ΔCO 2 sink. "Δ " is the time-dependent CO 2 emission equivalence of the surface albedo difference between SP and BAU; "ΔCO 2 sink" is the difference in the atmospheric CO 2 flux between SP and BAU. Vertical dashed green lines indicate the three time horizons (TH) of interest. 1 t/ha = 0.1 kg/m 2 . AFC, accelerated forest conversion; GSW, grasslands, settlements, wetlands; IFM, improved forest management (i.e., tCDR) between SP and BAU over time. For "Forests" and thus for IFM (Figure 5a), a CO 2 debt with a 45 year payback period is incurred as existing forests are first harvested resulting in large emissions from soil and living biomass pools in the short term (see and −780 (±301) t CO 2 -eq./ha at 2100 and 2150, respectively. In general, IFM gives lower mitigation benefit at all three THs relative to AFC.
For AFC (Figure 5b-d), the CO 2 debt payback periods are shorter because CO 2 emissions from soil and living biomass pools in the SP scenario are substantially lower in the short term ( Figures S5 and   S6). Because candidate areas in the AFC counterfactual scenario (i.e., BAU) are abandoned lands assumed to revert naturally to forest, negative surface albedo differences between SP and BAU over the medium to longer term are notably lower than would be expected F I G U R E 5 (a)-(d) Differences in mean atmospheric CO 2 -eq. stock between spruce planting (SP) and business-as-usual (BAU) scenarios per land use land use change and forestry category and hectare (±1SD [σ]). For NET, 1σ (gray shaded interval) is computed assuming independence between Δ and ΣCO 2 sink. "Δ " is the time-dependent CO 2 emission equivalence of the surface albedo difference between SP and BAU accumulated over time; "ΣΔCO 2 sink" is the difference in the atmospheric CO 2 flux between SP and BAU accumulated over time (i.e., the solution to Eq. (2)). 1 t/ha = 0.1 kg/m 2 . AFC, accelerated forest conversion; GSW, grasslands, settlements, wetlands; IFM, improved forest management otherwise had the reference (BAU) areas remained actively managed (i.e., had the original land use types at the start of the simulations been preserved or remained open; see Figure S7). As a result, the magnitude of the CO 2 -eq. impact of surface albedo change in the longer term is similar for all AFC land use types and on the same order as that for Forest (IFM)-or 40 (±25) and 34 (±24) t CO 2 -eq./ha at 2100 and 2150, respectively.
For 2100, largest tCDR potential is found for "Other" (Figure 5d; −659 ± 278 t CO 2 -eq./ha), with the total area-weighted tCDR potential of AFC estimated to be −621 (±304) and −937 (±334) t CO 2 -eq./ ha at 2100 and 2150, respectively. Because the total area classified as "Forest" greatly exceeds all others, the atmosphere would see a net increase in the stock of CO 2 by 151 (±122) Mt CO 2 -eq. at 2050 (Figure 6a) given our stand-level findings presented in Section 3.2. A net atmospheric stock decrease of −342 (±198) and −650 (±221) Mt CO 2 -eq. would, however, be realized at the 2100 and 2150 horizons for IFM (Figure 6a). For "Cropland" (Figure 6b) and "Other" (Figure 6d), the atmosphere may F I G U R E 6 (a)-(d) Means and SD of the estimated terrestrial carbon dioxide removal (tCDR) potential (atmospheric perspective) in Norway from AFC and IFM per land use land use change and forestry category, including surface albedo change radiative forcings (as ΣTDEE). The terrestrial CO 2 pool is disaggregated into living biomass ("LBC") and soils ("SOC"). Combined variance of "net tCDR" is estimated assuming independence between ΔLBC, ΔSOC, and Δα. Note the difference in scales for y-axes. AFC, accelerated forest conversion; GSW, grasslands, settlements, wetlands; IFM, improved forest management; TDEE, time-dependent emissions equivalent see either an increase (i.e., no tCDR) or decrease in the stock of CO 2 at 2050 depending largely on factors driving differences in primary productivity ( Figure 6 "ΔLBC," 1 − σ whiskers) between SP and BAU.

| tCDR potential of AFC and IFM in Norway
For "GSW" (Figure 6c), whether tCDR is realized at 2050 depends more on factors driving differences in SOC between SP and BAU.
The tCDR potential of AFC at 2100 and 2150 is between one and two orders of magnitude lower than that of IFM.
In general, although the net total tCDR potential (IFM + AFC) is driven by increases to CO 2 bound in living biomass (LBC; 365 and 601 Mt CO 2 at 2100 and 2150, respectively), an additional 99-214 Mt CO 2 is bound in soil (SOC) at 2100 and 2150, respectively.
On average, decreases to the surface albedo are only found to reduce the 2100 and 2150 tCDR potential by around 15% and 7%, respectively. Table 2 presents mean ∆NPV per metric ton tCDR computed only for those NFI plots where tCDR is realized-that is, where there is a net CO 2 -eq. removal from the atmosphere at a given TH. Although it is not considered realistic to harvest a Norway spruce stand after 30 years (i.e., at 2050), ∆NPV at 2050 is presented to be consistent with the 2050 mitigation potentials presented in Sections 3.1-3.3.

| Cost of tCDR
Mean ∆NPV per metric ton tCDR for both AFC and IFM is negative at 2050 (in 2019 USD) which can be interpreted as a mitigation cost of $14 and $25, respectively. Costs are greater for IFM given greater establishment costs and short-term emissions in SP. At 2100 and 2150, tCDR becomes profitable-that is-mean ∆NPV per metric ton tCDR is positive, although differences between IFM and AFC are minor. Differences between 2100 and 2150 are relatively minor as increases to ∆NPV and tCDR scale similarly.

| D ISCUSS I ON
We employed a bottom-up and spatially inexplicit method to estimate the tCDR potential of large-scale spruce forest planting in Norway, making use of local information surrounding viable candidate areas for IFM and AFC projects, forest carbon cycle dynamics, and surface albedo dynamics. Areas considered as candidates for IFM projects (spruce planting) comprised 86% of the total IFM + AFC area and were carefully selected to minimize adverse impacts to biodiversity while ensuring forest productivity gains. Areas identified as candidates for AFC were carefully selected to avoid competition with food and fodder production. That is, for arable lands, only recently abandoned croplands and pastures were considered. It is important to emphasize that what we have labeled "AFC" would technically be reported as an "aff-/reforestation" activity under the Kyoto Protocol.
We believe, however, that the AFC label is more appropriate since our counterfactual (or baseline) land use scenario (BAU) involves a process of natural forestation. While this makes our AFC results difficult to compare to those of many climate or integrated assessment modeling studies focused on aff-/reforestation where the "open" usage of the baseline is often preserved, we believe that the realism of our baseline land use scenario is a major strength of the study and that such realism would be difficult to replicate in studies larger in geographic scope. Given that recovering secondary vegetation (which includes old fields and recovering forests) covers nearly twice the area globally than croplands (Isbell, Tilman, Reich, & Clark, 2019), and given the need to maintain croplands to preserve food security (Smith et al., 2020), the use of a baseline in reforestation studies that ignores the recovering vegetation seems inappropriate and unjustified.
Our results are robust with respect to the general net outcome simulated for the second half of the 21st century-that spruce planting on abandoned lands currently undergoing a process of natural forestation results in a net atmospheric CO-eq. stock decrease (tCDR) by 2100. After controlling for simulated changes to site indices (site productivities) over this period, our mean stand-level living biomass estimates (and thus mean LBC stocks) in both scenarios agree well with observations ( Figure S9), as do our mean simulated rates of change (biomass increments; Figure S10)-giving us confidence in our LBC results which dominate the net (tCDR) outcomes ( Figure 6).
We are also confident in our albedo results, since the albedo predictions are based on empirical models calibrated within Norway for the same forested ecosystems, which have normalized prediction biases of <10% (see e.g., figure S9 of Bright and Astrup, 2019).
TA B L E 2 Mean ∆NPV (SD) in 2019 USD per metric ton tCDR (as CO 2 -eq.) to global-scale (Arora & Montenegro, 2011, Bala et al., 2007, Betts, 2000, Devaraju et al., 2015 studies that considered aff-/ reforestation projects at high latitudes, we find that the offsetting impact of surface albedo decreases is relatively low and negligible, which we attribute simply to our different assumptions surrounding the reference land usage (or cover type): here, the reference albedo is that of a growing forest, whereas in most previous studies,  Figure S11), whereas fine root litter inputs in SP (spruce forests) are also likely biased high ( Figure S12). These biases suggest that the simulated SOC sink difference between SP and BAU may be too high, although this is difficult to verify given the scarcity of litter turnover observations for the study region. According to Mayer et al. (2020), there is little evidence to support the notion that a choice in tree species alone can affect the magnitude and stability of soil C stocks, since this largely depends on both biotic and abiotic site properties that vary largely in space. Nevertheless, given the relative importance of the LBC sink, a 100% increase/decrease to the simulated difference in soil C stock between the SP and BAU scenarios-or 99 and 214 Mt CO 2 at 2100 and 2150, respectively ( Figure 6)-would not affect the main finding that the terrestrial CO 2 sink is enhanced resulting in reduced atmospheric CO 2 concentrations by the end of the 21st century.
Although this general overall finding is robust, the timing of the estimated tCDR benefits is uncertain, and, like any forecast, this uncertainty increases with increasing TH. While we have accounted for temperature and moisture change feedbacks on forest growth via application of site index models , these models do not capture the effects of changes to atmospheric CO 2 concentrations and water use efficiency (tree physiology and adaptability). There is some evidence from carbon isotropic discrimination studies that intrinsic water use efficiency has increased by ~20% over the 20th century for many boreal species (Saurer, Siegwolf, & Schweingruber, 2004;Saurer et al., 2014); however, there is also evidence suggesting that tree growth has not increased as expected due to overriding factors such as nutrient limitation or physiological acclimation to elevated CO 2 concentrations (Peñuelas, Canadell, & Ogaya, 2011).
Regarding nutrient limitation, the site index models we employed also do not explicitly take into account nutrient status and availability, only implicitly through a variable describing soil quality (see Section S1). Boreal forests are generally nitrogen (N) limited, and the availability of N generally constrains biomass production (Högberg, Näsholm, Franklin, & Högberg, 2017;Tamm, 1991). For instance, Van Sundert, Horemans, Stendahl, and Vicca (2018) found that C/N ratios as well as N% (SOC%) strongly explained most of the forest productivity variation across Sweden after factoring out climatic effects. Furthermore, N limitation has been found to curtail the response of forest biomass production to temperature changes (Sigurdsson, Medhurst, Wallin, Eggertsson, & Linder, 2013). There is some evidence suggesting that total N is influenced by N deposition (Olsson et al., 2009, Solberg et al., 2004, and since N deposition increases along a north-south gradient in Norway, nitrogen limitations and the associated uncertainty of our simulated growth responses may be larger for NFI plots located more in the north. See Section S6 for a more in-depth discussion on N-limitation. Although our study focuses on global radiative balance effects, it is important to note that forestry projects involving aff-/reforestation or tree species change can affect local and regional climate through its influence on energy and moisture budgets at or near the surface (Bonan, 2008;Devaraju, De Noblet-Ducoudré, Quesada, & Bala, 2018;Luyssaert et al., 2018). Forestation increases surface roughness and often total evapotranspiration, which can affect turbulent heat flux partitioning at the surface and both the transmissivity and emissivity of the lower atmosphere, in turn affecting air temperatures at or near the surface (Alkama & Cescatti, 2016;Lee et al., 2011;Novick & Katul, 2020). Forestation can also affect energy budgets of the surface and lower atmosphere via changes to biogenic aerosol regimes, which can affect atmospheric chemical and physical processes influencing both short-and long-wave radiation transfer (Spracklen, Bonn, & Carslaw, 2008;Unger, 2014 in Norway but has double the productivity (i.e., mean annual atmospheric CO 2 sequestration over the rotation period) to that of Norway spruce in the coastal regions of western Norway (Godal & Grønlund, 2014) where the largest total IFM + AFC candidate area potential is found. The reason behind the blacklisting of Sitka spruce is due to the potential threat on biodiversity outside plantations and further spreading, although some research has shown that such "invasion" threats can be greatly reduced if Sitka spruce plantations are given a 200 m buffer from protected areas (Nygaard & Øyen, 2017).

| CON CLUS IONS
In summary, we find that when spruce forests supplant existing broadleaved deciduous-dominant forests (so-called IFM), the atmospheric CO 2 stock decreases on average by −467 (±293) t CO 2 / ha at 2100 although this is offset by 43 (±23) t CO 2 -eq./ha (~12%) when taking into account the surface albedo difference between spruce and broadleaved forests. This reduction is not trivial as most of our NFI plots and the areas they represent lie in coastal regions of western Norway where climate conditions are favorable to Norway spruce production. We estimate even greater reductions at 2100 when spruce is planted on area presently classified at open (i.e., AFC projects)-or −661 (±302) t CO 2 /ha-which is offset on average by 40 (±25) t CO 2 -eq./ha due to the surface albedo change. However, it is important to keep in mind that the area suitable for AFC projects represents only 14% of the total area we considered as viable for IFM + AFC projects in Norway (Table 1).
Considering the total viable area (IFM + AFC), we estimate a tCDR potential (i.e., net atmospheric CO 2 -eq. stock change) of −447 Mt CO 2 -eq. at 2100, of which −340 Mt CO 2 -eq. is attributable to IFM. However, IFM + AR projects in Norway would increase atmospheric CO 2 by 152 Mt CO 2 -eq. at 2050 since the majority of tCDR projects involve IFM (tree species conversion) leading to large net CO 2 emissions from a reduced short-term primary productivity in forests and an increased CO 2 emission from soil and living biomass pools. The net atmospheric CO 2 -eq. stock change at 2050 from IFM projects alone would amount to 149 Mt CO 2 -eq. despite the counteracting effect of the short-term increase to the surface albedo from the initial clear-felling of existing forests. While this net emission is balanced by net removals in the second half of this century (i.e., −447 Mt CO 2 -eq. at 2100)-an outcome that is in line with the text of the Paris Agreement-it does not however align with the EU's policy target of becoming climate neutral by 2050 (European Commission, 2018). Most emission scenarios suggest that net global CO 2 emissions need to reach zero by about 2070 or 2050 to remain "likely below" a warming target of 2°C or 1.5°C, respectively (Rogelj et al., 2016), although the timing of zero emissions is a poor indicator of the likelihood of achieving a temperature target (in contrast to the total CO 2 emitted over time; Matthews, Zickfeld, Knutti, and Allen, 2018).
Although the estimated tCDR potential for Norway is even higher at 2150 relative to 2100 (Figure 6), the reversibility risk increases as the forest carbon stock becomes more vulnerable to changes in pest-or climate-driven disturbance regimes. As for fire regime changes, fire risk is inherently lower in Norway relative to boreal North America and central Asia given its wetter climate and different native tree species (Rogers, Soja, Goulden, & Randerson, 2015). Within Norway, fire activity has decreased since the 19th century due to effective fire suppression (more intensified forest management; Rolstad, Blanck, and Storaunet, 2017), and the future fire regime in Norway may even benefit from future regional climate changes (Flannigan, Bergeron, Engelmark, & Wotton, 1998).
In conclusion, Norway's total estimated net tCDR potential of −447 Mt CO 2 -eq. by 2100 equates to a global cooling of around −6 × 10 -4°C when applying an observationally constrained measure of the transient climate response to cumulative carbon emissions . Framed differently, this amount offsets roughly 8 years of Norway's production-based GHG emissions (Eurostat, 2017)-or a single year of CO 2 -eq. emissions embodied in Norway's oil and gas exports (Andrew, 2019). Although seemingly negligible when compared to Norway's current emission rate, tCDR via Norwegian AFC and IFM is profitable-with ∆NPV estimated here to be US$ 12 (±3) per ton CO 2 -eq. avoided at 2100. This is in contrast to most mitigation measures of the Norwegian land transport and process industry sectors which cost well over US$ 100 per ton CO 2 -eq. avoided (Faehn & Isaksen, 2016). While there is some uncertainty behind our ∆NPV estimates for AFC and IFM projects on steep terrain-such is often found on the west coast of Norway-even a doubling of costs would make these activities relatively cheap mitigation measures and perhaps even remain profitable. Regarding the implementation of large-scale IFM + AFC, the largest barrier is likely the gaining of social acceptance, which will require an open, honest, and scientifically informed dialogue about the trade-offs in ecosystem services (e.g., recreation, biodiversity) or other values provided by the current land usage.
Lastly, it is important to emphasize that we have not considered any additional longer term mitigation benefit that might be incurred upon forest harvest and utilization. For example, the increased supply of timber originating from AFC + IFM forests around the turn of the century could displace more emission-intensive materials in the global construction sectors (Kohlmaier, Weber, & Houghton, 1998;Valsta et al., 2017).

ACK N OWLED G EM ENTS
This research was supported by the Research Council of Norway, grant number 194051. The authors disclose no conflict of interest.
The authors thank Lars Vesterdal for making the original data behind a number of Danish litter flux experiments available to us.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.