Large‐scale forest restoration stabilizes carbon under climate change in Southwest United States

Abstract Higher tree density, more fuels, and a warmer, drier climate have caused an increase in the frequency, size, and severity of wildfires in western U.S. forests. There is an urgent need to restore forests across the western United States. To address this need, the U.S. Forest Service began the Four Forest Restoration Initiative (4FRI) to restore four national forests in Arizona. The objective of this study was to evaluate how restoration of ~400,000 ha under the 4FRI program and projected climate change would influence carbon dynamics and wildfire severity from 2010 to 2099. Specifically, we estimated forest carbon fluxes, carbon pools and wildfire severity under a moderate and fast 4FRI implementation schedule and compared those to status quo and no‐harvest scenarios using the LANDIS‐II simulation model and climate change projections. We found that the fast‐4FRI scenario showed early decreases in ecosystem carbon due to initial thinning/prescribed fire treatments, but total ecosystem carbon increased by 9–18% over no harvest by the end of the simulation. This increased carbon storage by 6.3–12.7 million metric tons, depending on the climate model, equating to removal of carbon emissions from 55,000 to 110,000 passenger vehicles per year until the end of the century. Nearly half of the additional carbon was stored in more stable soil pools. However, climate models with the largest predicted temperature increases showed declines by late century in ecosystem carbon despite restoration. Our study uses data from a real‐world, large‐scale restoration project and indicates that restoration is likely to stabilize carbon and the benefits are greater when the pace of restoration is faster.


