Seasonal and spatial patterns of erosivity in a tropical watershed of the Colombian Andes

Journal of Hydrology 314 (2005) 177–191 www.elsevier.com/locate/jhydrol

Seasonal and spatial patterns of erosivity in a tropical watershed of the Colombian Andes ´ lvaro Jaramillob,2 Natalia Hoyosa,1, Peter R. Waylena,*, A a

Department of Geography, University of Florida, P.O. Box 117315, Gainesville, FL 32611-7315, USA b Centro Nacional de Investigaciones de Cafe´, CENICAFE, Chinchina´, Colombia Received 3 May 2004; revised 16 March 2005; accepted 18 March 2005

Abstract The Dosquebradas Basin, in the central coffee growing region of Colombia, covers an area of 58 km2 between 1350 and 2150 m of elevation, with an annual precipitation of 2600–3200 mm. Seasonal erosivity (EI30), as defined by the Revised Universal Soil Loss Equation (RUSLE), was calculated for 11 years of record (1987–1997) from six pluviographic stations located within 21 km of the basin. Regression models for each station indicated that storm rainfall explained 61–70% of the variation in storm erosivity. Individual storms represented as much as 25% of the annual EI30 (10,409–15,975 MJ mm haK1 hK1 yrK1). At the seasonal scale, the explained variation increased to 75–86%. There was a significant difference between wet and dry seasons, with higher values and larger increases in erosivity per unit increase in rainfall during the wet seasons. Two pooled regression models, one for the wet and one for the dry seasons, were created and used to estimate seasonal erosivity for 10 stations with pluviometric data. Interpolation surfaces were created from seasonal values using the local polynomial algorithm. Spatial patterns of erosivity were related to (a) the regional elevation gradient, particularly important during the dry seasons, and (b) local topographic effects, particularly during the wet seasons. Our findings underscore the importance of using seasonal erosivity values and local rainfall intensity records in tropical mountainous regions characterized by marked rainfall seasonality and complex topography. q 2005 Elsevier B.V. All rights reserved. Keywords: Rainfall erosivity; RUSLE; Soil erosion; Andes; GIS

1. Introduction Soil erosion is defined as the detachment of soil particles from the surface by some erosive agent,

* Corresponding author. Fax: C1 352 392 8855. E-mail addresses: [email protected] (N. Hoyos), prwaylen@ geog.ufl.edu (P.R. Waylen), [email protected] (A Jaramillo). 1 Fax: C1 352 392 8855. 2 Fax: C57 6 850 47 23.

0022-1694/$ - see front matter q 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2005.03.014

and subsequent transportation to another location (Flanagan, 2002). The energy to carry out such work comes from sources like (a) physical energy including wind and water, (b) gravity through mass movements, (c) chemical reactions through solution and weathering, and (d) anthropogenic perturbations, for example through tillage (Lal, 2001). Specifically, soil water erosion encompasses the detachment and transportation of soil particles by raindrops and flowing water (Morgan, 1995). This is a complex process resulting from the

178

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

interaction of the soil itself, climate, relief, surface cover and land use practices. The influence of climate is referred to as erosivity, and is generally expressed in terms of rainfall amount and intensity, both of which determine its potential to cause erosion. Several erosivity indices have been developed to quantify the effect of raindrop impact and reflect the amount and rate of runoff likely to be associated with that rain. These include the modified Fournier index for Morocco (Arnoldus, 1977), the EI30 index for the US Midwest (Brown and Foster, 1987; Wischmeier and Smith, 1978; Wischmeier, 1959), the KEO25 index for southern Africa (Hudson, 1971) and the AIm index for Nigeria (Lal, 1976). From these, the EI30 index developed as part of the Universal Soil Loss Equation (USLE; Wischmeier and Smith, 1978) and later Revised USLE (RUSLE; Renard et al., 1997) has been the most commonly used. However, its calculation is elaborate and requires rainfall intensity records over a length of time frequently unavailable (O20 years recommended; Renard et al., 1997). To facilitate the calculation of this index, models to estimate it from other types of precipitation data (e.g. monthly or annual totals) have been developed (i.e. Angima et al., 2003; Mati et al., 2000; Yu and Rosewell, 1996a; Renard and Freimund, 1994). However, variability in model parameters indicates that the relationship between rainfall amount and erosivity is site specific, or in the best case, region specific (Yu, 1998; Lo et al., 1985). Our objective was to develop such relationships for a basin of the Colombian Andes, as a first step in the application of a soil erosion model using the RUSLE methodology. This basin is representative of other regions in tropical mountainous environments in terms of high spatial variability of rainfall, complex and steep topography, and limited availability of pluviographic data (Poveda et al., in press; El-Swaify, 1997). The specific objectives of this study were to (a) generate basic erosivity information from pluviographic rainfall stations located near the basin of interest, (b) develop a regional model to predict erosivity from more readily available types of rainfall data, (c) study temporal and spatial variability of erosivity, and (d) summarize these results as isoerodent maps.

