Introduction

Natural resources have been used to meet the needs of societies throughout the history of mankind, which has led to their overexploitation and eventual deterioration of extensive territories, which implies, in many cases, the unavoidable loss of biodiversity genetics. Reforestation, as a remediation strategy, is a continuous ecological restoration activity in Mexico.

With the creation of the Comisión Nacional Forestal (National Forestry Commission) (Conafor) in 2003, programs have been designed with the aim to exploit, conserve and restore forest resources through the granting of financial support to producers. The Forest Management Units (UMAFORES) were established on the basis of Article 12 of the General Law of Sustainable Forest Development, which considers the need to delimit the entire country under this scheme of ordination, with the purpose of exercising proper planning of activities and efficient management of forest resources.

A reforestation program aims to promote the revegetation of deforested areas by increasing the amount of biomass. In order to assess its relevance, several strategies are proposed, one of which involves considering the success of reforestation directly through the forest inventory of the supported areas, with the limitation of costs in time and resources, as well as logistics, as there are access difficulties. Another way of estimating this aspect is to evaluate indirectly, through the analysis of materials from remote sensing, before and after reforestation.

Calculating the amount of biomass through satellite imagery has favored the analysis of forest cover at different scales thanks to the wide range of spatial resolutions offered by sensors installed in satellites (^{Dube and Mutanga, 2015}; ^{Scottá and da Fonseca, 2015}; ^{Valdez et al., 2006}; ^{Yin et al., 2015}). Landsat satellite imagery has been used repeatedly for the estimation of forest variables (^{Barasa et al., 2010}, ^{Günlü et al., 2014}, ^{Xu and Guo, 2014}).

The methods used for the estimation of biomass by satellite imagery include supervised classification (^{Luong et al., 2015}, ^{Perumal and Bhaskaran, 2010}), unsupervised classification (^{Bharatkar and Patel, 2013}, ^{Wang et al., 2015}), as well as the generation of the Normalized Difference Vegetation Index (NDVI) (^{Ji et al., 2012}, ^{Sall et al., 2013}), which generates a numerical value indicating the vigor of vegetation by means of the relation between the amount of energy reflected in the near infrared and the one absorbed by vegetation in the of the red part (^{Tun et al., 2008}), which indicates the existing biomass in an indirect way, With the field information it is necessary to validate it by calculating the amount of biomass of the plants present in the field, using allometric equations (^{Aguilar et al., 2012}; ^{Balderas and Lovett, 2013}; ^{Vahedi et al., 2014}).

The results of the calculation of the biomass in the field and the generated with the algorithm of the NDVI have served as the basis for the development of mathematical models that allow to relate this information (^{Aguirre et al., 2012}; ^{Karlson et al., 2015}; ^{Lu et al., 2012}); once statistically validated, are used to estimate the aboveground biomass by extrapolating the values obtained from the image process to an area or region.

The verification of the results is done by comparing the data obtained by remote perception with respect to those calculated in the field (^{Cabral and Zamyatin, 2009}; ^{Seebach et al., 2011}). To analyze the dynamics of forest cover in a given period, the difference in the amount of biomass between two or more dates is estimated.

Assuming that reforestation actions should have a positive regional effect on biomass growth, the objective was to assess the impact of the reforestation program implemented in 2010 in the region covered by UMAFOR 2404.

Materials and Methods

The study area is located in the state of San Luis Potosí, located in the north central part of Mexico and comprises the administrative region of UMAFOR 2404. In this analysis, the properties that were supported in the reforestation program by the Mexican government in 2010, which are located between the coordinates 255 550 E, 2 509 507 N and 379 230 E, 2 382 393 N (Figure 1) and they have an area of 553 284 ha were given particular interest.

Figure 1 Geographical location of the Unidad de Manejo Forestal (Forest Management Unit) 2404. 

The climate varies between very dry (Bw) and subhumid temperate (Cw) (^{García, 1981}). There are nine types of vegetation in the area, three of which correspond to temperate forests with pure pine and oak masses, as well as mixed forests of both genera, shrubs formed by shrub oaks, shrub communities of Prosopis spp. and Acacia spp., submontane scrubland, desert scrubland, halophyte vegetation and grasslands (^{Rzedowski, 1965}).

