Carbon mapping in pine-oak stands under timber management in southern Mexico

Study area

The research was carried out in a forest under the management of San Juan Lachao, Pueblo Nuevo, Juquila, Oaxaca, Mexico, with coordinates of 16°09′30.26″N and 97°07′28.04″W and an average altitude of 1,900 m (Fig. 1). The predominant climate is warm subhumid (Cw), with an average annual temperature of 22 °C and an average rainfall of 2500 mm. The community of San Juan Lachao has a pine-oak forest of 3,201 ha, of which 573 ha are under forest management through the silvicultural development method (SDM) and 1,787 ha are under the Mexican Management Method for Uneven-aged Forests (MMOBI). This area contains species of timber importance, such as Pinus douglasiana Mtz., P. maximinoi HE Moore, P. devoniana Lindl., Quercus rugosa Née, Q. crassifolia Humb. & amp. Bonpl., and other broadleaved trees (Servicios Técnicos Forestales STF, 2011).

Tree inventory and laboratory analysis

In 2022, 12 circular sites of 400 m² with slope compensations were established in a targeted manner in stands harvested during 2013–2014 under the methodology used by Chávez-Pascual et al. (2017) and Miguel-Martínez et al. (2016). The sites were established with four silvicultural treatments: free thinning (FT), first thinning (T1), regeneration cut (RC) under the parent tree method and unmanaged area (UA), with three repetitions of each ST. The sites were geopositioned with a global positioning system (GPS) (Garmin eTrex 30, USA^®) MAPS 6 within each stand. Descriptive variables of the site were taken: altitude (m); average slope of the site (%), with a clinometer (Haglöf^®, EC II D-HS115, Spain); and exposure.

At the sites, a tree inventory was carried out, recording the variables: diameter at breast heigt (DAH, cm) with a diametric tape (Hartmann^®, 283D); total height (TH, m) with a clinometer (Haglöf^®, EC II D-HS115); and crown area (CA, m) with a tape measure. Trees were selected in proportion to the relative frequency of each species to obtain a sample (chip or slice); for Pinus spp., a Pressler drill (Haglöf^®, Sweden) was used to obtain a chip at a height of 1.30 m; in some species such as Quercus spp., Arbutus xalapensis Kunth, and Alnus acuminata Kunth, at least one tree of each species was cut per site to obtain a five cm thick slice (Rodríguez-Ortiz et al., 2019).

The samples were analyzed in the Agroecosystems Laboratory of the Technological Institute of the Valley of Oaxaca. Core increment and slices of all species were weighed on an analytical scale (Shimadzu^®, ATY224, ±1 mg) to determine their green weight (GW g). The green volume (GV, cm³) was determined with a digital Vernier caliper and by applying Newton’s formula (Romahn de la Vega & Ramírez, 2010). Subsequently, the samples were placed into a drying oven (Memmert^®, 100-800) at 102 °C until constant weight and dry weight (DW, g) were obtained. The specific gravity (SG, kg m⁻³) by species was obtained: SG = (DW/GV) × 100 (Bhardwaj et al., 2016). Subsequently, the samples were ground to determine their C content through the content of organic matter under the Walkey and Black method (SEMARNAT, 2000); 58% of the determined organic matter was C (Fig. 2).

The total tree volume with bark (TTVbark, m³) of all tree species found was determined with the forest biometric system (Vargas-Larreta et al., 2017); the product volume and SG (kg^m−3) provided the TB content (kg tree⁻¹) (Rodríguez-Ortíz et al., 2011).

Sentinel 2 spectral data

