Carbon dynamics of Western North American boreal forests in response to stand-replacing disturbances

North American boreal forests are known to be an important carbon pool in boreal ecosystems, but have experienced extensive tree mortality and carbon loss due to multiple agents of stand-replacing disturbances in recent decades. However, the impacts of these stand-replacing disturbances on forest dynamics are still unknown. We used a recently developed remote-sensing based stand-replacing disturbance product, coupled with aboveground biomass (AGB), gross primary productivity (GPP) and leaf area index (LAI) datasets to estimate the impacts of stand-replacing disturbances (e.g., fires, logging and insect outbreaks) on the carbon balance of western North American boreal forests during 2000 – 2012. Our results showed that fire, logging and insect outbreaks resulted in AGB losses of 23.4, 16.6, and 4.7 Tg/yr, respectively. In the post-disturbance periods, AGB did not recover to its pre-disturbed levels in the 10th year, which is longer than the recovery time of GPP and LAI. Furthermore, the losses of AGB, GPP and LAI in fire events were the dominant factors for forest recovery after stand-replacing fire. Vapor Pressure Deficit (VPD), soil clay content, temperature and precipitation were the important factors for forest recovery after stand-replacing insect outbreaks and stand-replacing logging. When removing the impact of environmental factors, our results showed a smaller magnitude of AGB, GPP and LAI loss relative to the results including these factors, although similar recovery trajectories were observed among the two results. The results have important implications for understanding the effects of stand-replacing disturbances on the carbon dynamics of boreal forests, which is required to adopt effective forest management strategies after disturbance.

North American boreal forests are known to be an important carbon pool in boreal ecosystems, but have experienced extensive tree mortality and carbon loss due to multiple agents of stand-replacing disturbances in recent decades.However, the impacts of these stand-replacing disturbances on forest dynamics are still unknown.We used a recently developed remote-sensing based stand-replacing disturbance product, coupled with aboveground biomass (AGB), gross primary productivity (GPP) and leaf area index (LAI) datasets to estimate the impacts of stand-replacing disturbances (e.g., fires, logging and insect outbreaks) on the carbon balance of western North American boreal forests during 2000-2012.Our results showed that fire, logging and insect outbreaks resulted in AGB losses of 23.4, 16.6, and 4.7 Tg/yr, respectively.In the post-disturbance periods, AGB did not recover to its pre-disturbed levels in the 10th year, which is longer than the recovery time of GPP and LAI.Furthermore, the losses of AGB, GPP and LAI in fire events were the dominant factors for forest recovery after stand-replacing fire.Vapor Pressure Deficit (VPD), soil clay content, temperature and precipitation were the important factors for forest recovery after stand-replacing insect outbreaks and stand-replacing logging.When removing the impact of environmental factors, our results showed a smaller magnitude of AGB, GPP and LAI loss relative to the results including these factors, although similar recovery trajectories were observed among the two results.The results have important implications for understanding the effects of stand-replacing disturbances on the carbon dynamics of boreal forests, which is required to adopt effective forest management strategies after disturbance.