2. Data and methods 2.1. Site description The Dosquebradas Basin is located in the coffeegrowing region of the Central Cordillera of Colombia. The basin covers an area of 58 km2 and has a semicircular shape with an open southern boundary formed by the Otu´n River scarp (Fig. 1). Elevations range from 1350 to 2150 m, with a relief characterized by a valley floor surrounded by slopes with moderate (up to 25%) to steep gradients (25% to more than 75%) to the west, north and east (Instituto Geogra´fico Agustı´n Codazzi [IGAC], 1988). Annual precipitation ranges between 2600 mm in the south and 3200 mm at the northern watershed divide, and has a bimodal distribution, with two wet (MAM, SON) and two dry (DJF, JJA) seasons (Guzma´n and Jaramillo, 1989). The predominant soil mapping unit in the basin is Chinchina´, which consists mostly of Melanudands (80%) derived from volcanic ash deposits (IGAC, 1988). About 16% (908 ha) of the basin has been urbanized, most of it concentrated on the valley floor. Major land uses in the rural area include coffee agriculture and pasture (62 and 18% of the rural area, respectively) (Corporacio´n Auto´noma Regional de Risaralda [CARDER], 1997). 2.2. Erosivity calculation for pluviographic stations We chose the EI30 index to calculate erosivity as it has been widely used, therefore providing the best opportunity for comparison with other locations. In addition, a study on soil erosion from runoff plots in southwest Colombia showed that this index had a better correlation with soil loss than others such as KE and single rainfall intensity values (Ruppenthal et al., 1996). Individual storm EI30 values were computed following the Revised Universal Soil Loss Equation RUSLE methodology (Renard et al., 1997), using rainfall intensity data from six gauging stations located within 21 km of the basin boundary (Fig. 2, Table 1). These data were provided by the National Coffee Research Center CENICAFE. A period of 11 years (1987–1997) was selected for analysis because it represented a reasonable compromise between the number of available stations and the amount of data to be analyzed. Within each station, storm EI30 values

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

179

Fig. 1. Location of the Dosquebradas basin, western flank of the central andean cordillera of colombia (modified from digital chart of the world Digital Chart of the World (DCW), 1992; Centro Internacional de Agricultura Tropical (CIAT), 1998 and Programa de las naciones unidas para el medio ambiente [PNUMA], 1998; CARDER, 1997).

180

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

2.3. Missing data All gauging stations had incomplete records for a period of 1–2 years between 1995 and 1997. During this period, total daily rainfall was available but intensity data were limited to those storms producing the highest rainfall and maximum intensity for that particular month. Records were completed by:

Legend

Fig. 2. Location of rainfall gauging stations relative to the study area, grouped by type of data available for this study. Stations operated by CENICAFE and IDEAM. Transverse mercator projection with origin at 4835 0 56 00 N and 77804 0 51 00 W, international spheroid 1924, bogota´ observatory datum.

were added on a seasonal basis and differences between seasons were determined through analysis of variance. Fig. 3 summarizes the methods followed to develop the regional erosivity model from pluviographic (intensity) and pluviometric (totals) data.

† Adding days when daily rainfall was equal to or greater than 12.7 mm. This threshold was chosen deliberately since it is part of the criteria used to define an erosive event according to the RUSLE methodology (Renard et al., 1997). In order to check that daily rainfall was a reasonable substitute for storm rainfall, correlation analyses were run for daily and storm rainfall in each station. † Developing a relationship between storm rainfall (independent) and storm EI30 (dependent) through regression analysis at each station separately. Only storms with total rainfall greater or equal to 12.7 mm were included. Both variables were log transformed and a linear regression analysis was performed separately for each station. The derived relationship was then used to complete missing EI30 values from observed daily rainfall records.

Table 1 Characteristics of available rainfall stations Station

Lat. N

Long. W

Elev. (m)

Period of record

Average annual rainfall (mm)

Years with missing data

Period used in this study

Cenicafe´a Catalinaa Naranjala Pta. Tratamientoa Jazmı´na Cedrala Bohemiab Combiab Ca´mbulosb Apto. Matecan˜ab Rosab Bosqueb Recuerdob Sta. Rosab Floridab Veracruzb

5800 0 4845 0 4858 0 4849 0 4855 0 4842 0 4853 0 4851 0 4849 0 4849 0 4850 0 4850 0 4857 0 4853 0 4854 0 4852 0

75836 0 75845 0 75842 0 75841 0 75838 0 75832 0 75855 0 75847 0 75850 0 75844 0 75842 0 75841 0 75845 0 75838 0 75840 0 75838 0

1310 1350 1400 1450 1600 2120 1020 1173 1240 1342 1440 1450 1590 1675 1660 1684

1941–2002 1987–1999 1950–2002 1971–1999 1961–1999 1977–1999 1963–2000 1981–1999 1969–1998 1947–2000 1978–1995 1978–1999 1970–2000 1978–1995 1978–1995 1977–2000

2524 2171 2715 2682 2559 2665 1871 2114 1802 2198 2654 2901 2780 2683 3390 2647

0 0 0 5 6 3 7 2 7 5 7 1 7 4 5 6

1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1997 1987–1995 1987–1997 1987–1997 1987–1995 1987–1995 1987–1997

a b

Pluviographic data used to calculate the EI30 index. Pluviometric data (monthly totals).

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

storms ≥ 12.7mm (n=577)

Storm EI30 values

+ storms ≥ 12.7mm daily rainfall for missing days Complete dataset (n=775)

Log storm EI30

Jazmín original records 1987-1997 (n=666 storms)

181

Log storm rain

Log seasonal EI30

+ seasonal models for other stations Comparisons among stations seasons Log seasonal rain

seasonal EI30 map seasonal rain pluviometric stations

Log dryseason EI30

Log wetseason EI30

Regional model

Log wetseason rain

Log dryseason rain

Fig. 3. Flowchart describing the development of a regional erosivity model and a seasonal erosivity map from pluviographic and pluviometric data. The example of Jazmı´n station is presented.

Several indices were used to assess the performance of regression models. These included the square root of the mean square error (RMSE), adjusted coefficient of determination (adj. R2), coefficient of efficiency (E) and modified coefficient of efficiency (E1). The coefficient of efficiency is defined as (Nash and Sutcliffe, 1970) PN ðOi K Pi Þ2 E Z 1:0 K PiZ1 N  2 iZ1 ðOi K OÞ where O refers to observed data, O is the observed mean and P is the model predicted data. A value of zero indicates that the observed mean is as good a predictor as the model itself, while negative values indicate that the observed mean is a better predictor than the model (Wilcox et al., 1990). Since values are squared, this coefficient is very sensitive to extreme values. This is overcome by the modified coefficient

