Forest Fires, Land Use Changes and Their Impact on Hydrological Balance in Temperate Forests of Central Mexico

Abstract

Temperate forests play a fundamental role in the provision, regulation, and support of hydrological environmental services, but they are subject to constant changes in land use (clearing, overgrazing, deforestation, and forest fires) that upset the hydrological balance. Through scenarios simulated with the Water Evaluation and Planning (WEAP) hydrological model, the present study analyzes the effects of forest fires and land use changes on the hydrological balance in the microwatersheds of central Mexico. The land use changes that took place between 1995 and 2021 were estimated, and projections based on the current scenario were made. Two trend scenarios were proposed for 2047: one with a positive trend (forest permanence) and the other with a negative trend (loss of cover from forest fires). The results show that with permanence or an increase in forest area, the surface runoff would decrease by 48.2%, increasing the base flow by 37% and the soil moisture by 2.3%. If forest is lost, surface runoff would increase up to 454%, and soil moisture would decrease by 27%. If the current forest decline trends continue, then there will be negative alterations in hydrological processes: a reduction in the interception of precipitation by the canopy and an increase in the velocity and flow of surface runoff, among others. The final result will be a lower amount of water being infiltrated into the soil and stored in the subsoil. The provision of hydrological environmental services depends on the maintenance of forest cover.

1. Introduction

Temperate forests provide and regulate environmental goods and services, which bring direct and indirect benefits to society to meet its needs and are essential for human well-being [1,2,3]. Ecosystem benefits are numerous and varied, and they are classified as support, provision, regulation, and cultural benefits [2,4,5]. Two of the most important ecosystem services that are provided by forests are provision and hydrological regulation. The capacity of a basin to provide hydrological services depends on topographic characteristics, vegetation cover, land use, and climatology [6,7]. However, in basins with dominant forest cover, it is the forest that ensures the integral, continuous, and stable flow of the natural elements that intervene in the hydrological cycle [8]. This is because the forest is the layer of the Earth’s surface that is responsible for capturing and buffering rain, controlling surface runoff processes, promoting water infiltration, and therefore influencing the recharge of aquifers to maintain stable levels, among other functions [9,10,11,12].

The current capacity of Mexican ecosystems to provide these services is deteriorating [4] since they have suffered constant alterations that have mainly taken place due to clearing, overgrazing, excessive logging, and forest fires [13]. The deterioration and loss of forest cover alters the hydrological cycle [14,15], generating a decrease in the recharge of aquifers, drying bodies of water, and generating sudden and uncontrollable runoff [16]. Changes in forestland use compromise the availability and quality of hydrological environmental services [17]. According to Roff et al. [5], many of the problems mentioned above result from the lack of investment in the protection and management of forests and other natural resources, leading to the depletion of natural vegetative cover and soils, the deterioration of watersheds, and the extinction of species. These effects often lead to considerable economic and social losses.

Therefore, it is necessary to determine the current and future conditions of the watershed behavior in Mexico because these watersheds are a source of hydrological resources at the local level [18]. One way to quantify the effect that an increase or decrease in temperate forest cover would have on water distribution is through modeling. Models are tools that allow us to design and project simulative processes based on the analysis of the behavior of the historical data of the variables involved [19,20,21]. In Mexico, studies that have applied predictive models have mostly been related to agricultural irrigation basins under climate change conditions and have presented results with acceptable effectiveness [18,22,23], but there is little evidence that highlights the importance of linking the natural processes that result from changes in temperate forest cover with the generation of hydrological environmental services.

The WEAP platform is a modeling tool that can be used for the planning and comprehensive distribution of water resources that was developed by the Stockholm Environment Institute (SEI), and it can be applied at different scales, from small catchment areas to large basins [24]. WEAP can provide an integrated assessment of the climate, hydrology, land use, irrigation facilities, water allocation, and water management priorities of the basin and uses a standard linear programming model to solve water allocation problems at any step, and its objective function is to maximize the percentage of the needs of the supply and demand centers with respect to the priority of supply and demand, the hydrological balance, and other limitations [25,26].

The WEAP model operates under the hydrological balance model to help with the management of water resources and the estimation of the environmental services of the water resource and its relationship with land use and climate change [5]. The model is based on multicriteria scenarios so that the researcher can carry out several types of study and can make specific and different types of hydrological: population growth or change, land use change, and climate change [18,24,27]. The scenarios allow us to make decisions, develop adaptation measures, and design regional policies based on the management and conservation of ecosystems [23,28]. WEAP can be adjusted according to the density of the available data; even in regions with scarce data, the model can support a complete hydrological representation [1]. Policymakers often lack the funding and/or experience to develop methods with which to evaluate the complex trade-offs that involve changes in land use, the management of forest areas, and the influence of both on environmental services [29].

This study analyzes the effects of forest fires and changes in temperate forest cover on the hydrological balance in the microwatersheds of central Mexico. We designed and projected scenarios using the WEAP hydrological model, which was assembled using geographic information systems (GIS), to determine the importance that these changes represent today and in the long term. The results can form a basis for new research on hydrological environmental services.

2. Materials and Methods

2.1. Study Area

The research was carried out in the three microbasins of the Texcoco, Chapingo, and San Bernardino Rivers in the municipality of Texcoco, Mexico. They are located between the 19°24′ and 19°28′ N parallels and the 98°43′ and 98°52′ W meridians. They are located at an altitude that ranges from 2243 to 4080 m asl and cover a total area of 77.40 km². These microwatersheds are part of the Texcoco aquifer, key 1507, which belongs to the Pánuco Hydrological Region, within the Hydrological-Administrative Region XIII, Aguas del Valle de México (Figure 1); the aquifer has been reported as being overexploited, and despite the fact that it is completely prohibited by national decrees, said regulatory instruments have not been sufficient to reverse problems related to overexploitation since there is a great demand for groundwater, mainly for urban public use in the region [30]. The temperature ranges between 5 and 16 °C, and the average temperature is 11 °C. Precipitation ranges between 500 and 800 mm per year, with an annual average of 652.5 mm. The predominant climate is temperate dry with cool and rainy summers, with thermal oscillation, interannual Ganges-type temperature variation, and the presence of heat waves (C(w0)bi’gw’’). There is also a semiarid climate with rainy summers (BS1w(w)ki’gw’’) that comprises 21.3% of the total area [31].

