Determinants of Above-Ground Biomass and Its Spatial Variability in a Temperate Forest Managed for Timber Production

The proper estimation of above-ground biomass (AGB) stocks of managed forests is a prerequisite to quantifying their role in climate change mitigation. The aim of this study was to analyze the spatial variability of AGB and its uncertainty between actively managed pine and unmanaged pine-oak reference forests in central Mexico. To investigate the determinants of AGB, we analyzed variables related to forest management, stand structure, topography, and climate. We developed linear (LM), generalized additive (GAM), and Random Forest (RF) empirical models to derive spatially explicit estimates and their uncertainty, and compared them. AGB was strongly influenced by forest management, as LiDAR-derived stand structure and stand age explained 80.9% to 89.8% of its spatial variability. The spatial heterogeneity of AGB varied positively with stand structural complexity and age in the managed forests. The type of predictive model had an impact on estimates of total AGB in our study site, which varied by as much as 19%. AGB densities varied from 0 to 492 ± 17 Mg ha−1 and the highest values were predicted by GAM. Uncertainty was not spatially homogeneously distributed and was higher with higher AGB values. Spatially explicit AGB estimates and their association with management and other variables in the study site can assist forest managers in planning thinning and harvesting schedules that would maximize carbon stocks on the landscape while continuing to provide timber and other ecosystem services. Our study represents an advancement toward the development of efficient strategies to spatially estimate AGB stocks and their uncertainty, as the GAM approach was used for the first time with improved results for such a purpose.


Introduction
Forests play an important role in mitigating climate change by fixing about 30% of fossil fuel CO 2 (C) emissions [1][2][3]. Temperate forests constitute about 25% of the world forest area and store some 14% of the total forest carbon [1,4]. At a global scale, 44% of the world forests are under some sort of management regime [5], and forest management has been a very significant factor contributing to the current terrestrial carbon sink of 2.4 Pg C, accounting for about 0.7 Pg C annually [1,3]. Managed forests are critically important to assisting with reducing atmospheric CO 2 concentration through initiatives like Reducing Emission from Deforestation and Forest Degradation, which includes the conservation and sustainable management of forests and enhancement of forest carbon stocks (REDD+) [6]. To implement the REDD+ scheme, assessing forest biomass stocks with a reliable accuracy is fundamental for project planning, and supports monitoring, reporting, and verification (MRV). Forest carbon stocks are usually calculated using AGB since typically 50% of AGB is carbon [7]. Incomplete information regarding the spatial variability (spatial distribution) of forest AGB introduces substantial uncertainty about current estimates of the carbon stocks present in managed forests and the global carbon budget [8,9]. Therefore, it is important to improve our knowledge of the spatial variability of forest biomass and its associated uncertainty in order to: (a) improve estimates of carbon stocks and understand carbon cycles [10]; (b) support future climate mitigation actions [11]; and (c) support sustainable forest management [12].
While timber production (usually evaluated as wood volume [13]) used to be the primary forest management's objective (about 1.53 billion m 3 per year in 2011) [14], some studies recognize the role of managed forests in contributing to carbon sequestration [15]. Forest management may enhance C sequestration in managed forest ecosystems [16,17] and optimize sink strength by using biofuels as substitutes for fossil fuels [3,18,19], while continuing to meet the timber demand of society [15]. There are mitigation opportunities involving forestry activities that could increase the gross terrestrial C uptake from roughly 4.0 to 6.2 Pg C annually [3]. This suggests the importance of assessing the AGB and timber production in managed forests to evaluate their contribution to climate change mitigation.
The AGB spatial variability and associated uncertainty, as well as the methods used to minimize it, have attracted considerable attention [20]. An uncertainty map at the pixel level is valuable for the interpretation of the AGB map [21,22] and for reporting emissions. The development of remote sensing technology combined with intensive site-based inventory approaches for estimating forest biomass stocks and forest area changes, are important for addressing monitoring problems and implementing REDD+ [23]. The combination of remote sensing and field measurements contributes to a robust and transparent MRV system.
The carbon pools in managed forests are affected by different drivers, mostly by forest management regimes and natural disturbances [24][25][26][27], the species composition of forests and forest type [28,29], plus the age structure [30,31]. AGB and carbon in forest ecosystems are also strongly affected by climate [24,32,33], as well as soil properties and topography [34][35][36]. However, it is difficult to clearly attribute the effect to forest management or other environmental factors [3]. Therefore, understanding the factors determining AGB spatial variability at the landscape level is a first step for properly estimating AGB through a prediction method.
The prediction method has a substantial effect on the accuracy of AGB estimates. Linear models are among the most popular choices for estimating AGB because of their simplicity and interpretability when compared to other models, such as random forest methods and Gaussian processes [37]. Random forest methods have performed well in many cases, while stepwise linear models are often the least effective, judging by r 2 and RMSE with cross-validation [38]. Other model approaches such as the generalized additive models (GAMs) may describe a complex relationship between the response variable and its predictors. This is especially useful in research fields such as ecology or biology in which simple models cannot capture the structure of the data and more complex models may be required [39]. However, this approach has not been used to assess the spatial variability of AGB. In addition, LiDAR technology has proven to be a valuable tool in mapping the spatial variability of AGB [21]. This type of high-resolution data offers great improvement for AGB estimation and can reduce its uncertainty in accordance with the Intergovernmental Panel on Climate Change (IPCC)'s Tier 3 for the forestry sector [20,38].
In this study, we used three types of models to characterize the spatial variability of AGB in managed forests in central Mexico. Mexico is covered by 65 million hectares of forests [5,40] with a variety of forest ecosystems (i.e., 51% temperate forests, 49% tropical forests) [41] due to its geographical location and topographic heterogeneity. Forests in Mexico are important for their current and potentially increasing accumulation of C in their biomass [42][43][44]. However, there is a lack of information on the spatial variability of AGB and the determinants that affect it, especially in managed forests (11% of total forested area within Mexico) [40,41]. Therefore, the objectives of this study were: (1) to analyze the determinants of AGB and its spatial variability in a managed forest landscape for timber production; (2) to evaluate the effect of different models (using the same input variables) for predicting AGB; and (3) to estimate the AGB estimation uncertainty across the forest landscape. We expected that AGB is driven by stand structure, climate, and topography; and that the prediction method has a substantial effect on the AGB estimates and uncertainty. This work is novel as it is the first known study that uses a GAM model to map AGB, which is compared to other commonly used models.

