Using matrix models to estimate aboveground forest biomass dynamics in the eastern USA through various combinations of LiDAR, Landsat, and forest inventory data

The ability to harmonize data sources with varying temporal, spatial, and ecosystem measurements (e.g. forest structure to soil organic carbon) for creation of terrestrial carbon baselines is paramount to refining the monitoring of terrestrial carbon stocks and stock changes. In this study, we developed and examined the short- (5 years) and long-term (30 years) performance of matrix models for incorporating light detection and ranging (LiDAR) strip samples and time-series Landsat surface reflectance high-level data products, with field inventory measurements to predict aboveground biomass (AGB) dynamics for study sites across the eastern USA—Minnesota (MN), Maine (ME), Pennsylvania-New Jersey (PANJ) and South Carolina (SC). The rows and columns of the matrix were stand density (i.e. number of trees per unit area) sorted by inventory plot and by species group and diameter class. Through model comparisons in the short-term, we found that average stand basal area (B) predicted by three matrix models all fell within the 95% confidence interval of observed values. The three matrix models were based on (i) only field inventory variables (inventory), (ii) LiDAR and Landsat-derived metrics combined with field inventory variables (LiDAR + Landsat + inventory), and (iii) only Landsat-derived metrics combined with field inventory variables (Landsat + inventory), respectively. In the long term, predicted AGB using LiDAR + Landsat + inventory and Landsat + inventory variables had similar AGB patterns (differences within 7.2 Mg ha−1) to those predicted by matrix models with only inventory variables from 2015–2045. When considering uncertainty derived from fuzzy sets all three matrix models had similar AGBs (differences within 7.6 Mg ha−1) by the year 2045. Therefore, the use of matrix models enabled various combinations of LiDAR, Landsat, and field data, especially Landsat data, to estimate large-scale AGB dynamics (i.e. central component of carbon stock monitoring) without loss of accuracy from only using variables from forest inventories. These findings suggest that the use of Landsat data alone incorporating elevation (E), plot slope (S) and aspect (A), and site productivity (C) could produce suitable estimation of AGB dynamics (ranging from 67.1–105.5 Mg ha−1 in 2045) to actual AGB dynamics using matrix models. Such a framework may afford refined monitoring and estimation of terrestrial carbon stocks and stock changes from spatially explicit to spatially explicit and spatially continuous estimates and also provide temporal flexibility and continuity with the Landsat time series.