The riverbeds in the study area experience ephemeral and intermittent runoff of a torrential nature, with short-term floods and dry streams during the dry season. The three riverbeds, in addition to rainwater, receive and conduct sewage, a situation by which, in many cases, they form parts of the drainage systems [30]. Sedimentary breccia rocks have the greatest distribution, and extrusive igneous rocks of the andesite type are abundant [32]. Hydrogeology indicates that the study area has a medium to high permeability [33]. The soil type with the greatest distribution is epipetric Phaeozem, followed by the epipetroduric Phaeozem [34]. The main types of land use and vegetation correspond to rainfed agriculture (29.6%), temperate forest (26.1%), and reforestation areas (14.8%). In the microbasins of the Texcoco River and Chapingo River, the land uses with the greatest surface areas are temperate forest and vegetation, at 37% and 24.2%, respectively. In the San Bernardino River microbasin, rainfed agriculture is the main land use, at 44.2% (

Table 1. Total and per-microwatershed land use and vegetation cover areas.

Land Use Class	Surface Total	%	Texcoco River	%	Chapingo River	%	San Bernardino River	%
Temperate forest	2021.4	26.1	1448.5	37.1	472.1	24.2	100.8	5.4
Reforestation	1145.6	14.8	199.9	5.1	395.4	20.3	550.3	29.2
Secondary vegetation	438.7	5.7	166.3	4.3	160.9	8.3	111.5	5.9
Grassland	52.2	0.7	48.6	1.2	3.6	0.2	0.1	0.0
Mine	212.2	2.7	39.4	1.0	87.0	4.5	85.9	4.6
Rainfed agriculture	2294.8	29.6	1042.8	26.7	420.4	21.6	831.6	44.2
Irrigated agriculture	540.3	7.0	354.6	9.1	89.1	4.6	96.7	5.1
Protected agriculture	80.8	1.0	60.2	1.5	18.7	1.0	1.9	0.1
Urban	949.9	12.3	545.9	14.0	300.0	15.4	104.0	5.5
Bodies of water	4.6	0.1	3.3	0.1	1.3	0.1	0.0	0.0
Total	7740.6	100.0	3909.4	100.0	1948.5	100.0	1882.6	100.0

).

2.2. Hydrological Simulation

To perform hydrological simulations, the WEAP used five methods. In this research, the rain runoff method (soil moisture method) was applied, which represents the basin through two soil layers: the root zone and deep zone (Figure 2). This method characterizes the vegetation cover and soil type, and through specific functions, the processes of evapotranspiration, surface runoff, subsurface runoff, percolation, and base flow are estimated [22,27,35]. The calculation of the hydrological model in the root zone and in the deep zone can be calculated with Equations (1) and (2) [27]. First,

$$S{W}_j\frac{d{z}_{i,j}}{dt}=P_e\left(t\right)-PET\left(t\right)K_{cj}\left(t\right)\left(\frac{5z_{1,j}-2z_{1,j}^2}{3}\right)-P_e\left(t\right)z_{1,j}^{\frac{LA{I}_j}{2}}-f_j{k}_s{z}_{1,j}^2-\left(1-f_j\right)k_j{z}_{1,j}^2$$

(1)

where the 1st term is the changes in the soil moisture; the 2nd factor is the effective precipitation (including irrigation); the 3rd term calculates the evapotranspiration; the 4th factor represents the surface runoff; the 5th term is the intermediate flow; and the 6th term represents percolation. Second,

$$D{W}_j\frac{d{z}_{2,j}}{dt}=\left(1-f_j\right)k_j{z}_{1,j}^2-k_2{z}_{2,j}^2$$

(2)

where the 1st term is the base flow, and the 2nd term is the percolation.

2.3. Data and Sources of Information

The monthly precipitation and the minimum and maximum temperature data from 15 nearby stations with at least 20 years of information were downloaded from the national database [36] (see Supplementary Materials: Table S1). Isotherms and isohyets were plotted according to the method of Gómez et al. [37]. Relative humidity and wind were averaged monthly, and the WEAP model was fed according to the location of the microwatersheds. The information was obtained from Gómez and Monterroso [38].

2.3.1. Hydrological Catchment Units

The subbasins were subdivided to the determine the value of the flow at a given time during the application of the model more precisely, either for calibration or for the simulation of the scenarios [24]. Following the method of Young et al. [39], the catchments were delimited using the QGIS software through the intersection of the microbasin maps and contour lines at 250 meters. The resulting catchment map was intersected with the land use and vegetation layer (map obtained in the estimation of land use changes) to specify the distribution of the different land use and vegetation types by hydrological unit. In total, 22 catchments were delimited, nine of which belonged to the Texcoco River microbasin, seven of which belonged to the Chapingo River, and six of which belonged to the San Bernardino River (see Supplementary Materials: Figure S1. Location of the Catchments).

2.3.2. Land Use Changes