Study Site
We conducted the present study in a Mexican temperate forest in the Zacualtipán forested region, in the state of Hidalgo (20 • Figure 1). This site is a member of the Mexican Network of Intensive Carbon Monitoring Sites (MEX-SMIC, its acronym in Spanish) [45,46], and the MexFlux network [47]. From the 1980 s, communities that own these forests started Forest Management Programs for timber production, generating even-aged stands of Pinus patula Schltdl. et Cham. through the seed-tree regeneration method. Natural regeneration is also combined with artificial reforestation (enrichment planting) when natural regeneration is insufficient. Thinning is the most common silvicultural treatment applied for producing wood for commercial purposes, decreasing competition between the remaining trees, and controlling tree species composition. Mountain cloud and pine-oak forest that is 80 years old, with no silviculture activities applied (designated as reference forest), represents 31% of the study area. The predominant climate is humid with summer rains averaging 2050 mm and an average annual temperature of 13.5 • C [48], both varying according to elevation and aspect. The typical phaeozem haplic soils (Hh) are rich in organic matter, and in the areas with steep slopes, calcaric regosol (Rc) soils are found [49]. variety of forest ecosystems (i.e., 51% temperate forests, 49% tropical forests) [41] due to its geographical location and topographic heterogeneity. Forests in Mexico are important for their current and potentially increasing accumulation of C in their biomass [42][43][44]. However, there is a lack of information on the spatial variability of AGB and the determinants that affect it, especially in managed forests (11% of total forested area within Mexico) [40,41]. Therefore, the objectives of this study were: (1) to analyze the determinants of AGB and its spatial variability in a managed forest landscape for timber production; (2) to evaluate the effect of different models (using the same input variables) for predicting AGB; and (3) to estimate the AGB estimation uncertainty across the forest landscape. We expected that AGB is driven by stand structure, climate, and topography; and that the prediction method has a substantial effect on the AGB estimates and uncertainty. This work is novel as it is the first known study that uses a GAM model to map AGB, which is compared to other commonly used models.

Study Site
We conducted the present study in a Mexican temperate forest in the Zacualtipán forested region, in the state of Hidalgo (20°37′4978′′ N and 98°37′51.01′′ W and 20°35′18.74′′ N and 98°35′23′′ W), central Mexico ( Figure 1). This site is a member of the Mexican Network of Intensive Carbon Monitoring Sites (MEX-SMIC, its acronym in Spanish) [45,46], and the MexFlux network [47]. From the 1980′s, communities that own these forests started Forest Management Programs for timber production, generating even-aged stands of Pinus patula Schltdl. et Cham. through the seed-tree regeneration method. Natural regeneration is also combined with artificial reforestation (enrichment planting) when natural regeneration is insufficient. Thinning is the most common silvicultural treatment applied for producing wood for commercial purposes, decreasing competition between the remaining trees, and controlling tree species composition. Mountain cloud and pine-oak forest that is 80 years old, with no silviculture activities applied (designated as reference forest), represents 31% of the study area. The predominant climate is humid with summer rains averaging 2050 mm and an average annual temperature of 13.5 °C [48], both varying according to elevation and aspect. The typical phaeozem haplic soils (Hh) are rich in organic matter, and in the areas with steep slopes, calcaric regosol (Rc) soils are found [49].

Permanent Sampling Plots
The intensive carbon monitoring site was established during 2012-2013 across a 3 × 3 km land area, adapting guidelines developed for similar monitoring sites in the U.S. [50]. A flux tower was instrumented in the center of the polygon in order to provide direct estimates of energy, carbon dioxide, and water fluxes. We established 24 permanent inventory plots following a systematic pattern and 27 plots randomly located following a gradient of harvesting years to spatially sample a forest recovery chronosequence, including a reference forest. The additional plots enabled a more effective representation of the variability of ecosystem carbon stocks (live biomass, dead organic matter, and soils) across the managed forest landscape. Sampling protocols were consistent with the design of the Mexican national forest inventory (Inventario Nacional Forestal y de Suelos-INFyS) [51], but included many more variables generally following the standard protocols described in Hoover [52]. We delimited 11.28 m radius subplots (400 m 2 area) to measure trees with a diameter at breast height (DBH) of ≥5 cm. Additionally, 5.04 m radius microplots (80 m 2 ) were employed to measure trees with 2.5 ≤ DBH < 5 cm and taller than 1.3 m. Understory biomass (herbs, shrubs, and trees < 2.5 cm in diameter) was destructively sampled at two circular microplots of 1 m 2 per subplot (two on each of the secondary sampling units 2, 3, and 4) ( Figure 2).