of efficiency, defined as (Legates and McCabe, 1999): PN jOi K Pi j E1 Z 1:0 K PiZ1 N  iZ1 jOi K Oj 2.4. Regional erosivity model The seasonal erosivity models from the pluviographic stations were used to build a regional model, so that rainfall data from pluviometric stations could be used to estimate seasonal erosivity. To determine if data from different stations should be pooled together, or if it was necessary to separate wet from dry seasons, the seasonal regression models were tested for (a) differences between wet and dry seasons within each station, (b) differences among stations, (c) differences among wet and dry seasons with all stations data pooled together. Comparisons between

182

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

models were done following the method explained by Chatterjee and Price (1977). Model performance was evaluated through the indices already mentioned. 2.5. Erosivity prediction for pluviometric stations Monthly rainfall data were available for 10 gauging stations located within 27 km from the center of the basin, operated by CENICAFE and the National Environmental Institute IDEAM (Fig. 2, Table 1). These data were grouped into seasons and used with the regional regression model to obtain seasonal EI30 values for each pluviometric station. Then, differences among seasonal EI30 values for each station were determined through analysis of variance. 2.6. Spatial patterns of seasonal erosivity Exploratory spatial data analysis was carried out to assess spatial variability while the relationship between erosivity and elevation was explored through scatter plots. The spatial trend of seasonal EI30 was modeled using ArcGIS Geostatistical Analyst (Johnston et al., 2001). Two interpolation methods were tested, inverse distance weighted and local polynomial. Universal polynomial was not used because it assumes very gradual change (Johnston et al., 2001). The final surface was chosen as that one combining the lowest prediction error and being physically meaningful.

the criteria used to calculate the EI30 index at each station. The regression models of storm rain and EI30 used to complete the data sets for each pluviographic station had very similar adjusted R2 and E values, while lower E1 values (Table 2). Two stations failed to comply with the assumption of equal variance of residuals (Naranjal and Pta. Tratamiento). However, they were still used because the failure was explained by the way in which storm EI30 values were calculated. That is, for a given rainfall amount the maximum EI30 would be obtained when all the rain falls within 30 min. Higher values could not be obtained regardless of all the rain falling within 20, 10 or 5 min. As a result, values of EI30 were limited in the upward direction as shown by the dashed line in Fig. 4, affecting the assumption of normality in the residuals. Correlation analyses between daily totals and storm rainfall yielded coefficients (Pearson’s r) ranging from 0.91 to 0.95, indicating that daily rainfall was a reasonable proxy for storm rainfall during years with missing storm data. In this way, data sets for each station were completed and seasonal EI30 values were obtained by adding storm EI30 values for each season. In general, seasonal EI30 values had a pattern similar to rainfall, with lower values in the dry seasons and higher values in the wet seasons (Table 3). 3.2. Regional erosivity model

3. Results 3.1. Erosivity calculation for pluviographic stations The intensity records for the 1987–1997 period yielded an average of 603 storms, which met

Seasonal regression models showed adjusted R2 and E values ranging from 0.75–0.86, and E1 ranging from 0.49 to 0.66 (Table 4). Comparisons of wet and dry season models within each station showed no difference among seasons except for Jazmı´n (Table 4). On the other hand, comparisons of regression models

Table 2 Regression model parameters and measurements of model performance for storm EI30 (log MJ mm haK1 hK1) versus storm rain (log mm) Station

n

Intercept

Slope

RMSE (log MJ mm haK1 hK1)

Adj. R2

E

E1

Catalina Jazmı´n Cenicafe´ Naranjal Pta. Tratamiento Cedral

441 577 573 596 471 523

K0.269 K0.375 K0.439 K0.407 K0.403 K0.292

1.693 1.769 1.819 1.795 1.798 1.664

0.26 0.26 0.25 0.25 0.26 0.28

0.62 0.68 0.70 0.69 0.69 0.61

0.62 0.69 0.70 0.69 0.69 0.61

0.38 0.42 0.47 0.44 0.43 0.37

4.0

4.0

3.5

3.5

3.0 2.5 2.0 1.5

y = -0.375+1.769x n = 577 Adj. R2 = 0.68

1.0

0.5 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 Log storm rain (a)

Log season EI30

Log storm E130

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

183

3.0 2.5 y = 0.063 + 1.234x n = 43 Adj. R2 = 0.75

2.0 1.5 2.0

(b)

2.2

2.4 2.6 2.8 3.0 Log season rain

3.2

Fig. 4. Relationship between rain (mm) and EI30 (MJ mm haK1 hK1) at the (a) storm and (b) seasonal scale. The dashed line represents the maximum possible EI30 value for each rainfall amount. The continuous line is the regression line calculated by the least squares method. The example of Jazmı´n station is presented.

among stations indicated that Cedral was always significantly different from the others. These comparisons were repeated after excluding one of Cedral’s data points that clearly stood out from the general trend. This point corresponded to the dry season of June–August 1991, with a extremely low EI30 value (51.5 MJ mm ha K1 hK1 ). However, excluding this point made no difference as the regression model for Cedral still came out as being different from the rest. Since this station clearly had a different relationship between seasonal rainfall and EI30, it was excluded from the regional model while data from all the others were pooled and tested again for differences among seasons. In this case, the relationship between seasonal rain and seasonal EI30 showed a significant difference between wet and dry seasons. As a result, two regional regression models were produced, one for the dry seasons and another one for the wet seasons (Fig. 5, Table 4).

Within each of these models, the explained variability (R2) could be inflated by the fact that measurements across stations were not independent from each other (spatial autocorrelation). However, the indices of performance for the regional models were not significantly higher than those for each individual station (Table 4). 3.3. Erosivity prediction for pluviometric stations Seasonal erosivity values at each pluviometric station had a pattern similar to rainfall amount (not shown), falling into one of the following categories: † Bimodal pattern, with significant differences among the wet and dry seasons (Apto. Matecan˜a, Bosque, Combia and Rosa). † Bimodal pattern with some seasons not significantly different (Bohemia, Veracruz, Sta. Rosa, Florida and Ca´mbulos).