Land use changes were estimated through a supervised classification using the GIS data in QGIS and the Semi-Automatic Classification Plugin (SCP) [40]. The supervised classification was based on the Normalized Difference Vegetation Index (NDVI), which highlights certain properties and allows for the spectral behavior of the forest vegetation and the soil to be differentiated. The NVDI is based on red and near-infrared reflectance, the difference of which increases as the density of green leaves increases and therefore increases with the concentration of canopy chlorophyll; therefore, it is a good indicator of the amount of green vegetation. Although the satellite images belong to the United States Geological Survey, they were downloaded directly from the SCP complement and corresponded to the following dates: 28 February 1995 (Landsat 5), 22 May 2008 (Landsat 5), and 23 March 2021 (Landsat 8). They had a spatial resolution of 30 meters per pixel and a percentage of cloudiness that was less than 20%. The original images were preprocessed starting with an assembly of false color bands: 4, 3, and 2 for the Landsat 5 images and 4 and 3 for the Landsat 8 image. The preprocessed images were used as the input data for the classification.

Ten classes of land cover were defined: urban (U), which includes buildings, streets, highways, and unused land; temperate forest (BT); reforestation (R), which refers to preferably deforested or degraded forest areas, where forest trees have been established by planting to recover forest cover; secondary vegetation (VS); grassland (P); mine (M); rainfed agriculture (AT); and irrigated agriculture (AR); protected agriculture (PA); and water bodies (CU). For each of the classes, at least 10 survey polygons were drawn manually. The supervised classification was performed with the c random forest classification complement [40], which is based on an algorithm that begins with a random selection of the predictor variables and results in a collection of identically tree-structured, independent, and distributed classifiers. Each individual tree casts a unitary vote for the most popular class, while the results of the classification are determined from the majority of the votes of each class [41]. The confidence of the results was analyzed through the evaluation report, and some errors were identified and corrected. With the maps obtained, the occupied area was estimated, and then the area that presented some type of land use change in two periods between 1995–2008 and 2008–2021 was determined (Figure 3). A land cover change tool that compares two classifications was applied to evaluate the changes in land cover. A confusion matrix was filled out to determine the changes in land use class as well as permanence.

2.3.3. Hydrological Parameters

The soil moisture method considers the root zone and the deep zone for the calculation of the hydrological balance [24]. The first includes six parameters: water storage capacity in the root zone (Sw), root zone conductivity (Ks), preferential flow direction (f), initial moisture of the upper layer (Z1), coefficient of crop (Kc), and leaf area index (LAI). The deep zone includes three parameters: water storage capacity in the deep zone (Dw), conductivity of the deep layer (Kd), and initial moisture of the lower layer (Z2). In The Supplementary Materials (Table S2), the initial values used for these parameters are shown.

2.4. Model Construction

In the WEAP platform, the following information was added: maps of the microwatersheds, riverbeds, and hydrometric stations and the centroid of the hydrological catchment units. The catchment units were used to establish the locations of the elements of the base scheme: river, hydrological unit, runoff/infiltration lines, and flow meter (Figure 4). The data structure within the catchments included the types of land use and vegetation and climate data. The values of the parameters used to calculate the water balance (Dw, Sw, f, Ks, Kd, Z1, Z2, Kc, LAI) were incorporated into the model using the Key Assumptions tool [42]. The parameters belonging to the upper soil zone were subdivided to incorporate a specific value for each type of land use and vegetation, while in the deep zone parameters, a single value per microbasin was used [42].

2.5. Calibration and Validation

The model was calibrated with information from the period of 1970–1999 and was based on the coincidence between the historical series of precipitation [36] and the time series of the flow measurement data. The latter were taken from the National Data Bank of Surface Waters (BANDAS) database [43]. The calibration period was 1980–1999 for the Chapingo River, 1980–1994 for the San Bernardino River, and 1970–1986 for the Texcoco River. These periods were selected because of the continuity in their data, and their extension provided variability in the climatic parameters (wet periods and dry periods). The temporal scale of the model was monthly because this was sufficient for the purposes of applying the model, along with the fact that much of the climate data was in this format.

In the calibration of hydrological balance parameters, the sensitivity analysis proposed by Jantzen et al. [44] mentions that the model is highly sensitive to the climate and land use parameters: precipitation, area, Dw, Kc, and LAI (Table 2). The parameters Dw, Kd, Sw, f, Z1, and Z2 were calibrated using the PEST tool, which automates the process of comparing the WEAP results with historical observations and modifying the parameters of the model to improve its accuracy [45].

The validation involved running the model using the values of the parameters that were determined during the calibration process to evaluate the predictive capacity of the model outside of the calibration period [46]. The validation period was 2001–2014 for the Chapingo River, 2001–2003 for the San Bernardino River, and 2011–2014 for the Texcoco River. As in the calibration process, the fit of the model was evaluated graphically and statistically with the measurement of the goodness of fit using the same estimators.

Measurement of Goodness of Fit

Lu and Chiang [47] mentioned that after the sensitivity analysis, the selection of the calibrated ranges and adjusted parameter values are identified based on the statistical criteria: (i) R2, the coefficient of determination, represents the linear correlation between the simulated and observed data. When the value of R2 is close to 1, it indicates a greater correlation. (ii) The Nash–Sutcliffe efficiency (NSE) represents the residuals of the measured data, and its value varies from −ꝏ to 1. An NSE value equal to 1 indicates that the simulation is the same as the observation, while an NSE > 0.5 is an acceptable value for the performance of a model [48]. (iii) Percentage bias (PBIAS) indicates whether the simulated data are overestimated or underestimated. When the PBIAS is greater than 0, the simulation is underestimated, and when the value is less than 0, the simulation is overestimated [49]. (iv) The ratio of the root mean square error to the standard deviation of the observations (RSR), which is defined as the ratio between the mean square error and the standard deviation. The smaller the RSR is, the better the simulation performance is [50]. The goodness of fit of the model was calculated using the RStudio program with the hydroGOF library [51]. The reference values that were used are shown in the Supplementary Materials (Table S3).