Permanent Sampling Plots
The intensive carbon monitoring site was established during 2012-2013 across a 3 × 3 km land area, adapting guidelines developed for similar monitoring sites in the U.S. [50]. A flux tower was instrumented in the center of the polygon in order to provide direct estimates of energy, carbon dioxide, and water fluxes. We established 24 permanent inventory plots following a systematic pattern and 27 plots randomly located following a gradient of harvesting years to spatially sample a forest recovery chronosequence, including a reference forest. The additional plots enabled a more effective representation of the variability of ecosystem carbon stocks (live biomass, dead organic matter, and soils) across the managed forest landscape. Sampling protocols were consistent with the design of the Mexican national forest inventory (Inventario Nacional Forestal y de Suelos-INFyS) [51], but included many more variables generally following the standard protocols described in Hoover [52]. We delimited 11.28 m radius subplots (400 m 2 area) to measure trees with a diameter at breast height (DBH) of ≥5 cm. Additionally, 5.04 m radius microplots (80 m 2 ) were employed to measure trees with 2.5 ≤ DBH < 5 cm and taller than 1.3 m. Understory biomass (herbs, shrubs, and trees < 2.5 cm in diameter) was destructively sampled at two circular microplots of 1 m 2 per subplot (two on each of the secondary sampling units 2, 3, and 4) ( Figure 2).

Figure 2.
Plot design used for forest measurement. Plot design is composed of four grouped circular subplots. Subplot 1 is the center of the cluster, with the other three plots located 45.15 m away at azimuths of 360°, 120°, and 240°, respectively. Understory biomass was sampled at 1 m 2 microplots, tress with DBH ≤ 2.5 were sampled in a central circular plot of 80 m 2 , and trees with DBH ≥ 5 cm were measured in the circular plot of 400 m 2 . 1 m 2 microplots are not to scale.

Field Measurement and AGB Estimates
Field data consisted of individual tree information recorded during 2013. We used 204 subplots of 400 m 2 , 204 microplots of 80 m 2 , and 240 microplots of 1 m 2 in order to quantify the AGB defined as all living vegetation, both woody and herbaceous, above the soil [20]. All plots were geo-referenced on the ground with a GPS Garmin GPSmap ® 62 s [53]. We tagged all trees, identified species, and measured the DBH at a 1.3 m height of live trees ≥ 5 cm diameter in the 400 m 2 plots, and trees with 5 > DBH > 2.5 cm and taller than 1.3 m in the 80 m 2 subplots. When irregularities (e.g., bole irregularities, buttress roots) were present at the sampling height, we took DBH measurements 5 cm

Field Measurement and AGB Estimates
Field data consisted of individual tree information recorded during 2013. We used 204 subplots of 400 m 2 , 204 microplots of 80 m 2 , and 240 microplots of 1 m 2 in order to quantify the AGB defined as all living vegetation, both woody and herbaceous, above the soil [20]. All plots were geo-referenced on the ground with a GPS Garmin GPSmap ® 62 s [53]. We tagged all trees, identified species, and measured the DBH at a 1.3 m height of live trees ≥ 5 cm diameter in the 400 m 2 plots, and trees with 5 > DBH > 2.5 cm and taller than 1.3 m in the 80 m 2 subplots. When irregularities (e.g., bole irregularities, buttress roots) were present at the sampling height, we took DBH measurements 5 cm above the irregularity [54]. We also recorded the total height of each tree. The total number of plots amounted to 204, where we examined and measured 11,015 trees.
We measured understory biomass (herbs, shrubs, and trees < 2.5 cm) in 240 microplots of 1 m 2 by destructive sampling, took sub-samples to the laboratory, and oven-dried them at 70 • for dry weight determination. The chronosequence in managed forest covered 31 years. More than 31 units (stands) of different ages were identified from one to 31 years old. The stand age was defined as the time since harvesting. It was obtained during the field survey and from forest management plans. GIS-based forest inventory polygons were used to extract the age for each stand in the complete area. The reference forest showed no recent harvesting or disturbance, and has an estimated age of 80 years [55]. We used four to 20 plots per stand age according to the stand area ( Figure 1), except for the 31-year stand, where only one plot was used. We estimated AGB (oven dry mass) by applying local species-specific allometric equations [55] (Table 1). For each plot, we computed the AGB per plot by summing all sample tree biomass estimates. Next, we obtained the biomass values in megagrams per hectare (Mg ha −1 ) considering the size that each plot or subplot represented (0.04 and 0.008 ha). We estimated understory biomass per 1 m 2 microplot by using the ratio dry weight:total fresh weight and extrapolated to biomass per hectare (ha). Additionally, for each plot, we computed the basal area per hectare (m 2 ha −1 ) and tree density (N ha −1 ).

