Description of the study place

This research was carried out in the El Malacate micro-basin with an extension of 149.2 ha. The micro-basin is part of the Pátzcuaro Lake closed watershed, municipality of Tzintzuntzan, state of Michoacán de Ocampo, with a mean altitude of 2300 m, between 19° 36’ 50’’ and 19° 38’ 30’’ N, and 101° 35’ 10’’ and 101° 36’ 20’’ W. According to the Köppen climate classification, modified by ^{García (2004}), the climate in the area is temperate sub-humid C (w2) with annual mean temperature between 12 °C and 18 °C, temperature of the coldest month between -3 °C and 18 °C, and temperature of the warmest month between 20°C and 30°C, and an annual precipitation of 200 to 1800 mm with a mean value of 920.3 mm (^{IMTA, 2009}).

Experimental data

The cartographic information used was the digital elevation model (DEM), map of soil and vegetation uses, and map of the type of soil by INEGI at a scale of 1:50 000. The map of land and vegetation uses was verified in the field where a predominance of forest use was identified. To verify the type of soil, a visit was carried out along the watershed, and the type of soil nitic Luvisol was found, based on the FAO-UNESCO classification, made up of four layers, with a superficial layer of 25 cm of loam clay sand texture (32.5 of clay, 51.3 of loam, and 6.09 % of very fine sand), and fine granular structure with 4.98 % of organic matter.

The daily precipitation in 2013 was recorded, appraised at intervals of 0.2 mm with a HOBO RG3-M pluviograph, with an annual total of 769.6 mm. The hydrographs were also used (Figure 1), as well as the amount of sediments at the exit of the basin, obtained from 14 runoff events recorded in 2013. The hydrographs were obtained through a long-throated flume with an ultrasonic sensor to measure the water depth, and an automatized system to record the information. The sediment production was obtained by quantifying the amount of solids in suspension from samples taken at the exit of the basin during the runoffs.

Figure 1 Hydrographs of the runoff events in 2013. 

Estimation models in the production of sediments

USLE Model The Universal Soil Loss Equation (USLE) was proposed to compute longtime average soil losses from sheet and rill erosion under specified conditions (^{Wischmeier and Smith, 1978}), and is defined as:

(1)

where, according to the international metric system of units (^{Foster et al. 1981}), A is the annual erosion rate per unit of area (Mg ha^-1 year^-1), R is the rainfall and runoff erosivity factor (MJ mm ha^-1 h^-1 year^-1), K is the soil erodibility factor (Mg ha^-1 MJ^-1 mm^-1), L is the slope length factor (dimensionless), S is slope steepness factor (dimensionless), C is the cropping management factor (dimensionless), and P is the support practice factor (dimensionless).

MUSLE model

The MUSLE model (^{Williams, 1975}) is a modification of the USLE model, and it consists in replacing the R factor of erosivity of the rain because of superficial runoff and the peak runoff rate of a storm, with the aim of calculating the sediment yield per event and expressed as:

(2)

where Y is the sediment yield in a specific event (Mg event^-1), Q_s is the runoff volume (m³) and q_p is the peak runoff rate of the event (m³ s^-1).

R factor

This factor represents the potential capacity that rain drops have to cause erosion and it was calculated in two ways: through the EI₃₀ index (^{Wischmeier and Smith, 1978}), and with an equation adapted by ^{Figueroa et al. (1991}), who consider a quadratic functional relationship with the mean annual precipitation registered.

The EI₃₀ method was used because there is daily information about precipitation, registered at intervals of minutes. This method is defined as the product of the total kinetic energy of a rain event by the maximum intensity in 30 min, so that the total value of the factor of the rain erosivity R is the sum of values of erosivity of individual storms, calculated as follows: (^{Wischmeier and Smith, 1978}):

(3)

where E is the kinetic energy of the rain from each storm (MJ ha^-1 mm^-1), I₃₀ is the maximum intensity of rain in any 30 min from each storm (mm h^-1), and n is the number of storms in the year.

The kinetic energy of each storm is the eddition of the kinetic energies of the rain intervals, differentiated by their rain intensities, and it was calculated with the equation proposed by ^{Wischmeier and Smith (1978}):

(4)

(5)

where e_j is the kinetic energy of the interval of rain j of a storm (MJ ha^-1 mm^-1), and I_j is the intensity of the rain from interval j (mm h^-1), and pp_j is the precipitation of the interval j (mm). For intensities greater than 76 mm h^-1 an e value of 0.283 was established.

According to the criterion to calculate EI₃₀ proposed by ^{Wischmeier and Smith (1978}), the storms with precipitations lower than 12.5 mm with rain periods separated more than 6 h were not considered, unless 6.2 mm or more had fallen in 15 min (I₁₅≥24.8 mm h^-1). The existence of rain periods in long storms was defined when the sheet was less than 1.2 mm in a lapse of 6 h (^{Renard et al., 1997}).

