Methods

For this work were gathered topographic and soil charts, from which was generated the information to feed the soil loss prediction model.

With 25 vector topographic charts with 10 m equidistance contours to offset curves, and with the Envi 4.5 software was generated the Digital Elevation Model (MDE) of the study area with spatial resolution of 1 m. The procedure for get the subbasins consisted of using the Arc Hydro tool in Arc GIS 10.3.

The Modified Universal Soil Loss Equation (MUSLE) was developed to empirically estimate soil loss for small basins by changing the erosivity factors and the transport efficiency of sediments produced by the rain kinetic energy Universal Soil Loss Equation (USLE) (^{Wischmeier & Smith, 1978}) for a function of the storm's runoff volume product from peak runoff (^{Williams, 1975}). This equation is applied for average annual soil loss or to estimate daily loss, for this study the annual loss was estimated, the equation of the MUSLE is as follows:

Where:

The procedure for applying the MUSLE model in each of the subbasins is as indicated in the diagram of Figure 1, and the procedure for obtaining the maps of each of the factors in the equation is as described below.

It was used the procedure for assessing the K factor using the Food and Agriculture Organization of the United Nations (FAO) method to determine the soil units within the study area and their texture according to the classification of the National Territory Studies Commission (CETENAL), then with soil order and texture data, the FAO´s table was used to obtain K.

To obtain the L factor it was used the Equation (2) recommended by the methodology of the Soil Loss Prediction Manual (^{Figueroa et al., 1991}):

Where:

To obtain the S factor, which involves the slope, it is based on the figure S_t out with the Equation (3):

Where:

The S factor is based on S _t and is calculated using the Equation (4):

Generally factor C is given in terms of the average annual value for a combination of culture systems, management and rainfall. In this case the relative soil losses due to factor C were obtained from tables as recommended by ^{Wischmeier and Smith (1978)}.

To calculate de runoff was used the Curve Number method, which was developed by the United States Soil Conservation Service (USSCS).

In this method each type of soil is assigned a hydrological group, which are based on the infiltration capacity, which relates to the hydraulic conductivity of the surface of the saturated soil.

From the soil charts 1:250, 000 and the information of the soil profiles was generated the hydrological group map that is based on the type of soil, texture and other physical properties of the soil.

The land use map was obtained from the supervised classification of satellite images, for which the maximum likelihood algorithm was used as a decision rule for the classification of images in the two periods.

This algorithm has its basis on assuming that samples from each class and each band are normally distributed, so the decision rule for an unknown X pixel is:

Assign X to class c only if probability is met;

Where:

Where:

Thus, to assign the X vector of an unknown pixel within a class, the decision rule first calculates the average class pixel, then assigns it to the one with the largest probability value (^{ITT, 2008}).

The runoff was estimated by the Curve Number (NC) method with the Equation (6):

Where:

To use this equation was generated an annual precipitation map of each period with information from the ERIC III Rapid Weather Information Extractor (^{IMTA, 2013}), and with the information about hydrological group and soil use derived from the classification of images, the NC map was generated, and with that the map of S. The value of Q was transformed to m³ to get drained volume. The Equation (7) was used to calculate peak flow on each subbasin:

Where:

Thus, when calculating the parameters of Equation (7) and resolving it in each subbasin the peak flow q_p raster map was generated.

With the parameters of Equation (1) represented through raster maps, the MUSLE model was implemented in the Raster Calculator tool in ArcMap 10.3, this tool integrates the map algebra technique so that it generated the erosion map in the study periods.

The original model of the MUSLE is y=11.8(Q*qp)0.56KCSL and the way it was implemented in the map algebra module was like: Yp = 11.8*PowQ_m3_ras*qp_m3_ras,0.56*K_ras*C_ras.*SL_ras[/p]

Digital elevation model

The resulting Digital Elevation Model (DEM) has the spatial precision of one meter, it was generated in the ENVI 4.5 software. The DEM accuracy is important because boundaries are generated to define the hydrological subbasins. The maximum elevations of 4,680 meters are located approximately to the center of the DEM that correspond to those of the volcano crater and the northwestern limits, corresponding to the Sierra de las Cruces; while the lower elevations are to the southeast (Guerrero State) and southwest (Morelos State).

Hydrologic subbasins

The subbasins were generated using the DEM in a radius of 60 km within the influence area of the Nevado de Toluca, 52 of them resulted from this process as shown in Figure 2, which were characterized geometrically as defined by the relief of the region corresponding to mountains with dissection, mountain ranges and recent of them with shields. In subbasins the average area is 183 km², the average slope is 12%, the drainage network is third order and the predominant shape is elongated.