Table 3 Seasonal EI30 values for pluviographic stations Station

n

Seasonal EI30 (MJ mm haK1 hK1 seasonK1)a DJF

Catalina Jazmı´n Cenicafe´ Naranjal Pta. Tratamiento Cedral

43 43 43 43 43 42

1947.2 3004.5 3318.9 3006.1 2895.5 2920.0

MAM a ab a a a

3220.1 4528.4 4952.2 5182.6 4907.7 3384.4

JJA b a b b a

1853.7 2242.2 3131.5 3409.5 2565.9 818.6

SON a b ab a b

3345.2 4310.1 4669.6 4350.1 4673.2 4298.3

b a ab b a

Values followed by different letters within the same line are significantly different from each other (aZ0.05, Tukey–Kramer test). a Seasons: DJF, December–February; MAM, March–May; JJA, June–August; SON, September–November.

184

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

Table 4 Regression model parameters and adjusted R2 for seasonal EI30 (log MJ mm haK1 hK1 seasonK1) versus seasonal rain (log mm) Station Catalina Jazmı´nc Cenicafe´ Naranjal Pta. Tratamiento Cedral Wet season Dry season

c

Slope

RMSE (log MJ mm haK1 hK1 seasonK1)

Adj. R2

E

E1

Different from

43 43 43

K0.169 0.063 K0.883

1.314 1.234 1.591

0.10 0.10 0.11

0.81 0.75 0.86

0.81 0.75 0.86

0.62 0.49 0.57

43 43 42 110 105

K0.250 K0.240 K1.338 K0.967 K0.365

1.354 1.353 1.690 1.601 1.399

0.08 0.11 0.17 0.08 0.12

0.81 0.81 0.86 0.79 0.76

0.81 0.81 0.85 0.79 0.76

0.62 0.57 0.66 0.56 0.54

Cenicafea Cedralb Cenicafeb Cedralb Catalinaa Jazmı´nb Cedralb Cedralb Cedralb Allb

Models different at aZ0.05. Models different at aZ0.01. Wet and dry models different within the station at aZ0.05.

† No significant (Recuerdo).

difference

among

seasons

3.4. Spatial patterns of seasonal erosivity

Log season EI30

Exploratory data analysis indicated increased erosivity from west to east, and from south to north (Fig. 6). In both seasons, the west–east trend was of polynomial character, with a slight decrease towards the easternmost and highest station (Cedral). By comparison, the south–north trend had a linear increase, with a steeper slope during the dry season, when highest erosivity values were centered around the northeast stations (Naranjal, Cenicafe´, Florida and Sta. Rosa), and the lowest values to the southwest (Bohemia, Ca´mbulos) and southeast (Cedral). During

4.0

4.0

3.5

3.5

3.0 2.5 y = -0.967+1.601x n = 110 Adj. R2 = 0.79

2.0 1.5 2.0

(a)

the wet season, the highest values were at the center (Florida and Bosque) and lowest to the southwest (Bohemia and Ca´mbulos). These trends suggested a relationship between elevation and erosivity, with the latter increasing to about 1750 m and then decreasing to the lower values at Cedral (2120 m; Fig. 6). The interpolation surfaces reflected the noted trends. The local polynomial interpolation method was selected because it had lower mean prediction errors and comparable root mean square prediction errors relative to the inverse distance weighted method (Fig. 7, Table 5). In addition, surfaces calculated with the inverse distance weighted method had anomalous contours. All 16 stations were included to calculate the interpolation surfaces.

2.2

2.4 2.6 2.8 3.0 Log season rain

Log season EI30

a b

Intercept

n

3.0 2.5

1.5 2.0

3.2

(b)

y = -0.365+1.399x n = 105 Adj. R2 = 0.76

2.0 2.2

2.4 2.6 2.8 Log season rain

3.0

3.2

Fig. 5. Regional erosivity models for the (a) wet and (b) dry seasons of seasonal rain (mm) and EI30 (MJ mm haK1 hK1 seasonK1), with all pluviographic stations data included except cedral.

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

185

Z

(a)

2200

N E X

Elevation (m)

Cedral 1900 1600

1300 1000 1500

3500

5500

6500

EI30 (MJ mm ha-1 h-1 season-1) Z

(b)

2200

Elevation (m)

Cedral

Z

N

E

1900

1600

1300

1000 1500

3500

5500

6500

EI30 (MJ mm ha-1 h-1 season-1)

Fig. 6. Trend analysis plots and variation of seasonal EI30 (MJ mm haK1 hK1 seasonK1) with elevation for the (a) wet and (b) dry seasons. In trend plots, each dot represents a rainfall station with elevation above the horizontal plane being proportional to the seasonal erosivity value.

4. Discussion 4.1. Seasonal patterns The regression models developed in this study between erosivity and rain amount revealed that, at the storm level, rainfall explained between 61 and 70% of the variability in EI30. Model performance was very similar for all stations except Catalina and Cedral, where it was less satisfactory. When both variables were grouped at the seasonal level, the explained variability increased to 75–86%. Yu and Rosewell (1996b) also noted a decrease in explained variability when using a daily erosivity model for New South Wales, Australia, to predict single storm EI30 (i.e. underestimation of EI30 for severe storms). Model performance at the seasonal level was similar for all stations but Jazmı´n. This may be related to the significant differences among regression models for the wet and dry seasons found only at this station. On the other hand, differences in model parameters