Introduction
North American boreal forests play a vital role in the global carbon cycle (Dieleman et al. 2020).Yet, the carbon balance of these forests is threatened by increasingly frequent stand-replacing disturbances, including forest fires (Mack et al. 2021), insect outbreaks (Bright et al. 2020), and timber harvests (Masek et al. 2011).Thus, understanding the carbon dynamics of North American boreal forests in response to standreplacing disturbances is crucial for comprehending the mechanisms that control the global carbon cycle.
Stand-replacing disturbances are defined as events that result in the complete death of the living tree biomass, which means the complete removal of tree cover or a change in land cover from forest to non-forest (Hansen et al. 2013;Pugh et al. 2019;Zhang et al. 2022).Fire is the dominant stand-replacing disturbance in North American boreal forests (Bond-Lamberty et al. 2007;Johnstone et al. 2010).It impacts carbon dynamics through the combustion of aboveground biomass, since regional tree species are susceptible to crown fires (Rogers et al. 2015).In recent decades, North American boreal forests have experienced an increasing trend in burned area, with an increase of 3.40 × 10 4 ha per year since 1959, especially in western Canada (Hanes et al. 2019).Insect outbreaks have been increasing severe in recent years, affecting up to 1.76 × 10 7 ha of forests (Natural Resources Canada 2017).For example, the mountain pine bark beetle has caused widespread tree mortality, particularly in British Columbia (Kurz et al. 2008).Also, timber harvest activities occupied approximately one quarter of Canada's boreal forest (Schroeder et al. 2011).These three stand-replacing disturbance agents (e.g., fires, logging and insect outbreaks) impact the carbon balance of North American boreal forests together (Turner et al. 2015).Therefore, clarifying the contribution of each agent to North American boreal forest disturbances is important for understanding the carbon dynamics of North American boreal forests.
The carbon loss resulting from various disturbance agents in North American boreal forests has been estimated using field surveys (Amiro et al. 2010), remote sensing (Wang et al. 2021;Wulder et al. 2020), and ecological models (Kurz et al. 2008;McGuire et al. 2010).For example, based on field surveys, the carbon loss in North American forests following disturbances ranged from 200 g C m − 2 yr − 1 to 1270 g C m − 2 yr − 1 (Amiro et al. 2010).The cumulative effect of the beetle outbreak resulted in more than 540 Tg of carbon loss during 2000-2020 (Kurz et al. 2008).A recent remote-sensed study estimated the carbon loss of 800 Tg and 74 Tg from fires and timber harvest, respectively (Wang et al. 2021).These disturbances have a large impact on the carbon balance of North American boreal forests.Compared to field surveys, remote sensing and ecosystem model-based methods are more suitable for monitoring carbon loss resulting from large-scale forest disturbances.However, the existing remote sensing-based methods are often constrained by the inclusion of disturbance datasets from various sources, which encompass both stand-replacing and non-stand-replacing disturbances, resulting in the inability to make an inter-comparison between different studies.Especially, the lack of available disturbance history products, that represent multiple stand-replacing disturbances simultaneously, limits the understanding of the response of forest carbon dynamics to multiple stand-replacing disturbances (e.g., fires, logging, and insect outbreaks) (National Forestry Database 2018;Senf et al. 2017;Wees et al. 2021;White et al. 2017).
A recently released long-term forest disturbance product potentially provides new insights into understanding the response of North American boreal forests to stand-replacing disturbances.This product provides annual data on individual stand-replacing disturbances (i.e., fires, logging, and insect outbreaks) at 30 m resolution from 1987 to 2012, respectively (Zhang et al. 2022).This dataset provides unprecedented information to assess recent trends in disturbance and recovery.
Furthermore, knowledge of post-disturbance recovery is particularly important, considering the impact of the increasing frequency of disturbances on forest carbon dynamics (Goetz et al. 2012;Turner 2010).
Numerous indicators associated with forest dynamics have been used to characterize recovery after disturbances (Chu and Guo 2013), such as the Normalized Difference Vegetation Index (NDVI) (Wang et al. 2022), Enhanced Vegetation Index (EVI) (Jin et al. 2012) and net primary productivity (NPP) (Amiro et al. 2000;Bond-Lamberty et al. 2004).However, these indicators mainly represent the recovery of vegetation photosynthetic activity after disturbances, rather than the recovery of aboveground biomass.Note that the recovery of photosynthetic activity does not mean the recovery of aboveground biomass (Fan et al. 2022;Myers-Smith et al. 2020).Moreover, most studies on post-disturbance forest recovery rely on a single indicator, limiting the ability to capture the multidimensional characteristics (e.g., vegetation cover, carbon storage, biomass) of forest recovery after disturbances.A recent annual aboveground biomass (AGB) product over the North American boreal forests provides an opportunity to estimate biomass loss and recovery after disturbances (Wang et al. 2021).Also, the interannual changes of LAI (Zhao et al. 2021) and gross primary productivity (GPP) (Goulden et al. 2011) characterize the degree of vegetation cover and the amount of carbon fixed by forests, respectively.Thus, the combined information from LAI, GPP, and AGB can be used to characterize the biomass loss and the multidimensional characteristics (e.g., vegetation cover, carbon storage, biomass) of forest recovery after disturbances.
The purpose of this study is to quantitatively estimate the impact of stand-replacing disturbances (i.e., fires, logging and insect outbreaks) on carbon dynamics in western North American boreal forests from 2000 to 2012, using time series of remotely sensed LAI, GPP and AGB.This work was carried out by (1) assessing the effects of forest stand-replacing disturbances (e.g., fires, logging, and insect outbreaks) on AGB; (2) investigating the recovery trajectories of LAI, GPP and AGB; (3) analyzing the importance of the environmental factors (e.g., climatic factors, soil properties, radiation) that affect forest recovery in postdisturbance periods.

Study region
The study region is located in the western North American boreal forests, mainly focusing on the Arctic Boreal Vulnerability Experiment (ABoVE) core domain, (100 degreesW-168 degreesW, 52 degreesN-74 degreesN) (Loboda et al. 2019) (Fig. 1), which covers Alaska, Northwest Territories, Yukon, British Columbia, Alberta, Saskatchewan and Nunavut.The major tree species consist of spruce, poplar, pine, fir, birch, hemlock and aspen (United States Department of Agriculture 2009).

Forest disturbance dataset
The historical forest disturbance dataset from Zhang et al. (2022) was used to identify stand-replacing disturbances.This dataset provides annual maps of three stand-replacing disturbance agents (i.e., fires, logging and insect outbreaks) with 30 m spatial resolution over the study region from 1987 to 2012 (Zhang et al. 2022).Stand-replacing disturbances were defined as events that caused a land cover change from forest to non-forest, i.e. stand-replacing fire, stand-replacing logging and stand-replacing insect outbreaks (Zhang et al. 2022).The disturbance agents were derived using the continuous change detection and classification algorithm (CCDC) based on the Landsat time series of Thematic Mapper (TM) and Enhanced Thematic Mapper (ETM+) observations.The overall accuracy of this dataset was 96.7%.