LiDAR Data
We acquired LiDAR returns and feature heights in ASPRS LAS file format (American Society for Photogrammetry and Remote Sensing, LASer file format) and the Digital Terrain Model, as LiDAR-derived bare earth elevation (m, EGM96 geoid) in raster data (GeoTIFF) at a nominal 1 m spatial resolution, through the Goddard's Lidar, Hyperspectral, & Thermal imager (G-LiHT) webpage [57]. LiDAR data was collected in April 2013 in an attempt to avoid clouds associated with the rainy season [58]. We processed LiDAR data using FUSION software version 3.60+ [59]. We calculated LiDAR metrics [59] (Table 2), the mean, and the standard deviation values from the cloud of points within all the study sites in the 3 × 3 km landscape at a 5 m resolution. We used a threshold of 2 m as a minimum height above ground to reduce the noise within the near-ground cloud of returns caused by low vegetation and imperfections of the ground. A canopy threshold height (height break above the ground) of 3.0 m was used to compute LiDAR canopy cover metrics according to the specific characteristics of the study site [53]. We used the LiDAR metrics as predictor variables in the models for estimating the spatial variability of AGB.

Climatic and Topographic Data
Since the sampled plots did not have climate data at a fine spatial resolution, we gathered meteorological data from recording stations in the region. First, we resampled the 1 m DTM into a 5 × 5 m pixel size grid by using the downscaling method of bicubic spline [60] implemented in the System for Automated Geoscientific Analysis (SAGA) software v2.2 [61]. Then, we used R programming language [62] to fit annual mean temperature (AT) and annual total precipitation (AP) models by using a LiDAR-derived digital terrain model (DTM). Temperature was strongly correlated with latitude and elevation, while precipitation was better correlated with latitude and longitude.
We estimated the AT and AP by applying the fit equations at the pixel level (using latitude, longitude, and elevation) for the entire study area. We obtained topography parameters by using the LiDAR-derived DTM (1 m resolution) [57,61]. SAGA software was used for the analysis of spatial data, and to derive the primary topographic metrics [63]: Analytical Hillshading, Slope, Aspect, Plan Curvature, Profile Curvature, Convergence Index, Closed Depressions, Total Catchment Area, Topographic Wetness Index, LS-Factor, Channel Network Base Level, Channel Network Distance, Valley Depth, and Relative Slope Position, according to Wilson [64] (Table 2).

Statistical Models and Analysis
All data were harmonized in a geographic information system and we applied a clipping process from the gridded information (stand age, LiDAR metrics, climatic and topographic parameters) at the 198 plots (400 m 2 ) with SAGA [61].
We developed several kinds of models to fit our data for predicting the AGB density (Mg ha −1 ) across the landscape (spatial variability of AGB) from the set of predictor variables ( Table 2): empirical statistical models, a multiple linear model (LM), generalized additive models (GAM) [65], and Random Forest (RF) [66] using 10-fold repeated cross-validation [67]. All the analyses were conducted in the R programming language [62]. We analyzed all possible combinations between variables to understand which predictors had more weight in the explanation of the spatial variability of AGB and thereby to interpret the main significant effects.
We used LM to identify highly-correlated variables and determine the best-fit. In order to find a model with the greatest explanatory power, a backward stepwise multiple linear regression was performed to automate the selection of the best explanatory variables. Multiple regressions are often affected by overfitting and co-linearity among variables. Multicollinearity was tested by the variance inflation factor (VIF) [68]. We removed the variables with a VIF above 5 from the model. Bootstrap samples were used from the data to implement the 10-fold cross-validation analysis.
We fitted GAM, a semi-parametric approach [65], for predicting the nonlinear responses of biomass to the set of predictor variables. GAMs are regression models that generalize the family of generalized linear models (GLMs), by replacing the linear functional form with a sum of smooth functions [69,70]. GAMs have been strongly accepted in several domains as a flexible modeling technique, suited for capturing nonlinear, unspecified relationships between predictor variables and the response variable [71,72]. GAM has been used to explore the AGB [73], but not to map its spatial variability. In this study, we performed GAM regression assuming a Gaussian error distribution. GAM regression was carried out using the R packages mgcv [74], and gamclass [75] to implement the 10-fold cross-validation. Stepwise selection using the Akaike information criterion (AIC) [76] was used to control the complexity of the GAM model.
The RF algorithm [66] was used as the other prediction model. RF is a machine learning algorithm that has been proven to reduce bias and overfitting. It works well with random inputs, and in some cases, tends to be more accurate than simple regression techniques for biomass estimation [37]. We used the rfcv function within the randomForest package [77] for the 10-fold cross validation analyses. The search of the covariates was pursued with the tune.randomForest function in the R e1071 [78] and VSURF [79] packages, with 1000 iterations of bootstrapping to improve performance estimates and control over-fitting.
For accuracy assessment, we used the 10-fold repeated cross-validation [67] (each model was re-fitted 10 times using 90% of the data and predictions derived from the fitted models were compared with observations of the remaining 10%) [80]. We derived the coefficient of determination (R 2 ) and root mean squared error (RMSE) to quantify the goodness of fit of the models. Those with the lowest AIC [76] root mean square error (RMSE) and the highest R 2 (adjusted for number of predictors in the model) were considered [81]. Finally, we applied spatial prediction models using tiled raster stacks (covariates) to map the AGB by using the R raster package [82].