MUSLE equation Factors

The maps resulting from the K, LS, C factors have spatial resolution of 10 meters, which corresponds to the píxel size of the satellite image in map algebra.

The K factor map shows higher values (0.079) in the subbbasins southwest of the Nevado de Toluca, which correspond to the maximum slope values, while there were the minim values in the subbasins that descend from towards Toluca valley.

The LS factor by inclination and slope length is 0.995 in areas where the topography is rugged and match to the subbasins that descend from the Volcano and to the southern area, where the surface is characterized by its rugged relief; while in Toluca valley there are values of 0.34 that become to the subbasins to the north.

Factor C due to vegetal cover and soil management is related to the soil use map. The C-factor map has values ranging from 0 to 1 having lower values in places where there is well covered forest as in subwatersheds; Río del valle, Toluca, Rio Verde and Progreso, which descend from the Nevado de Toluca and the Sierra de las cruses (northeast of the study area), while in the valley and towards Temazcalcingo the C values are close to or equal to 1 due to the limited plant coverage , where are located the subbasins of Progreso-Huautla, Temascaltepec, Valle de Bravo, San Felipe and Medio Ixtapan.

Soil use and hydrological group

Figure 3 shows the distribution of soil hydrological groups, where 4 groups were identified that vary according to the physical characteristics of the soil. The dominant group is B with 66.9% of the area, which corresponds to medium textured soil with moderate permeability in wet state, followed by group D with 14.7% corresponding to very impervious soils that cause runoff and therefore erosion.

The result of the digital classification of satellite images shows 6 types of land use, which are observed in Figure 4 for the two periods analyzed, corresponding to: 1) freshwater; 2) forest; 3) matorral; 4) bare soil; 5) agriculture, and 6) urban areas. The classes were defined according to the INEGI land use nomenclature.

The Figure 4 (right) shows the change in land use in the southern and southeastern site of the study area, which is due to deforestation and the Rosa Tumba y Quema (R-T-Q) system for agricultural activity; it is also associated with fires caused by high temperatures, since the ^{Diakite study (Diakite, 2008)} found significant increase in temperature in the area surrounding the volcano, It also coincides with studies where the significant increase in temperature is the factor that detonates fires as indicated by ^{González, Lara, Urrutia and Bosnich (2011)}; ^{Yong (2003)}; ^{Bonebrake et al. (2014)}; ^{Guariguata (2009)}; ^{Laucirica, Mancino, Uboldi and Michalijos (2011)}, and ^{Huerta-Martínez e Ibarra-Montoya (2016)}.

Runoff

For estimate the runoff by the Soil Conservation Service (SCS) method, the Equation (6) was applied for that the Curve Number (CN) maps were obtained as shown in Figure 5. The larger CN values favor runoff by the hydrological condition of the soil and its very poor coverage. The variation in the CN value from one period to another depends basically on the change in land use, such as plant cover, as the physical characteristics of the soil were not changed over time.

The runoff generated with the Curve Number method for each period shows spatial variation due to the hydrological group factors and plant coverage involved in its calculation, this variability ranges from 52 to 108 m³/ha in 1989 and from 41 to 131.3 m³/ha in 2014. The variation of runoff is reflected in erosion; the areas of least runoff are located in the subbasins that descend from the Nevado de Toluca and the Sierra de las Cruces, and the areas of greatest runoff are located in the subbasins of Progreso-Huautla south of the study area.

Erosion map

Erosion maps for each period are shown in Figure 6; the 2014 map has higher erosion values in t/ha; for 52 subbasins in the studied area 36 showed a positive increase in erosion compared to 1989.

Subbasins with negative values are due to a favorable change in their vegetal cover conditions, i.e. they were regenerated some forests or grasslands that were used for agriculture.

Based on erosion levels (^{FAO, 1980}) and geodynamic classification, the subbasins studied are grouped into two units:

The increase in erosion occurred mainly in the valley subbasins, as well as those descending from the crater where the plant cover area was reduced leaving the soil unprotected; meanwhile, the reduction in the rate of soil loss was identified to the northwest and southwestern of the volcano.

According to the soil loss values, there is mild erosion, as they occur 0 to 10 t/ha/year and is characterized by being laminar type with few gutters, presents up to 25% of the horizon "A" eroded so it corresponds to the stable geodynamic environment.

Moderate erosion is also present in areas where 10 to 50 t/ha/year occur, is gutter type, presenting 25% to 50% of the thickness of the eroded "A" horizon corresponding to the penestable geodynamic medium.