Vegetation parameters
The LAI dataset was obtained from MODIS observations (MOD15-A2HGF) (Myneni et al. 2015), which is available from January onwards.The temporal and spatial resolutions of LAI are 8 days and L. Yu et al. m, respectively.The LAI time series was quality-screened using the quality assessment layers included in the product, and only LAI observations with good quality were used in this study.This LAI dataset was combined into yearly composites by calculating the mean of all the 8day LAI observations with good quality within each year.The yearly composite values were then aggregated to 0.05 degree spatial resolution using a simple spatial averaging method.
GPP was used to define the recovery of vegetation productivity in response to stand-replacing disturbances.The VPM GPP V20 dataset provided the yearly GPP product at 0.05 degree spatial resolution from 2000 to 2016 (Zhang et al. 2017), which was retrieved by simulations of the improved light use efficiency theory, MODIS satellite data and NCEP Reanalysis II climate data.Additionally, the low-quality or missing observations were filled using the most advanced gap-filling and smoothing algorithms.This dataset has a satisfactory performance (R 2 = 0.74, RMSE = 2.08 g C m − 2 day − 1 ) when assessed against in situ GPP data from flux towers (Zhang et al. 2017).
The AGB dataset developed by Wang et al. (2021) provides annual AGB density for live woody species for the western North American boreal forests at 30 m spatial resolution over the time period 1984-2014.The annual AGB density was predicted using Machines machine learning algorithm by applying the GLAS LiDAR data to allometric equations.The accuracy of the AGB model was evaluated using an independent validation dataset with the coefficient of determination (R 2 ) and Root Mean Square Error (RMSE) was 0.8 and 20.77 Mg ha − 1 , respectively.This AGB dataset was spatially aggregated to 0.05 degree resolution using a simple averaging method.
Accounting for the availability of LAI, GPP, AGB and forest disturbance data, we mainly focused on the study period from 2000 to 2012.

Environmental parameters
Six environmental parameters were used to evaluate the impact of environmental factors on forest recovery in the post-disturbance period (Table S1).The air temperature at 2 m (Temp), precipitation (Prec), and soil moisture (SM) at the depth of 0-7 cm were obtained from the ERA5 monthly average dataset, which has a spatial resolution of 0.25 degree (Hersbach et al. 2020).Vapor Pressure Deficit (VPD) was calculated using the dew point temperature (Td), Temp and the surface pressure according to the method described in Yuan et al. (2019).Td and the surface pressure were also obtained from the ERA5 monthly average dataset (Hersbach et al. 2020).Photosynthetically Active Radiation (PAR) at 1 degree spatial resolution was derived from the CER-ES_SYN1deg monthly product, version Ed4A (Loeb et al. 2013).Soil clay content at the depth of 0-5 cm was provided by SoilGrids 2.0 at a spatial resolution of 250 m (Poggio et al. 2021).
The monthly Temp, Prec, SM, VPD and PAR data were combined into yearly products using the simple average method, and all environmental parameters were resampled to 0.05 degree spatial resolution using the nearest neighbor interpolation method.

The impact of stand-replacing disturbances on AGB
To investigate the impact of stand-replacing disturbances on AGB, we used the method proposed by Qin et al. (2021) to calculate the AGB loss caused by fires, logging, and insect outbreaks, respectively.Taking the calculation of AGB loss caused by fire as an example.As shown in Eq.
(1), first, using the 30 m forest disturbance dataset (Zhang et al. 2022), we calculated the gross fire area (A fire ) from 2001 to 2012 at pixel scale.Second, we multiplied the gross fire area (A fire ) by the 30 m AGB density in 2000 (AGB 2000 ) (Wang et al. 2021) to calculate the AGB loss at 30 m grid cell.Third, The AGB losses were summed up as the total AGB loss (AGBloss fire ).Correspondingly, the AGB loss caused by logging and insect outbreaks were calculated using Eq. ( 2) and Eq. ( 3), respectively. ) (2) where AGBloss fire , AGBloss logging , and AGBloss insect represent AGB loss caused by fires, logging and insect outbreaks, respectively.A fire , A logging and A insect represent the gross area affected by fires, logging and insect outbreaks form 2001-2012, respectively.AGB 2000 represents the AGB density in 2000.