Uncertainty Quantification
To estimate the potential uncertainty of AGB predictions, we used the difference of 95% prediction intervals (0.925-0.025) to calculate the absolute values of AGB uncertainty (Mg ha −1 ) per pixel by LM and GAM methods. For a normal distribution, this interval would be about ± 2 standard deviations about the mean. This uncertainty is limited to prediction error and serves as a reference to which the uncertainty due to model selection can be compared. This uncertainty represents the amount of AGB that cannot be explained for the predictor variables.
To estimate the uncertainty of the RF AGB predictions, we used the Quantile Regression Forests (QRF) approach to infer a pixel-based full conditional distribution [83] of the AGB response to predictors given available data. The statistical distribution of predicted values per pixel, estimated as the 95% prediction intervals, was obtained using the quantreg Forest package in R [84]. This technique allows a non-parametric and alternative way of estimating the reliability of AGB predictions.

AGB and Effects of Forest Structure
AGB of trees >5cm varied from 0.49 ± 0.66 to 166.6 ± 42.4 Mg ha −1 in stands ranging from three to 28 years old (Figure 3). Some recently harvested stands have up to 30 Mg ha −1 because of the seed trees remaining which are removed after regeneration is established. Biomass in trees of diameter class 2.5-5cm DBH ranged from 0 to 4.6 ± 1.3 Mg ha −1 , increasing with stand age, but not significantly. Understory biomass varied from 1.1 ± 0.7 to 12.8 ± 7.4 Mg ha −1 , but did not show a trend with increasing stand age. The reference forest had an average AGB of 159.5 ± 67.4 Mg ha −1 , which was not statistically different (p > 0.05) from managed older stands (20-28 years). However, it is noteworthy that the old-growth stand reached a maximum AGB of 287.3 compared to 223.4 Mg ha −1 in some managed older stands (Figure 3). Tree density and basal area were positively correlated with the AGB (p < 0.0001). They explained up to 90% of its variability. However, these forest management variables are not usually mapped.
LiDAR-derived forest metrics were strongly correlated with AGB. A canopy height metric (the 95th percentile height for all returns above 2 m) (eP95), and a canopy density metric (All returns above 3 m/first returns) × 100) (a1a3), were the two most significant predictors of AGB (p < 0.0001). Both forest structure variables showed positive effects in all the models (p < 0.0001) (Figure 4). managed older stands (Figure 3). Tree density and basal area were positively correlated with the AGB (p < 0.0001). They explained up to 90% of its variability. However, these forest management variables are not usually mapped. LiDAR-derived forest metrics were strongly correlated with AGB. A canopy height metric (the 95th percentile height for all returns above 2 m) (eP95), and a canopy density metric (All returns above 3 m/first returns) × 100) (a1a3), were the two most significant predictors of AGB (p < 0.0001). Both forest structure variables showed positive effects in all the models (p < 0.0001) (Figure 4).

Climate and Topographic Variables
There were spatial patterns in the distribution of climatic values: precipitation increased from the southwest to the northeast, while temperature showed the opposite pattern. These variables did not individually influence the AGB, but the interaction between stand age and precipitation had a positive effect (p < 0.0001). Topographic factors, which mostly influenced temperature and  managed older stands (Figure 3). Tree density and basal area were positively correlated with the AGB (p < 0.0001). They explained up to 90% of its variability. However, these forest management variables are not usually mapped. LiDAR-derived forest metrics were strongly correlated with AGB. A canopy height metric (the 95th percentile height for all returns above 2 m) (eP95), and a canopy density metric (All returns above 3 m/first returns) × 100) (a1a3), were the two most significant predictors of AGB (p < 0.0001). Both forest structure variables showed positive effects in all the models (p < 0.0001) (Figure 4).  . Effect for each factor on AGB by the generalized additive model (GAM). Canopy height (the 95th percentile height for all returns above 2 m), and canopy density (All returns above 3 m × 100/first returns) are metrics derived from LiDAR.

Climate and Topographic Variables
There were spatial patterns in the distribution of climatic values: precipitation increased from the southwest to the northeast, while temperature showed the opposite pattern. These variables did not individually influence the AGB, but the interaction between stand age and precipitation had a positive effect (p < 0.0001). Topographic factors, which mostly influenced temperature and

Climate and Topographic Variables
There were spatial patterns in the distribution of climatic values: precipitation increased from the southwest to the northeast, while temperature showed the opposite pattern. These variables did not individually influence the AGB, but the interaction between stand age and precipitation had a positive effect (p < 0.0001). Topographic factors, which mostly influenced temperature and precipitation patterns, had a limited effect on AGB. Slope and relative slope position as determined by LiDAR-derived DTM had a negative effect, while valley depth had a positive effect, on AGB (VIF < 1.2), explaining up to 11.4% of AGB variability (p < 0.0001).

Spatial Prediction Models
The linear model explained 80.9% (AIC = 1782, RMSE = 25.7) of the spatial variability of AGB (Table 3). LIDAR-derived variables related to forest structure, which are the 95th percentile height for all the returns above 2 m (eP95), canopy density (all the returns above 3 m/first returns) × 100) (a1a3), and stand age (management variable), were the three most significant predictors of AGB (p < 0.0001). The optimal GAM model explained 89.8% (AIC = 1712, RMSE = 17) of the variability in AGB, having forest structure variables (height, density, and age) as the most significant predictors. This semiparametric approach confirmed the importance of forest management variables (eP95, a1a3, and stand age) ( Table 3) because the management actions directly define them. The RF model also consistently showed forest management variables as the most important, explaining 88.7% (RMSE 19.6) of the AGB spatial variability.