Table 2. Parameters, resolution, and sensitivity of land use and climatological variables.

Parameters	Abbreviation	Resolution	Sensitivity *
Area	A	Land use	High
Water storage capacity in the deep zone	Dw	Basin	High
Deep layer conductivity	Kd	Basin	Moderate
Initial moisture of the lower layer	Z2	Basin	No influence
Water storage capacity in the root zone	Sw	Soil type	Moderate
Root zone conductivity	Ks	Soil type	Moderate
Preferential flow direction	f	Soil type	Moderate
Initial moisture of the upper layer	Z1	Soil type	No influence
Crop coefficient	Kc	Land use	High
Leaf area index	LAI	Land use	High
Precipitation	P	Basin	High
Temperature	T	Basin	Moderate
Wind	V	Basin	Low
Relative humidity	HR	Basin	Low
Latitude	L	Basin	Low
Cloud fraction	FN	Basin	Low

* Sensitivity is the influence of a parameter (model input) on the rate of change of the model outputs [52]. It allows us to identify the key parameters and their required precision in the calibration [53].

Table 2. Parameters, resolution, and sensitivity of land use and climatological variables.

Parameters	Abbreviation	Resolution	Sensitivity *
Area	A	Land use	High
Water storage capacity in the deep zone	Dw	Basin	High
Deep layer conductivity	Kd	Basin	Moderate
Initial moisture of the lower layer	Z2	Basin	No influence
Water storage capacity in the root zone	Sw	Soil type	Moderate
Root zone conductivity	Ks	Soil type	Moderate
Preferential flow direction	f	Soil type	Moderate
Initial moisture of the upper layer	Z1	Soil type	No influence
Crop coefficient	Kc	Land use	High
Leaf area index	LAI	Land use	High
Precipitation	P	Basin	High
Temperature	T	Basin	Moderate
Wind	V	Basin	Low
Relative humidity	HR	Basin	Low
Latitude	L	Basin	Low
Cloud fraction	FN	Basin	Low

* Sensitivity is the influence of a parameter (model input) on the rate of change of the model outputs [52]. It allows us to identify the key parameters and their required precision in the calibration [53].

2.6. Construction of Scenarios

The scenario design was carried out by considering that the main cause of land use change is the high frequency of forest fires. Forest fire events have occurred every five years, with magnitudes ranging between 200 and 300 hectares of affectation. Thus, three scenarios were designed (

Table 3. Main characteristics of the land use change scenarios.

Scenario	Year	Trend	Characteristics
1	2021	Current	Decrease in temperate forest area and increase in reforestation areas with respect to the areas occupied in 1995.
2	2047	Positive	Increase in temperate forest area, recovery of degraded areas, and without any new forest fires.
3	2047	Negative	Decrease in temperate forest area and restoration areas; increase in urban areas, rainfed agriculture, and secondary vegetation with respect to current land use areas.

), two of which are projected to the year 2047 since this time span is the same number of years over which the land use changes were evaluated (26 years, between 1995 and 2021).

Scenario 1 was called the current scenario and refers to 2021. It incorporates the land use trends observed in the period of 1995–2021. Scenario 2 assumed that there would be no forest fires, that there would be a considerable ecosystem restoration, both by reforestation and by natural regeneration, and that there would be a decline in the agricultural surface. Therefore, it was called a positive trend scenario and was projected to 2047. Scenario 3 considered the ensuing situation of forest fires affecting between 200 and 300 hectares continuing to occur every five years. These fires would not allow the establishment of new reforestation or the natural regeneration of the forest. Considerable advances in agricultural areas would be expected to replace forested areas. Therefore, it was determined to be the negative scenario and was projected to be the scenario representing 2047. The last two scenarios included the new surface values of the land covered by each of the land use and vegetation classes.

2.7. Hydrological Balance

The scenarios were used to determine the hydrological balances in the microwatersheds. In the current scenario, the inputs (precipitation, water stored in the soil from the previous year) and the output flows (evapotranspiration, water storage in the soil, base flow, interflow, surface runoff) were analyzed. Then, land uses were replaced according to the scenario to identify the impact of these land use changes at the microwatershed level on the water balance. In this way, it was possible to quantify how these changes would influence the availability of hydrological environmental services.

3. Results and Discussion

3.1. Changes in Land Use

We must know the historical and current trends that can be observed in the spatiotemporal fluctuation of the changes in the surface of the dominant cover classes [54] to determine the effects that these changes will have on the behavior of the hydrological environmental service balance. The results of the supervised classifications for 1995, 2008, and 2021 are shown in Figure 3. In the period of 1995–2008, the land use changed in 1246 ha, which represents 16% of the total area (Figure 3d). For the period of 2008–2021, changes in land use affected 1906 ha, representing 24% of the total area (Figure 3e). Between 1995 and 2008, the area of temperate forest decreased by 57.6 ha, but in the period of 2008–2021, 323.9 ha were lost, and the sum of the changes in both periods amounted to a decrease of 16.3% (

Table 4. Area (hectares) of the land use and vegetation classes in 1995, 2008, and 2021 and the trends of change between 1995–2008 and 2008–2021.

Land Use Class	1995	%	2008	%	2021	%	1995–2008	2008–2021	Total Change	Change ha/Year
Temperate forest	2403.0	31.0	2345.4	30.3	2021.4	26.1	−57.6	−323.9	−381.6	−14.7
Reforestation	968.3	12.5	1104.6	14.3	1145.6	14.8	136.4	41.0	177.4	6.8
Secondary vegetation	137.0	1.8	300.7	3.9	438.7	5.7	163.7	138.0	301.7	11.6
Grassland	74.0	1.0	102.6	1.3	52.2	0.7	28.6	−50.3	−21.8	−0.8
Mine	147.8	1.9	124.2	1.6	212.2	2.7	−23.6	88.0	64.4	2.5
Rainfed agriculture	2530.0	32.7	2399.8	31.0	2294.8	29.6	−130.3	−105.0	−235.3	−9.0
Irrigated agriculture	664.4	8.6	483.6	6.2	540.3	7.0	−180.8	56.8	−124.1	−4.8
Protected agriculture	20.0	0.3	32.2	0.4	80.8	1.0	12.2	48.6	60.8	2.3
Urban	791.7	10.2	844.3	10.9	949.9	12.3	52.6	105.6	158.2	6.1
Bodies of water	4.3	0.1	3.2	0.0	4.6	0.1	−1.1	1.4	0.3	0.0

), mainly because the surface of the lower temperate forest was affected by two forest fires that occurred in 2012 and 2017 [55,56] (see Supplementary Materials: Figure S2). The changes in land use as a consequence of a high frequency of forest fires coincide with the results of [57]. Highly severe forest fires have the potential to interrupt a wide variety of hydrological processes and functions in forested watersheds, such as interception, infiltration, evapotranspiration, and storage [58,59], and their effects can result in an increase in surface runoff and erosion and an increase in the sediment on the riverbeds [60,61,62].