between Cedral and all other stations may be pointing to the need of using longer data series to capture the effect of regional climatic anomalies at the other stations. For instance, during El Nin˜o phase of ENSO (El Nin˜o/Southern Oscillation ENSO) decreased annual rainfall and soil moisture have been reported for this region, particularly severe during the dry season (Poveda et al., 2002; 2000; Guzma´n and Baldio´n, 1997a). Erosivity showed differences among wet and dry seasons, not only in terms of its magnitude, but also in its relationship to rainfall. The seasonal regression models had a significantly higher slope for the wet seasons than for the dry seasons. Therefore, for every unit increase in rainfall, there was a larger increase in erosivity during the wet seasons. We believe this was related to a higher frequency of large and intense storms during the wet seasons. An analysis of the 20 highest I30 (maximum 30-min intensity) values for each of the pluviographic stations revealed that 71% of them occurred during the wet seasons. Similar results

186

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

(a)

+

+

+ +

+

(b)

+

+

+ +

++

+ Fig. 7. Erosivity surface (MJ mm haK1 hK1 seasonK1) and Voronoi map of prediction error (%) for the (a) wet and (b) dry seasons. Transverse mercator projection with origin at 4835 0 56 00 N and 77804 0 51 00 W, international spheroid 1924, bogota´ observatory datum.

have been found by other local and regional studies (Poveda et al., 2002; Sua´rez de Castro and Rodrı´guez, 1962). In addition, the number of rainfall events for any given rainfall amount was higher during the wet season, as expected. As a result, and considering that the rainfall generating mechanism for both seasons is the same (i.e. convection), the likelihood of having a large storm resulting in a high EI30 value, would be greater during the wet season. Differences in erosivity between seasons confirmed the importance of considering seasonality for erosivity calculations, particularly in this region where maximum erosivity values coincide with the two coffee harvesting seasons (minor harvest in April through June, and major harvest in October through December; Chalarca´, 1998). Agricultural practices in coffee plantations, including fertilization and harvesting, have been identified as critical for soil erosion (Ataroff and Monasterio, 1997).

4.2. Spatial patterns The spatial patterns of erosivity showed large variability within relatively short distances. For instance, during the wet season, observed erosivity Table 5 Comparison of different interpolation methods Method

Power

Prediction errors (MJ mm haK1 hK1 seasonK1) Mean

Wet season Inverse distance weighted Local polynomial Dry season Inverse distance weighted Local polynomial

RMSE

2

274.7

809.8

1

K85.7

705.1

2

149.1

475.0

1

41.4

424.2

Local polynomial was selected for both seasons.

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