Forest recovery post stand-replacing disturbances
Forest recovery post stand-replacing disturbances, as measured by LAI, GPP and AGB, was studied by analyzing pixels that were disturbed only once during 2000-2012.Taking the pixels that were disturbed only once by fire as an example, the pixels were determined using the area fraction of fire within a 0.05 degree grid cell using the dataset from Zhang et al. (2022): first, the annual fire area fraction was calculated at the resolution of 0.05 degree as the proportion of the summed areas of forest fires within each 0.05 grid cell.Second, the pixels with a fire area fraction greater than or equal to 10% in the fire year, and less than 5% in the remaining years corresponded to the pixels that were disturbed only once by fire.Correspondingly, the pixels that were disturbed only once by logging and insect outbreaks were calculated using the same methods.Based on the pixels that were disturbed only once, forest recovery after fire and logging was analyzed by building the time series of LAI, GPP and AGB recovery.Taking the recovery of LAI after fire as an example, the LAI time series were first shifted to align with the fire years of all the fire events considered (Fig. 2).Year zero (e.g., "0 year") represents the year of fire, negative values represent the years before fire (e. g., "-1 year" means one year before fire) and positive values represent the years after fire (e.g., "1 year" means the 1st year after fire).
Then, at ith year after fire, the recovery of LAI (Recovery LAI,i ) was calculated as the difference of LAI between ith year in the post-fire period (LAI i ) and its values for the year prior to fire (LAI − 1 ) using Eq. ( 4).Correspondingly, the recovery of GPP (Recovery GPP,i ) and the recovery of AGB (Recovery AGB,i ) were calculated using Eq. ( 5) and Eq. ( 6), respectively.
where Recovery LAI,i , Recovery GPP,i and Recovery AGB,i represent the recovery of LAI, GPP and AGB at ith year, respectively.LAI i , GPP i , and AGB i are the LAI, GPP and AGB from the year before disturbed year to the 12th year after the disturbance, respectively.i ranges from − 1 to 12. AGB − 1 , GPP − 1 , and LAI − 1 are the AGB, GPP and LAI in the year prior to disturbance, respectively.Correspondingly, the recovery of LAI, GPP and AGB after logging were calculated using the above methods.
Since insect outbreaks mainly occurred in 2005-2009, the length of the recovery time series built using the above method was limited.To better understand the response of carbon dynamics to insect outbreaks, we took the first year (the year 2000) of our study period as the benchmark year.Taking LAI as an example, LAI recovery (Recovery LAI,j ) after insect outbreaks was calculated as the difference of LAI between jth year (LAI j ) and the year 2000 (LAI 2000 ) using Eq. ( 7).Correspondingly, GPP and AGB recovery (Recovery GPP,j and Recovery AGB,j ) after insect outbreaks were calculated using Eq. ( 8) and Eq. ( 9): Recovery GPP,j = GPP j − GPP 2000 (8) where Recovery LAI,j , Recovery GPP,j and Recovery AGB,j represent the recovery of LAI, GPP and AGB at jth year, respectively.LAI j , GPP j , AGB j are the LAI, GPP and AGB from 2000 to 2012, respectively.j ranges from 2000 to 2012.LAI 2000 , GPP 2000 and AGB 2000 are the LAI, GPP and AGB in 2000, respectively.

The effects of environmental factors on forest recovery
Environmental factors (e.g., Temp, Prec and SM) may affect forest recovery after disturbances (Danneyrolles et al. 2023;Dobor et al. 2018).The method proposed by Cuevas-González et al. (2009) was used to remove the influence of environment factors on the recovery of LAI, GPP and AGB after stand-replacing disturbances.This method assumes that the environmental factors have the same effects on the forest changes (e.g., the changes in LAI, GPP and AGB) over the disturbed areas and those paired adjacent undisturbed areas.Thus, the influence of environmental factors on forest changes can be removed by calculating the difference of LAI, GPP and AGB between the disturbed areas and those paired adjacent undisturbed areas.
We took LAI recovery after stand-replacing fire as an example to illustrate the method proposed by Cuevas-González et al. (2009).As shown in Fig. 3, the central pixel (P5) was disturbed by stand-replacing fire, defined as the disturbed pixel.The adjacent undisturbed pixels (P1, P3, P7 and P9) were selected in the 3 × 3 pixels around the disturbed pixel.Then, we calculated the differences between the disturbed LAI values (P5) and the averaged LAI values over the undisturbed pixels (P1, P3, P7 and P9) to represent the LAI recovery without the effects of environmental factors.Note that if there are no adjacent undisturbed pixels in the 3 × 3 pixels, adjacent undisturbed pixels will be found over a wider area, e.g., 5 × 5 or 7 × 7 pixels around the central disturbed pixel.
Following the method above to remove the effect of environmental factors, we calculated the forest recovery (indicated by LAI, GPP and AGB) in response to the stand-replacing fire, logging and insect outbreaks, respectively.
Furthermore, to investigate the importance of environmental factors in the recovery of forests after stand-replacing disturbances, we considered nine variables to explore the drivers of forest recovery in the post-disturbance period, including Temp, Prec, VPD, SM, PAR, soil clay content, the losses of LAI, GPP and AGB in fire and logging events (i.e., ΔLAI, ΔGPP and ΔAGB) (Table S1).ΔLAI, ΔGPP and ΔAGB were calculated as the difference between the values in the disturbance year and one year before the disturbance year.Note that ΔLAI, ΔGPP and Fig. 2. Example of how the time series was built based on the forest disturbance dataset from 2000 to 2012."0 year" on the time scale of the time series represents the year of fire, "-1 year" represents the one year before fire, positive values represent the years after fire (e.g., "1 year" means the 1st year after fire).
L. Yu et al.ΔAGB were not used in the predictive variables to investigate the impact of insect outbreaks on forest recovery dynamics.The increase of LAI, GPP and AGB in the post-disturbance period was used as the target variables, respectively.
The Random Forest model (RF) was used to quantify the importance of each environmental factor in determining the recovery of LAI, GPP and AGB after stand-replacing disturbances.The RF model employs an ensemble learning approach for classification and regression (Breiman 2001).This model can effectively handle both categorical and continuous predictors, prevent overfitting, and accurately measure the importance of variables (Iverson et al. 2008).The 'TreeBagger' function in the RF model based on MATLAB 2019b was used to estimate the importance of the predictive variables.To increase the robustness of the results, we conducted 100 runs and took the mean value of all runs as the output result.The variables with the highest importance were considered as the key drivers for the recovery of LAI, GPP and AGB after fires, logging and insect outbreaks.