Introduction
In recent decades, the western U.S. has seen an increase in the frequency, size, and severity of wildfires (Dillon et al. 2011, Dennison et al. 2014. A century of fire suppression, along with a warmer climate, below-average winter precipitation, and earlier spring snowmelt have contributed to this increase (Westerling 2016, Abatzoglou andWilliams 2016). This phenomenon is particularly pronounced in the southwestern U.S., where the trend of warmer temperatures and extended droughts are projected to increase in the future. Climate-induced changes are expected to lead to increased tree mortality due to drought stress, insect

Accepted Article
This article is protected by copyright. All rights reserved. outbreaks, and larger and more severe wildfires (Williams et al. 2010, McDowell et al. 2015. In addition to the economic costs, destruction of habitat, and air and water quality (Western Forestry Leadership Council 2010, Combrink et al. 2013, large wildfires can reduce carbon storage for decades after fire Allen 2002, Hurteau et al. 2008).
There is an urgent need to accelerate the pace and scale of restoration in forests across the western U.S. to reduce vulnerability to catastrophic wildfire and improve resilience to climate Restoration of forest structure (mechanical thinning) and processes (prescribed fire) is one of the few adaptation strategies available to forest managers. Because 4FRI's goal is to accelerate the pace of restoration across tens of thousands of hectares, 4FRI presents a rare opportunity to evaluate the extent to which restoration efforts aimed at restoring fire regimes can also be applied as mitigation strategies to stabilize carbon under climate change. Previous studies have projected that restoration under climate change will result in long-term increases in forest carbon due to reductions in wildfire activity (Loudermilk et al. 2016, Krofcheck et

Accepted Article
This article is protected by copyright. All rights reserved.
al. 2017, Hurteau 2017, Liang et al. 2018. Mediated in part by low regeneration rates under warmer temperatures, other studies have projected that extreme climate change will result in forest type changes, biomass loss, and deforestation despite restoration treatments (Tarancón et al. 2014, Flatley and Fulé 2016, Loehman et al. 2018. The objective of this study was to explore how the pace of restoration in ponderosa pine forests under a large-scale restoration program would influence carbon dynamics and wildfire severity given projected changes in climate. Specifically, we estimated forest carbon fluxes and pools and wildfire severity under a moderate and rapid 4FRI implementation schedule and compared those to status quo restoration and no harvest using the LANDIS-II simulation model. This study had several unique features that add value to previous efforts. It was one of the few that was conducted at a large-scale and under an accelerated pace of restoration. Additionally, we incorporated U.S. Forest Service's (USFS) 4FRI implementation plan, with initial thinning prescriptions followed by maintenance fire into the LANDIS-II model to match, as closely as possible, the actual restoration underway in the forests for the next several decades. Finally, we estimated the fate of carbon from harvested products and the resulting storage/emissions. In comparison to status quo restoration and no harvest, we hypothesized that ecosystem carbon would initially decline under accelerated thinning and prescribed fire but subsequently increase by reducing fuel loads available for wildfires under warmer and drier conditions caused by climate change.

Accepted Article
This article is protected by copyright. All rights reserved.

Study Area
The study area was the initial phase of 4FRI, approximately 400,000 ha of the Coconino and Kaibab National Forests including the city of Flagstaff in northern Arizona, U.S.A. (Figure 1). We included the entire footprint of 4FRI Phase 1 in our study but the actual 4FRI EIS did not include some areas within that footprint, as they already had separate, previously approved restoration plans that were similar in scope to 4FRI. We simulated restoration of the entire project area following the methods used by USFS for the 4FRI restoration.
The elevation of the study area ranges from 1780 to 3850 m with a mean elevation of 2190 m. The average temperature is 10.2° C and average annual precipitation is ~460 mm.
The 4FRI landscape is dominated by contiguous ponderosa pine tree cover intermixed with lesser amounts of dry and wet mixed conifer forests types and oak. Ponderosa pine forests in this area have shifted from naturally open conditions to high densities of small diameter trees in the last century, due to fire suppression, harvesting, and episodic regeneration events (Covington and Moore 1994).

No Harvest/Control:
In this scenario, we assumed no harvest or prescribed fire occurred in the project area for the entire length of the simulation. Wildfire was still simulated.

Status quo restoration:
We estimated the current pace of restoration by averaging the area across which thinning and prescribed fire occurred in years 2010-2016. This included 1,200 ha of thinning and ~3,600 ha of prescribed fire

Accepted Article
This article is protected by copyright. All rights reserved.
per year for a period of 20 years followed by ~3,600 ha of prescribed fire per year for the remaining 70 years of the simulation.
3. Moderate pace of restoration: This included 12,000 ha of thinning and 7,000 ha of prescribed fire per year, such that entire project area was treated in a period of 20 years, followed by 7,000 ha of prescribed fire per year for the remaining 70 years of the simulation.
4. Fast pace of restoration: This was an accelerated pace of restoration that closely matched 4FRI's target restoration pace. It included restoring the same area as the moderate pace of restoration but in half the time; 24,000 ha of thinning and 14,000 ha of prescribed fire per year for the first 10 years followed by 14,000 ha of prescribed fire per year for the remaining 80 years of the simulation.

LANDIS-II Model
The LANDIS-II model is a spatially explicit landscape model that simulates forest succession and disturbances, and tracks tree growth using species-age cohorts rather than individual trees (Scheller et al. 2007). Life history attributes are species-specific model input parameters that inform how trees grow and compete. Species-age cohort biomass is tracked in each active cell across the landscape. To simulate succession, we used LANDIS-II v. 6.2 and the Net Ecosystem Carbon and Nitrogen (NECN) Succession extension (v. 4.2.4; formerly the Century Succession extension) that simulates monthly growth, mortality, and reproduction and also tracks carbon, nitrogen, and water cycling both above-and belowground (Scheller et al. 2011

Accepted Article
This article is protected by copyright. All rights reserved.

GCM selection and climate inputs
We used daily precipitation and temperature data from randomly selected years from 1981 to 2010, from the Parameter-elevation Relationships on Independent Slopes Model (PRISM; PRISM Climate Group, Oregon State University 2016) to run model spin-up. Future climate projections were based on 4 Global Circulation Models (GCMs) that were selected on their correlation with historical climate data and their representation of climate extremes in the future. To do this, we calculated correlation coefficients of monthly PRISM vs. modeled historical  data from each GCM using four climate variables that previous studies have found to be important to ponderosa pine growth and recruitment (Feddema et al. 2013, Kolb et al. 2013, Petrie et al. 2016: winter precipitation, summer average temperature, spring average temperature, and summer relative humidity. Using those same four climate variables plus annual precipitation, we also calculated the difference between the averages from 1981-1999 and 2081-2099. We then choose 4 GCMs based on those that were highly correlated with the historical period and fell into each of these categories: below and above average increases in modeled temperature and below and above average changes in modeled precipitation. Chosen GCMs were INMCM4 (warm/dry) , BCC-CSM1-1 (warmer/wetter), IPSL CM5A LR (hot/drier), MIROC ESM CHEM (hotter/wet) and were obtained from Multivariate Adaptive Constructed Analogs datasets (Appendix S1:Figure S1 and S2; Abatzoglou and Brown 2012). RCP 8.5 models were used in all cases. Monthly precipitation (mm), maximum temperature (°C), minimum temperature (°C) and the standard deviation of each of those variables were input into the NECN succession extension. Daily precipitation (mm), average temperature (°C), wind speed (km h-1 ), wind direction (degrees), and relative humidity (%) were used to calculate a fire weather index, fine fuel moisture code, and buildup index for use in the Dynamic Fire extension. Wind speed and wind direction were not expected to change with future climate change (Eichelberger et al. 2008), so historical daily

Accepted Article
This article is protected by copyright. All rights reserved.
wind direction and wind speed (Western Regional Climate Center 2018) were randomized and used for future climates. Future fire weather was input as different files for each decade so that fire weather would change into the future by decade.

Initial communities
LANDIS-II requires an initial communities input file that specifies the species and age cohorts for every cell on the landscape. We included ten species in our model: ponderosa pine (Pinus ponderosa), aspen (Populus tremuloides), pinyon pine (Pinus edulis), juniper (we used 4 similar species in the genus Juniperus that all occur in the region), Gambel oak (Quercus gambelii), white fir (Abies concolor), subalpine fir (Abies lasiocarpa), Engelmann spruce (Picea engelmannii), blue spruce (Picea pungens), and Douglas fir (Pseudotsuga menziesii). We created a raster layer with 100-m grid cells and used USFS inventory data from the 4FRI project (USFS unpublished data) to assign an initial community to each cell.
To do this, we developed diameter to age regressions for each of our species, based on USFS Forest Inventory and Analysis (FIA) data (U.S. Department of Agriculture, Forest Service 2017). We used the regressions to convert tree-level 4FRI inventory data to ages and binned ages into 10-year bins. Unique species-age cohort combinations were created and assigned to every cell on the landscape. Our study area was divided into 14 ecoregions based on a combination of soil texture/depth information derived from U.S. Forest Service soil surveys (U.S. Department of Agriculture, Forest Service 1991, 1995) and four precipitation classes (PRISM Climate Group, Oregon State University 2016; Appendix S1: Figure S3).

NECN Succession
The NECN Succession extension simulates regeneration, growth, and mortality of trees, wood and litter decomposition, soil accumulation and decomposition, and available soil water (Scheller et al. 2011(Scheller et al. , 2012. NECN tracks all pools and fluxes of carbon and nitrogen as well as the biomass in each species/age cohort in all cells. NECN uses species and

Accepted Article
This article is protected by copyright. All rights reserved. ecoregion parameters to define how cohort growth is affected by soil and climate (DataS1).
Maximum biomass and maximum aboveground net primary productivity (ANPP) were defined for each species and each ecoregion. In addition to species-specific parameters which outline many physiological characteristics (DataS1), NECN succession further organizes species into functional groups. This classification helps dictate tree growth responses to environmental conditions, with many parameters inherited from the Century model (Parton et al. 1993). Probability of establishment, which dictates the likelihood of successful regeneration, is the product of species-specific traits and environmental and edaphic conditions. As cohorts grow, they compete for resources within each cell and can disperse across grid cells. Cohorts die when they reach species longevity values or are killed by disturbance (i.e. harvest or fire). Some species, functional group, and NECN parameters used in this study were taken from similar studies that used the same or similar species (Scheller et al. 2011, Loudermilk et al. 2013, Gustafson et al. 2015, Dymond et al. 2016, Hurteau et al. 2016).

Wildfire simulation
We used the Dynamic Fuel extension to assign fuel types to the landscape based on dominant species/age structure in each 100-m grid cell. Fuel types influence the spread and severity of wildfire on the landscape. Fuel types can change with every time step and are reassigned if species composition and ages change (via growth or disturbance) in a cell from the previous time step. We used 12 pre-disturbance fuel types based on species and ages that were expected to burn similarly (DataS1). We also used 5 post-disturbance fuel types that had increased crown base height, and thus lower fire severity, because we expected that treatment would alter the vertical structure of the trees on the landscape (Graham et al. 1999, Jain et al. 2007

Accepted Article
This article is protected by copyright. All rights reserved.
We used the Dynamic Fire extension to simulate wildfire on the landscape. We divided the landscape into three fire regions based on the probability of fire ignitions. Fire ignition probability was determined by reproducing an analysis from Dickson et al. (2006) using topographic roughness, presence of ponderosa pine, and precipitation as independent variables and historical fires greater than 1 ha (Short 2017) as the dependent variable. Fire size distribution, seasonal foliar moisture content (FMC), and number of ignitions were input for each fire region (DataS1). Fire size distribution and number of ignitions for each fire region were obtained from Short (2017). FMC was calculated for each season (spring/summer/fall) using the latitude/longitude and average elevation of each ecoregion (Forestry Canada 1992), and Julian dates for the beginning and ends of each season based on local fire staff knowledge and fire database information (Short 2017). Fire weather index was also calculated at a daily time scale for the entirety of the simulation, using daily climate projections. Fire weather index was binned into severity classes (1-5), to provide specific weather information for future simulated fires. When running LANDIS-II, fire occurrence is determined based on the number of ignitions for each fire region and the probability of fire initiation. We used the number of fires > 1 ha from the fire database (Short 2017) as the number of ignitions in each fire region and the probability of fire initiation was determined by the dominant vegetation type. Fire spread in LANDIS-II is determined by fire weather, topography, and fuel conditions (Sturtevant et al. 2009). Once a fire is initiated, the fire spreads until it reaches the pre-determined fire size as selected from a lognormal distribution from user-input µ (mean fire size) and σ 2 (standard deviation of fire size). Future fire weather was randomly selected from a distribution of fire weather days from each decadal fire weather file. Likelihood of cohort mortality is based on the severity of the fire and the age class of the cohort; this is parameterized in the Dynamic Fire input text file (DataS1). For each age cohort that is killed in a fire, a user-defined proportion remains on site as dead wood

Accepted Article
This article is protected by copyright. All rights reserved. and a proportion is volatilized. These proportions for wood and litter vary by fire severity and are parameterized in the fire reduction table in the NECN Succession text file (DataS1).

Management
We used the Biomass Harvest extension to simulate the 4FRI thinning prescriptions and prescribed fire prescriptions. Biomass Harvest requires the landscape to be divided into 'management units', within which targets are defined for the amount of area treated. We divided the landscape into 6 management units based on the prescriptions that were to be placed there. We developed 15 LANDIS-II prescriptions based on the prescriptions used by USFS (unpublished data) in the 4FRI restoration project: grassland prescriptions that removed all trees; savanna prescriptions that removed most trees; nine prescriptions that removed ponderosa pine at various levels; and four prescribed fire prescriptions that included initial prescribed fire and maintenance prescribed fire for grasslands, savannas, and the remaining project area (DataS1). To realistically simulate current and anticipated management across the 4FRI landscape, we designed our management prescriptions to mimic the thinning practices used by managers, including age and diameter targets for individual species. Because LANDIS-II is a cohort-based model and cannot simulate individual tree removal, we then used expected pre-and post-treatment condition 4FRI data from the USFS (unpublished data) to estimate the biomass removed from each cohort in each prescription in LANDIS-II (Table 2). The Dynamic Fire extension is best suited for simulating wildfire effects, therefore we used the Biomass Harvest extension to simulate prescribed fires. As is designed for 4FRI, we implemented prescribed maintenance fires to continue after initial thinning and prescribed fire (10 or 20 years, depending on scenario).

Accepted Article
This article is protected by copyright. All rights reserved.

Model calibration/validation
LANDIS-II runs a model spin-up for a length of time that is equal to the age of the oldest cohort, placing cohorts on the landscape in each time step corresponding to their current age. The oldest cohort in our initial communities layer was estimated to be 800 years old, resulting in a model spin-up time period of 800 years. After model spin-up, at year 0, we compared the biomass in each ecoregion from LANDIS-II to the biomass at each ecoregion (Appendix S1: Figure S4) from USFS data. LANDIS-II biomass values were less variable than inventory data but the median LANDIS-II biomass values in each ecoregion matched inventory data well, reliably falling within the interquartile range. As is common with LANDIS-II, initial biomass values were slightly overestimated because spin-up does not include disturbances . The average biomass from the inventory data (11,396 g/m 2 ) and the average biomass from LANDIS-II calibrations at year 0 (14,931 g/m 2 ) were within expected bounds (Appendix S1: Figure S4).
To calibrate the Dynamic Fire extension, we ran the simulations for 30 years with historical climate input for growth and fire weather and compared the mean fire size, standard deviation of the fire size, fire return interval, and area burned per year for each fire region to fires from a fire database (Short 2017; Appendix S1: Table S1). We prioritized area burned per year over the other metrics in the calibration to ensure that it was similar to historical data. The average area burned across the project area in the historic database was 3,387 ha/yr and in the LANDIS-II calibrations was 3,651 ha/yr (Appendix S1: Table S1). To calibrate the Biomass Harvest extension, we compared the biomass removed overall and from each prescription to USFS data (Appendix S1: Figure S5).

Accepted Article
This article is protected by copyright. All rights reserved.

Simulations and analysis
We ran 10 replicate simulations for each scenario and each GCM from 2010 -2099 using an annual time step for all extensions. We conducted sensitivity analyses by evaluating aboveground biomass output averaged for 2, 5 and 10 replicates and found little difference between averages of 5 and 10 replicates; thus, we ran 10 replicates of each scenario. For each scenario/GCM, we evaluated total ecosystem carbon (TEC), individual pools of carbon (i.e. live carbon, dead carbon, and soil carbon), net ecosystem production (NEP), wildfire emissions, wildfire severity and net ecosystem carbon balance (NECB) at each time step and averaged across replicates. NEP is gross primary productivity minus ecosystem respiration with positive values indicating the landscape is acting as a sink for carbon and negative values indicating a source of carbon to the atmosphere. NECB is NEP minus carbon lost from fire and prescribed fire, with positive values indicating the landscape is acting as a sink for carbon and negative values indicating a source of carbon to the atmosphere. NECB also often includes carbon lost from thinning but, in this study, a separate analysis outside of LANDIS-II was used to determine the fate of carbon post-thinning and thus, the carbon storage/emissions from thinning (see below) so thinning was not included in our NECB estimates. TEC values reported are those output from LANDIS-II and do not include carbon that remains sequestered in harvested material (see below).
Fire severity class is created for each cell in each time step and ranges from 1 -5 with 5 representing the most severe wildfire. We refer to fire severity classes as: 1= low; 2=lowmoderate; 3=moderate; 4=moderate-high; and 5=high. Fire severity was evaluated across time steps for each cell and compared across scenarios. Fire sizes were summed for each time step and averaged across replicates to compare among scenarios. We used ArcMap 10

Accepted Article
This article is protected by copyright. All rights reserved.

Harvested-Products Analysis
In addition to the models in LANDIS-II, we developed and applied a simple flow analysis to determine the carbon storage and emissions during and following harvest through local product creation. Although not a comprehensive Life Cycle Assessment (Oneil et al. 2017), we included: emissions from harvesting equipment, emissions from transportation of wood (stem wood and biomass) to wood products facilities, emissions from biomass products that decompose in our analysis timeframe, and storage from wood products. To calculate carbon emissions from harvesting activities, we followed the methodology outlined by Markewitz (2006) combined with local operator and machine productivity rates for restoration thinning. To determine carbon emissions from transportation, we used local knowledge of individual truck load weights (stem wood and residual biomass separately), average distance traveled to product facilities (wood products and chipped wood separately), average miles per gallon (mpg) of transport vehicles, and U.S. Energy Information Administration's carbon dioxide emissions coefficients (U.S. Energy Information Administration 2016). To account for carbon storage and emissions from resulting forest products and waste associated with this production, we made several key assumptions: 1) the approximate amounts of wood products were as follows: pallets -55%; single family homes (post-1980) -20%; manufacturing -5%; railroad ties -15%; and furniture -5% (pers comm Stephen Horner, Campbell Global) and each product had a carbon half-life as indicated in Skog and Nicholson (1998) , 2) chipped wood/biomass used for soil amendments was an emission for the length of analysis assuming decomposition rates for that material of ~10%/year (Harmon et al. 1996); and 3) any residual biomass material created from processing wood products was an emission for the length of our analysis. We used local industry knowledge of proportions of wood versus soil amendment products, and the

Accepted Article
This article is protected by copyright. All rights reserved.
proportions of waste created when making wood products. Additional information can be found in Appendix S2.

Results
Despite initial decreases in carbon in the first two decades due to accelerated harvest and prescribed fire, the moderate-and fast-4FRI scenarios resulted in greater carbon storage by the end of the century than the status-quo and no-harvest scenarios and that pattern remained consistent among the climate models ( Figure 2). Depending on the climate model, the overall increases in TEC for the fast-4FRI scenario was 9%-16% higher than the noharvest scenario, equating to an increase of 6.3 million -12.7 million metric tons of carbon across the 4FRI project area. Climate models that predicted the largest temperature increases by late century showed greater decreases in TEC in all restoration scenarios by the end of the century (Figure 2 c,d).
Given that the relative changes in carbon pools and fluxes in the four climate models were similar, we present the remaining results only for the hot/drier climate model here because it resulted in intermediate carbon estimates. Results from additional models can be found in the supplementary materials. Under the hot/drier climate model, live carbon pools accounted for 52%, soil carbon pools accounted for about 40%, and dead carbon pools accounted for about 8% of the overall TEC differences between the fast-4FRI and the noharvest scenario at the end of the simulation (Table 3; Figure 3; Appendix S1: Figure S6). By 2099, the no-harvest scenario had 1.5 times more wildfire emissions and 33% less production than the fast-4FRI scenario (Figure 4; Appendix S1: Figure S7, S8). In terms of carbon fluxes, ~70% of the TEC increase in the fast-4FRI relative to the no-harvest scenario was explained by decreased wildfire emissions while greater production explained ~30%. In the fast-4FRI scenario, ~22% of the area burned was burned at a high severity (moderate-high or high) compared to ~66% in the no-harvest scenario ( Figure 5; Appendix S1: Figure S9). There were

Accepted Article
This article is protected by copyright. All rights reserved. no significant differences in annual area burned among the scenarios (Appendix S1: Figure   S10). NECB exhibited high interannual variability but was consistently higher in the fast-4FRI scenario than in the no-harvest scenario for approximately 50 years after initial treatment, indicating that an accelerated harvest schedule can increase carbon sink strength ( Figure 6; Appendix S1: Figure S11). Annual NECB values between the two scenarios converged in the last 40 years when perhaps warmer temperatures and drier conditions in the hot/drier climate model (Appendix S1: Figure S1) neutralized the effects of restoration.
The harvested-products analysis showed that ~5% of the harvested carbon remained sequestered in wood products at the end of our study, adding an additional ~53 g C/m 2 or 207,000 metric tons of carbon in the fast 4FRI scenario. In contrast, approximately 137,000 metric tons of carbon was emitted by the harvesting equipment and vehicles transporting the material off site. Thus, there was only a net gain of 70,000 metric tons of carbon from the harvested material, which represents only a fraction (~0.5%) of the carbon gains modeled in LANDIS-II from the fast-4FRI scenario (Table 3). Because these gains were so small, we did not include them in the overall total carbon storage for the landscape.

Discussion
Our results demonstrate that large-scale forest restoration can increase the potential for carbon storage and stability and those benefits could increase as the pace of restoration accelerates. The fast-4FRI scenario could increase total carbon storage by ~12.7 million metric tons over the no-harvest scenario in the hot/drier climate model, equating to the removal of the annual emissions from 110,000 passenger vehicles or the electricity consumed in over 90,000 homes per year until the end of the century. Increases in ecosystem carbon were greater when the 4FRI program was completed in 10 years versus 20.

Accepted Article
This article is protected by copyright. All rights reserved.
Potential gains in TEC in the fast-4FRI scenario compared to the no-harvest scenario were the result of greater productivity, explaining ~30% of the increases, and a reduction in wildfire emissions, explaining ~70% of the increases. Production increases were apparent in the 4FRI scenarios subsequent to initial thinning and prescribed fire (10 or 20 years, depending on the scenario) and growth continued to remain higher in the 4FRI scenarios for ~50 years. Additionally, smaller proportions of the landscape burned at high severity as the pace of restoration increased, resulting in lower tree mortality ( Figure 5). In the case of ponderosa pine, for example, the low and low-moderate fire severity categories that occurred more frequently in the 4FRI scenarios killed only trees <10 years old, whereas the moderatehigh and high severity categories, that occurred more frequently in the no-harvest scenario, killed all trees <80 years old and <340 years old (nearly all the trees in most sites), respectively.
Climate change effects on forest carbon were most apparent in the last half of the century and these effects were more pronounced in the climate models that predicted hotter temperatures (Figure 2). In the hot/drier and hotter/wet climate models, ecosystem carbon decreased in the last half of the century likely because higher temperatures and/or less precipitation led to lower rates of growth and establishment (Appendix S1: Figure S12) that could not keep pace with increased mortality from wildfire (Appendix S1: Figure S13), maintenance prescribed fire, and tree senescence. This effect was greatest in the hotter/wet climate model that predicted a 9.7°C increase in summer temperature by the end of the century, suggesting that a ~10°C increase in summer temperature could represent a threshold where the effects from restoration are diminished or eliminated. This is demonstrated in Figure 4b, where the slopes of the lines representing cumulative NEP for 4FRI scenarios and the no-harvest scenario become similar in the last 30 years of the century, and in Figure 6, where annual NECB values become similar and overall negative during the same time frame.

Accepted Article
This article is protected by copyright. All rights reserved.
While no other studies have incorporated the level of detail to represent an actual landscape-scale forest restoration project as we did, other studies have found similar results indicating that forest restoration can increase carbon storage and decrease wildfire severity (Loudermilk et al. 2014, Hurteau et al. 2016, Krofcheck et al. 2017, Hurteau 2017, Liang et al. 2018. Because 4FRI prescriptions were developed to achieve heterogeneity in size-class distributions across a large landscape, the prescriptions in this study were more varied than the basic fuels reduction treatments in previous studies (Table   2). Previous studies have removed as little as 30% of live biomass focused on the smaller size classes (Hurteau et al. 2016, Krofcheck et al. 2017, e.g., Hurteau 2017, but in some cases removed up to 95% of post-settlement trees (Flatley and Fulé 2016) or examined different thinning intensities and objectives (Loudermilk et al.,2014).
Still other studies found that, despite restoration efforts, extreme climate change conditions could result in significant biomass losses and deforestation (Tarancón et al. 2014, Flatley and Fulé 2016, Loehman et al. 2018. We also saw decreases in TEC and aboveground biomass in the no-harvest and status-quo scenarios implying that without restoration there is a strong likelihood that forest biomass will decline with climate change due to high severity fires and low productivity. Additionally, our study found lower recruitment across all scenarios with extreme increases in temperature (Appendix S1: Figure   S12), implying that, despite restoration efforts, biomass loss is still possible. We found similar vegetation type changes as Loehman (2018) in a similar time frame but Flatley and Fulé (2016) found that vegetation type changes are possible further into the future.
However, when comparing our study to other studies, even those done with LANDIS-II in the same region, we found widely varying assumptions and inputs and it is unknown what effect these had on the results. Without sensitivity analyses to evaluate the effect of different inputs, it is difficult to compare across studies leaving systematic review or meta-

Accepted Article
This article is protected by copyright. All rights reserved.
analyses difficult (see James et al. 2018). Future research would benefit from more transparency in model assumptions and limitations and in reporting a common set of metrics for all studies.
We conducted a separate analysis evaluating the fate of carbon from the harvested products in this region. LANDIS-II assumes all harvested material is removed from the ecosystem and we believe that assumption is valid in this region with the current wood product types. Our results suggest that there is little net gain in carbon from harvested material because the amount of carbon sequestered in wood products at the end of the study timeframe was roughly equivalent to the emissions from the equipment and transport vehicles used during the harvest. Most of the wood products in this region were pallets (55%) with a carbon half-life of only six years but if the harvested biomass were converted to longerlasting wood products, we would expect that this number would substantially increase (Finkral and Evans 2008). We did not include harvest in our estimates of NECB, but in accounting for the emissions from harvest, we would expect the NECB values for the fast 4FRI landscape would decrease by ~45 g C/m 2 /yr for the first 20 years, including immediate emissions (material burned at processing) and decomposition over time (chipped and residual biomass). This indicates that the 4FRI landscapes were likely a greater source of carbon to the atmosphere during the first 20-30 years than shown in the LANDIS-II results ( Figure 6).

Model limitations
We found that ~40% of the increases in carbon storage occurred in soil carbon, a stable long-term carbon pool that is less vulnerable to losses from forest disturbances (Scheller et al. 2011). In LANDIS-II, thinning, prescribed fire, and wildfire all add biomass to the dead pools of carbon on the landscape. Dead-pool carbon quickly decomposes to become soil organic carbon and remains in soil carbon pools long term. However, there are limitations to the harvesting process used in LANDIS-II because equal amounts of

Accepted Article
This article is protected by copyright. All rights reserved.
aboveground and belowground live carbon are removed from the system. Instead of transferring the belowground live carbon (roots) to the dead pools and eventually soil pools, the model removes the belowground pools from the ecosystem in the same way as aboveground live pools. This leads to fewer inputs into the dead pools and thus, fewer inputs into soil organic carbon, and a conservative estimate of the amount of carbon being stored in the soil. Additionally, black carbon is not accounted for in LANDIS-II so it is not known how this would vary by scenario.
While the pace of restoration influenced wildfire severity in this study, neither restoration or climate change significantly affected area burned per year (Appendix S1: Figure   S10). We believe this is an artefact of how the LANDIS-II Dynamic Fire extension fits wildfire sizes to a lognormal distribution that cannot be altered with changing climatic conditions. This distribution shortens the right-side tail of the fire-size distribution in such a way that large wildfires that have already occurred and larger fires that are predicted to occur (Kitzberger et al. 2017, Keyser andWesterling 2019) are not simulated. However, despite the limitations of LANDIS-II to model it, warming-induced increases in wildfire size will likely occur (Westerling 2016) and result in higher carbon emissions. Restoration treatments, including prescribed fire, will remain important for reducing fire severity, but may not consistently reduce wildfire size (Cochrane et al. 2012, Price et al. 2015. The purpose of 4FRI is to reestablish forest structure to improve forest health and the resiliency of forests to climate change, fire and other disturbances (U.S. Department of Agriculture, Forest Service 2013). Other studies have found that large-scale forest restoration projects are likely to have co-benefits that include maintaining biodiversity (Hurteau et al. 2014), improving water supply (Robles et al. 2014), improving insect resistance (Fettig et al. 2013, Kolb et al. 2016, and increasing soil function (Sánchez Meador et al. 2017). This study illustrates that large-scale forest restorations are also likely to have potential co-benefits of

Accepted Article
This article is protected by copyright. All rights reserved.
increased carbon storage, but climate change could diminish the effects of restoration if higher predicted temperatures are realized. We conclude that the 4FRI project will increase potential carbon stores in ponderosa pine forests by reducing tree mortality from severe wildfires and sustaining productivity under adverse climatic conditions. Our results also suggest that carbon increases can result from initial thinning followed by maintenance fires.
Prescribed fire is less expensive and easier to implement than thinning and can provide benefits not achieved with thinning alone (North et al. 2012, Hood et al. 2015 and both are required to achieve the desired results. 4FRI is only one of 16 western US landscapes in which the CFLR program has funded accelerated restoration of dry forests across 1.1 million hectares (U.S. Department of Agriculture, Forest Service 2018). However, the scale of the fire risk is much larger -of the estimated 25 million ha of frequent-fire conifer forests on federal lands in the western US (Ager et al. 2013), fire regimes have been significantly altered and restoration is needed on more than 11 million ha (Menakis et al. 2004). All 11 million ha are likely not suitable for mechanical thinning. However, even if only half are suitable and accelerated restoration in those landscapes has similar effects on the carbon cycle as 4FRI, substantial increases in carbon sequestration and stabilization at the subcontinental scale could be realized. While accelerated forest restoration on its own is unlikely to be a sufficient long-term response to climate change, carbon stabilization from accelerated restoration would buy time to retain forest cover, understand climate effects and develop management strategies that reduce the loss of forest cover from climate-induced impacts.
Short of enhancing the current CFLR program or developing new programs that address the urgent need to accelerate restoration, this study and others project potentially dramatic changes in forest carbon pools and forest cover in the Western US.

Accepted Article
This article is protected by copyright. All rights reserved.

Acknowledgments
We would like to thank the Anne Ray Foundation for funding this project. We would like to acknowledge two anonymous reviewers for helping to improve this manuscript. We would like to thank Mark Nigrelli and the U.S. Forest Service for generously providing unpublished data related to 4FRI and this project. We would like to thank Stephen Horner for providing estimates of the wood products coming from harvest materials in the region. We      carbon (g C/m 2 ) includes aboveground (wood and leaves) and belowground (coarse and fine roots) live carbon; B -Total soil carbon (g C/m 2 ) includes all soil organic matter pools; C -Total dead carbon (g C/m 2 ) includes dead wood, dead leaves, and dead roots. Shaded areas