Land Use Changes, Disturbances, and Their Interactions on Future Forest Aboveground Biomass Dynamics in the Northern US

Land use change (LUC), disturbances, and their interactions play an important role in regional forest carbon (C) dynamics. Here we quantified how these activities and events may influence future aboveground biomass (AGB) dynamics in forests using national forest inventory (NFI) and Landsat time series data in the Northern United States (US). Total forest AGB predictions were based on simulations of diameter growth, mortality, and recruitment using matrix growth models under varying levels of LUC and disturbance severity (low (L), medium (M), and high (H)) every five years from 2018 to 2098. Land use change included the integrated effects of deforestation and reforestation/afforestation (forest [F]→agriculture [A], settlements [S, urbanization/other], and A&S→F), specifically, conversion from F→A, F→S, F→A&S, A→F, S→F, and A&S→F. Disturbances included natural and anthropogenic disturbances such as wildfire, weather, insects and disease, and forest harvesting. Results revealed that, when simultaneously considering both medium LUC and disturbances, total forest AGB predictions of LUC + fire, LUC + weather, LUC + insect & disease, and LUC + harvest indicated substantial increases in regional C stocks (± standard deviation) from 1.88 (±0.13) to 3.29 (±0.28), 3.10 (±0.24), 2.91 (±0.19), and 2.68 (±0.17) Pg C, respectively, from 2018 to 2098. An uncertainty analysis with fuzzy sets suggested that medium LUC under disturbances would lead to greater forest AGB C uptake than undisturbed forest C uptake with high certainty, except for LUC + harvest. The matrix models in this study were parameterized using NFI and Landsat data from the past few decades. Thus, our results imply that if recent trends persist, LUC will remain an important driver of forest C uptake, while disturbances may result in C emissions rather than undisturbed forest C uptake by 2098. The combined effects of LUC and disturbances may serve as an important driver of C uptake and emissions in the Northern US well into the 21st century.


Introduction
The accelerating rise of carbon (C) emissions and resulting net increases in atmospheric carbon dioxide (CO 2 ) is a substantial forcing factor in global warming [1]. The dynamics of the terrestrial C cycle have come under increasing focus because of climate and land use change (LUC, deforestation and reforestation/afforestation). Such changes, when combined with disturbances, are the primary drivers of global changes in forest C stocks [2,3]. Although recent studies and reports suggest the C uptake in forests has partially offset net C emissions in the United States (US) [2,4], regional assessments of silvicultural systems including, but not limited to, the partial harvesting of large trees for forest management. In order to explore the sensitivity of LUC and disturbances, LUC and all disturbance events were simulated with low (L), medium (M), and high (H) severity/intensity. To investigate this we first used matrix growth models and the FVS to predict forest AGB C dynamics over the past 20 years (1998-2018) as well as the future 80 years . Secondly, we quantified the separate effects of LUC and disturbances with low, medium, and high severity, and the comprehensive effects of both LUC and disturbances on forest AGB dynamics. Finally we adopted fuzzy sets to represent variability in forest AGB predictions resulting from uncertainty in the simulated LUC, disturbances, and the interactive effects of these factors on forest AGB across the Northern US.

Research Region
The study domain covered eleven states in the Northern US including Minnesota, Wisconsin, Illinois, Indiana, Michigan, Ohio, Pennsylvania, New York, Vermont, New Hampshire, and Maine ( Figure 1). The major land-cover types included forest (F), agriculture (A), and settlements/other (S, Figure 2). Forests were the dominant land use across nearly half (50.86%) of the Northern US. Forest land use was greatest (82.65%) in the upper Great Lakes and New England. Agriculture land use was most dominant in the Great Plains extending into Illinois, Indiana, and Ohio. Settlements/other was the dominant land use (17.28%) across portions of Ohio and New York ( Figure 2). Across the Northern US, forestland areas of deforestation for F→A, F→S, and F→A&S decreased 331,711 ha, 408,482 ha, and 740,193 ha, respectively, during the remeasurement period of five years (one measurement during 2008-2012 and remeasured five years later during 2013-2018). Meanwhile, forestland areas of reforestation/afforestation for A→F, S→F, and A&S→F increased 946,968 ha, 323,456 ha, and 1,270,424 ha, respectively, during the five years (Table 1 and Figure 2). Forestland areas under medium LUC (F→A&S + A&S→F) increased 530,231 ha across the study area. This study region was ideal for examining LUC and disturbance impacts on forest AGB given the abundant forest cover in these areas and the wide range of land use and disturbance types affecting these systems.    agriculture, and settlements/other) summarized for discrete hexagons across the Northern US during the remeasurement period of five years. Land use change was estimated for each hexagon that had at least eight sample plots within it, with change calculated as the difference in the percent land use between time one and two by land use category [5].