According to the Sistema Nacional de Información Forestal (National Forestry Information System), in 2010, support was given for the reforestation of an area of 2 893 ha for UMAFOR 2404, of which 97 % was used for reforestation with plants obtained by vegetative propagation and 3 % to plants coming from seed.

Nine Landsat images corresponding to the orbit / line 27/45, 28/44 and 28/45, of the years 2000, 2010 and 2015 were used; the first formed the baseline of comparison (Table 1). Landsat 5 TM images were downloaded from http://glcfapp.glcf.umd.edu:8080/esdi/ and those of 2015, from http://earthexplorer.usgs.gov/. Images were processed by the Imagine 201 Erdas TM program (^{ERDAS, 2010}) and were geometrically fixed to make their register and thus allow their comparison.

Table 1 Characteristics of the Landsat satellite images used for temporal analysis of biomass in forest plantations. 

To reduce the volume of data, subpictures were created that cover only the selected areas. A digital cover in vector format was used that contained the limits of the farms with works of reforestation in 2010 provided by Conafor. The resulting subimages were applied the NDVI generation algorithm with the following relation:

For the creation and validation of the model, the information recorded at the Inventario Nacional Forestal y de Suelos (INFyS) (National Forest and Soil Inventory (INFyS)) sampling sites from 2005 to 2007 was taken; data from the adult woodland record collected in 55 conglomerates were used, 44 for the development of the models and 1 for their validation (Figure 3). In order to calculate the aboveground biomass (trunk, branches and foliage) allometric equations were developed at the genus level or, if they exist at the species level, which have as input the values of the diameter and height variables, which were generated by several authors (^{Návar, 2009}; ^{Rodríguez et al., 2006}; ^{Rodríguez et al., 2007}; ^{Rojas et al., 2015}) for similar regions. The total biomass value per unit area was obtained with the sum of the total aboveground biomass of the trees recorded by conglomerate and extrapolated to values per hectare.

In order to establish these prediction models, the estimated biomass information for the 44 field sampled clusters was estimated with the NDVI values for each date, and a sample of 25 pixels was taken from the images corresponding to the values of the conglomerate.

An analysis was performed to determine the degree of correlation between the field data and those obtained from the satellite images, for which the Pearson correlation index “r” was calculated.

With the results obtained in the previous processes, linear, exponential and polynomial prediction models were developed.

In order to choose the model, the values of the root mean square error (RMSE) statistics and the R² coefficient of determination were considered.

To compare the variance in the three-year estimates, a variance analysis (ANOVA) was performed with a 95 % confidence interval (α = 0.05).

Once the prediction models were selected, they were applied to each of the images from which came a digital chart with the amount of biomass per unit area for the study area. The total biomass value of the selected area was obtained with the sum of the total aboveground biomass per unit area.

Results and Discussion

The distribution of the vegetation index defined by the NDVI values is shown in Figure 2. The gray tones represent eight ranges of vegetation cover, the most obscure ones refer to the low numbers and the clear ones to the high ones.

As in the actual study, ^{Alatorre et al. (2014)} found a tendency of increase in the NDVI numbers from the same sensor in Chihuahua State in a multitemporal analysis (1986-2011).

a) 2000, b) 2010 and c) 2015.

Figure 2 Distribution of the NDVI values 

The spatial distribution of the sampling clusters used in both, modeling (dot) as in validation (triangle) are shown in Figure 3. In them, 28 genera and 66 tree species were recorded: Quercus with 22 species was the best represented, followed by Pinus with nine species; Acacia, Arbutus, Crataegus, Juglans, Juniperus, Prunusa and Yucca with two species; Acer, Arctostaphylos, Buddleia, Carya, Casimiroa, Eysenhardtia, Fraxinus, Gochnatia, Jacaranda, Lysiloma, Nolina, Persea, Pithecellobium, Platanus, Prosopis, Pyrus, Rhus, Stenocereus and Zanthoxylum with one species.

Figure 3 Spatial distribution of the clusters of the Inventario Nacional Forestal y de Suelos (Conafor, 2007) in the study área 