Spatial Variability of AGB
AGB varied across the forested landscape, driven mainly by forest structure and stand age ( Figure 5). The three models overestimated AGB in the plots where measured biomass was very low or near zero, but they tended to underestimate AGB for plots with very high biomass measurements ( Figure 6). LiDAR metrics performed well in predicting AGB variability. The prediction methods have an impact on AGB estimates. GAM predicted the highest values of AGB, while RF predicted more conservative values (Figure 7), with a 19% difference between methods.
Forests 2018, 9, x FOR PEER REVIEW 9 of 20 precipitation patterns, had a limited effect on AGB. Slope and relative slope position as determined by LiDAR-derived DTM had a negative effect, while valley depth had a positive effect, on AGB (VIF < 1.2), explaining up to 11.4% of AGB variability (p < 0.0001).

Spatial Prediction Models
The linear model explained 80.9% (AIC = 1782, RMSE = 25.7) of the spatial variability of AGB (Table 3). LIDAR-derived variables related to forest structure, which are the 95th percentile height for all the returns above 2 m (eP95), canopy density (all the returns above 3 m/first returns) × 100) (a1a3), and stand age (management variable), were the three most significant predictors of AGB (p < 0.0001). The optimal GAM model explained 89.8% (AIC = 1712, RMSE = 17) of the variability in AGB, having forest structure variables (height, density, and age) as the most significant predictors. This semiparametric approach confirmed the importance of forest management variables (eP95, a1a3, and stand age) ( Table 3) because the management actions directly define them. The RF model also consistently showed forest management variables as the most important, explaining 88.7% (RMSE 19.6) of the AGB spatial variability.

Spatial Variability of AGB
AGB varied across the forested landscape, driven mainly by forest structure and stand age ( Figure 5). The three models overestimated AGB in the plots where measured biomass was very low or near zero, but they tended to underestimate AGB for plots with very high biomass measurements ( Figure 6). LiDAR metrics performed well in predicting AGB variability. The prediction methods have an impact on AGB estimates. GAM predicted the highest values of AGB, while RF predicted more conservative values (Figure 7), with a 19% difference between methods.

Mapping AGB Uncertainty
The highest uncertainty predicted (expressed by absolute values of AGB) by each prediction method was observed in areas associated with high AGB values (Figure 8). The uncertainty of the estimates, expressed as the range of the 95% confidence intervals divided by their mean (in standardized form) [85], was larger in young stands and open areas (Figure 9). The prediction method affected the AGB uncertainties, with the GAM method produced the lowest uncertainty values.   The width of the plot shows the distribution's density of estimated AGB across the landscape, highly concentrated around the median.

Mapping AGB Uncertainty
The highest uncertainty predicted (expressed by absolute values of AGB) by each prediction method was observed in areas associated with high AGB values (Figure 8). The uncertainty of the estimates, expressed as the range of the 95% confidence intervals divided by their mean (in standardized form) [85], was larger in young stands and open areas (Figure 9). The prediction method affected the AGB uncertainties, with the GAM method produced the lowest uncertainty values.

Mapping AGB Uncertainty
The highest uncertainty predicted (expressed by absolute values of AGB) by each prediction method was observed in areas associated with high AGB values (Figure 8). The uncertainty of the estimates, expressed as the range of the 95% confidence intervals divided by their mean (in standardized form) [85], was larger in young stands and open areas (Figure 9). The prediction method affected the AGB uncertainties, with the GAM method produced the lowest uncertainty values.

Discussion
We identified determinants of AGB and derived spatially explicit estimates by three empirical methods in a managed forest for timber production. AGB stocks were compared with an 80-year-old reference forest. Our results highlight the importance of including historic anthropogenic disturbances (i.e., forest management) for a better estimation of landscape-scale carbon dynamics and ultimately, the Mexican carbon budget as affected by climate mitigation activities. We emphasized the estimation of the AGB uncertainty by building a map that provides a detailed spatially explicit characterization of AGB uncertainty, which is a novel approach for better characterizing its variation across the landscape, identifying the most uncertain areas, and developing strategies to reduce it.

Discussion
We identified determinants of AGB and derived spatially explicit estimates by three empirical methods in a managed forest for timber production. AGB stocks were compared with an 80-year-old reference forest. Our results highlight the importance of including historic anthropogenic disturbances (i.e., forest management) for a better estimation of landscape-scale carbon dynamics and ultimately, the Mexican carbon budget as affected by climate mitigation activities. We emphasized the estimation of the AGB uncertainty by building a map that provides a detailed spatially explicit characterization of AGB uncertainty, which is a novel approach for better characterizing its variation across the landscape, identifying the most uncertain areas, and developing strategies to reduce it.

Discussion
We identified determinants of AGB and derived spatially explicit estimates by three empirical methods in a managed forest for timber production. AGB stocks were compared with an 80-year-old reference forest. Our results highlight the importance of including historic anthropogenic disturbances (i.e., forest management) for a better estimation of landscape-scale carbon dynamics and ultimately, the Mexican carbon budget as affected by climate mitigation activities. We emphasized the estimation of the AGB uncertainty by building a map that provides a detailed spatially explicit characterization of AGB uncertainty, which is a novel approach for better characterizing its variation across the landscape, identifying the most uncertain areas, and developing strategies to reduce it.