Statistical analysis
To analyze the forest area loss and AGB loss caused by stand-replacing disturbances, zonal statistical methods were used in this study.Firstly, we calculated the forest area loss and AGB loss resulting from fire, logging and insect outbreaks from 2000 to 2012 for both the study region and each province, respectively.Secondly, we computed the annual forest area loss caused by fire, logging and insect outbreaks from 2000 to 2012, respectively.These statistical analyses were performed in Matlab 2019b and Arcgis 10.5.

Carbon loss caused by stand-replacing disturbances
Our results showed that from 2000 to 2012, fire affected a total forest area of 4.7 × 10 6 ha, followed by logging (2.3 × 10 6 ha) and insect outbreaks (6.9 × 10 5 ha) (Fig. 4).This indicates that fire was the major stand-replacing disturbance in the western North American boreal forests.The largest forest fire area (9.2 × 10 5 ha) was observed in (Fig. 5a), accounting for 20% of the total forest fire area during the study period (Fig. 4).The maximum forest logging area was found in (2.4 × 10 5 ha) (Fig. 5b).The period between 2006 and 2009 had the largest forest area (3.6 × 10 5 ha) disturbed by insect outbreaks (Fig. 5c), accounting for 52% of the total forest area disturbed by insect outbreaks between 2000 and 2012 (Fig. 4).
Spatially, forest fire was mainly observed in Alaska, Alberta, Saskatchewan and Yukon (Fig. 6a, Fig. 7a), accounting for almost 47% of the total forest area disturbed by fire.Among them, the largest forest fire area was observed in Alaska (Fig. 7a).The outbreak of forest insect in southern of the study region was an order of magnitude greater in area than the rest of the study region (Fig. 6e).Especially, British Columbia had both the largest forest insect outbreak (3.4 × 10 5 ha) and logging area (9.8 × 10 5 ha), accounting for 49% and 43% of the total insect outbreak and logging disturbance area, respectively (Fig. 6c, Fig. 6e, Fig. 7a).In addition, the forest area in Alberta was strongly disturbed by fire and logging, contributing 18% and 50% of the total area disturbed by fire and logging, respectively (Fig. 6a, Fig. 6c, Fig. 7a).
At the provincial scale, British Columbia had the highest amount of AGB loss resulting from the three disturbance agents (13.2 Tg/yr), of which logging and insect outbreaks contributed 65% and 22%, respectively (Fig. 7b).Sequentially, Alberta had an AGB loss derived from disturbances of 12.5 Tg/yr, of which fire and logging contributed 36% Fig. 3. Example of 3 £ 3 matrix with the disturbed pixels and the adjacent undisturbed pixels.Fire, logging and insect represent that the pixel was disturbed by fire, logging and insect outbreaks, respectively.Non disturbance represents the undisturbed pixel.P1 ~ P9 represents the numbering of the pixels, respectively.and 59%, respectively, followed by Alaska (8.9 Tg/yr), where fire accounted for 90% of the carbon loss (8.06 Tg/yr) (Fig. 7b).In sum, the AGB loss in these three provinces accounted for 78% of the total AGB loss in the study region.
The spatial distribution of AGB loss caused by stand-replacing disturbances (Fig. 6b; Fig. 6d; Fig. 6f) follows the forest area loss map (Fig. 6a; Fig. 6c; Fig. 6e).For example, AGB loss caused by fire mainly occurred in Alaska, Alberta, Saskatchewan and Yukon (Fig. 6b), while AGB loss caused by logging (Fig. 6d) and insect outbreaks (Fig. 6f) are more prominent in southern areas of the study domain.