National Forest Inventory Data
Comprehensive NFI data were extracted from the US Department of Agriculture Forest Service, Forest Inventory and Analysis (FIA) database (https://apps.fs.usda.gov/fia/datamart/datamart.html) in 11 states of the Northern US for a total of 78,458 plots ( Figure 3, with one measurement between 2008 and 2012 remeasured five years later between 2013 and 2018). Each permanent ground plot was composed of four smaller fixed-radius (7.32 m) subplots spaced 36.6 m apart in a triangular Figure 2. Percent land use (a,b,c) and land use change (LUC) (d,e,f) by general land use classes (forest, agriculture, and settlements/other) summarized for discrete hexagons across the Northern US during the remeasurement period of five years. Land use change was estimated for each hexagon that had at least eight sample plots within it, with change calculated as the difference in the percent land use between time one and two by land use category [5].

National Forest Inventory Data
Comprehensive NFI data were extracted from the US Department of Agriculture Forest Service, Forest Inventory and Analysis (FIA) database (https://apps.fs.usda.gov/fia/datamart/datamart.html) in 11 states of the Northern US for a total of 78,458 plots ( Figure 3, with one measurement between 2008 and 2012 remeasured five years later between 2013 and 2018). Each permanent ground plot was composed of four smaller fixed-radius (7.32 m) subplots spaced 36.6 m apart in a triangular arrangement around a central subplot. For each permanent ground plot, plot-level estimates of forest AGB (Figure 3), trees per hectare, basal area per hectare, and measurements of slope and elevation were used (Table 2). Land use was characterized by a hexagonal grid from Woodall et al. [5] to facilitate analysis at spatial scales optimized for plot sampling intensity of the NFI. Land use change was estimated for each hexagon that had at least eight sample plots within it, with change calculated as the difference in the percent land use between time one and two by land use category. Detailed information about land use characterizations can be found in Woodall et al. [5]. arrangement around a central subplot. For each permanent ground plot, plot-level estimates of forest AGB (Figure 3), trees per hectare, basal area per hectare, and measurements of slope and elevation were used (Table 2). Land use was characterized by a hexagonal grid from Woodall et al. [5] to facilitate analysis at spatial scales optimized for plot sampling intensity of the NFI. Land use change was estimated for each hexagon that had at least eight sample plots within it, with change calculated as the difference in the percent land use between time one and two by land use category. Detailed information about land use characterizations can be found in Woodall et al. [5].

Landsat Data
Landsat time-series data were obtained from the Google Earth Engine (https://earthexplorer.usgs.gov/). Landsat time-series data were acquired by Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), Landsat 5 Thematic Mapper (TM), and Landsat 5 Multispectral Scanner (MSS) instruments from 2008 to 2018. Landsat 7 ETM+ and Landsat 8 OLI/TIRS data were processed upon download but used predicted ephemeris, initial bumper mode parameters, or initial TIRS line-of-sight model parameters. The data were reprocessed with definitive ephemeris, updated bumper mode

Landsat Data
Landsat time-series data were obtained from the Google Earth Engine (https://earthexplorer.usgs. gov/). Landsat time-series data were acquired by Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), Landsat 5 Thematic Mapper (TM), and Landsat 5 Multispectral Scanner (MSS) instruments from 2008 to 2018. Landsat 7 ETM+ and Landsat 8 OLI/TIRS data were processed upon download but used predicted ephemeris, initial bumper mode parameters, or initial TIRS line-of-sight model parameters. The data were reprocessed with definitive ephemeris, updated bumper mode parameters and refined TIRS parameters, and the products were transitioned to either Tier 1 or Tier 2 and removed from the Real-Time tier (https://landsat.usgs.gov/landsat-collections).

Description of the Matrix Models
Matrix growth models were applied to project AGB dynamics of remaining forests controlling for diameter growth, mortality, and recruitment as follows [35,36]: in which y t and y t+1 are number of live trees per hectare at time t and t + 1, respectively; G t is a transition matrix describing how trees grow or die between t and t + 1; R t is a recruitment vector representing the trees naturally recruited in the smallest diameter class between t and t + 1; ε is a random error; and G t and G it are state-and time-dependent transition matrices. G it is a submatrix of G t , where in which a ijt represents the probability that a tree of species group i and diameter class j stays alive in the same size class between t and t + 1. The transition probability of growth, b ijt , represents a tree of species i, and diameter class j stays alive and grows into diameter class j + 1 between t and t + 1, assuming that trees were evenly distributed within a diameter class. b ijt was estimated as the annual tree diameter growth, g ijt , divided by the width of the diameter class. a ijt and b ijt are related by where m ijt is the tree mortality of species i and size class j between t and t + 1. R t represents the number of trees naturally recruited into the smallest size class of each species group, between t and t + 1. R it is a sub vector of R t representing recruitment of species group i at time t, where The annual diameter growth of the tree of species group i and size class j from t and t + 1 is represented by the following model ( [27]; all notations defined in Table 2): (5) in which α i 's are parameters to be estimated with the generalized least squares (GLSs, [37]) for diameter growth of species group i and size class j. θ is an error term.
Tree mortality of species group i and diameter class j at time t, m ijt = P(M ijt = 1|x) is estimated with a Probit function [38]. M ijtk is a binary variable representing whether a tree of species i and diameter class j died (M ijtk = 1) or not (M ijtk = 0): where Φ is the standard normal cumulative function, δ i 's are parameters estimated by maximum likelihood, and ξ is an error term.
For the recruitment of species group i, R i is estimated with a Tobit model [39]: with where Φ is the standard normal cumulative distribution function, ϕ is the standard normal probability density function, and µ is an error term. The Tobit model explicitly accounts for unobserved recruitment values that are left-censored at the preset diameter limit.
The forest AGB C density is estimated with the following model [27]: in which ψ i is a vector of forest AGB C density, ω i 's are parameters to be estimated with the GLS, ω is an error term, and ψ i is calculated with Equation (9). For the future simulations using matrix growth models, a certain percentage of trees were removed under disturbances and added under LUC (F→A&S + A&S→F) every five years. In order to backcast forest AGB dynamics using matrix growth models from 2018 to 1998, we applied similar annual rates of mortality and recruitment, minus the annual rate of tree growth, without considering LUC and disturbances. Stand density was set to zero if stand age was less than 20 years.

Model Calibration and Validation
We adopted the same variables and functional form for the underlying diameter growth, mortality, recruitment, and AGB models as those established in previously published matrix models [27] to maintain accuracy and avoid overfitting and multicollinearity between predicted variables. A variety of stand level and Landsat variables selected for two species groups (deciduous and coniferous) were used in the matrix growth models (Table 3). Physiographic variables, elevation (E) and slope (S), were used to control for site productivity [40]. Landsat variables, tasseled cap brightness (TCB), disturbance index (DI), enhanced vegetation index (EVI), shortwave infrared surface reflectance (SWIR), tasseled cap greenness (TCG), soil adjusted vegetation index (SAVI), and tasseled cap angle (TCA) were used to replace total stand basal area (B) as key predictors, due to the correlation between the Landsat predictors and forest stand characteristics [27]. For model calibration purposes, 62,766 permanent ground plots (80%) and plot-corresponded Landsat data were used to estimate parameters of the matrix growth model. In addition, model validation (the rest: 15,692 permanent ground plots, 20%) followed Ma et al. [27], where matrix models were used to enable various combinations of LiDAR, Landsat, and field data, especially Landsat data to estimate large-scale AGB dynamics (i.e., the central component of carbon stock monitoring) without loss of accuracy from only using variables from forest inventories ( Figure 4). Finally, the 15,692 permanent ground plots and plot-corresponded Landsat data were applied to make the future predictions. Table 3. Estimated parameters of the diameter growth model, mortality model, recruitment model, and aboveground biomass model in the matrix growth models.

Land Use Change and Disturbance Scenarios
Recent studies have suggested that LUC and forest disturbances (e.g., harvest, wildfire, insects and diseases, and weather) simultaneously affect regional forest C uptake [2,5,6,24]. Thus we used three general land use categories following Woodall et al. [5]: forest (F), agriculture (A), and settlements (S, urbanization/other). In addition, we explored effects of disturbances such as harvest, wildfire, insects and diseases, and weather with low, medium, and high severity/intensity (Table 4). Annual medium disturbance rates (percent of the forest area influenced per year) were obtained and bounded from Wear and Coulston [6] for the Northern US based on empirical information. Moreover, average annual LUC and disturbance rates increased by 1% from low-to-medium to medium-to-high severity for sensitivity analysis. Based on the average annual disturbance rates, we assumed 50% removal (i.e., mortality and subsequent C emission) of all trees subject to fire, insects and diseases, and weather, and only removal of large trees (greater than 40 cm) under harvesting. In order to study the effects of LUC, disturbances, and their comprehensive effects on forest C dynamics, we predicted future total forest AGB dynamics considering scenarios of LUC and disturbances every five years (same period with historical LUC) from 2018 to 2098. The scenarios were defined as follows: (a) LUC included F→A, F→S, F→A&S, A→F, S→F, A&S→F, and F→A&S + A&S→F with low, medium, and high intensity; (b) Disturbances included fire, weather, insect and disease, and harvest with low, medium, and high severity; (c) Medium LUC + disturbances included LUC + fire, LUC + weather, LUC + insect and disease, and LUC + harvest. Table 4. Average annual disturbance rates with low (L), medium (M), and high (H) severity in the study region.

Land Use Change and Disturbance Scenarios
Recent studies have suggested that LUC and forest disturbances (e.g., harvest, wildfire, insects and diseases, and weather) simultaneously affect regional forest C uptake [2,5,6,24]. Thus we used three general land use categories following Woodall et al. [5]: forest (F), agriculture (A), and settlements (S, urbanization/other). In addition, we explored effects of disturbances such as harvest, wildfire, insects and diseases, and weather with low, medium, and high severity/intensity (Table 4). Annual medium disturbance rates (percent of the forest area influenced per year) were obtained and bounded from Wear and Coulston [6] for the Northern US based on empirical information. Moreover, average annual LUC and disturbance rates increased by 1% from low-to-medium to medium-to-high severity for sensitivity analysis. Based on the average annual disturbance rates, we assumed 50% removal (i.e., mortality and subsequent C emission) of all trees subject to fire, insects and diseases, and weather, and only removal of large trees (greater than 40 cm) under harvesting. In order to study the effects of LUC, disturbances, and their comprehensive effects on forest C dynamics, we predicted future total forest AGB dynamics considering scenarios of LUC and disturbances every five years (same period with historical LUC) from 2018 to 2098. The scenarios were defined as follows: (a) LUC included F→A, F→S, F→A&S, A→F, S→F, A&S→F, and F→A&S + A&S→F with low, medium, and high intensity; (b) Disturbances included fire, weather, insect and disease, and harvest with low, medium, and high severity; (c) Medium LUC + disturbances included LUC + fire, LUC + weather, LUC + insect and disease, and LUC + harvest. Effects of LUC and disturbances were added in the forest AGB dynamics every five years from 2018 to 2098. Forest AGB C changes induced from LUC and disturbances were estimated as follows: forest AGB C changes from LUC and disturbances = forestland area changes from LUC and disturbances × forest AGB C density.

Fuzzy Sets representing LUC and Disturbance Uncertainty
Uncertainty analysis is a necessary procedure in model simulation [16]. Uncertain LUC and disturbances in future have led to high variability in predicted values of forest AGB. The averages of predicted forest AGB are important point estimates, but understanding the associated risk, ranges, or sets indicating prediction uncertainty is essential. To more accurately quantify this uncertainty, we used fuzzy sets, which involved defining membership functions that determined the level of uncertainty [41]. The theory of fuzzy promulgates the concept of a possibility distribution as a fuzzy restriction acting as a flexible constraint on the average values that may be assigned to a variable [41]. A trapezoidal fuzzy set was used, mathematically expressed as f (x; a, b, c, d) = max (min (x − ab − a, 1, d − xd − c), 0). [a,b] and [c,d] are the uncertainty intervals with membership degrees ranging from 0 to 1. [b,c] represents the certainty interval for which the membership degree is 1. [a,d] is a measure of the total range of uncertainty arising from LUC and disturbance occurrences. Given the average value of predicted forest AGB (X) and its relative standard deviation (S r ) from simulations, a, b, c, and d values can be calculated as follows [42]:

Forest Vegetation Simulator Comparisons
In order to more thoroughly evaluate our matrix model predictions of AGB accumulation, forest inventory plot attributes (e.g., tree abundance by species and diameter class, and AGB) were projected using a forest stand vegetation simulator (FVS). Similar to the matrix growth models, FVS is a simulation tool that relies on empirical growth and yield models that have been extensively used in the US [43]. FVS is used broadly to predict forest stand dynamics, and is an approved quantification tool by the American C Registry [44]. Stand condition by the year 2018, which was the same as the matrix models, was executed as the starting point. The TREEBIO report from the FVS was used to provide plot-level AGB outputs for the same research area. Then, area-wide total forest AGB was projected for the future 80 years at decadal intervals until 2098, without considering LUC and disturbances.

Results
Based on the 80-year future predictions (2018-2098), total forest AGB of the Northern US would increase over time from 2018 (1.88 Pg C) and converge to~2.60 Pg C by 2098 without considering LUC and disturbances (scenario none, Figure 5). Total forest AGB would increase 38.30% under the condition of natural growth over the next 80 years. Comparing matrix model predictions of AGB with FVS predictions, our models predicted higher forest AGB over time, with a 24.13% increase in forest AGB by the end of period. In addition, we also made backward validation of total forest AGB from 2018 to 1998 without considering LUC and disturbances ( Figure 5). Forest AGB increased from 1.67 to 1.88 Pg C, with an increase of 12.57% during 20 years from 1998 to 2018.
Considering only disturbances of low, medium, and high severity, total forest AGB subject to all disturbance types would decrease by the year 2098 compared with natural growth without considering LUC and disturbances (scenario none) (Table 5 and Figure 5). With medium severity, total forest AGB due to fire, weather, insect and disease, and harvest would become 2.42, 2.24, 2.11, and 1.90 Pg C, respectively, at the end of period. From low-to-medium and medium-to-high severity (1% annual disturbance rate increase), the total forest AGB subject to fire would decrease by 0.11 and 0.12 Pg C by the year 2098, respectively. In addition, forest AGB subject to fire, weather, insect and disease, and harvest with high severity would decrease by 11.54%, 18.08%, 23.08%, and 28.85%, respectively, from 2018 to 2098, in comparison with a scenario without any affecting factors (referred to hereafter as "scenario none") ( Figure 5). When only considering LUC, predicted total forest AGB subject to conversions of F→A, F→S, A→F, S→F, A&S→F, and medium F→A&S + A&S→F is projected to increase over time from 2018 to 2098 (Figure 6a), while F→A&S would result in a decrease in forest AGB. The forest AGB subject to conversions of F→A&S, A&S→F, and medium F→A&S + A&S→F would be 1.86, 4.93, and 3.81 Pg C, respectively, by the year 2098 (Figure 6a). Deforestation including F→A and F→S would have forest AGB increases of 10.11% and 25.00%, respectively, while F→A&S would decrease 1.06% by the year 2098 compared to the year 2018. Reforestation and afforestation including conversions of A→F, S→F, and A&S→F would increase forest AGB by 137.77%, 71.28%, and 162.23%, respectively, by the end of simulations. LUC (F→A&S + A&S→F) with low, medium, and high intensity (i.e., 1% increase of annual LUC rates from low-to-medium and medium-to-high) would result in 79.52%, 102.66%, 119.26% increases in total forest AGB, respectively, over the next 80 years (Figure 7). Additionally, predicted forest AGB subject to conversions of F→A, F→S, and F→A&S were less than scenario none (i.e., no land use change or disturbances) while A→F, S→F, A&S→F, and all intensities of F→A&S + A&S→F were greater than scenario none (Figures 6a and 7).
When simultaneously considering both LUC and disturbances in medium severity, total forest AGB predictions of LUC + fire, LUC + weather, LUC + insect and disease, and LUC + harvest increased over time from 2018 to 2098 (Figures 6b and 7). Total forest AGB of LUC + fire, LUC + weather, LUC + insect and disease, and LUC + harvest would reach to 3.29, 3.10, 2.91, and 2.68 Pg C, respectively, by the year 2098. With medium severity, LUC + fire (annual rate 1.5%) and LUC + harvest (annual rate 7%) would increase 75.00% and 42.55% of forest AGB over the next 80 years ( Figure 6). Aboveground biomass on forestland subject to LUC + harvest and LUC + insect and disease would reduce AGB by 18.54% and 11.55% relative to AGB accumulation subject to LUC + fire by the year 2098. LUC + fire had the smallest increase, while the LUC + harvest had the largest increase for all predictions over the next 80 years.
To account for variability induced by medium LUC and disturbances, fuzzy sets were constructed for predicting forest AGB in the year 2098 ( Figure 8). Considering only LUC, A&S→F clearly outperformed the other scenarios by storing the highest total forest AGB with high certainty. Considering only disturbances, fire would lead to the highest total forest AGB under all disturbance scenarios. Weather and insect and disease would have a similar total forest AGB in the year 2098. When considering both medium LUC and disturbances, most LUC + disturbances would in general lead to higher forest AGB than undisturbed forest growth with high certainty.      To account for variability induced by medium LUC and disturbances, fuzzy sets were constructed for predicting forest AGB in the year 2098 ( Figure 8). Considering only LUC, A&S→F clearly outperformed the other scenarios by storing the highest total forest AGB with high certainty. Considering only disturbances, fire would lead to the highest total forest AGB under all disturbance scenarios. Weather and insect and disease would have a similar total forest AGB in the year 2098. When considering both medium LUC and disturbances, most LUC + disturbances would in general lead to higher forest AGB than undisturbed forest growth with high certainty.

Matrix Growth Models for Predicting Forest AGB Dynamics from 2018 to 2098
Matrix growth models are useful for predicting forest population dynamics and associated forest C dynamics under varying natural disturbances, management, and climate scenarios through consideration of diameter growth, mortality, and recruitment [15,16,26,27,35,36,45]. Moreover, matrix growth models are highly accurate in characterizing forest C dynamics over large-scale forest areas based on empirical forest inventories and Landsat data [27]. Applying selected Landsat variables, we revalidated the strength of Landsat dependent matrix growth models in their applications for predicting large-scale forest AGB dynamics [27]. Such an assertion is supported by Ma et al. [27], who suggested that Landsat variables can be used to replace stand basal areas to predict AGB at landscape scales. Moreover, we found that predicted forest AGB increased over time and converged to ~2.60 Pg C by 2098 through natural accretion of biomass (i.e., absent of disturbance and LUC). This finding is supported by prior studies which suggested that total forest AGB accumulation should increase towards a theoretical maximum within a certain period of time [26,27,[45][46][47][48][49]. Over the next 80 years, we found that total forest AGB would increase 38.30% across the study area without considering LUC and disturbances. This is consistent with other studies that have estimated potential , and A&S→F) included specific changes from F→A, F→S, F→A&S, A→F, S→F, and A&S→F. Disturbance included fire, weather, insect and disease, and harvest. All LUC and disturbance scenarios use medium intensity/severity. None represents forests that naturally grow without considering LUC and disturbances.

Matrix Growth Models for Predicting Forest AGB Dynamics from 2018 to 2098
Matrix growth models are useful for predicting forest population dynamics and associated forest C dynamics under varying natural disturbances, management, and climate scenarios through consideration of diameter growth, mortality, and recruitment [15,16,26,27,35,36,45]. Moreover, matrix growth models are highly accurate in characterizing forest C dynamics over large-scale forest areas based on empirical forest inventories and Landsat data [27]. Applying selected Landsat variables, we revalidated the strength of Landsat dependent matrix growth models in their applications for predicting large-scale forest AGB dynamics [27]. Such an assertion is supported by Ma et al. [27], who suggested that Landsat variables can be used to replace stand basal areas to predict AGB at landscape scales. Moreover, we found that predicted forest AGB increased over time and converged to~2.60 Pg C by 2098 through natural accretion of biomass (i.e., absent of disturbance and LUC). This finding is supported by prior studies which suggested that total forest AGB accumulation should increase towards a theoretical maximum within a certain period of time [26,27,[45][46][47][48][49]. Over the next 80 years, we found that total forest AGB would increase 38.30% across the study area without considering LUC and disturbances. This is consistent with other studies that have estimated potential regional forest AGB C stock increases of 15%-50% over the next 50-100 years [46,50,51]. Besides the standard FVS outputs (stem density, basal area, volume, etc.), tailored C output provided estimates of total forest AGB of the Northern US from 2018 to 2098. In comparison with the FVS, the matrix growth models projected greater total forest AGB over time and predicted more total forest AGB (24.13%) by the year 2098. The possible explanation is that FVS projects recruitment and regeneration if (and only if) the stand density is below a predesignated threshold value [30,52], and such underestimation of recruitment in FVS predictions would systematically affect the accuracy of the entire model. As a substantial portion of tree regeneration in the study region is derived from natural regeneration [53], restrictions on tree recruitment during modeling exercises such as the FVS might reduce future estimates of tree stocking and hence AGB accumulation.

Backward Validation of Matrix Models from 2018 to 1998
As signatories to the United Nations Framework Convention on Climate Change (UNFCCC), the US has reported economy-wide greenhouse gas emissions and removals each year from 1990 to the near present [4]. In addition, the Intergovernmental Panel on Climate Change (IPCC) requires that emissions and removals associated with changes in land use be attributed to specific activities. In order to harmonize IPCC reporting requirements with scenario testing for determining national future commitments, we made backward validations of forest AGB from 2018 to 1998 and found that forest AGB increased from 1.66 to 1.88 Pg (13.25%) during this time period across the Northern US. This provides a baseline for reporting C emissions and removals from 1990 to the present. Moreover, such backward validations of large-area AGB could be further applied to provide spatially explicit contemporary baseline assessments as a means of refining US greenhouse gas (GHG) estimates.

Effects of LUC, Disturbances, and Their Interactions on Future Forest AGB Dynamics
Land use change may substantially influence future C dynamics over time [54,55]. Quantifying forest C dynamics associated with LUC is important for forest C monitoring and UNFCCC reporting [5]. Conversions from agriculture and settlements/other to forestland have been occurring, which substantially offsets net C emissions across the globe. In addition, post potential forest AGB C uptake induced by decreased forest land area could be offset by increases in forest AGB via elevated forest growth [6]. Medium LUC increased total forest AGB by 1.21 Pg C over the next 80 years. This finding is consistent with prior studies that have strongly suggested LUC as a major driver of the forest C sink in the US [5,10,11,56]. Nevertheless, Birdsey et al. [57] estimated that the settlement and subsequent deforestation and/or intensive harvest of US forests emitted over 4.0 Pg C between 1715 and 1935. Most of those C emissions would not be compensated by future northern forest C accumulation via LUC and natural forest growth. The strong C uptake estimated in privately owned forestland implies that positive land use planning and forest management might be one opportunity to maintain C sink strength [5,53]. Although maintaining forestland via conservation measures will be one policy consideration in the future, additional approaches such as reforestation and afforestation concomitant with tightly coupled forest management and long-lived forest product industries might be needed to ensure maximum terrestrial C accumulation within diverse land uses in a context of global change (i.e., invasive, browse, CO 2 enrichment).
Natural and human disturbances such as wildfire, insect and disease outbreaks, weather, and timber harvesting could cause forest C uptake losses that potentially alter the directionality of C fluxes in forests [19,58,59]. Disturbance events may decrease forest AGB accumulation if annual disturbance rates increase every 1% across the Northern US (Tables 4 and 5). Except for growth enhancements, forest C stocks accrete mainly due to the net effect of disturbances and regrowth. Forest C stock accumulation under disturbances would be less than C stored from natural growth if disturbances led to additional C stock losses than the resulting offset of forest recovery and growth enhancements (e.g., density reduction of overstocked stands). Another possible explanation for such a decrease is that higher severity of natural disturbances could increase tree mortality and suppress tree growth in stands, which in turn may reduce forest AGB growth [16,60,61]. Moreover, some of the predicted C stock losses may be captured in harvested wood products (0.58 Pg C). Although natural disturbances may initiate a regeneration legacy that contributes considerably to the increased rate of forest C stocks [58], it might be more than offset by C sequestration potential losses induced by greater mortality compared with forest natural growth. However, disturbance scenarios with low/medium fire and low weather may lead into higher total forest AGB than no disturbances during some periods in the future simulations due to the low and/or mixed severity nature of such disturbances, which can create mortality in the suppressed cohorts of many stands. Such mortality can lead to increased growth of survivors in combination with the regeneration of new, potentially rapidly growing cohorts [62]. The sensitivity analysis of disturbances suggested that disturbances with higher severity would lead to less forest AGB sequestration potential across the study region. This is consistent with other studies that imply that the majority of forest C stocks decrease with greater harvest levels and increased disturbance rates [58,63]. Maintaining or increasing forest C stocks will be very challenging if disturbance percentages increase substantially over the next several decades [63,64]. Moreover, disturbances would not be the only process influencing forest AGB dynamics, as other global change factors, including increased atmospheric CO 2 and N deposition, may also be other factors to consider. Without accounting explicitly for these factors in this study, we propose that tremendous uncertainty surrounds the true risk of reductions in C sequestration rates. Additionally, although changes in the frequency and severity of disturbances are forecast to occur under global climate change, we did not study the influences of climate change and interactions between disturbance events and a warmer climate. In all, our study has addressed the effects of most types of disturbances on future forest AGB dynamics in the Northern US.
Although LUC and disturbances have collectively affected terrestrial C dynamics over time, few studies have attempted to characterize the interactive impact of simultaneous changes in LUC and disturbance processes on future forest AGB dynamics [7,20,22,34,54,65]. Our findings indicate that contemporary LUC substantially contributes to forest C sink, and forest C sequestration rates may increase as well under interactions of LUC and disturbances in Northern US forests. Combined medium LUC and disturbances could keep the region's forests as a C sink, which suggests that positive effects of medium LUC would not be offset by negative influences of medium disturbances. Under a 7.0% annual harvest rate (including partial cutting of large trees), we estimate a decrease of 4.26% of total forest AGB in the LUC + harvest scenario compared with a natural growth (scenario none). This is consistent with other studies that have focused on the potential for active LUC planning and altered forest management practices, which has been viewed as a potential opportunity to maintain regional C sink strength [5,63].

Uncertainty Estimation
Taking into consideration uncertainty arising from LUC and disturbances is useful and important in large-scale predictions of forest C dynamics [66] and is a requirement in UNFCCC reporting [26]. Long-term forest C simulations were subject to the bias caused by the changes in LUC (i.e., from non-forest to forest or from forest to non-forest). Disturbances could affect C sequestration, as C could rapidly move from living biomass to dead organic matter or to the atmosphere via combustion [67]. In this study, we only focused on parameter uncertainty in LUC and disturbances, which mostly stems from measurement errors, missing data, or unrepresentative data [68]. We used fuzzy sets to estimate the parameter uncertainty in terms of membership degree, as fuzzy set theory can provide a useful point of departure for the construction of a conceptual framework, which may prove to have a much wider scope of applicability [41]. To handle uncertainty of parameters in the forest AGB predictions, a two-step approach can be adopted. First, Monte Carlo simulations [69] could be used to produce random realizations of possible parameters, thus generating the average and relative standard deviations of outputs. Then, uncertainty in forest AGB predictions of LUC and disturbance scenarios with fuzzy sets could determine the level of output uncertainty caused by uncertain parameter inputs. The uncertainty range of outcomes related to LUC and disturbance scenarios showed LUC and disturbances to be perhaps the most important uncertainty regarding the strength of the forest C sink. Nevertheless, hindered by the complexity of modeling fundamental disturbance processes in forests, we only applied annual disturbance rates to explore impacts of disturbances on forest AGB dynamics every five years. More importantly, we only examined the possible ranges in forest AGB predictions as an index of uncertainty while not explicitly examining detailed regimes of harvest, weather, wildfire, insect and disease, or other disturbance regimes. Admittedly, this study only dealt with very limited sources of uncertainty. Exploring more detailed and deeper sources of uncertainty is appropriate to estimate the uncertainty associated with forest AGB dynamic predictions. As such, caution should be extended to applying this study's models in the context of individual drivers of C flux.

How Might We Increase the Forest AGB C Sink in the Northern US?
Modeling forest AGB C dynamics under LUC and disturbances is needed to help design policies increase forest C uptake in the US. In order to enhance regional AGB C sinks, development of long-term land use planning policy should explicitly recognize the importance of maintaining forestland use and maximizing reforestation/afforestation activities. Land use and forest dynamics should be simultaneously considered when projecting forest AGB C uptake across landscapes [5]. Forest AGB C operations should take into account warmer and more variable climates, as well as associated changes of disturbance severity and frequency [70]. Moreover, tremendous forest AGB C stocks still exist that need to be maintained, even as forest AGB C uptake may decrease in the future. Research efforts should focus on exploring and understanding the influence of changing disturbance regimes on post-disturbance forest conditions. Monitoring efforts should reach a clear identification of forest AGB C dynamics for different LUC and disturbance categories.

Future Work
Although we used matrix models combining plot-level physiographic variables (elevation and slope) with Landsat data to project future total forest AGB dynamics influenced by LUC and disturbance scenarios, we did not produce spatially explicit estimates regarding where risks and opportunities lie across the landscape for C accumulation. Additionally, we only considered forest AGB C gain/loss inducements from increases/decreases of forestland. In the future, we can directly obtain the physiographic variables from the digital elevation model (DEM) of the US Geological Survey (USGS). Dynamics of larger area AGB and other C pools (belowground biomass, soil C, dead wood C, and litter C) could be mapped using matrix models to allow for the inclusion of DEM and Landsat data to improve pixel-level prediction under LUC and disturbances. Transfers of C pools to and from each land use category could be considered as well. Such a development would allow forest C estimation to move from just being spatially explicit to spatially explicit and continuous while incorporating the temporal flexibility that comes with Landsat time series considering LUC.

Conclusions
There is a confluence of factors that may drive future trajectories of AGB C uptake in forests, especially when disturbances and LUC are considered. In the Northern US, when only considering LUC, total forest AGB of A&S→F and medium F→A&S + A&S→F would increase from 1.88 to 4.93 and 3.81 Pg C, respectively, by the end of the 21st century (year 2098). When only medium disturbances were considered, predicted forest AGB associated with forests subject to fire and harvest would increase from 1.88 to 2.42 and 1.90 Pg C, respectively, over the next 80 years. When both medium LUC and disturbances are considered, forest AGB predictions of LUC + fire, LUC + weather, LUC + insect and disease, and LUC + harvest would increase from 1.88 to 3.29, 3.10, 2.91, and 2.68 Pg C, respectively. At a medium level, an uncertainty analysis with fuzzy sets indicated that LUC + fire would lead to the highest forest AGB, and LUC + harvest would lead to the lowest forest AGB by the end of period. This study confirms that LUC is an important driver of the forest C sink in the Northern US, while disturbances may lend themselves toward C emissions compared with natural growth (scenario none), exclusive of impending global change threats such as invasives or CO 2 enrichment. Landsat data coupled with empirical field observations of forest C pools were used in this study, and the opportunity exists to more tightly link national C reporting requirements (i.e., 1990-baseline year reporting) with the need for spatially explicit estimates of future forest C dynamics as affected by LUC and disturbance scenarios. Finally, as several caveats (i.e., C emission rates per specific disturbance events) were applied in this study for the development of initial C trajectories (~80 years into the future), there are opportunities to explore alternative emission and removal factor estimates with emerging science.