Due to the degradation of the temperate forest by fires, the secondary vegetation showed an increase of 301.7 ha between 1995–2021 (Table 4) since it can sprout quickly after a fire, and if fires are frequent, this type of vegetation can thrive for a long time, and it is unlikely to be successfully replaced [63]. A strong incidence of forest fires has led to forest areas in Mexico where secondary vegetation predominates [64,65].

There was a positive response in the recovery of the degraded area, with a total increase of 184.2 ha for reforestation areas, which will greatly encourage the infiltration levels to be much higher since by reducing surface runoff, the infiltration capacity of rain within the mineral soil will increase [66,67].

Table 5. Matrix of the dynamics of land use changes (hectares) from 1995 to 2008 and from 2008 to 2021.

Land Use Class	BT 1	R 2	VS 3	P 4	M 5	AT 6	AI 7	AP 8	U 9	CA 10
Reference class 1995	Current class 2008
Temperate forest	-	32.6	49.3	31.7	-	2.7	-	0.8	0.4	-
Reforestation	18.2	-	100.2	0.1	5.4	49.3	8.9	0.4	35.0	-
Secondary vegetation	34.7	4.3	-	0.9	-	-	2.7	-	-	-
Grassland	3.5	-	1.2	0.0	-	-	-	-	-	-
Reference class 2008	Current class 2021
Temperate forest	-	131.5	280.9	18.3	-	4.1	-	-	-	-
Reforestation	14.4	-	31.9	0.1	8.2	197.1	53.8	3.4	42.6	0.2
Secondary vegetation	26.6	106.2	-	0.9	-	42.2	0.2	0.6	0.6	0.1
Grassland	68.9	0.1	0.5	-	-	0.2	-	-	-	-

¹ temperate forest, ² reforestation, ³ secondary vegetation, ⁴ grassland, ⁵ mine, ⁶ rainfed agriculture, ⁷ irrigated agriculture, ⁸ protected agriculture, ⁹ urban, ¹⁰ water bodies.

shows the dynamics of the land use classes. In the of period 1995–2008, the temperate forest changed to reforestation and secondary vegetation at amounts totaling 32.6 and 49.3 ha, respectively, and the reforestation areas mainly changed to secondary vegetation and rainfed agriculture at the amounts of 100.2 and 49.3 hectares, respectively. For the period of 2008–2021, 280.9 ha of the temperate forest changed to secondary vegetation, and 131.5 ha of it changed to reforestation areas. The microwatershed results are presented in the Supplementary Materials (Tables S4 and S5).

3.2. Future Land Use Scenarios

Land use scenarios can consider the aspects of changes in vegetation cover created by future trends such as decreases or increases in forest areas or changes in the types of crops due to economic trends. All of these considerations should be studied so that when modeling the scenario, it is clearly known which variables and functions will be taken into account when defining the scenario [24].

Table 6. Area (%) of land use in 1995 and 2008 and scenarios 2021 and 2047.

Land Use Class	Area (%)
	1995	2008	Current_2021	Positive_2047	Negative_2047
Temperate forest	31.0	30.3	26.0	48.5	13.6
Reforestation	12.5	14.3	14.9	9.6	2.4
Secondary vegetation	1.8	3.9	5.7	0.3	4.9
Grassland	1.0	1.3	0.7	0	2.8
Mine	1.9	1.6	2.7	2.8	3.6
Rainfed agriculture	32.6	31.0	29.6	18.2	48.7
Irrigated agriculture	8.6	6.3	6.9	3.6	5.6
Protected agriculture	0.3	0.4	1.1	2.6	2.1
Urban	10.2	10.9	12.2	14.3	16.3
Bodies of water	0.1	0.0	0.1	0.1	0.1

shows the percentages of the land cover surfaces with which the proposed scenarios were designed and projected, taking the area percentages of the land use classes in 1995 and 2008 as a reference. Scenario 1 (Current_2021) mainly consists of 2021.4 ha of temperate forest, which is equivalent to 26% of the total area, with 1145.6 ha of reforestation, which corresponds to 14.9% and 2294.8 ha of rainfed agriculture land, which continues to be the land use with the most coverage, even after it decreased by 235.2 ha (most of which become urban areas). Scenario 2 (Positive_2047) is projected to have the largest area of temperate forest (3757 ha), which is equivalent to 48.5% of the total area, with 743 ha of reforestation (9.6%), and with 18.2% of the surface being rainfed agriculture. Scenario 3 (Negative_2047) is projected to have the smallest forest area, and its coverage would only cover 13.6% of the total area (1050 ha), but it would have the highest percentage of urban areas (16.3%) and 3767 ha of rainfed agriculture, which is equivalent to 48.7% of the total area. The Supplementary Materials (Table S6) present additional information on the evolution of land uses.