Spatial Variability of AGB
Plot-level AGB. Our plot-level AGB estimates across the chronosequence demonstrated how tree biomass increases with stand age. This pattern corroborates previous studies conducted in temperate forests [86,87], confirming that stand age is an important variable for estimating biomass in managed forests for timber production. Estimated AGB in our site was similar to other managed forests in the region [48,88,89], and higher than in some other pine-oak managed forests in México [90]. There was no statistical difference in AGB (>0.05) between the oldest managed stands (20-31-years-old) and the reference forest (Figure 3), revealing that forests managed for timber production in Mexico have a high potential to accumulate biomass and achieve high carbon stocks in a relatively short period of time. However, it is still necessary to integrate other carbon pools (i.e., below-ground biomass, soil organic carbon, etc.) into the analysis of managed forest carbon dynamics to fully quantify the impacts of management.
We noticed that AGB varied among stands of the same age, with smaller variation in young stands (sd = 2.2 Mg ha −1 ) but greater variation in old stands (sd = 42 Mg ha −1 ), and the greatest variation in the reference forest (sd = 64 Mg ha −1 ). Basal area and tree density controlled by management practices [15] such as thinning in our site [88] affected the AGB (p < 0.0001). Age-related differences in AGB among plots are also likely due to heterogeneous site conditions and species composition through variations in wood specific gravity [91], soils, microclimate, and human disturbances over time of forest development. The average AGB in the reference forests in our site was low compared with the range of temperate forests in North America (106-366 Mg ha -1 ) [4], although we did observe AGB > 280 Mg ha -1 in some plots (Figure 3). Most of the mature forests in our site are mixed pine-oak with a relatively low tree density. They appear to be nearing the old-growth stage with large gaps and more dead trees, which typically results in a lower average and a wide range of AGB (Figure 3). This suggests an opportunity to assess the role that biodiversity plays in determining AGB. Additionally, unmanaged natural forests have experienced natural disturbances, such as fire and storms affecting the AGB.

Estimation of AGB across the Landscape
Our results support using LiDAR-based estimates of biomass to improve AGB mapping and support management decisions. LiDAR information is well-suited for inferring forest height, as well as vertical and horizontal forest structure [37] such as canopy cover [92], which successfully estimated AGB across the landscape. Our findings support previous studies that show that AGB can be estimated using LiDAR and ground measurements with a higher accuracy than traditional assessments based on RMSE [93][94][95]. We found that the LiDAR-based AGB varied from values near 0 (recently harvested stand) to the highest values of 302 ± 25.7 and 492 ± 17 Mg ha −1 on the temperate forest landscape due to changes in forest structure and stand age, mainly attributed to forest management practices, that produced a mosaic of different age classes, tree heights and densities, and thus biomasses (Figures 3  and 4). These AGB estimates were higher than the observed data from field sampling alone, suggesting the need to increase the sampling effort to capture a larger variability of forest structure and thus AGB. For forest and land managers, to have this detailed knowledge of the spatial distribution of AGB is helpful for designing silvicultural practices and harvesting schedules to increase carbon sequestration on their lands, in light of REDD+ and carbon mitigation programs.

Effects of Climate and Topography
Topography affects the AGB variability in forest landscapes [9,96]; however, in the present study, this influence was overshadowed by the various structural conditions existing in the stands that make up the forest, a product of the intensive management activities (silvicultural regime) implemented in the last 40 years. Slope and relative slope position had a negative effect on AGB. Valley depth was positively correlated with AGB, representing the vegetation remaining along rivers or streams to protect the water. While topography and climate explained 57% of the AGB in the reference forest, they only explained 11.4% of it in the managed forest, likely due to the conditions within stands derived from the planned silvicultural activities. However, this study provides critical information on disturbances, climatic, and topography which can be used to test and parameterize process-based models in temperate managed forest for assessing spatially explicit carbon dynamics.