The trajectories of forest recovery
In the post-disturbance period, our results showed that AGB did not recover to its pre-disturbed levels in the 10th year after the fire (Fig. 8).Over the study region, fire caused an overall decrease of 0.04 Tg in AGB relative to the pre-fire value (Fig. 8a), followed by a negative AGB trend over the next 8 years, reaching a minimum value of AGB in year 8 (-0.07 Tg).Although there was an increasing trend in the 9th and 10th years after the fire, AGB did not fully recover to its pre-fire level.Relative to the declining AGB trend, the average recovery time for GPP after fire disturbance was 8 years (Fig. 8a).The minimum GPP values were observed in the first year after the fire (-169 g C/m 2 yr − 1 ), followed by a positive recovery trajectory over the next 10 years, with GPP exceeding the pre-fire level in year 8.This suggests that GPP has the ability to recover quickly after fire disturbance.Similar rapid post-fire recovery was observed in LAI (Fig. 8a), it took about 8 years for LAI to fully recover after fire disturbance.
The results also showed that the recovery trajectory of AGB after logging disturbance was similar to that observed after fire disturbance, i. e., AGB did not recover to the pre-logging levels in the 10th year after logging (Fig. 8c).It took 3 years for GPP and LAI to fully recover after logging disturbance (Fig. 8c), which was considerably faster than the recovery of GPP and LAI from fire disturbance (8 years) (Fig. 8a).For insect outbreaks, AGB decreased rapidly after insect outbreaks in 2005, and there were no signs of recovery in 2012 (7 years after insect  outbreaks) (Fig. 8e).In contrast, GPP and LAI recovered to their predisturbed levels after 7 years post insect outbreaks (Fig. 8e).
Note that environmental factors may affect forest recovery after stand-replacing disturbances.We investigated forest recovery by removing the influence of environmental factors.Over the fire disturbed regions, our results showed a smaller magnitude of AGB loss (0.02 Tg) when the effect of environmental factors was removed (Fig. 8b) compared to the results with the effect of environmental factors (0.04 Tg) in the fire year (Fig. 8a), although similar recovery trajectories were observed for the two results.Similar results were also observed over the logging disturbance and insect outbreak regions that the results by removing environmental factors had a smaller AGB loss, i.e., AGB loss was 0.02 Tg (Fig. 8d) over the logging region in logging year and 0.03 Tg (Fig. 8f) over the insect outbreak region in the most serious year without the effect of environmental factors, the corresponding AGB loss was 0.03 Tg (Fig. 8c) and 0.06 Tg (Fig. 8e) with the effect of environmental factors.These results suggest that environmental factors could potentially enhance the AGB loss during the disturbance events, but did not obviously impose an extra influence on the AGB recovery.This conclusion was also supported by the GPP and LAI results over the fire and logging disturbance regions (Fig. 8a, Fig. 8b, Fig. 8c, Fig. 8d): they both had similar recovery trajectories between the results with and without the effect of environmental factors.
Note that, over the insect outbreak regions, the recovery trajectories of GPP and LAI were obviously different when environmental factors were considered compared to those trajectories without environmental factors (Fig. 8e, Fig. 8f): GPP and LAI recovered to its pre-disturbed levels in year 7 from the results when including environmental factors, but continuous GPP and LAI decreasing trends were observed from the results without environmental factors 7 years after insect outbreaks.This suggests that environmental factors potentially dominated the recovery trajectories of GPP and LAI over the insect outbreak regions.

Environmental factors influencing forest recovery
The analysis of environmental factors on forest recovery showed that over fire disturbance regions, the losses of AGB in fire events (ΔAGB) is the most important factor for AGB recovery (Fig. 9a).Similar results can also be inferred from the recovery of GPP and LAI (Fig. 9b; Fig. 9c): the losses of GPP and LAI in fire events (ΔGPP and ΔLAI) were the dominant factors for forest recovery after the fire.
Over logging disturbance regions, VPD and soil clay content are the important factors for AGB and GPP recovery (Fig. 9d, Fig. 9e).In addition, water availability (e.g., precipitation and soil moisture), temperature and radiation availability (e.g., PAR) played important roles in forest recovery (Fig. 9f).Similar results were found from the insect outbreak regions (Fig. 9g, Fig. 9h, Fig. 9i).This indicates that climatic  factors have a vital role in the recovery of forests after insect outbreaks.

Carbon losses caused by stand-replacing disturbances
During 2000-2012, we estimated a forest loss (7.45 × 10 6 ha) was considerably lower than the results (2.62 × 10 7 ha) reported by Hansen et al. (2013) (Supplementary Fig. S1).This is not surprising considering the different definitions of 'forest loss' between ours and Hansen et al. (2013).The term 'forest loss' is defined here as a land cover change from forest to non-forest (Zhang et al. 2022), while 'forest loss' in Hansen et al. (2013) was defined as the temporary or permanent loss of tree cover.The forest areas that had a temporary or permanent loss of tree cover but did not show a change in the land cover type were not recognized as 'forest loss' by the product of Zhang et al. (2022).In  c) and (e) showed the forest recovery trajectories from fires, logging and insect outbreaks, respectively.Correspondingly, (b), (d), and (f) showed the forest recovery trajectories derived by removing the influence of environmental factors, respectively.Colored shading for each section denotes ± 1 standard deviation.For fire-and loggingdisturbance, the recovery of AGB, GPP and LAI were calculated as the difference of AGB, GPP and LAI between each year in the 14-year time series and its value for the year prior to disturbance.The results only show the recovery trajectories for 10 years after disturbances.For insect outbreaks, the recovery of AGB, GPP and LAI were calculated as the difference of AGB, GPP and LAI between each year from 2000 to 2012 and the year 2000.2022) product only accounted for three disturbance agents (e.g., fires, logging, and insect outbreaks), while the disturbance agents in Hansen et al. (2013) included all natural or anthropogenic factors (e.g., fire, logging, windstorm) that result in tree cover removal.Note that, the CCDC algorithm used in Zhang et al. (2022) requires five consecutive observations to detect a disturbance, thus the disturbance areas in the last two years of the study period were underestimated (Zhang et al. 2022), which could be another reason for the smaller areas in Zhang et al. (2022) relative to Hansen et al. (2013).
Stand-replacing fire, as the dominant factor in forest loss over the North American boreal forests, resulted in an AGB loss of 5.03 Mg/ha yr − 1 , being 8 times higher than the results reported by Wang et al. (2021) (0.64 Mg/ha yr − 1 ), confirming that stand-replacing fire is the main driver of carbon loss in these regions.The higher estimated carbon loss may be attributed to stand-replacing fire was considered in this study, instead of the wildfire including both stand-replacing fire and the additional inclusion of non-stand replacing fire in Wang et al. (2021).Stand-replacing fire completely burns the aboveground tree, shrub, and herbaceous biomass (Zhang et al. 2016), whereas non-stand-replacing fire burns incompletely and leaves substantial residual vegetation at a site, resulting in a relatively smaller amount of AGB loss (Wulder et al. 2020).Similarly, insect outbreaks resulted in a carbon loss of 6.63 Mg/ ha yr − 1 , which was much higher than the estimate from Kurz et al. (2008) of 0.72 Mg/ha yr − 1 .This is partly due to the fact that standreplacing insect outbreaks were considered in our study, whereas Kurz et al. (2008) considered non-stand replacing insect outbreaks.Also, stand-replacing logging resulted in a carbon loss of 7.66 Mg/ha yr − 1 , which was 9 times higher than the results from Wang et al. (2021) (0.84 Mg/ha yr − 1 ), suggesting that the effect of logging on the carbon loss over the North American boreal forest is substantial.