3.3. Model Calibration

Arnold et al. [46] said that calibration is an effort to better parameterize a model for a given set of local conditions, thus reducing the uncertainty of the prediction. Droogers and Immerzeel [21] held that calibration can be considered parameter estimation or, more generally, an optimization.

Texcoco River: Figure 5a compares the average monthly flow data observed and those simulated by the model. According to the goodness of fit, the model is very good at predicting the hydrological response of the watershed (estimators NSE = 0.98, R2 = 0.93, RSR = 0.15, PBIAS = 5.3). The model fails to simulate the peak flow rates, and the simulated rates are closer to the base flow rates. This behavior of the model differs from the results presented by Ingol-Blanco and McKinney [68], who reported that WEAP better simulates the peak runoff of the base flows in the Conchos River basin.

Chapingo River: In the microbasin, the goodness of fit indicates that the model has a very good ability to predict the hydrological response of the microbasin (NSE = 0.93, R2 = 0.81, RSR = 0.26). The simulation overestimates the flow due to its PBIAS of −13.1 since in most years, the peak flows that are predicted are above those observed. However, this percentage is considered good, and the simulated curves of the base flows are observed closer to the records reported in this period (Figure 5b).

San Bernardino River: In this microwatershed, the goodness of fit indicates that our model is very good (NSE = 0.84, R2 = 0.76, RSR = 0.38) and good according to the PBIAS of −13.5. The simulated flow shows cyclical behavior because the WEAP program allows the calculation to be performed even when there are no precipitation records in a given time. This is possible because WEAP provides the ability to model the flow with full accuracy if the data are available or to calculate them if they are not available [69]. When the precipitation data entered are the values that were calculated using isohyets, then the model fails to simulate the variations in the peak flows, and the simulation is closer to the variation in the base flows (Figure 5c).

Goodness of fit. According to the criteria set by Moriasi et al. [52], the ranges of calibration values are generally classified as being very good by all four estimators, which indicates that the model correctly simulates the hydrological responses of the three microwatersheds. Previous studies have confirmed the capacity of the WEAP model to reproduce the hydrological processes of hydrographic basins in different parts of the world. Among these is the study conducted by Abera Abdi and Ayenew [70], who reported R2, NSE, and RSR values of 0.82, 0.8, and 0.44, respectively. Asghar et al. [71] calibrated their model with an NSE and R2 of 0.85 and 0.86, respectively. Al-Mukhtar and Mutar [72] reported values of R2 and PBIAS during the calibration of 0.70 and 2.52%, respectively. Nevárez et al. [73] obtained values of 0.84, 0.73, and −15.92 for the R2, NSE, and PBIAS estimators, respectively. The validation information and the goodness of fit of the validation are presented in the Supplementary Materials (Figure S3).

3.4. Hydrological Balance

The evaluation of the water resources in a basin requires a correct estimation of the hydrological balance, that is, understanding the cycle in its different phases: the way in which the water that is received from precipitation or mist is distributed between the evapotranspiration process, runoff, and infiltration [74].

The effects of current and evolutionary trends in the land use changes of the temperate forest cover on the hydrological balance were analyzed by evaluating the behavior of the inflow (precipitation, water stored in the forest in the previous year) and the outflows (evapotranspiration, water storage in the soil, base flow, interflow, surface runoff) of the hydrological balance emitted by the model. The flows were taken as indicators to determine how the hydrological resources were distributed in the microwatersheds of the study area. The values were compared to the historical reference data to obtain the percentage change or evolution (Ev%) under the conditions of the projected scenarios in 2021 and 2047 (

Table 7. Inflow and outflow of the hydrological balance of the Chapingo River, Texcoco River, and San Bernardino River.

Scenario		Reference (mm/Year)	Current (mm/Year)	Ev (%) Current	Positive (mm/Year)	Ev (%) Positive	Negative (mm/Year)	Ev (%) Negative
Chapingo River
Inflows	Precipitation	631.3	635.8	0.7	624.5	−1.1	624.5	−1.1
	Water stored in the soil the previous year	111.0	108.5	−2.3	126.3	13.8	74.5	−32.9
Outflows	Evapotranspiration	532.2	534.7	0.5	536.0	0.7	497.7	−6.5
	Water stored in the soil	159.3	153.3	−3.8	165.8	4.1	127.6	−19.9
	Base flow	1.5	4.2	180.0	9.3	520.0	9.1	506.7
	Inter flow	10.9	12.3	12.8	14.0	28.4	13.1	20.2
	Surface runoff	38.4	40.0	4.2	25.9	−32.6	64.6	68.2
Texcoco River
Inflows	Precipitation	635.7	635.3	−0.1	632.8	−0.5	632.8	−0.5
	Water stored in the soil the previous year	128.5	130.7	1.7	139.4	8.5	125.3	−2.5
Outflows	Evapotranspiration	532.3	533.7	0.3	538.7	1.2	513.1	−3.6
	Water stored in the soil	204.2	199.8	−2.2	204.2	0.0	186.4	−8.7
	Base flow	0.8	4.1	412.5	10.8	1250.0	9.4	1075.0
	Inter flow	2.3	2.7	17.4	3.1	34.8	3.7	60.9
	Surface runoff	24.6	25.9	5.3	15.2	−38.2	45.6	85.4
San Bernardino River
Inflows	Precipitation	645.0	645.0	0.0	645.0	0.0	645.0	0.0
	Water stored in the soil the previous year	72.4	70.7	−2.3	102.4	41.4	54.1	−25.3
Outflows	Evapotranspiration	544.6	544.7	0.02	558.5	2.6	518.5	−4.8
	Water stored in the soil	114.8	108.1	−5.8	140	22.0	91.8	−20.0
	Base flow	1.1	2.8	154.5	6.8	518.2	6.7	509.1
	Inter flow	12.7	13.5	6.3	15.3	20.5	15.9	25.2
	Surface runoff	44.3	46.8	5.6	27.1	−38.8	66.4	49.9

).