The R factor is calculated with the equation adapted by ^{Figueroa et al. (1991}) as a methodological alternative to the watersheds where it is not possible to calculate the EI₃₀ for lack of detailed information from the precipitation throughout the year. The watershed of study is located in the hydrological zone V (^{Bravo et al., 2009a}), with the following equation:

(6)

where P is the annual precipitation of 2013 in millimeters.

K factor

The term erodibility of the soil is used to indicate the susceptibility of a soil to erosion. It is defined as the rate of soil loss for each additional unit of EI₃₀ when the slope and its length, the plant cover, and the conservation practices of the soil remain constant and its values are equal to one (^{Wischmeier and Smith, 1978}).

As the sum of the loam particles plus the very fine sand was less than 70 % of the total of soil particles, the K factor was calculated with the equation proposed by ^{Wischmeier and Smith (1978}):

(7)

where M is a parameter of the particle size, MO is the percentage of organic matter (%), C _estr is a code of soil structure, C_perm is a code of soil permeability, and 0.313 is a conversion factor of the system of English units to the international metric system (^{Foster et al, 1981}). The parameter of the particle size was calculated as:

(8)

where m_l is the percentage of loam content (diameter of particles of 0.002 to 0.05 mm), m_amf is the percentage of the content of very fine sand (diameter of 0.05 to 0.10 mm), and m_a is the percentage of the content of clay (diameter of particles smaller than 0.002 mm).

The value of the soil permeability was determined in function of their texture using the functional relations described by ^{Saxton and Rawls (2006}).

Another method was that of ^{FAO (1980)}, which is used in México, since the K factor is determined only in function of the unit to which the soil belongs in the FAO-UNESCO classification and of the texture of the superficial layer.

LS factor

This factor represents the effect of the topography on soil erosion. The erosion increases as the length (L) of the plot increases in the direction of the slope, and the inclination of the surface (S) increases, calculated as (^{Wischmeier and Smith, 1978}):

(9)

where:

(10)

(11)

where λ is the plot slope length (m), s is the slope of the plot (%), and the m exponent of the factor of the slope length considers the proportion of erosion caused by the erosion between drains, due to the impact of the drops, and that caused in drains by the dragging force of the flow (^{Renard et al., 1997}). The value of m varies in function of the slope of the plot and was calculated with the equation proposed by McCool (^{Renard et al., 1997}):

(12)

where θ is the angle of slope.

The classical equation to calculate the LS factor, proposed by ^{Wischmeier and Smith (1978}), does not represent well the effect of the topography on the soil erosion for slopes greater than 16 % (^{Pongsai et al., 2010}), so that the S factor was calculated with the equation proposed by McCool (^{Renard et al., 1997}) for slopes greater than 9 %:

(13)

where θ is the angle of the slope.

Evaluation of the model’s adjustment