Introduction
Forest ecosystems have the largest terrestrial carbon (C) stock and represent a majority (∼80%) of all aboveground C (Pacala et al 2001, Houghton et al 2009, Pan et al 2011. In concert with international efforts to reduce greenhouse gas (GHG) emissions, nations such as the US have been monitoring and mapping forest aboveground biomass (AGB) dynamics using data from the national forest inventory (NFI) conducted by the US Department of Agriculture (USDA) Forest Service, Forest Inventory and Analysis (FIA) program in the US (e.g. Heath et al 2011, Woodall et al 2015a. Recent work has highlighted that the accuracy and precision of AGB dynamics over a range of spatial scales could be further improved by incorporating vegetation indices as predictors derived from remotely sensed data that are temporally consistent and spatially continuous (Deo et al 2017a). Providing spatially continuous and temporally consistent estimates of AGB dynamics in forests at large-scales is crucial for evaluating existing land use policies and land management practices intended to help mitigate GHG emissions (Woodall et al 2015b, Deo et al 2017b. Antiquated approaches to predict AGB dynamics have been extensively used to provide robust estimates relying solely on forest inventory plots (Hall et al 2006, Ma et al 2018, but they are inefficient and difficult to implement over large spatial and temporal scales. Indeed the technologies of light detection and ranging (LiDAR) and timeseries Landsat surface reflectance high-level data products provide not only relatively cost-effective means, but also accuracy and spatial resolution to predict AGB dynamics over large domains of space and time (Powell et al 2010, Wulder et al 2012, Babcock et al 2018.
Generally, tree-and stand-level attributes from forest inventory data have been recognized as the most accurate method to provide information on biomass prediction at local or regional scales (Fang et al 2001, Brown 2002. However, these methods are usually labor-intensive and costly when applied to large regions (Houghton 2005) and are often logistically prohibitive in remote areas. The approach of linking ground measurements with remotely sensed data has greatly improved efficiency and cost-effectiveness of AGB estimation in forests, and can directly provide AGB estimates based on regression models (Running et al 1999, Hu et al 2016. Such a process has become increasingly popular with the advent of LiDAR and synthetic aperture radar (Lefsky et al 2002, Goetz and  . However, approaches for quantitatively predicting large-scale AGB dynamics in the longterm through incorporation of remotely sensed data, especially Landsat surface reflectance high-level data products and forest inventory plots are limited. To quantitatively evaluate AGB dynamics at national scales, it is crucial to develop a relationship between field inventory variables related to ABG, such as live tree basal area, and co-located spatial predictors extracted from LiDAR and Landsat time-series data (Margolis et al 2015). The relationship between inventory plot estimates and co-located auxiliary variables can then be exploited to achieve a reduction in error variance (when compared to an absence of auxiliary data) in the estimates of population dynamics (Powell et al 2010, Gregoire et al 2016. Because Landsat spectral data can characterize forest cover types (Hall et al 2006) and LiDAR provides highly visible three-dimensional profile of canopy structures, improved estimation of forest AGB can be attained when LiDAR data are combined with Landsat-derived fine resolution metrics (Hudak et al 2002, White et al 2016, Deo et al 2017b. However, it is not easy to implement more suitable AGB dynamic predictions to actual AGB dynamics for any model-assisted estimator (e.g. Forest Vegetation Simulator) that depends on field plots and remote sensing data being established as control variables (Babcock et al 2018).
Substantial improvement in computational environments and techniques have led to enhancements and flexibility in prediction models. For example, matrix models apply transition matrices to estimate the dynamics of ecological populations based on three main components including forest growth, mortality, and recruitment (e.g. Solomon et al 1987, Caswell 2001, Fieberg and Ellner 2001, Liang and Picard 2013. Since the 1940s (Lewis 1942, Leslie 1945, researchers have widely employed matrix models to study forest ecosystem dynamics (Usher 1969). Although matrix models vary in their degrees of complexity related to capturing forest growth, recruitment, and mortality, they have been applied to estimate forest population dynamics and associated forest C dynamics under different disturbance, management, and climate scenarios (Solomon et al 2000, Liang and Zhou 2014, Wear and Coulston 2015, Ma et al 2016, Ma and Zhou 2017. Therefore matrix models afford an opportunity to predict AGB dynamics through incorporation of forest inventory attributes with co-located remotely sensed variables. As signatories to the United Nations Framework Convention on Climate Change (UNFCCC), the United States (US) reports economy-wide GHG emissions and removals, which includes the land sector, each year from 1990 to near present. This report is intended to be a tool to evaluate the US contributions to global emissions but also to evaluate the effectiveness of policies and practices within the US (Hockstad and Hanel 2018). Given that most land use polices and land management practices occur at local or regional scales, having spatially and temporally resolved information is critical to evaluate existing policies and practices that can be used to inform future activities. The objective of this study was to evaluate whether remotely sensed variables, especially Landsat time series, can be used in place of forest inventory estimates (i.e. stand basal area [B]) in matrix models to predict AGB dynamics without substantial loss of accuracy and precision under uncertainty. The row and column vectors of the matrix were number of trees per unit area classified by field plot and by species group and diameter class. To test this we (1) developed matrix models using parameters from LiDAR strip samples, Landsat time-series, and forest inventory plots (i.e. stand basal area [B], elevation [E], plot slope [S] and aspect [A], and site productivity [C]) for specific study sites in the eastern US -Maine (ME), Minnesota (MN), Pennsylvania and New Jersey (PANJ), and South Carolina (SC), (2) predicted short-term stand basal area and long-term AGB changes among matrix models at each site and across sites, (3) adopted fuzzy sets to represent variability in predictions caused by uncertainty, and (4) compared the accuracy and precision of the predictions from each matrix model to evaluate whether remotely sensed variables could be used in place of field variables when such plots are missing across spatial and temporal domains of interest.

Research region
The research region is located in the eastern USA and contains a wide range of climatic conditions, ecoregions, and species composition, and covers four Landsat scenes corresponding to the WRS-2 paths/rows of 12/28, 27/27, 14/32 and 16/37. In which scenes 12/28, 27/27 and 16/37 lie completely within the ME, MN and SC while scene 14/32 covers parts of PANJ, respectively (figure 1).

Forest data
We included 322 remeasured permanent ground plots from the USDA, Forest Service, FIA database (previous measurement period was 2007-2010 with the same plots re-measured in 2012-2015, which resulted in an inconsistent time gap between plot measurements, USDA Forest Service 2017). The number of FIA plots was 96 in ME, 109 in MN, 42 in PANJ, and 75 in SC. The FIA plots were selected such that each represents a single condition and stand basal area is above zero. Each permanent ground plot is made up of four smaller fixed-radius (7.32 m) subplots spaced 36.6 m apart in a triangular arrangement with one subplot in the center (USDA Forest Service 2015). For model calibration and validation purposes, 258 permanent ground plots (80% of each state, 77 in ME, 87 in MN, 34 in PANJ, and 60 in SC) were used to establish the matrix model and the rest 64 permanent ground plots (20% of each state, 19 in ME, 22 in MN, 8 in PANJ, and 15 in SC) were randomly selected to test the matrix model's performance.
For each permanent ground plot, plot-level attributes included site productivity, trees per hectare, basal area per hectare, slope, aspect, and elevation. Tree-level data, including species, diameter at breast height (DBH), and status (recruitment, live, or dead), was also collected on each plot. Tree height data was not included in the study as only certain trees' height were collected in the FIA database. The four states (ME, MN, PANJ, and SC) were largely dominated by several species in terms of stand basal area. These species included: balsam fir (Abies balsamea), red spruce (Picea rubens), and northern white-cedar (Thuja occidentalis) in ME; quaking aspen (Populus tremuloides), balsam fir (Abies balsamea), and black spruce (Picea mariana) in MN; pitch pine (Pinus rigida), red maple (Acer rubrum), and white ash (Fraxinus americana) in PANJ; and loblolly pine (Pinus taeda), sweetgum (Liquidambar styraciflua), and red maple (Acer rubrum) in SC (table S1 is available online at stacks.iop.org/ ERL/13/125004/mmedia). For simplicity and computational efficiency, we classified all tree species into two species groups from the FIA program for each state: Deciduous and Coniferous (table S1). Within each species group, trees were further categorized into 15 DBH classes of 5 cm increments, except for the first class (2.54-7 cm) and the last (72 cm and above).

LiDAR and landsat data
LiDAR strip samples were acquired in summer (28 June-3 September) 2014 at the four sites and colocated Landsat time-series data, captured between 1984-2015, were obtained from the Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on-demand interface (http://espa.cr.usgs.gov/). The LiDAR strips covered about 2590, 3056, 4454 and 3571 square kilometers in 24, 24, 35 and 29 flight lines for ME, MN, PANJ, and SC, respectively. The average width of the strips was about 1.2 km. The LiDAR datasets were summarized to obtain the grid metrics for canopy cover and vertical strata density as defined in Deo et al (2017aDeo et al ( , 2017b. In addition, the Landsat time-series data were acquired by Landsat-5 Thematic Mapper (TM) and Landsat-8 Operational Land Imager (OLI) sensors at the peak of the growing seasons (July-September) from 1984-2015. The high-level products included surface reflectance of individual bands (visible and infrared), some derived spectral indices, and cloud masks for each image (https://landsat.usgs.gov/cdr-ecv) at 30 m spatial resolution (Deo et al 2017b).

Description of the matrix models
We used matrix models to predict AGB dynamics in this study. Recently, matrix models have been extended to account for forest management, disturbances and climate change due to their accuracy and robustness in depicting forest populations (see Liang and Picard 2013 and references therein). In this study, matrix models were used to control for tree growth, mortality and recruitment as follows (Liang et al 2005): in which y t and y t+1 are stand density (i.e. number of trees per unit area) vectors by species group and diameter class at time t and t +1, respectively. R t is a regeneration vector. ε is an error term. G t and G it are state-and time-dependent transition matrices describing how trees grow or die to model dynamics of forest population, 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 diameter class between t and t+1. The transition probability of upgrowth, 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 probability of tree mortality between t and t+1. R t is a state-and time-dependent recruitment vector representing the number of trees naturally recruited in the smallest diameter 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 diameter class j from t and t+1 is represented by the following model (Liang et al 2007; all notations defined in table 1): where Φ is the standard normal cumulative function, δ i 's are parameters estimated by maximum likelihood for mortality of species group i and diameter class j. ξ is an error term. Recruitment of species group i at time t, R it is estimated with a Tobit model (Tobin 1958): where Φ is the standard normal cumulative distribution function and j is the standard normal probability density function. The Tobit model explicitly accounts for unobserved recruitment values that are leftcensored at the preset diameter limit (2.54 cm). β i 's are parameters estimated for recruitment of species group i. μ is an error term. The AGB is estimated with the following model (all notations defined in table 1, Ma et al 2018): in which i y is a column vector of AGB for all plots, ω i 's are parameters to be estimated with the GLS for AGB. θ is an error term and i y is calculated with equation (9).
In order to explore whether remotely sensed variables, especially Landsat variables, can be used in place of forest inventory estimates in matrix models to predict AGB dynamics, LiDAR and Landsat variables were used to replace stand basal area in the above models, resulting in three types of matrix models for this study (all variables defined in table 1). The first is a matrix model using only forest inventory variables (matrix1), the second matrix model used LiDAR and Landsat variables to replace basal area of matrix1 (matrix2), and the third matrix model used only Landsat variables to replace basal area of matrix1 (matrix3) in ME, MN, PANJ and SC. The accuracy of all matrix models was examined using FIA inventory data. To avoid compromised type I error rates, we used hierarchical partitioning (HP) for the screening of explanatory variables. Explanatory variables were selected by the average independent contribution of each variable to the overall goodness-of-fit (Chevan andSutherland 1991, Mac Nally 2000). The HP method is useful since it can be applied to all regression methods such as ordinary least squares, probit, logistic, and loglinear regression (Chevan and Sutherland 1991). The HP analysis was conducted with the hier.part package of the R program (Nally and Walsh 2004).

Fuzzy sets representing uncertainty
Uncertainty inducing from disturbances led to high variability in predicted values of AGB. The averages of predicted AGB are important point estimations but to understand the associated risk, ranges or sets indicating uncertainty in predictions are essential. Here we used fuzzy sets which involved defining membership functions that determined the level of uncertainty (Zadeh 1965). Membership degree of an element, which may fall anywhere in the interval [0, 1], is expressed by a membership function. 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) (Zadeh 1965). a, b, c, d were four threshold values representing range of uncertainty. [b, c] represented the certainty interval for which the membership degree is 1. [a, b) and (c, d] were the uncertainty intervals with membership degrees ranging from 0 to 1. [a, d] was a measure of total range of uncertainty arising from disturbances. Following Weckenmann and Schwan (2001), given the average value of predicted AGB X ( ) and its relative standard deviation (S r ) from simulations, a, b, c, d values can be calculated as follows Zhou 2017, Ma et al 2018 Total returns in the vertical height interval of 1-4 m divided by 3 DenStra2 Total returns in the vertical height interval of 4-8 m divided by 4 DenStra3 Total returns in the vertical height interval of 8-16 m divided by 8 DenStra4 Total returns in the vertical height interval of 16-32 m divided by 16 DenStra5 Total returns in the vertical height interval of 32-64 m divided by 32 Stratum0 Total return proportion in the vertical height interval 0-1 m Stratum1 Total return proportion in the vertical height interval 1-4 m Stratum2 Total return proportion in the vertical height interval 4-8 m Stratum3 Total return proportion in the vertical height interval 8-16 m Stratum4 Total return proportion in the vertical height interval 16-32 m Stratum5 Total return proportion in the vertical height interval 32-64 m CovMean Percentage of all returns above mean (all returns above mean×100/total count of all returns) CovMode Percentage of all returns above mode per pixel (all returns above mode×100/total count of all returns) CovDBH Percentage of all returns above DBH per pixel (all returns above DBH×100/total count of all returns) ElevAAD Average absolute deviation of elevations of all returns above DBH ElevAv Average of elevations of all returns above DBH CRR Canopy relief ratio ((mean-min)/(max-min)) of elevations of all returns ElevCM Cubic mean of elevations of all returns ElevCV Coefficient of variation of elevations of all returns above DBH ElevIQ Interquartile range of elevations of all returns above DBH First L-moment of elevations of all returns above DBH ElevL2 Second L-moment of elevations of all returns above DBH ElevL3 Third L-moment of elevations of all returns above DBH ElevL4 Fourth L-moment of elevations of all returns above DBH ElevLCV L-moment coefficient of variation of elevations of all returns above DBH ElevLkurt L-moment kurtosis of elevations of all returns above DBH ElevLskew L-moment skewness of elevations of all returns above DBH EMADmed Median of the absolute deviations from the overall median of elevations EMADmod Mode of the absolute deviations from the overall mode of elevations ElevMax Maximum of elevations of all returns above DBH ElevMin Minimum of elevations of all returns above DBH ElevMod Mode of elevations of all returns above DBH

Integrated framework
It is important to comprehensively assess forest AGB dynamics using LiDAR, Landsat and forest inventory data under uncertainty using matrix models. The predictions of large-scale AGB dynamics may be further improved by including LiDAR and Landsat variables that are temporally consistent and spatially continuous facilitating estimation over large spatial domains. In this study, we developed an integrated matrix-modeling framework to project the forest AGB dynamics of eastern US forests under uncertainty with (1) the forest inventory measurements, (2) the LiDAR strip samples, (3) the Landsat time-series. We synchronously coupled the three datasets to predict stand basal area in the short-term and used remotely sensed variables, especially Landsat variables, to replace stand basal area to project AGB in the long-term with fuzzy sets under uncertainty (figure 2).   smallest DBH (D) while coniferous species in SC had the highest annual growth rate (g). Coniferous species in MN had the lowest annual growth rate (g) (table S3). The average interval between two inventories was five years.

Parameters of the three matrix models
The dependent variables (i.e. the average annual rates of growth, mortality, recruitment, and AGB) were estimated from field, LiDAR, and Landsat attributes using repeated measurements of 258 permanent ground plots in the four states ME, MN, PANJ, and SC (tables 2-9). The explanatory variables were selected based on statistical and biological significance. The primary control variables for the matrix1, stand basal area (B), slope (S), and site productivity (C) were all significant (α0.05) except for mortality models. The other variables elevation (E), and aspect (cos A and sin A) were mostly significant (α0.05) (tables 4, 6, 8). The rate of tree growth was significantly higher for lower stand basal areas, higher site productivity, lower elevation, and lesser slopes across all species groups in all states (table 4). Mortality rate declined significantly with site productivity, and increased significantly with stand basal area, elevation, and slope (table 6). Recruitment strongly declined with increasing stand basal area (table 8). In AGB models, greater stand basal area, lower elevation, and lesser slope had greater potential for AGB sequestration (table 10). Control variables from LiDAR and Landsat were used to replace stand basal area in the matrix1 including diameter growth, mortality, recruitment, and AGB models (matrix2 and matrix3). In the matrix2 and matrix3, starting from a large amount of attributes from both LiDAR and Landsat and only Landsat (table 1), the predictors were selected by applying stepwise regression for pruning collinear variables and multiple linear regressions at significance level α<0.05. Based on those selected variables (table 2), we further chose a subset in the final models of the matrix2 and matrix3 with adjusted R 2 (0.78-0.88 and 0.37-0.51, table 3) based on statistical and biological significance and contribution to the goodness-of-fit as determined through HP. In the matrix2, we selected the final predictors from both LiDAR and Landsat to replace stand basal area, such as canopy height model (CHM), total return proportion in the vertical height interval 1-4 m (Stratum1), and total return proportion in the vertical height interval 8-16 m (Stratum3) in the ME; 5th percentile of elevations of all returns above

Predictions and errors
For the 64 validation plots of the four states, the average stand basal area by species and diameter class predicted by the matrix1, matrix2, and matrix3 in the short term all fell within the 95% confidence interval of the observed values in all species-diameter classes, demonstrating high accuracy of the matrix models ( figure 3). Predicted basal area were well aligned with the mean observed values. Based on the root mean square error (RMSE) of basal area, the three matrix models with diameter class (matrix1, matrix2, and matrix3) had a high accuracy as well (figure 3). The matrix1 had the lowest RMSE (7.53-10.58 m 2 ha −1 ), while the matrix3 had the highest (9.34-11.68 m 2 ha −1 ). Based on the 30 years prediction from 2015 onward, average predicted AGB of the three matrix models increased over time and converged to ranges of 84-90, 65-69, 99-106, and 78-85 Mg ha −1 in the ME, MN, PANJ, and SC, respectively, at the year 2045 (figure 4). Average AGB increase per year would reach to ranges of 0.93-1.13, 1.03-1.18, 1.02-1.12, and 1.01-1.19 Mg ha −1 in the ME, MN, PANJ, and SC, respectively. The AGB predicted by the matrix2 and matrix3 all fell within the 95% confidence interval of the AGB predicted by the matrix1 ( figure 4).

Uncertainty analysis
To account for variability in the simulation results, fuzzy sets were constructed for three matrix models based on equation (10) (figure 5). Matrix1, matrix2, and matrix3 could lead to similar AGBs in the ME, MN, PANJ and SC, given the amount of overlapping among the fuzzy sets. In sum, the existing overlaps among three matrix models in the ME, MN, PANJ and SC (figure 5) suggested the possibility of similar disturbance effects on AGB for the three matrix models across study region when taking account of uncertainty.

Discussion
It is critically important to improve predictions and projections of C estimates over large forest areas by developing matrix models that accurately characterize C dynamics based on empirical forest inventories and remotely sensed data or only Landsat data. Importantly, such work can serve as a key refinement of the US National GHG Inventory for improvement of the consistency and accuracy of GHG emission and removal estimates from the forestland category of UNFCCC reporting. In this study, we developed and tested the performance of a matrix-modeling framework that can predict AGB dynamics with forest field data in combination with LiDAR strip samples and Landsat time-series, especially Landsat data (all fell within the 95% confidence interval of observed Table 4. Estimated parameters of the diameter growth models using variables from FIA plots.
Diameter growth models Note. Significance levels: * <0.05; ** <0.01; *** <0.001. Table 5. Estimated parameters of the diameter growth models using remotely sensed variables to replace basal area from both LiDAR and Landsat predictors and only Landsat predictors.
Diameter growth models values). In the matrix models, our study suggested that forestland with lower stand basal area, higher site productivity, flatter slopes, and lower elevation may have lower mortality rates, greater diameter growth rates, and greater recruitment generating a higher potential for forest carbon sequestration across our study sites. In the AGB simulation models, AGB density was found to be positively correlated with stand basal area, supporting the conclusion that forest AGB increases with forest stocking. Additionally, we also demonstrated the strength of LiDAR and Landsat dependent matrix models in their application to predicting regional AGB dynamics. Among the selected variables from both LiDAR and Landsat data, LiDAR metrics were much more influential to stand basal area and associated AGB than Landsat variables. This finding is consistent with previous studies (Deo et al 2017a(Deo et al , 2017b. Nevertheless, we found that Landsat variables alone can also be used to replace stand basal area to predict AGB in the study region although with higher RMSE (9.34-12 Mg ha −1 ). This finding may improve the ability to predict AGB dynamics across large scales with the Landsat time series. Model accuracy has been significantly influenced by spatial and temporal matching of field sampling and remote sensing data given they are site-and species-specific and scale-dependent (McRoberts 2010, Huang et al 2015, Deo et al 2017b. Furthermore, precise short-and long-term prediction is possible with these accurate models. Prior to this study, previous work had estimated AGB employing LiDAR , Zolkos et al 2013, Babcock et al 2015 and Landsat (Zheng et al 2004, Lu 2005, Hall et al 2006. However, more accurate predictions of site-and species-specific AGB dynamics that combined both LiDAR and Landsat with forest inventory were rare (Babcock et al 2018, Deo et al 2017b). Stand basal area was a statistically significant and frequently used predictor in the matrix model for diameter growth, mortality, and recruitment in the previous studies (Liang et al 2005, 2011, Liang 2010, Ma et al 2016, Ma and Zhou 2017. Hence, it was logical to assume that the matrix1 would perform better than the matrix2 and matrix3. This assumption was consistent with our results and the differences of short-and long-term prediction (differences within 7.2 Mg ha −1 ) and RMSE (8.34-12.12 Mg ha −1 ) among the three models were not obvious. This result suggested that stand basal area can be replaced by selected variables from both LiDAR and Landsat and only Landsat in the matrix models while the accuracy of AGB prediction was not highly impacted. This implies that the matrix models can be formulated by integrating field and remotely sensed data of closely matching resolutions, and subsequently provide accurate large-scale predictions of AGB dynamics. Therefore, it revealed that plot-level (e.g. E, S, and A) and remotely sensed information or only Landsat data might be used to project large-scale forest AGB dynamics without the detailed basal area calculations that would be more costly in terms of measurement and processing time. The physiographic variables (E, S, and A) can be directly obtained from the digital elevation model of the US Geological Survey (USGS). In the future, the dynamics of larger area AGB and other carbon pools (belowground biomass, soil carbon, dead wood carbon, litter carbon) should be investigated using matrix models to allow for the inclusion of more field and remote sensing samples, especially Landsat time series, to improve pixel-level prediction under land-use change and disturbances. Such refinements would enable large-scale GHG monitoring to move from only being spatially explicit to both spatially explicit and spatially continuous with the addedtemporal flexibility enabled by the Landsat time series. Table 6. Estimated parameters of the mortality models using variables from FIA plots.
Considering uncertainty is imperative in large-scale predictions of forest carbon dynamics (Miehle et al 2006) and is a requirement in UNFCCC reporting. Disturbances could affect carbon sequestration through rapid flux of carbon from living biomass to dead organic matters or via combustion to the atmosphere (Hicke et al 2012). In our study, we only examined the possible ranges in AGB predictions as an index of uncertainty (AGB differences of three matrix models within 7.6 Mg ha −1 at the year 2045) while not explicitly examining individual sources of uncertainty resulting from climate change, wildfire, wind damage, insect, disease, or other natural disturbances. As such, caution should be extended to applying this study's models in the context of individual drivers of carbon flux. How to combine all these sources of uncertainty into AGB predictions largely remains an open issue. Large-scale ecological modeling studies can gain credibility in the prediction of AGB dynamics following equally dynamic disturbance regimes through the use of suitable uncertainty analysis techniques. This study's matrix models made it possible to provide explicit AGB dynamics of eastern US forests. Prediction of AGB using LiDAR, Landsat, and forest inventory plots, which has been seldom studied but is key to local, regional, and national forest C monitoring programs, can now be analyzed by matrix models. Due to their simplicity and accuracy, the matrix models are also applicable to estimate forest population dynamics under climate change (Liang et al 2011), assessments of fire regime responses to changing climate (Ma et al 2016) and assessments of harvesting regimes under climate and fire uncertainty (Ma and Zhou 2017). In addition, the matrix models may also be incorporated with adaptive management of timber and biomass resources (Zhou et al 2008), land-use change (Woodall et al 2015b), and detailed regimes of natural disturbances (Wear and Coulston 2015). The presented model illustrates the benefit of including physiographic variables (E, S, and A) in the matrix models which can control for large-scale spatial variations in regions while still being able to achieve adequate prediction accuracy (Peterson et al 2014).
Matrix models maybe a convenient tool to predict forest C dynamics (Ma et al 2017). Not unlike other models, the matrix models have their limitations, and our results should be interpreted in appropriate contexts (Liang et al 2011). One weakness of our study's matrix models is that their geographic coverage was limited to the study region in the eastern US. Another limitation is that the models were unable to reflect land-use change (i.e. from non-forest to forest or from forest to non-forest) and major natural disturbances such as wildfires, insects, hurricanes, and widespread diseases. Therefore, their long-term simulations were subject to the bias caused by the changes in land use, if there was any. However, the matrix models can provide a useful baseline tool for the estimation of future AGB dynamics across the entire reporting period from 1990 to the present. Future work should focus on addressing the impacts of land-use change and detailed regimes of natural disturbances on AGB dynamics through incorporating the present models with a transition matrix of land use change (Woodall et al 2015b) and stochastic elements of disturbances Buongiorno 2006, Wear andCoulston 2015) at a national scale. More importantly, the matrix models can be used in the future to map spatially continuous and temporally explicit forest population dynamics across scales with a geospatial matrix over long-term periods (Liang 2012, Liang andZhou 2014).   Table 9. Estimated parameters of the recruitment models using remotely sensed variables to replace basal area from both LiDAR and Landsat predictors and only Landsat predictors.

Conclusion
In the short-term, average stand conditions predicted by three matrix models considered in this study were well aligned with the mean observed values in the eastern USA-ME, MN, PANJ, and SC and consistent with the broad differences in forest ecological conditions between these regions (e.g. greatest diameter growth in SC, highest recruitment in ME). Based on the 30 years prediction, predicted AGB using matrix models incorporating both LiDAR and Landsat attributes and only Landsat attributes to replace stand basal area in the field sample variables (LiDAR+Landsa-t+inventory and Landsat+inventory) had similar patterns of the AGB predicted by the matrix models with all forest inventory variables (inventory). After making comparisons of predictions of AGB dynamics among inventory, LiDAR+Landsat+inventory, and Landsat+inventory, we found that stand basal area of the matrix models can be replaced by selected variables from LiDAR and Landsat in terms of predicting AGB dynamics over short-and long-term periods without loss of accuracy. Similar predicted AGBs were found in the fuzzy sets using inventory, LiDAR + Landsat + inventory, and Landsat + inventory when taking account of uncertainty. Therefore we concluded that remotely sensed attributes, especially for Landsat variables, can be used to replace stand basal area in the matrix models over our study's time period and domain. In addition, Landsat data alone incorporating elevation, plot slope and aspect, and site productivity is able to produce useful estimates of AGB dynamics using matrix models. More importantly, the ability to predict AGB dynamics over large scales was improved using Landsat time series which suggests the ability to move beyond the spatially explicit NFI plots to all pixels within the study regions. Such a refinement may move GHG monitoring from being only spatially explicit to being spatially explicit and spatially continuous estimates while also providing temporal flexibility with the Landsat time series.
Technical advances such as this may pave the way for researchers to start quantifying biomass dynamics beyond that observed on plots at a particular measurement date towards assessing the interaction between disturbances, land use change, and other activities may influence those dynamics (i.e. attribution) across space and time. Given the availability of Landsat variables, the opportunity exists to further refine quantification of AGB dynamics over large spatial domains in a relatively cost-effective manner in consideration of disturbances and land use change while enabling valid comparisons with domains that have limited field data.