3.4.1. Current Hydrological Balance

The results of Scenario 1 (Current_2021) indicate that the decrease in the the temperate forest cover (−15.9%) from 2008 as a result of a frequent incidence of forest fires influenced the increase in the surface runoff in the three microwatersheds: 4.2%, 5.3%, and 5.6%. Thus, the microbasin of the San Bernardino River (smallest of the three microbasins studied) presents the greatest increase in runoff at present, in line with the fact that it only has 5.5% temperate forest. coverage There was also an increase of 153.3% in the base flow in the Chapingo River, 412.5% in the Texcoco River, and 154.5% in the San Bernardino River. This increase is related to the increase in the reforestation areas (18.3%) estimated for this scenario. It is worth mentioning reforestation takes more than 5 years, which is the minimum time necessary for a basin to reach a new equilibrium [75]. Therefore, these changes are compensation for the lost forest area. According to the study by [76], there are increases in the base flow that are related to larger forest cover areas due to greater the infiltration and recharge of the underground storage.

3.4.2. Future Projections

The simulation of the hydrological balance in Scenario 2 (Positive_2047) mainly implies an increase in the temperate forest cover (85.9%) over what it currently covers. The increase assumes that there will be no forest fires with significant effects in the next 26 years, a condition that would allow considerable restoration of the ecosystem through both reforestation and natural regeneration. It assumes a decrease in the rainfed agricultural surfaces (−38.7%) that are near or within the preferably forested areas. The results indicate that under these conditions, a decrease in surface runoff is expected in all three microwatersheds: Chapingo River −32.6%, Texcoco River −38.2%, and San Bernardino River −38.8%. These decreases indicate that there would be lower velocity runoff, a condition that occurs when the surfaces are protected by temperate forest cover [67,77].

The projected inflows suggest that there would be an increase in the previous-year soil water storage of 13.8%, 8.5%, and 41.4% in the three respective microbasins, which would increase the water available for infiltration and for the recharge of aquifers [78]. The highest percentages of base flow occur in this scenario: 520% in the Chapingo River, 1250% in the Texcoco River, and 518% in the San Bernardino River microbasin. The increase in base flow corroborates the studies by Price and Jackson [79] and Price et al. [76], who evaluated the base flow of 30 streams in the highlands of the Appalachians, and their results indicate a positive relationship between the forest coverage of the basin and the discharge of the base flow. Increases in the evapotranspiration of 0.7%, 1.2%, and 2.6% would also be expected, and the latter outflow is directly related to the increase in forest cover [80]. However, according to the study reported by Qiu et al. [81], although the increase in forest cover, which is mainly due to the establishment of reforestation, has positive responses with the increase in humidity in the superficial layer of the soil, in the subsoil and in the deep layers of the soil, they led to a significant reduction in the humidity; this negative effect is due to the fact that large-scale restorations consume much more water from the deep soil than natural regeneration; therefore, it will be difficult for local precipitation to replenish the decrease in humidity in the deep zones, and, in turn, this will have negative effects on plant growth and water resources.

The hydrological balance scenario, Scenario 3 (Negative_2047) assumes a 50.1% decrease in temperate forest cover due to the persistence of highly severe forest fires that would result in an 81.7% decrease in reforestation, preventing the establishment of the natural regeneration of the forest. Therefore, in degraded areas, there would be a considerable advance in rainfed agricultural areas (64.2%), which would take the place of what would preferably be forested areas. Under these conditions, the model indicates that surface runoff would increase [54,57] in the three microbasins: Chapingo River (68.2%), Texcoco River (85.4%), and San Bernardino River (49.9%). The elimination of vegetation cover and soil organic matter will lead to changes in the hydrological processes by reducing the interception of precipitation and modifying the structure of the surface soil [14,15,82]. The microbasin of the Texcoco River presents the greatest increase in surface runoff because it would be the one with the greatest loss of forest (687.5 ha). A decrease in the water stored in the soil from the previous year of 32.9%, 2.5%, and 25.3%, respectively, would be expected. The alteration of these flows indicates that the water that is available for hydrological environmental services for the recharge of aquifers will be compromised throughout the study area [83].

The decrease in evapotranspiration of 6.5% in the Chapingo River microbasin, 3.6% in the Texcoco River, and 4.8% in the San Bernardino River is also related to the alteration of the hydrological balance [57,84]. The results coincide with the effects of land use changes on the hydrological components in the Chongwe River basin presented by Tena et al. [85], showing that the actual annual evapotranspiration decreased from 840.6 mm to 796.3 mm due to a decrease of 41.11% in forest cover.

3.4.3. Environmental Hydrological Service in Temperate Forests

The analysis of the hydrological balance according to the catchment groupings where the temperate forest predominates allowed us to more clearly identify the variations in and the distribution of the inflow and outflow of the hydrological balance in the three microbasins related to the changes in the temperate forest cover.

Table 8. Inflow and outflow of the hydrological balance in catchments where the cover of the temperate forest predominates in the Chapingo River, Texcoco River, and San Bernardino River.

Scenario		Reference (mm/Year)	Current (mm/Year)	Ev (%) Current	Positive (mm/Year)	Ev (%) Positive	Negative (mm/Year)	Ev (%) Negative
Chapingo River
Inflows	Precipitation	769.6	769.6	0.0	769.6	0.0	769.6	0.0
	Water stored in the soil the previous year	157.7	152.4	−3.3	159.0	0.8	124.7	−20.9