The correspondence between the series of values observed and those estimated are expressed in terms of four indexes: Model’s efficiency (EF) proposed by ^{Nash and Sutcliffe (1971}) and used by ^{Besteiro and Gaspari (2012}), and ^{Muche et al. (2013}); standard error (ES) (^{Pongsai et al. (2010}); coefficient of determination (CD) (^{Rosenmund et al., 2005}); and percentage of bias (PBIAS) (^{Moriasi et al., 2007}). The equations of the indicators are given by:

(14)

(15)

(16)

(17)

where O_i are the values observed, P_i are the values predicted, Ō is the average of the values observed, and n is the number of data evaluated.

The model increases the precision of the prediction of the values observed as the values of the indicators EF, ES, CD and PBIAS approach one, zero, one, and zero, respectively.

The correspondence between the pairs of data observed and estimated was made through the relative estimation error, calculated as (^{Pongsai et al., 2010}; ^{Arekhi et al., 2012}):

(18)

The optimal value of ER is zero; positive values indicate the percentage of over-estimating of the model, while negative values quantify the percentage of sub-estimating of the model.

Calibration of the MUSLE model

Figures 2 and 3 show a potential relationship between the production of sediments and the runoff volumes and peak runoff rate of each storm, as ^{Williams (1975}) points out. The adjustment of these relationships was R²>0.95, but the standard deviations were high; when relating the sediment production with the runoff volume, the standard error was 8.2 and 6.6 Mg when it depended on the peak expenditure, as ^{Sadeghi et al. (2007}) reported.

When relating the production of sediments with the product of the runoff volume by the peak runoff rate, term by which ^{Williams (1975}) substituted the erosivity factor of R rain of the USLE, there is good adjustment in a potential-type model (R²=0.983; Figure 4). However, the value of the power (0.703) is higher than the one found by Williams (1975) in the US (0.56), and the standard error was 4.4 Mg. A model of linear type also presented an adjustment R²=0.945 with a lower standard error to the potential type (3.1 Mg), but with the disadvantage of not representing the phenomenon for products Q_sxq_p close to zero (Figure 4).

Figure 2 Sediment production in function of the runoff volume. 

Figure 3 Sediment production in function of the peak runoff rate. 

Figure 4 Sediment production in function of the runoff volume and the peak runoff rate. 

The errors made by the potential relationship between the sediment production and the product Q_s xq_p, added to the estimation errors of the K, LS and C factors, led to strong errors between the values observed and predicted by the MUSLE model. In fact, the results from Table 3 and the statistical indicators of columns 3 and 4 from Table 4 show the need for adjustment of the factors from the MUSLE model.

Table 3 Relative errors in the estimation of sediment production with the MUSLE model. 

^† The K factor was calculated with Equation (7). ^¶ The K factor was calculated with the FAO methodology.

Table 4 Statistical indicators of the MUSLE model predictions before and after calibration. 

According to the total results from sediment production (Table 3), it was found that the MUSLE1 model (K factor calculated with Equation 7) underestimated the results in 12 of the 13 rain events, and the MUSLE 2 model (K factor calculated with the methodology by FAO) underestimated the results in all the events. The results have the tendencies reported by ^{Sadeghi et al. (2007}), who found underestimations by the MUSLE model for short-lasting small storms. This is explained because the MUSLE model was developed for grassland watersheds with larger storms and slopes smaller than those of the watershed in study (^{Williams, 1975}).

Both MUSLE models underestimated the results, but it is notable that the FAO model is less precise and the relative estimation error calculated with Equation (18) for the annual data reported in Table 3 was -68.2 %, while with the other model it was -87.8 %. These results imply that the FAO model did not represent adequately the K factor and it is more advisable to use Equation (7) because, in addition to considering the soil texture, it considers the structure, permeability and organic matter content.

For the calibration of the model, only one of the factors is adjusted and it tends to be the one for which there is little information, or in which the conditions of the site of study are not contemplated in the tables, nomograms and equations suggested in publications about the subject. This is indicated by ^{Pongsai et al. (2010}), ^{Besteiro and Gaspari (2012}) and ^{Sadeghi et al. (2014}).

Since the P factor is well-defined with a value of one because there are no management practices, the topography changes very little in time; the characteristics of the land and vegetation uses remain similar throughout the year because it is a watershed that is mostly forest and the growth of plant species is slow. In our study, we set out to adjust the K factor due to the dynamic behavior of the soil, which during the year acquires different degrees of humidity and vulnerability to erosion (^{Kinell, 2015}). Proof of the difficulty to represent the value of the K factor is the great difference in the estimations between the model that included the K value calculated with Equation (7), and the one obtained with the FAO method.

The value of the calibrated K factor was 0.107 and was obtained through the technique of least squares regression. In México there are no experimental values of K from similar watersheds for their comparison. In micro-basins of the Pátzcuaro Lake watershed and other neighboring ones, studies were performed in runoff plots of some crops with different types of farming, but in most of the cases the calibrated factor corresponds to crop C (^{Tapia et al., 2002}; ^{Bravo et al., 2009b}) and it is assumed that the other factors of the USLE model are the correct ones and calculated with the methodology proposed by ^{Wischmeier and Smith (1978}). Table 4 shows that the values of three (EF, CD and PBIAS) of the four statistical indicators of adjustment are quite close to their optimums. The value next to the unit of the Nash and Sutcliffe (EF) indicator establishes that the calibrated model is efficient. The ES value implies that there was a standard error of 4.78 Mg in the production of sediments and given the 149.2 ha of the watershed, it translates into a standard error of 0.03 Mg ha^-1. The CD value next to the unit indicates that the values predicted suggested the trend of those measured, compensating the errors; thus, in some events the results were overestimated and in others underestimated (Figure 5). The final balance resulted in an average underestimation of 0.17 Mg (PBIAS).

Figure 5 Sediment production measured and simulated with the calibrated MUSLE model. 

The calibrated MUSLE model can be used to estimate in an approximate way the production of sediments in other watersheds of similar conditions to those of our study, and where there are no gauge data as observed by ^{Prabhanjan et al. (2015}) with the application of the SWAT module, which is supported by the MUSLE model (^{Ayana et al., 2012}). There are similar basins, as is verified in the results obtained in a discharge basin into the Tzurumútaro drainage in the state of Michoacán, also mostly forest, and the production of sediments was 0.8 Mg ha^-1 year^-1 from 2007 to 2010 (^{Huerto and Vargas, 2014}), which is similar to the 1.2 Mg ha^-1 year^-1 measured in the micro-basin of our study.

Calibration of the USLE model

For the review of the USLE model results, the values of the K, LS, C and P factors of the calibrated MUSLE model were taken up again, and the measured value of the R factor, since the MUSLE model predicts the soil erosion better (^{Muche et al., 2013}). In the methodology established by ^{Wischmeier and Smith (1978}) and in the revision of the USLE model by ^{Renard et al. (1997}), it is not specified clearly whether the criterion to calculate the EI₃₀erosivity index must be applied only for the storms where there are runoffs or for all the storms; it is also not specified whether the erosivity of storms that do not satisfy such a criterion for the calculation of EI₃₀ but which did produce runoffs should be contemplated.

If the EI₃₀ of the storms with presence of runoffs from Table 2 (1154.9) and without the presence of runoffs (992.1) are considered, resulting from adding the values outside the parentheses of Table 5, the total annual value in 2013 was 2147.0 MJ mm ha^-1 h^-1, quite near to the value of 2163.2 obtained with Equation (7), proposed by ^{Figueroa et al. (1991}). The good representation of the R factor of the estimative method is evident for the year of study, with an error of 0.75 % with regards to the total measure of events with or without runoffs (2147.0 MJ mm ha^-1 h^-1). However, if the mean annual precipitation of 920.3 mm (^{IMTA, 2009}) is used to calculate R, according to ^{Figueroa et al. (1991}) instead of using the one recorded in 2013, the R value is 2464.7 MJ mm ha^-1 h^-1 and higher than 14.8 % to the measured one (2147.0 MJ mm ha^-1 h^-1).

Table 5 Monthy values of EI₃₀ (MJ mm ha^-1 h^-1) from storms without runoffs. 

The numbers between parentheses indicate the of the month when the event took place.

According to the magnitude of the production of sediments in some storms that did not fulfill the criteria for the inclusion of the EI₃₀ in the R factor, it is necessary to consider the storms where there were runoffs and sediment production, even when they do not comply with the requirements fixed, as it happened in the event of September 24^th when there were two rain periods (Table 2), and the sediment production in that storm was significant (14.6 Mg) compared to the production of the other runoff events (Table 3). Similarly, in the runoff event on May 23^rd, with considerable sediment production (21.8 Mg) (Table 3), the erosivity index was not considered in one of the two rain periods (Table 2). In addition to this event, there were no storms immediately before as possible influence, since the nearest one took place on May 14^th with a precipitation of 1 mm in 2.7 h.

With values of 10.1, 0.062 and 1.0, respectively, for the LS, C and P factors, the calibrated value of K (0.1067) for the MUSLE model, and the measured value of R (2147) used, the USLE model resulted in an erosion rate of 143.0 Mg ha^-1 year^-1. The sediment production at the exit of a basin ranges from 33 to 100 % of soil losses, from sheet and rill erosion, represented by the USLE model, in the US (^{NEH, 1983}; ^{Calhoun and Fletcher III, 1999}). For the analysis of the result, if it is considered that the sediment production in the study basin was due only from sheet and rill erosion, and it is assumed that 37 % of its value reached the discharge point of the watershed obtained from the national engineering manual from the US (^{NEH, 1983}) in function of the area of the micro-basin, a sediment production of 52.9 Mg ha^-1 year^-1 would be expected, which is higher than the one measured in the study basin (1.2 Mg ha^-1 year^-1) and even the one observed in basins with high erosion rates such as Las Huertitas, which discharges 15 Mg ha^-1 year^-1 to the Cointzio Dam in the state of Michoacán (^{Duvert et al., 2010}).

The USLE model estimates the average annual soil loss of many years in the terrain, and, an erosion of 3.25 Mg ha^-1 year^-1 would need to take place in order to satisfy the soil delivery rate to the watershed’s exit of 37 % (1.2 Mg ha^-1 year^-1). Pero this amount represents large variability compared to the 143.0 Mg ha^-1 year^-1 estimated for the year 2013, and implies that for some years there would need to be much higher erosion than this value. These results show that the delivery rate considered does not represent the characteristics of the study basin, and that the cultivation factor considered was not the adequate one.

The results from the USLE model suggest the benefit of adapting the methodologies used for: 1) obtaining the R factor, and the shape of its representation through the Q_sxq_p product in the MUSLE model, since the value of the power of this term (0.7031) was higher than the 0.56 proposed by ^{Williams (1975}), as is shown in the potential functional relation in Figure 4; 2) calculating the factor of the crop C and the soil K; and, 3) calculating the delivery rate to the basin’s exit where the contributions of sediments from erosion in gullies, streambed, streambanks, roadbanks and mass landslides would need to be considered (^{NEH, 1983}; ^{Calhoun and Fletcher III, 1999}). This adaptation would allow a compatibility of results between the USLE and MUSLE models, when applied to the study basin.