In order to determine the correspondence between the estimated field data and those obtained from the satellite image processing, the Pearson correlation index was calculated, resulting in values of 0.69, 0.65 and 0.64 for 2000, 2010 and 2015, respectively, which indicate a high correlation. In working similar problems, ^{Gunawardena et al. (2015)} in Sri Lanka and ^{Zhu and Liu (2015)} in Peru made similar correlations by imaging the sensors that were used in the present study.

From the aerial biomass prediction models developed with the NDVI values and the field estimates, the polynomial type showed the best determination coefficients (R²) and Mean Squared Error values (Table 2); in them the adjusted models are shown, as well as the values of R² and RMSE for each one.

B = Biomass; NDVI = NDVI value.

Table 2 Fitted models 

The polynomial model was selected for three years. In all cases, the selected model yielded the highest value of R² and the lowest RMSE. As far as the R² values are concerned, those calculated in this work are inferior to those presented by ^{Ruiz et al. (2014)}, who recorded values between 0.9765 and 0.9914 in the aboveground biomass calculation for Quercus laurina Bonpl. and Quercus crassifolia Bonpl. in a pine-oak forest in Oaxaca State.

Also, the RMSE of the present study is smaller than that described by ^{Martínez et al. (2016)}, whose values are distributed in a range of 0.68 to 0.99 for biomass by allometric equations for tree components (stem, branches and leaves) in the state of Durango. Those in the present study are similar to the lower value reported by ^{Powell et al. (2013)}, which was 22 Mg ha^-1, and the maximum 221 Mg ha^-1, which surpasses it. In contrast, the own results exceed most of those determined by ^{Main et al. (2011)}, who calculated 10.22 and 12.64 Mg ha^-1 in two models, while a third model showed a value of 40.64 Mg ha^-1. It is important to consider that non-linear models can overestimate calculated values (^{Aguirre et al., 2012}).

The selected models were applied to the image containing the NDVI values for the corresponding year to estimate the amount of biomass. Only the positive NDVI numbers representing surface biomass were taken. Figure 4 shows images of biomass estimation maps in the three years. The values presented in this study are similar to those of ^{Aguirre et al. (2012)}, who reached productions between 5 and 150 Mg ha^-1 per unit area in this study area depending on their vegetation cover.

In the images it can be observed that the increase of the biomass has occurred in different zones. For example, in 2000 the highest amounts are concentrated to the center and northwest; in 2010, to the center, east and southeast; and in 2015, in the center-south and northwest.

a) 2000, b) 2010 and c) 2015

Figure 4 Biomass spatial distribution images 

With the resulting values from the calculation of biomass and through the polynomial model the variance for each one of the years was analyzed (Table 3).

Table 3 Analysis of variance of biomass values and polynomial model. 

According to the information in Table 4, there was no significant difference between the values of the variance for the three years, which reinforces that the value of F was lower than the critical value in them. These results surpass those of ^{Trigueros et al. (2014)}, who registered a value of 2.712 in the calculation of normal diameter, height and basal area; however, there was agreement with this analysis, since the values of F were below the critical value in the three periods.

Total aboveground biomass for each of the analyzed years is shown in Table 4.

Table 4 Biomass estimated values. 

From 2000 to 2015 there was a considerable increase in the amount of biomass (423 %); From 2000 to 2010 was 234 %. In addition, the average biomass per ha had a similar behavior when going from 0.26 to 1.10 Mg ha^-1. These values are lower than those reported by ^{Velásquez et al. (2012)}, who calculated an average annual increase of 2.9 Mg ha^-1 in a region of the Colombian Andes. Also, the biomass results per hectare of the present study were lower than those of ^{Restrepo et al. (2016)}, who record an average of 4 Mg ha^-1 in a forest of Colombia.

Finally, the Pearson correlation coefficient was calculated in the clusters used for validation, which showed 0.69, 0.66 and 0.73 for 2000, 2010 and 2015, respectively. These values indicate a high correlation.

Conclusions

From the total estimated values, there was a 423 % increase in the amount of biomass, from 2000 to 2015, which reveals that the reforestation program has had the expected effects from its original design.

The models that more accurately estimated the biomass were chosen based on the values of R² and RMSE, but the values obtained with the nonlinear models could be overestimated given their own characteristics