Outflows	Evapotranspiration	657.4	653.8	−0.5	660.6	0.5	620.5	−5.6
	Water stored in the soil	245.7	237.3	−3.4	249.2	1.4	208.7	−15.1
	Base flow	0.91	0.92	1.4	1.0	10.4	1.0	6.0
	Inter flow	10.8	10.4	−3.6	10.9	1.0	8.4	−22.4
	Surface runoff	12.5	19.7	57.1	6.9	−45.2	55.8	346.3
Texcoco River
Inflows	Precipitation	769.6	769.6	0.0	769.6	0.0	769.6	0.0
	Water stored in the soil the previous year	138.0	137.9	−0.05	142.3	3.1	119.7	−13.3
Outflows	Evapotranspiration	634.7	634.5	−0.04	643.1	1.3	597.0	−6.0
	Water stored in the soil	254.5	246.9	−3.0	258.2	1.4	223.1	−12.3
	Base flow	1.2	2.9	136.3	2.0	67.3	1.8	43.5
	Inter flow	0.0	0.0	0.0	0.0	0.0	1.1	0.0
	Surface runoff	17.2	23.3	35.9	8.6	−49.6	66.5	287.4
San Bernardino River
Inflows	Precipitation	769.6	769.6	0.0	769.6	0.0	769.6	0.0
	Water stored in the soil the previous year	155.9	143.9	−7.67	160.6	3.0	102.5	−34.2
Outflows	Evapotranspiration	646.7	638.5	−1.27	651.7	0.8	583.0	−9.8
	Water stored in the soil	245.5	225.0	−8.34	249.3	1.6	179.1	−27.0
	Base flow	2.1	2.1	1.39	2.8	34.0	2.4	15.0
	Inter flow	20.6	18.7	−9.05	21.1	2.4	18.1	−12.1
	Surface runoff	10.8	29.4	172.24	5.4	−49.8	89.7	729.6

shows the distributions of the inflow and outflow of the hydrological balance by microbasin as well as the percentage of the evolution (Ev%) of the change in the flows under the three proposed scenarios. The percentage of evolution was generated based on the flows that existed in the reference scenario.

In the inflow of the three microwatersheds, the decrease in the water stored in the soil in the previous year under the conditions of Scenario 3 (negative) stands out at −20.9% in the microwatershed of the Chapingo River, −13.3% in the Texcoco River, and −34.2% in the San Bernardino River because changes in the soil cover strongly affect the distribution of soil moisture and the hydraulic properties of the soil [86,87]. In all of the outflows of the scenarios of the three microwatersheds, there is variation in the distribution of water resources, but the decrease in water storage in the soil stands out: −15.1% in the Chapingo River, −12.3% in the Texcoco River, and −27% in the San Bernardino River. The capacity of the soil to store and regulate the flow of water largely depends on its infiltration rate and depth, but disturbances such as changes in land use and fires can alter these soil capacities [88]. These alterations are caused because without forest cover to protect the soil, the soil is exposed to intense rainfall that can cause water erosion, compaction [89], and alterations in their structure (higher apparent density, lower field capacity, disappearance of organic matter and microfauna) [88]. The decrease in soil moisture coincides with the statistical results of a study along the karst slopes of Southwest China, which showed that changes in soil cover strongly affect the distribution of soil moisture. When compared to bare soil areas, the forest, shrub, and grass areas had 30.5, 20.1, and 10.2% higher soil moisture values, respectively [86].

Regarding the evolution of the changes that the surface runoff presents and would present, it can be observed that currently (Scenario 1), this flow has increased by 57.1% in the Chapingo River, 35.9% in the Texcoco River, and 172.2% in the San Bernardino River microwatersheds. However, under the conditions of the negative scenario, it is expected that surface runoff will increase 346.3% in the Chapingo River, 287.4% in the Texcoco River, and 729.6% in the San Bernardino River. These values would decrease under the conditions of Scenario 2 (positive): −45.2%, −49.6%, and −49.8%, respectively. The decrease in surface runoff is directly associated with the rain that is intercepted by the canopy because it is an important process in the water balance of forests [90]. The temperate forest canopy is also an important factor for soil conservation due to its role in reducing the erosive impact of rainfall [91]. Another factor that influences the reduction of surface runoff is the incorporation of organic matter into the soil by the temperate forest since its presence at high values reduces the risks of soil water erosion [92] and promotes water infiltration to the root zone and the deep zone. Therefore, the hydrological functions, as highlighted by Esse et al. [93], are more efficient in watersheds with forest cover and that are made up of native or exotic species than in agricultural cover basins with annual crop rotations.

4. Conclusions

In the study area, the temperate forest land area decreased by 16% from 1995–2021. The incidence of forest fires has been an important cause of this decrease. The presence of secondary vegetation is an indicator of the degradation and recovery of the ecosystem, as it can sprout quickly after a forest fire. Changes in land use and the presence of forest fires are included in the WEAP software to measure the hydrological response of watersheds. The WEAP model is a tool that can be incorporated into the forestry sector to determine the current and future behavior of hydrological resources and to support decision-making for the integral management and valuation of hydrological environmental services. Doing so would ensure the integral, continuous, and stable flow of environmental services provided by temperate forests.

The effect of maintaining a healthy forest without forest fires translates into lesser surface runoff, lower runoff velocity, and more water being available in the soil for infiltration and aquifer recharge. That is, it maintains the hydrological environmental service. Forest fires reduce rain interception and increase the velocity and flow of surface runoff, which modify the structure and composition of soils, thereby compromising hydrological environmental services. New studies are suggested to investigate the modifications of the water balance during forest fires.