The spectral image of the study area was downloaded, divided into 13 bands, on the Sentinels Scientific Data Hub platform (https://scihub.copernicus.eu/). The discharge was at level 2A, which involves preprocessing at the atmospheric correction level, (Table 1).

The image was processed with QGIS 2.16.16^® software, where atmospheric correction was first applied to each band. To visualize the image in true color, a set of bands 2, 3 and 4 was constructed, and all three had a 10 m resolution per pixel. Subsequently, a set of bands 3, 4 and 8 (10 m pixel⁻¹) was constructed to visualize the vegetation in infrared bands (Fig. 3).

Bands 4 (red) and 8 (near infrared (NIR)) were used to calculate the NDVI per pixel. The following formula was applied: NDVI = (NIR –red)/(NIR + red), where red and NIR are the spectral reflectance measurements acquired in the red (visible) and near infrared regions, respectively. These spectral reflectances are proportions of the radiation reflected on each individual spectral band. The NDVI, real visualization and infrared vegetation layers were cut from the area under use and from the 4 stands where the sampling sites were located (Fig. 4).

Other vegetation indices were also tested in this study:

$GNDVI=\left(\mathrm{NIR}-\text{Green}\right)/\left(\mathrm{NIR}+\text{Green}\right)$, where: GNDVI = Green Normalized Difference Vegetation Index, NIR =near infrared band, Green =green band.

$EVI=2.5\times \left(\left(NIR-\mathrm{Re}d\right)/\left(\left(NIR\right)+\left(C1\times \mathrm{Re}d\right)-\left(C2\times Blue\right)+L\right)\right)$, where: EVI = Enhanced Vegetation Index, NIR =near infrared band, Red = red band, Blue = blue band, C1 and C2 = Atmospheric correction coefficients (6 and 7.5), L = soil influence correction factor.

$SAVI=\left(\left(NIR-\mathrm{Re}d\right)/\left(NIR+\mathrm{Re}d+L\right)\right)\times \left(1+L\right)$, where: SAVI = Soil Adjusted Vegetation Index, NIR = near infrared band, Red = red band L = soil luminosity correction factor (0.428).

$AVI={\left[\mathrm{NIR}\times \left(1-\mathrm{Red}\right)\times \left(\mathrm{NIR}-\mathrm{Red}\right)\right]}^{1/3}$ , where: AVI =Advanced Vegetation Index, NIR =near infrared band, Red =red band.

Data analysis

For the TB and C variables, the assumptions of normality and homogeneity of variances were verified (Shapiro–Wilk and Bartlett test, respectively, α = 0.05). Once the NDVI values per pixel (10 m resolution) were obtained, the average value of the NDVI per sampling site was determined using a 3  × 3 pixel window surrounding the central point. For this part of the process, a validation was carried out through “raster statistics for points”. Allometric models were adjusted to estimate C and TB (t pixel⁻¹ and t ha⁻¹) based on the NDVI (Table 2). Subsequently, with the models, the increase in TB and carbon was determined by site, stand and management area through the NDVI value of each pixel.

The STs were differentiated using a generalized linear model (PROC GLM) and the Tukey means test (α = 0.05). All analyses were performed in the Statistical Analysis System (SAS) program (SAS Institute, Inc., 2017).

Model adjustment to estimate biomass and carbon

The linear models without an intercept (β₀) were those that best adjusted for the estimation of TB and C as a function of the NDVI, being highly significant (p ≤ 0.01). The NDVI explained 92.2% and 91.9% of the existing variation in TB and C, respectively. The coefficients of variation (CVs) ≤ 32.2% were moderate-high due to the variation in tree density and coverage between the applied STs. The standard deviation of the models was low, which indicates that the estimated TB and C values were slightly removed from the trend of the mean (Table 2).

The model was validated through a Student’s t test (α = 0.05) comparing TB and C per stand (t), estimated with field data vs. estimated with spectral data. No differences were found between both types of data (p > 0.05) (Figs. 5A and 5B). Therefore, the models were efficient at predicting TB and C through the NDVI since they generated reliable results when estimating at a large scale (estimation).

Biomass and carbon estimations

For both variables, TB and CC, the stands under FT (21.4 ha) and T1 (26.4 ha) presented similar average estimated values of 140 t TB ha⁻¹ and 63 t C ha⁻¹ and 137 t TB ha⁻¹ and 61 t C ha⁻¹, respectively, such that the total amounts per the FT and T1 stands were 1,343 t C and 1,620 t C, respectively. The minimum and maximum values between the areas under the two types of STs were similar because their function was to decrease tree density when they are in the pole stage (Fig. 6).

The stand with RC (14.4 ha) presented values of 139 t TB ha⁻¹ and 62 t C ha ⁻¹, with minimum and maximum values for C between 41 and 81 t ha⁻¹, showing a good accumulation of C in response to the regeneration of the parent tree treatment. The UA (17 ha) presented average values of 129 t TB ha⁻¹ and 58 t C ha⁻¹, with a range of variation between the minimum and maximum values of 37 to 79 t C ha⁻¹. The stands subjected to the STs (FT, T1, and RC) presented the same maximum value of carbon accumulation (81 t ha⁻¹), which was 2% higher than that in the UA (Fig. 6).

The FT stand presented 8.6% more carbon (t ha⁻¹) than the UA stand, the FT and UA stands had tree densities of 1166 trees ha⁻¹ and 825 trees ha⁻¹, respectively. This small difference in carbon values occurred because the FT stands maintained a contemporaneous structure with similar ND, crown diameter (CD) and TH, while the UA stands were highly heterogeneous, presenting variation in their area.

In the area managed by San Juan Lachao (3201 ha), a total of 404,048.08 t TB and 181,049.56 t C were estimated. An average accumulation of 126 t TB h⁻¹ and 57 t C ha⁻¹ was estimated, with minimum and maximum values of 27 to 209 t TB ha⁻¹ and 12 to 94 t C ha⁻¹, respectively. The areas with high accumulations of TB and C were observed in an intense red color, unlike the green areas, which showed the least accumulation (Fig. 7).

All STs resulted in similar amounts (p > 0.05) of TB and C by genus (Pinus and Quercus), as well as total TB amounts (inferred in field sampling). On the other hand, the tree C estimated with the NDVI in the spectral image was higher (p ≤ 0.01) in the RC (66.6 t ha⁻¹) than in the UA (57.9 t ha⁻¹), which showed greater heterogeneity; the same behavior was found for TB (Table 3).

This difference was visible because the UA had a high heterogeneity, so making an inference with field data was not feasible, unlike the estimation with the satellite image, which was more accurate in terms of irregular stand structure (Table 3). In addition, the areas with STs had a more homogeneous trend with a systematically distributed canopy, and the opposite occurred for the ASM, where heterogeneity led to an irregular structure.

Model adjustment

The stands with STs (FT, T1, and RC) showed similar TB and C values but varied from those of the UA due to the age of the stand, the crown area, the size of the trees, edaphic factors and altitude (Rajput, Bhardwaj & Pala, 2017; Ruiz-Díaz et al., 2014). The areas with T1 and FT thinnings presented an average of 138.5 t TB ha⁻¹ and 62 t ha⁻¹ of carbon. The area with RC presented 139 t TB ha⁻¹ and 62 t C ha⁻¹ (Fig. 6). Regarding the above information, Aguirre-Salado et al. (2009) reported in contemporary stands of P. patula Schl. et Cham. under harvested conditions, an estimate determined by remote sensing was 55 t C ha⁻¹.

The best fitted models in this work were ”simple linear” as a function of the NDVI (10 m pixel⁻¹), showing significance (p ≤ 0.01) and adjusted R² = 0.92, with which the total managed area (3201 ha) having 126 t TB ha⁻¹ (values from 27 to 209 t TB ha⁻¹) and 57 t C ha⁻¹ (values from 12 to 94 t C ha⁻¹), which shows potential values of TB and C (Fig. 7). Thurner et al. (2014) adjusted simple linear models to estimate carbon (R² = 0.70–0.90) in mixed forests of North America, Europe and Asia, where they found an average of 58 ± 22.1 t C ha⁻¹. Reyes-Cárdenas et al. (2019) estimated the forest TB in northern Mexico, adjusting an exponential model based on the NDVI, and found lower values (0.85 to 157 t ha⁻¹) than those in this study.

Choosing the best model depends on the type of vegetation that is being evaluated and the type of optical satellite that is being used; thus, the evaluator is responsible for choosing the best model; for example, Aguirre-Salado et al. (2009) adjusted a multiple regression model using the NDVI and the water stress index as independent variables (R² = 0.70). D’Oliveira et al. (2020) estimated the biomass of 10 forest plots in the southwestern Brazilian Amazon, adjusting multiple linear models based on LiDAR variables, (R² = 0.90, RMSE = 13.23 t ha⁻¹), presenting values from 11.1 to 273 t ha⁻¹. On the other hand, Rex et al. (2020) estimated biomass with LiDAR images (229.10 t ha⁻¹) in tropical forests with low intensity logging in Pará, Brazil, although in this case using ordinary least squares (OLS) regression and a selection of variables through principal component analysis (R² = 0.35 - R² = 0.53), and the model had less significance than that of other.

The model was validated through Student’s t test (α = 0.05) comparing TB and C per stand (t), and the values were estimated with field data vs. estimated with spectral data, showing no significant differences (p > 0.05); correlating field values with spectral data provides effective and accurate information (Vaghela et al., 2021). Reyes-Cárdenas et al. (2017) correlated field data with spectral data and obtained similar values. Aguirre-Salado et al. (2012) adjusted a linear model to estimate aerial biomass based on the NDVI and continuous vegetation fields (CVFs) (adjusted R² = 0.77 CME = 26.00 t ha⁻¹), showing high validation correlation coefficients (r = 0.87**); Verly et al. (2023) estimated the C of the Atlantic Forest with Sentinel 2 images, showing a high correlation coefficient between the estimated and observed carbon averages (r = 0.84).

Biomass and carbon estimation

The UA presented TB and C values of 129 and 58 t ha⁻¹ (Fig. 6), respectively, which indicated that it is necessary to use satellite images for estimations since they provide values that can be extrapolated to the area of study. Perea-Ardila, Andrade-Castañeda & Segura-Madrigal (2021) found similar TB and tree C values as a function of the NDVI in a high Andean forest of Boyacá, Colombia, determining 168.0 ± 11.2 t TB ha⁻¹ and 84.0 ± 5.61 t C ha⁻¹. Similarly, Clerici et al. (2016) reported 180.7 ± 23.8 t ha⁻¹ of biomass in 400 m² plots in Andean forests in Cundinamarca, and Yepes-Quintero et al. (2011) reported TB and CC of 102.38 ± 25.22 t ha⁻¹ and 51.19 ± 12.61 t ha⁻¹, respectively, in highland forests of Antioquia. On the other hand, Bhardwaj et al. (2016) reported an estimated biomass with the NDVI of 169.05 to 265.83 t ha⁻¹ in subtropical forests of the northwestern Himalayas. Therefore, satellite images are a potential tool for differentiating managed and unmanaged areas (Avogadro & Padró, 2019).

After testing various vegetation indices in this work, the NDVI was best adjusted (Table 2) to estimate TB and CC (0.09–0.69), providing positive values and indicating areas with healthy and vigorous vegetation (Fig. 3); however, although in some areas the NDVI values were low, the values of TB and CC were positive. This result was in contrast to those of Aguirre-Salado et al. (2012), where the estimates tended to be negative when the pixels had low NDVI values. The NDVI is an excellent indicator because it had three visible red bands that increased accuracy and had the near infrared band (Wang et al., 2018). Similarly, their the estimates of C for areas under timber management were highly accurate given the methodology that was selected was correctly followed (Yan et al., 2016).

The NDVI to estimate TB and C in this work was crucial, since the NIR spectrum contains information about the chemical composition (infrared light energy) through the photosynthetic pigment “chlorophyll”; chemical bond information such as C-H, N-H, S-H, C = O and O-H (Guerrero-Maestre et al., 2008). For this reason, the NDVI was adapted to the type of vegetation evaluated, where it has a high reflectance in the near infrared and low in the red band (Jensen, 2015).

Although there are other methods to evaluate the vegetation, such as the use of LiDAR (Light Detection and Ranging) images, in this study the most feasible was the use of satellite images, due to the scale to be evaluated; and it is suitably viable since it provides us with information on areas with and without timber forest management. Using Sentinel 2 satellite images in vegetated areas significantly increases estimate accuracy (Polat, Akcay & Balik Sanli, 2022).

Through remote sensing, it is possible to extract information about the vegetation in pine-oak forests that are unmanaged and classify the vegetation more easily that through observation methods (Ancira-Sánchez & Treviño Garza, 2015; Polat, Akcay & Balik Sanli, 2022). However, factors such as topography, altitude, slope, precipitation, and temperature must always be considered, as indicated in Olthoff, Martinez-Ruiz & Alday (2016). By determining the TB and CC of a large UA, this area can be classified into smaller areas, and suitable silvicultural management methods and treatments can be determined through a management plan (Agus et al., 2004; Puletti, Chianucci & Castaldi, 2018), in addition to promoting C cycle mechanisms (Herold et al., 2011).