Comparing Prediction Methods
The three AGB fitted models performed well when predicting AGB and associated spatial distributions (Figures 5-7). Comparatively, GAM was better for representing ABG variations across the landscape (Figures 5 and 7 All fitted AGB models tended to overestimate AGB in sites with lower field values, and underestimate AGB in sites where measured AGB is high ( Figure 6). This pattern is similar to previously reported studies [21,94,100], perhaps due to increasing forest structure differentiation with forest development, in contrast with the models representing more average conditions. Compared to LM, the RF method resulted in a better model performance, similar to previous studies [37,101]. However, in our case, the GAM method performed better than LM and RF based on RSME and R 2 . This method was previously applied with success in ecological and vegetation modeling [102,103]. Although the GAM method has only been used as an exploratory tool in biomass studies [8], in our study, the GAM method improved AGB predictions more than the other two methods. It showed high potential for use in similar analyses. GAMs may describe a complex relationship between the response and predictor variables. The resulting model suggests that AGB is primarily a function of stand age across the region. Therefore, generating a spatially-explicit age map is a necessary step for AGB estimation across larger regions where LiDAR data is not available.

Spatial Uncertainty
Several recent studies have incorporated more thorough statistical techniques to estimate the uncertainty of various biomass estimators across landscapes [104]; however, few studies have reported uncertainties in the spatial variability of AGB. In this study, we mapped the estimated uncertainty of AGB and identified the areas with high absolute values (Figure 8), which could assist land managers in better planning sampling strategies and management practices. Higher absolute values of uncertainty are associated with higher predicted values, showing that forest complexity may affect biomass predictions. A similar pattern was found by Su, et al. [21] in different forest types in China. The uncertainties of the estimates, expressed as the standardized range of the 95% confidence intervals (Figure 9), were generally small when observations were dense, and the highest in young stands and open areas, probably associated with the lack of mapping understory and small tree size (<5 cm DBH) biomass, which contributed up to 4.6 and 15.9 Mg ha −1 , respectively. The variability in AGB due to remaining seed trees left after harvesting also could contribute to these uncertainties.
Additionally, we assumed the stand age based on the years after harvesting; however, the first successional stage of managed forest could last for several years, and forest structure variables related to age (i.e., basal area) may not represent the real stand age since some large scattered trees are left after harvesting due to the used regeneration method. Although useful, the lack of continuous spatial coverage of the predictors themselves (i.e., stand age) is a problem when one is attempting to use interpolation methods to spatially predict AGB. Therefore, it makes sense to use forest structure variables to better predict stand age in managed forest across larger areas in order to reduce AGB uncertainties.
This estimated uncertainty can result in different correlations between field-measured AGB and predictors. This suggests an opportunity to increase efforts to quantify the uncertainties associated with different sources of error that interact and propagate, such as measurement error, sampling error, the AGB allometric equation [105,106], plot location uncertainty, and the edge effect [21], to improve the estimates of AGB and carbon stocks. It is necessary to quantify the map uncertainty at different spatial scales including different sources of error. Nonetheless, the current study represents an important approach to reporting uncertainty for improving estimates of carbon stocks to support sustainable forest management, as envisioned by REDD + programs.

Improving AGB Estimates in Managed Forests
The implications of our study for estimating biomass carbon stocks and improving forest management are significant. Over the forested study site (∼823 ha), the AGB estimated by the LM method was 78,928.0 (±242) Mg, by the RF method was 86,356 (±220) Mg, and the GAM method was 94,654.0 (±190) Mg. These results are higher than those obtained by Torres-Vivar, et al. [89] using different approaches. In that study, the authors estimated 80,431 Mg of AGB by inventory using a regression estimator, and 77,765 Mg using a ratio estimator. The results of this study suggest that it is important to evaluate the prediction method in its application to management decisions, since the method could lead to the significant underestimation of carbon storage. Particularly, the GAM model appears to be the preferred choice among those tested for estimating AGB, because of its better explicative power and the minimum RMSE, which does a better job simulating nonlinear relationships between predictors and observed biomass than LM. Most biomass prediction studies have the issue of truncating biomass distribution by underpredicting in areas of high biomass and overpredicting in areas of low biomass [107]. As demonstrated by this study, the GAM method is more robust than the other two types of models (LM and RF), and able to improve biomass predictions. In addition, the validation of the highest AGB values predicted by three methods using independent data from five sampling plots further illustrated that GAM estimates were solid, consistent, and more accurate than the other methods for the study site. Compared to other data-based methods, GAM also has the advantage of providing spatial estimates and their uncertainties across the landscape. In sum, these modeling techniques and detailed information (maps, uncertainty) are very important for increasing interoperability among data sources [108] and developing carbon management strategies of climate change mitigation that facilitate the implementation of REDD+ programs and actions.

Conclusions
Detailed knowledge of the spatial distribution of above-ground forest biomass is important for improving estimates of carbon sources and sinks in forests. In this study, we applied three models that were able to explain from 80.9 to 89.8% of the variability in the above-ground tree biomass relating LiDAR metrics of tree height and canopy density, and stand age. The analysis showed that spatial heterogeneity in AGB increased with increasing stand structural complexity and age in the managed forests, varying from 0 to 302-492 Mg ha −1 . The older managed stands had a similar average AGB (166.6 ± 42.4 Mg ha −1 ) to the reference stand (159.5 ± 67.4 Mg ha −1 ). This reveals that managed forests in Mexico have a high potential to efficiently accumulate biomass and carbon in a relatively short period of time, which has important implications for supporting carbon mitigation programs such as REDD+, and to contribute to the sustainable management of forests. Although climate and topography variables contribute to AGB variability, they were overshadowed by the effects of management practices for predicting AGB. The statistical protocol developed in this work is able to show the spatial pattern of uncertainty associated with each AGB prediction, which could be implemented in other studies of AGB. Uncertainty was affected by the prediction method, was not spatially homogeneously distributed, and was concentrated in areas with higher AGB values, indicating that those areas may need more intensive sampling in order to effectively manage their biomass stocks, as well as to identify the most uncertain areas and develop working hypotheses to reduce uncertainty. The empirical models generated in this study (i.e., GAM) could be useful to fill in spatial and temporal gaps in observations, since GAM estimates were consistent and more accurate than the other two methods tested in our study. This approach allowed us to identify areas with high AGB values that were not sampled.
Our study results represent an advancement toward the development of efficient methods to identify and enhance carbon stocks in managed forests, as a component of mitigation activities.