Forest recovery in the post-disturbance period
Our results demonstrated that the recovery of GPP and LAI was faster than that of AGB in the post-disturbance period, partly attributed to the long-term legacy of disturbance on forest AGB, and it took many years for AGB to fully recover relative to the rapid recovery of GPP and LAI (Nguyen et al. 2020).The differential recovery of GPP, LAI and AGB can be explained by the post-disturbance recovery process of vegetation, which can be divided into four stages: understory grass, shrubs, opencanopy forest and closed-canopy forest (Goulden et al. 2011).Understory grasses and shrubs, as pioneer species in the post-disturbance period, can reestablish LAI and GPP rapidly.In contrast, the accumulation of AGB in understory grass and shrubs is slow, and will increase more rapidly as trees establish and grow.In addition, moss was the first ground species established in the post-disturbance period (Schimmel and Granstrom 1996).The contribution of moss to the satellite greenness signal is large in young stands growing after disturbances, which may result in an overestimate of LAI and GPP (Yuan et al. 2014).
The results also showed that the recovery of GPP and LAI from logging was significantly faster than the recovery from fire.This could be because the increase of GPP and LAI from the residual vegetation accelerates the forest recovery over the logging regions, as the growth of residual vegetation is enhanced due to both more light reaching the understory and the alleviated competition for water after logging (McDowell et al. 2008).In addition, the fast recovery over the logging region could be because more productive sites in North American boreal forests are targeted for harvest (Wulder et al. 2020), and if harvested areas are replanted they will recover quickly (Fredeen et al. 2007).

Factors influencing forest recovery in the post-disturbance period
ΔAGB, ΔGPP and ΔLAI have been identified as indicators of the severity of disturbance (Yang et al. 2022), and these indicators were also found to be the most important factors for forest recovery in fire disturbance areas in this study.This finding is supported by previous studies, which have shown that more severe burns tend to result in higher vegetation recovery (Cai et al. 2013;Shvetsov et al. 2019).For instance, an analysis of forest recovery in the Central Siberia boreal forests after fire found that NDVI had the highest recovery rates in the most severely burned forests (Cuevas-González et al. 2009).
Note that forest recovery in logging disturbance areas is sensitive to environmental factors (e.g., VPD, soil clay content, and precipitation).For example, the increased VPD may induce stomatal closure in plants to prevent excessive water loss because of the high evaporative demand of the air, which may lead to a negative impact on forest recovery (Yuan et al. 2019).Soil clay content is also a significant factor for forest recovery in logging regions, as soil clay content partly controls the soil moisture available for plant roots, e.g., plants cannot make effective use of soil moisture in clay soils as they have more negative water potential values (Fensham et al. 2015).Also, clay soils with higher nutrients may have harmful effects, forests in high nutrient clay soils may be adapted to higher nutrient availability, which may lead to a higher embolism risk in the event of water stress (Gessler et al. 2016).
Climate factors such as temperature and precipitation are the most important factors for the forest recovery over the insect outbreak region, partly because climate factors can affect insect population dynamics by influencing fecundity, survival and dispersal (Kurz et al. 2008).Early spring temperature was shown to have positive and significant influence on insect outbreaks, because the higher larval survival rates associated with warmer spring may trigger the outbreak of insects (Chen et al. 2018).Also, temperature in winter and spring can impact the survival of hatched eggs and larvae of defoliator insects in Canada (Uelmen et al. 2016).Moreover, physiological stress triggered by drought (Speer et al. 2001) has been demonstrated to reduce the vitality of trees, making them more susceptible to insect attacks (Guarín and Taylor 2005).