at Pta. Tratamiento was 46% higher than at Catalina (4790 vs. 3283 MJ mm haK1 hK1 seasonK1) although these stations were only 10 km apart. Several factors explain this variability, such as elevation, local topography, and environmental conditions controlling the soil-atmosphere interactions (Poveda et al., in press; Poveda, 2004; Pe´rez, 1983). Seasonal rainfall in this area increases with elevation up to about 1750 m and then decreases. Regional studies also show increasing rainfall with elevation up to a level that varies between 1200 and 2200 (Universidad Nacional de Colombia, 1997; Guzma´n and Jaramillo, 1989; Pe´rez, 1983). Accordingly, elevation effects on rainfall amount partially explain our increasing erosivity values from the southwest to the northeast. In addition, rainfall intensity also seemed affected by elevation, with consistently lower I30 values at the highest station (Cedral, 2120 m). For example, a comparison of the upper quartile values among stations showed that Cedral’s average (45.3 mm hK1) was significantly lower than all other stations, except Catalina’s (47.7 mm hK1). Furthermore, the maximum I30 value at Cedral was the lowest (85 mm hK1 in 30 min) of all stations. Studies in other tropical mountainous regions show contradictory results. In southwest Mexico, Millward and Mersey (1999) found that erosivity and precipitation increased with elevation, while in Honduras Mikhailova et al. (1997) found that erosivity decreased from 200 to 1200 m of elevation. The latter study examined datasets from Costa Rica, Sri Lanka and southeastern US and found that in all cases there was an inverse relationship between erosivity and elevation, which was attributed to fewer large drops formed by accretion and coalescense at higher elevations (Mikhailova et al., 1997). The erosivity values from our study were similar to those found at other locations in southwest (Ruppenthal et al., 1996) and central Colombia (Rivera and Go´mez, 1991; Table 6). Comparable studies in other tropical regions show a wide range of values, while in temperate regions numbers are generally lower and found to decrease with increasing latitude (Table 6). The difference between tropical and temperate regions is related to differences in the amount of precipitation, intensity and kinetic energy of rain, the latter being a function of raindrop size and terminal velocity (Hudson, 1995; Lal, 1990;

187

Table 6 Comparison of annual EI30 values found in this and other studies Location

Years of record

Maximum I30 (mm hK1)

Average annual EI30a (MJ mm haK1 hK1 yK1)

This study

11

85–112

10409–15975 (6)

15

67–85

2

75–84

20

1

5065–13,083 (9) 6345–10,060 (2) 2237–12,881

15 O22

N/A N/A

12

N/A

18

N/A

N/A

N/A

N/A

Other tropical sites Central Colombia Southwest Colombia Northeast Brazil Honduras Southwest Mexico Central Kenya Malaysia Tropical Australia Hawaii Temperate sites Southern Australia Southern Portugal Southeastern US Northeastern US

Referenceb

(a) (b) (c) (d) (e)

N/A

3385–7297 (8) 4300–9600 (30) 3308–13,566 (1) 13,600–21,600 (2) 1079–33,481 (41) 1700–23,828

O30

N/A

250–500 (99)

(j)

28

N/A

697–3742 (32)

(k)

O22

N/A

4680–10212

(i)

O22

N/A

1020–3400

(i)

(f) (g) (h) (i)

a

In parenthesis, number of stations used. (a) Rivera and Go´mez (1991); (b) Ruppenthal et al. (1996); (c) Dias and Silva (2003); (d) Mikhailova et al. (1997); (e) Millward and Mersey (1999); (f) Angima et al. (2003); (g) Yu et al. (2001); (h) Yu (1998); (i) Renard et al. (1997); (j) Yu and Rosewell (1996a); (k) Santos and Azevedo (2001. b

El-Swaify and Dangler, 1982). The main rainfall generating mechanism in most tropical regions is convection. As a result, the tropics receive more rain at higher intensities than the temperate regions, dominated by midlatitude cyclones (El-Swaify and Dangler, 1982; Sheng, 1982). Consequently, larger EI30 values are observed in the tropics. Hudson (1971) proposed a threshold intensity between erosive and non-erosive rains of 25 mm hK1, and indicated

188

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

that in temperate regions 5% of the rains exceeded this value, while the figure was 40% in tropical regions. Although Hudson did not specify whether his values were derived from 5 to 30 min, etc. readings, an inspection of the 30-min intensity data from our study is instructive. The percentage of the rains with an intensity greater than 25 mm hK1 ranged from 38% (Cedral) to 52% (Cenicafe´). This is in line with Hudson’s estimate and with the 40–49% figure from southwest Colombia (Ruppenthal et al., 1996). An additional issue that needs to be addressed when comparing erosivity values from tropical and temperate regions is the uncertainty associated with the energy–intensity relationship used in the EI30 calculations. The relationship between storm kinetic energy and intensity used in the EI30 index is derived from drop size measurement studies made in Mississippi, USA, and partially from other studies in the United States (Washington DC, Miami), Zimbabwe and Australia (Van Dijk et al., 2002; Brown and Foster, 1987). Van Dijk et al. (2002) evaluated other studies on this relationship and proposed a general energy– intensity equation based on data from tropical, temperate and desert conditions. However, differences between the values obtained with their general equation and the one from the EI30 index were modest (wK8%), particularly when compared to natural variations between storms, and the uncertainty associated with the estimation of other RUSLE model factors (Van Dijk et al., 2002). Still, the extent to which either equation (Van Dijk’s or EI30’s) applies to tropical mountainous conditions is still questionable (Van Dijk et al., 2002). More research in these environments, particularly on the effect of altitude on drop size distribution and terminal velocity, is needed to elucidate the nature of the kinetic energy–intensity relationship (Van Dijk et al., 2002; Ruppenthal et al., 1996). Finally, it is important to assess the contribution of large, infrequent storms to erosivity, as they have been found to cause most of the erosion. Runoff plots in the central coffee region of Colombia (Sua´rez de Castro and Rodrı´guez, 1962) between 1949 and 1956, showed that 57% of all rainfall events generated runoff, but a mere 10% of them accounted for 89% of the soil loss. In southwest Colombia, Ruppenthal et al. (1996) reported that 61% of the annual soil loss from a runoff plot was caused by three large rainstorms,

Table 7 Contribution of the largest storm EI30 in record (1987–1997) to the annual EI30 value Station

Maximum storm EI30 (MJ mm haK1 hK1) Year

Value

Catalina Jazmı´n Cenicafe´ Naranjal Pta. Tratamiento Cedral

1992 1988 1987 1991 1996 1994

Annual EI30 (MJ mm haK1 hK1 yK1)

Contribution of storm to annual EI30 (%)

2191 2440 2683 1918 3826

8789 18991 19514 14286 22095

25 13 14 13 17

1692

12631

13

which in turn accounted for 30% of the annual EI30. Although soil loss was not measured in this study, the contribution of single events to the annual EI30 is illustrative of their influence. Evaluation of the largest EI30 value for each station indicated that it accounted for 13 to 25% of the respective year’s EI30 (Table 7). 4.3. Interpolation surfaces and prediction errors The interpolation of erosivity values from the 16 rainfall stations resulted in two isoerodent surfaces, one for the wet seasons and one for the dry seasons (Fig. 7). Variability in the dry seasons was lowest in the NW–SE direction, while for the wet seasons this direction shifted to the NNW–SSE. The erosivity gradient was steeper in the wet season than in the dry season. The spatial distribution and magnitude of prediction errors varied between seasons (Fig. 7). No spatial autocorrelation of residuals was found in either surface (Moran’s I). In the wet season, the model greatly underpredicted values at Bohemia (K43%) and Florida (K29%). For the dry season, the model overpredicted at Cedral (C67%) and Jazmı´n (C24%). Bohemia and Cedral are located at the western and eastern boundaries of the interpolation surface, and as such, larger errors from the interpolation procedure can be expected. They may also be influenced by different local processes, as Bohemia, with the lowest elevation, is located in the Cauca Valley, while Cedral, with the highest elevation, is located in the upper Otu´n River canyon. These stations also differ in their daily rainfall pattern, with

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

Bohemia’s rainfall occurring mostly at night (80%), while Cedral’s occurring mostly during the day (70%) (Guzma´n and Baldio´n, 1997b). Stations at elevations in between have approximately equal proportions of rain falling during the day and night. The overestimation at Cedral could also result from the unusually low measured EI30 values, largely influenced by the two dry seasons (June–August) of 1991 and 1993. In 1991 no events (according to the RUSLE criteria) were recorded for June and August, resulting in a very low EI30 (51.5 MJ mm haK1 hK1 seasonK1). Similarly, in 1993, no events were recorded in June, and only one event was recorded for July and one for August. The EI30 value for this season was also unusually low (126 MJ mm haK1 hK1 seasonK1). Prediction errors at these two stations were not likely to have a large impact on the basin as they were relatively far. This was not the case for Jazmı´n and Florida. Overestimation at Jazmı´n during the dry season may have resulted from the fact that all surrounding stations had higher values. By comparison, underestimation of the wet seasonal values at Florida was related to this station having unusually high EI30 relative to the surrounding ones. To assess the effect of removing this station from the prediction, an interpolation surface omitting Florida was created and subtracted from the original one. The results showed a mean difference in seasonal EI30 of 123 MJ mm haK1 hK1, with the largest difference at the northern boundary of the basin (300–400 MJ mm haK1 hK1). However, we decided to leave this station as part of the analysis because similar annual rainfall values (above 3000 mm) have been recorded at other stations in the region.

5. Conclusions Several findings from this research are relevant to the study of soil erosion potential in the central coffee region of Colombia as well as other areas in tropical mountains. First, the significant differences in erosivity between wet and dry seasons should discourage the use of annual averages. Second, the spatial variability within each season points to the necessity of calibrating existing relationships between rainfall totals and erosivity with local rainfall intensity data. Thirdly, the use of longer time series may reveal

189

unusually low or high values that can affect the regression models developed. Our calculated annual average EI30 values fall within the range observed in other tropical regions, while they are higher than the ones referenced for temperate regions, confirming that in general, higher amounts of rainfall and frequency of intense storms occur in the tropics. Finally, the spatial patterns of erosivity in this region show a general trend that seems to be associated with elevation, superimposed on a finer scale local variability. The latter one emphasizes the importance of local processes on precipitation amount and intensity in this topographically complex area.

Acknowledgements We would like to thank CENICAFE for providing the pluviographic data. The College of Liberal Arts and Sciences from the University of Florida provided financial support to travel to Colombia and gather these data.

References Angima, S.D., Stott, D.E., O’Neill, M.K., Ong, C.K., Weesies, G.A., 2003. Soil erosion prediction using RUSLE for central Kenyan highland conditions. Agriculture, Ecosystems and Environment 97, 295–308. Arnoldus, H.M.J., 1977. Methodology used to determine the maximum potential average annual soil loss due to sheet and rill erosion in Morocco. FAO Soils Bulletin 34, 39–51. Ataroff, M., Monasterio, M., 1997. Soil erosion under different management of coffee plantations in the Venezuelan Andes. Soil Technology 11, 95–108. Brown, L.C., Foster, G.R., 1987. Storm erosivity using idealized intensity distributions. Transactions of the American Society of Agricultural Engineers 30, 379–386. Chalarca´, J.E., 1998. Vida y hechos del cafe´ en Colombia. Comu´n Presencia Editores, Bogota´. Chatterjee, S., Price, B., 1977. Regression Analysis by Example. Wiley, New York. Centro Internacional de Agricultura Tropical (CIAT), 1998.Anon., 1998. Programa de las Naciones Unidas para el Medio Ambiente (PNUMA), Atlas de indicadores ambientales y de sustentabilidad para Ame´rica Latina y el Caribe. Cali, Colombia 1998. Corporacio´n Auto´noma Regional de Risaralda (CARDER), 1997.Anon., 1997. Informacio´n digital Municipio de Dosquebradas, escala 1:25,000. Pereira 1997.

190

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191

Digital Chart of the World (DCW), 1992. Digital chart of the world layer interface for Colombia. Penn State University Libraries. http://www.maproom.psu.edu/dcw/. Accessed in August 2004. Dias, A.S., Silva, J.R.C., 2003. Rainfall erosivity in Fortaleza, Ceara´ state, Brazil. I. Distribution, Probability of Occurrence and Return Period–1st Approximation, Revista Brasileira Ciencias Solo 27 2003 p. 335–345. El-Swaify, S.A., 1997. Factors affecting soil erosion hazards and conservation need for tropical steeplands. Soil Technology 11, 3–16. El-Swaify, S.A., Dangler, E.W., 1982. Rainfall erosion in the tropics: a state of the art. In: Kussow, W., El-Swaify, S.A., Mannering, J. (Eds.), Soil Erosion and Conservation in the Tropics. American Society of Agronomy/Soil Science Society of America Special Publication 43, Madison, pp. 1–25. Flanagan, D.C., 2002. In: Lal, R. (Ed.), Erosion Encyclopedia of Soil Science. Marcel Dekker, New York, pp. 395–398. Guzma´n, O., Baldio´n, J.V., 1997a. El evento ca´lido del Pacı´fico en la zona cafetera Colombiana. Cenicafe´ 48, 141–155. Guzma´n, O., Baldio´n, J.V., 1997b. Regionalizacio´n clima´tica de una zona montan˜osa tropical empleando el ana´lisis multivariado. Cenicafe´ 48, 260–274. Guzma´n, O., Jaramillo, A., 1989. Estudio clima´tico de Risaralda y Quindı´o. Boletı´n Te´cnico Cenicafe´ 15, 1–64. Hudson, N., 1971. Soil Conservation. Cornell University Press, Ithaca. Hudson, N., 1995. Soil Conservation, third ed. Iowa State University Press, Ames, IA. Instituto Geogra´fico Agustı´n Codazzi (IGAC), 1988.Anon., 1988. Estudio general de suelos del Departamento de Risaralda. Tomos I y II. Bogota´ 1988. Johnston, K., Ver Hoef, J.M., Krivoruchko, K., Lucas, N., 2001. Using ArcGIS Geostatistical Analyst. ESRI Press, Redlands. Lal, R., 1976. Soil erosion on alfisols in Western Nigeria. III— Effects of rainfall characteristics. Geoderma 16, 389–401. Lal, R., 1990. Soil Erosion in the Tropics: Principles and Management. McGraw-Hill, New York. Lal, R., 2001. Soil degradation by erosion. Land Degradation and Development 12, 519–539. Legates, D.R., McCabe, G.J., 1999. Evaluating the use of goodnessof-fit measures in hydrologic and hydroclimatic model validation. Water Resources Research 35, 233–241. Lo, A., El-Swaify, S.A., Dangler, E.W., Shinshiro, L., 1985. In: ElSwaify, S.A., Moldenhauer, W.C., Lo, A. (Eds.), Effectiveness of EI30 as an erosivity index in Hawaii Soil Erosion and Conservation. Soil Conservation Society of America, Ankeny, pp. 384–392. Mati, B.M., Morgan, R.P.C., Gichuki, F.N., Quinton, J.N., Brewer, T.R., Liniger, H.P., 2000. Assessment of erosion hazard with the USLE and GIS: a case study of the Upper Ewaso Ng’iro North basin of Kenya. International Journal of Applied Earth Observation and Geoinformation 2, 78–86. Mikhailova, E.A., Bryant, R.B., Schwager, S.J., Smith, S.D., 1997. Predicting rainfall erosivity in Honduras. Soil Science Society of America Journal 61, 273–279.

Millward, A.A., Mersey, J.E., 1999. Adapting the RUSLE to model soil erosion potential in a mountainous tropical watershed. Catena 38, 109–129. Morgan, R.P.C., 1995. Soil Erosion and Conservation. Wiley, Essex. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models. I. A discussion of principles. Journal of Hydrology 10, 282–290. Pe´rez, A., 1983. In: Van der Hammen, T., Pe´rez, A., Pinto, P. (Eds.), Algunos aspectos del clima La cordillera central Colombiana, transecto Parque de los Nevados, vol. 1. Studies on tropical Andean ecosystems, Cramer, Vaduz, pp. 38–47. Poveda, G., 2004. La hidroclimatologı´a de Colombia: una sı´ntesis desde la escala inter-decadal hasta la escala diurna. Revista de la Academia Colombiana de Ciencias Exactas, Fı´sicas y Naturales 107, 201–222. Poveda, G., Jaramillo, A., Mantilla, R., 2000. Influencia del evento ca´lido del Pacı´fico en la humedad del suelo y el ´ındice normalizado de vegetacio´n en Colombia. Cenicafe´ 51, 263– 271. ´ lvarez, J.F., Arias, P., Poveda, G., Mesa, O., Toro, V., Agudelo, P., A Moreno, H., Salazar, L., Vieira, S., 2002. Diagno´stico del ciclo annual y efectos del ENSO sobre la intensidad ma´xima de lluvias de duracio´n entre 1 y 24 horas en los Andes de Colombia. Meteorologı´a Colombiana 5, 67–74. Poveda, G., Mesa, O., Salazar, L.F., Arias, P.A., Moreno, H.A., Vieira, S.C., Agudelo, P.A., Toro, V.G., Alvarez, J.F., in press. The diurnal cycle of precipitation in the tropical Andes of Colombia. Monthly Weather Review 2005, 133 (1) 228–240. Renard, K.G., Freimund, J.R., 1994. Using monthly precipitation data to estimate the R-factor in the revised USLE. Journal of Hydrology 157, 287–306. Renard, K.G., Foster, G.R., Weesies, G.A., McCool, D.K., Yoder, D.C., 1997. Predicting soil erosion by water: a guide to conservation planning with the revised universal soil loss equation (RUSLE). USDA Handbook 703, Washington, DC. Rivera, J.H., Go´mez, A., 1991. Erosividad de las lluvias de la zona cafetera central Colombiana (Caldas. Quindı´o y Risaralda). Cenicafe´ 42, 37–52. Ruppenthal, M., Leihner, D.E., Hilger, T.H., Castillo, J.A., 1996. Rainfall erosivity and erodibility of Inceptisols in the southwest Colombian Andes. Experimental Agriculture 32, 91–101. Santos, N., Azevedo, M., 2001. A new procedure to estimate the RUSLE EI30 index, based on monthly rainfall data and applied to the Algarve region, Portugal. Journal of Hydrology 250, 12–18. Sheng, T.C., 1982. In: Kussow, W., El-Swaify, S.A., Mannering, J. (Eds.), Erosion problems associated with cultivation in humid tropical hilly regions Soil Erosion and Conservation in the Tropics. American Society of Agronomy/Soil Science Society of America Special Publication 43, Madison, pp. 27–39. Sua´rez de Castro, F., Rodrı´guez, A., 1962. Investigaciones sobre la erosio´n y la conservacio´n de suelos en Colombia. Federacio´n Nacional de Cafeteros, Bogota´.

N. Hoyos et al. / Journal of Hydrology 314 (2005) 177–191 Universidad Nacional de Colombia., 1997. Estudio de oferta y demanda hı´drica en la subregio´n No. 1 del Departamento de Risaralda (municipios de Pereira, Dosquebradas, Santa Rosa de Cabal y Marsella). Posgrado en aprovechamiento de recursos hidra´ulicos, Medellı´n. Van Dijk, A.I.J.M., Bruijnzeel, L.A., Rosewell, C.J., 2002. Rainfall intensity-kinetic energy relationships: a critical literature appraisal. Journal of Hydrology 261, 1–23. Wilcox, B.P., Rawls, W.J., Brakensiek, D.L., Wight, J.R., 1990. Predicting runoff from rangeland catchments: a comparison of two models. Water Resources Research 26, 2401–2410. Wischmeier, W.H., 1959. A rainfall erosion index for a universal soil-loss equation. Soil Science Society of America Proceedings 23, 246–249.

191

Wischmeier, W.H., Smith, D.D., 1978. Predicting rainfall erosion losses: a guide to conservation planning. USDA Handbook 537, Washington, DC. Yu, B., 1998. Rainfall erosivity and its estimation for Australia’s tropics. Australian Journal of Soil Research 36, 143–165. Yu, B., Rosewell, C.J., 1996a. Rainfall erosivity estimation using daily rainfall amounts for south Australia. Australian Journal of Soil Research 34, 721–733. Yu, B., Rosewell, C.J., 1996b. An assessment of a daily rainfall erosivity model for New South Wales. Australian Journal of Soil Research 34, 139–152. Yu, B., Hashim, G.M., Eusof, Z., 2001. Estimating the R-factor with limited rainfall data: a case study from peninsular Malaysia. Journal of Soil and Water Conservation 56, 101–105.