Limitations of the study
Some uncertainty remains in understanding the response of North American boreal forests to stand-replacing disturbances: (1) Landsatbased AGB data reflects mostly changes in the upper canopy instead of the whole above ground tree, because Landsat data are largely sensitive to changes in the vegetation canopy, not in tree trunks.The recovery of vegetation canopy could be much faster than the recovery of trunks.(2) The disturbance product, AGB, and LAI were aggregated to 0.05 degree resolution in order to explore forest recovery in response to standreplacing disturbances, yet the coarse spatial resolution product failed to separate pixel-scale carbon loss and carbon gain.This may result in the rapid recovery of forests.Furthermore, there may be other agents of disturbance in a 0.05 degree pixel (e.g., drought, windstorm), which are not considered in our study.(3) The study period of this study was from 2000 to 2012, which is not long enough to investigate the recovery of AGB.Previous studies have shown that the recovery of AGB in boreal forests following stand-replacing disturbances may take decades (Goulden et al. 2011;Wang et al. 2021).Thus, future studies should consider long time series of datasets to accurately estimate the recovery of AGB after stand-replacing disturbances.(4) Land use change will have an impact on the carbon dynamics in boreal forests, which were not accounted for in this study.Future studies should consider the impact of land use change on forest recovery in North American boreal forests, considering that previous studies often ignored the effects of land use change (Cuevas-González et al. 2009;Hicke et al. 2003;White et al. 2017).

Conclusions
The response of boreal forests to stand-replacing disturbances in western North America was studied using 12-year of AGB, GPP, and LAI datasets.By coupling data from disturbed and adjacent undisturbed control pixels, we were able to separate interannual variations due to environmental factors from those caused by forest disturbances.We found that stand-replacing fire was the major disturbance agent in western North American boreal forests, followed by stand-replacing logging and stand-replacing insect outbreaks.The recovery trajectories of AGB, GPP, and LAI are different, suggesting the different sensitivity of these vegetation indicators to the different parts of the forest canopies (trunks, stems, leaves) and their different regrowth time.By removing the impact of environmental factors, our results showed a smaller magnitude of AGB, GPP and LAI loss relative to the results including the effect of environmental factors.We also found that the dominant factors for forest recovery in response to different disturbances are heterogeneous.The synergistic use of the recently released disturbance and AGB products, with GPP and LAI data provides a new insight into the responses of North American boreal forests to multiple stand-replacing disturbances.The asynchronous recovery of LAI, GPP and AGB in the post disturbance period highlights the importance of considering the multidimensional characteristics of forests (such as vegetation cover, carbon storage, and biomass) for future management of carbon dynamics in North American boreal forests.Considering that the recovery of AGB in the post disturbance period takes more than 10 years, it is crucial to implement timely management strategies that focus on preventing disturbances and enhancing the potential for forest restoration and reforestation.This approach is pivotal in preserving the North American boreal forests as a long-lasting carbon sink.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The study region.The study region (red borders) is located in the western North American boreal forests (100 degreesW-168 degreesW, 52 degreesN-74 degreesN), which spans Alaska and western Canada.The background shows the Europe Space Agency (ESA) Climate Change Initiative (CCI) land cover map in 2012 (ESA 2012), which was aggregated to 0.05 degree resolution by the dominant class.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. The disturbance area and AGB loss of the three disturbance agents from 2000 to 2012.The blue bar represents the disturbance area from fire, logging, insect outbreaks and the total disturbance area of the three agents.The orange bar represents the AGB loss from fire, logging, insect outbreaks and the total AGB loss of the three agents.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Disturbance area of the three disturbance agents at an annual scale from 2000 to 2012.(a) fire area; (b) logging area; (c) insect outbreak area.Note that the y-axis is different for each plot.

Fig. 6 .
Fig. 6.Spatial patterns of disturbance agents and spatial patterns of AGB loss from 2000 to 2012.The left panel represents the spatial patterns of disturbance agents.(a) fire; (c) logging; (e) insect outbreak.The 30 m grid cells were aggregated to 0.05 degree by recording the disturbed area fraction of the disturbance dataset within 0.05 degree grid cells.Correspondingly, the right panel represents the spatial patterns of AGB loss by fire (b), logging (d), and insect outbreak (f), respectively.The 30 m grid cells were aggregated to 0.05 degree spatial resolution by recording the sum of AGB loss.The unit of AGB loss is Tg.

Fig. 7 .
Fig. 7. Disturbance area and AGB loss in each province from 2000 to 2012.(a) disturbance area.(b) AGB loss.The orange bar represents fire, the yellow bar represents insect outbreaks, and the blue bar represents logging.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Trajectories of forest recovery in response to disturbances.(a), (c) and (e) showed the forest recovery trajectories from fires, logging and insect outbreaks, respectively.Correspondingly, (b), (d), and (f) showed the forest recovery trajectories derived by removing the influence of environmental factors, respectively.Colored shading for each section denotes ± 1 standard deviation.For fire-and loggingdisturbance, the recovery of AGB, GPP and LAI were calculated as the difference of AGB, GPP and LAI between each year in the 14-year time series and its value for the year prior to disturbance.The results only show the recovery trajectories for 10 years after disturbances.For insect outbreaks, the recovery of AGB, GPP and LAI were calculated as the difference of AGB, GPP and LAI between each year from 2000 to 2012 and the year 2000.

Fig. 9 .
Fig. 9.The relative importance of environmental variables to forest recovery.(a), (b) and (c) showed the importance of environmental variables for AGB, GPP and LAI recovery in fire disturbance regions, respectively.(d), (e) and (f) showed the importance of variables for AGB, GPP and LAI recovery in logging disturbance regions, respectively.(g), (h) and (i) showed the importance of variables for AGB, GPP and LAI recovery in insect outbreak regions, respectively.The change in color from green to yellow indicates a decrease in importance.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1
AGB loss caused by three disturbance agents.