Determining Rainfall Erosivity in Costa Rica: A Practical Approach

Abstract Rainfall erosivity (R-factor) is an important variable used in soil erosion estimation models. In Costa Rica, the R-factor was computed for 106 stations across the country by Wilhelm-Günther Vahrson in 1990. These results provided the main input information for this study, which used them to estimate the R-factor for Costa Rica. Regression equations were computed to estimate the R-factor for the Caribbean slope, the Pacific slope, and the country as a whole. Forward stepwise analysis and multiple regression analysis were employed to determine the regression coefficients for each developed equation. Elevation and monthly rainfall had a strong influence on the definition of the R-factor equations. The Modified Fournier Index variable was included only in the national-scale equation as a good proxy for the mean annual precipitation effect at each site. Inclusion of elevation in all equations reflects the importance of the transitional effect of high-intensity convective rainfall in the lowlands and low-intensity orographic rainfall in the highlands. This study provides an easy way to estimate the R-factor using regression equations that require only simple and readily available geophysical information. The use of these equations in conjunction with soil and land-use maps as well as digital elevation models will allow the estimation and evaluation of soil erosion on a watershed scale in Costa Rica. This will also improve the application of other hydrological models that require soil erosion as an input variable to estimate sediment yields.


Introduction
Soil erosion and sediment yield reduce water storage capacity in lakes and reservoirs, modify fluvial and coastal geomorphology, pollute water resources, and affect the hydrology and ecology of downstream ecosystems (Summerfield 1991;Wood and Armitage 1997;Syvitski et al 2005;Zhang et al 2006;Huggett 2007;Bilotta and Brazier 2008). Highlands and steep mountains are very susceptible to soil degradation; this process is triggered by the steep slopes, land-cover change, and high rainfall erosivity capacity (Zhou et al 2008;García-Ruiz 2010). Therefore, information on soil erosion estimations in mountainous regions is highly important in efforts to prevent suspended sediment yield and loss of soil. However, in mountainous regions, particularly in the tropics, there is a lack of physical data and practical hydrological models to help in the design and evaluation of soil conservation programs. On the other hand, it is well known that the quality and geographical coverage of hydrometeorological data in poor countries are obstacles to hydrological studies and water resource planning (Kim et al 2005). This is an important worldwide issue for watershed management, land-use planning, and environmental impact analysis, since soil erosion estimation is an important way to evaluate the impact of soil conservation practices and land-use change.
Soil erosion prediction and conservation have been an important issue worldwide, and the development of empirical methods to estimate annual rates has facilitated estimates of soil loss from specific land units (eg plots, catchments, and watersheds). Soil erosion estimates from regions with scarce physical data have employed the universal soil loss equation (USLE) (Wischmeier and Smith 1978) or the revised version (RUSLE) (Renard et al 1996a;Calvo-Alvarado and Gregory 1997). However, these methods had not been extensively used in tropical regions due to the reduced availability of physical data to calibrate the parameters, such as rainfall erosivity, soil erodibility, land cover, topography, and soil conservation (Calvo 1998).
Rainfall erosivity (the R-factor) expresses the potential capacity of raindrops to cause soil detachment, a process that is strongly dependent on rainfall intensity (Ferro et al 1999;Angulo-Martínez et al 2009). The R-factor is expressed in megajoules (MJ) mm ha 21 h 21 y 21 and is defined as the mean annual sum of individual storm erosion index values: Total storm kinetic energy (MJ) multiplied by the maximum rainfall intensity in 30 minutes (Wischmeier and Smith 1978). The R-factor is particularly difficult to calibrate due to the data set that is necessary and to the required geographical sampling distribution over time and space; such requirements are impossible to fulfill in poor tropical nations. In tropical regions, rainfall intensities of 40 mm h 21 to more than 160 mm h 21 can occur in time intervals of 10 seconds to 20 minutes (Ajayi and Ofoche 1984). Consequently, Rfactor estimation requires extensive records of rainfall events at small time intervals (30 minutes or less) to extract needed information such as rainfall duration, intensity, and total volume. In addition, strong rainstorms in the tropics are related to convective precipitation mechanisms, producing from 1% to 8% of the total annual rainfall in the Neotropics (Rickenbach et al 2012); however, in the highlands, the orographic rainfall increases as the result of the condensation of air humidity as temperature decreases with altitude, which has an impact on the definition of the R-factor that has not been fully studied in tropical regions.
Recently, the widespread use of geographical information systems (GIS) in combination with remotesensing technologies has prompted the use of the USLE or RUSLE to simulate the impact of land-use change on total watershed soil erosion and sediment yields (Calvo-Alvarado and Gregory 1997;Calvo 1998;Jain and Das 2010;Chen et al 2011;Prasannakumar et al 2011;Lin et al 2013). As a consequence, a possibility has arisen of developing alternative methods of estimating the R-factor as a function of geographical variables that are easy to obtain in data-scarce locations. In Honduras, R-factor determination was carried out using linear regressions, obtaining good results based on elevation and mean annual precipitation (Mikhailova et al 1997). However, information related to methods for estimating the Rfactor and determining soil erosion in tropical regions is still scarce. Hence, this article provides a simplified regression relationship between the R-factor (dependent variable) and easily obtained geographical data (independent variables) at all locations across Costa Rica.

Methodology Site description
Costa Rica is located in the Neotropical region (8u-11uN; 82u-86uW), with an annual precipitation of 1300-8500 mm y 21 . Costa Rica has special hydrological characteristics due to a combination of several factors. First, the Talamanca, Central, and Guanacaste Mountain Ranges divide the country longitudinally into the Caribbean and Pacific slopes. Because of the steep topography and the prevailing trade winds, orographic precipitation has a strong influence. The northwestsoutheast ridges of the Caribbean slope are rainy almost all year-round. The northeast trade winds are active from November through March. From May to October, the entire country is influenced by the passage of the Intertropical Convergence Zone (ITCZ) and by the southwest trade winds. As a result of these influences, only the Pacific slope experiences a dry season between December and April, because the humidity coming from the northeast trade winds is retained on the Caribbean slope. Temporal and spatial distributions of rainfall in Costa Rica are also modulated by changes in oceanic and weather phenomena such as El Niñ o (warm phase) and its opposite, La Niñ a (cold phase), which are considered the highest expressions of climatic variability (Calvo 1990;Vahrson 1990;Manso et al 2005;Guzmá n and Calvo-Alvarado 2013). Rainwater chemical analysis has detected a strong marine influence on Costa Rica's rains (Hendry et al 1984;Eklund et al 1997;Clark et al 1998), as the result of water vapor transport by the wind from the Caribbean Sea and Pacific Ocean to inland locations.

Data
Rainfall erosivity values (the dependent variable) were obtained from Vahrson (1990) for 115 meteorological stations across Costa Rica ( Figure 1). The R-factor was originally reported in 100 foot ton inch acre 21 h 21 , and these values were translated into MJ mm ha 21 h 21 y 21 by multiplying each value by 17.02 (Renard et al 1996b). The mean annual rainfall erosivity estimated by Vahrson (1990) made use of Equation 1, proposed by Woodward (1975). This equation calculated rainfall erosivity (R) based on maximum annual rainfall (MAR) events of 6 hours duration and 2 years return period (MAR 2;6 ). Ferná ndez and Salas (1999) define the return period ''as the average number of trials (usually years) to the first occurrence of an event of magnitude greater than a predefined critical event.'' The MAR 2;6 records were obtained from pluviograph data with 10 years of records from each of the selected meteorological stations. As a consequence of Costa Rica's position on the Central American isthmus, extreme events are very common due to the frequent impact of tropical storms and convective storms associated with the passage of the ITCZ. Consequently, extreme rainfall events analysis was performed using the Gumbel (1945) method, which is well applicable to the climatic conditions of Costa Rica (Vahrson and Fallas 1988).
Data from the meteorological stations included mean monthly precipitation, mean annual precipitation (MAP), elevation, and the Modified Fournier Index (F index). The F index as proposed by Arnoldus (1980), was used in this analysis as a substitute for precipitation in the R-factor estimation. Ferro et al (1999) concluded that this index is linearly correlated to the R-factor. The F index depends on mean monthly rainfall (p i ) and mean annual precipitation (P) as shown in Equation 2.
Based on data availability, 106 of the 115 meteorological stations were selected for this analysis. The complete data set is presented in the Supplemental data, Table S1 (http://dx.doi.org/10.1659/MRD-JOURNAL-D-13-00062.S1).

Data analysis
The data were organized into three groups to account for regional effects, as previously discussed by Calvo-Alvarado et al (2009): data from the Pacific slope, the Caribbean slope, and a combination of the two to allow analysis on a national scale. For each group, the R-factor was determined through multiple regression analysis. The selected independent variables are those that are easy to obtain for each site, and hence variables such as rainfall intensity were omitted in this analysis. The nonselection of rainfall intensities is due to the fact that in Costa Rica-and generally in countries in the tropics-rainfall networks with extended and complete records (more than 10 years) usually collect only daily rainfall. The main purpose of our work was to develop equations to estimate the R-factor for any site in Costa Rica that lacks this information, and this can only be accomplished using variables that are easy to obtain.
Forward stepwise analysis was applied to select the best independent variables to estimate the R-factor for each station. This analysis selected the independent variables used to perform the multiple regression analysis obtaining the predictor coefficients. The software STATISTICA 6.0 (StatSoft Inc 2003) was employed in all the analyses, to which a probability value of 0.05 was applied. All the assumptions for normal distribution were considered and tested to verify the validity of the analysis.

Results
The precipitation (P) from the 106 stations is 3077.3 mm y 21 . March is the driest month (84.6 mm mo 21 ), and October is the wettest (408.1 mm mo 21 ). The Pacific slope has a well-defined dry season from December to March, with precipitation rates lower than 100 mm mo 21 , and another short dry season between June and July. The Caribbean slope lacks a dry season due to the constant influence of northeast trade winds responsible for precipitation rates .100 mm mo 21 , but rainfall diminishes during March ( Figure 2). P is 2725.2 mm y 21 for the Pacific slope (60 stations) and 3536.7 mm y 21 for the Caribbean slope (46 stations). This difference of more than 800 mm y 21 is the result of the humidity carried by the trade winds from the Caribbean Sea. Vahrson (1990) identified a decreasing trend of Rfactor with elevation. The Costa Rican lowlands show higher erosive rains than the uplands. The lowlands have a variety of rains as a consequence of local convective and frontal systems, which lead to a wide range of R-factor values. On the other hand, locations higher than 2000 m above sea level (masl) reported R-factor values lower than 5000 MJ mm ha 21 h 21 y 21 (Figure 3). The characteristic rain of the uplands results from the orographic effect and the presence of relatively few local convective systems. Rfactor values decrease at higher elevations on both slopes, the Pacific and the Caribbean.
Increasing P magnifies the R-factor on both slopes ( Figure 4A). A similar trend is observed with the F index, due to its linear relation with the P value ( Figure 4B); however, a wide variation among F index and R-factor values can be observed along the value ranges ( Figure 4C).
Forward stepwise analysis underlines the effect at the national scale of monthly precipitation registered in January, April, and November, plus the F index and elevation (Equation 3). The choice of monthly precipitation depicts the importance of precipitation seasonality, while elevation introduces topography as a key variable that indirectly considers the effect of orographic rainfall in the R-factor definition. The Rfactor on the Caribbean slope is affected by April and July precipitation, as well as elevation (Equation 4), while for the Pacific slope, it is defined by September precipitation and elevation (Equation 5). In these equations, R represents the R-factor; F is the Modified Fournier Index (mm); E is elevation (masl); p with subscripts denotes monthly precipitation (mm) (mm month 21 ) for January (Jan), April (Apr), July (Jul), September (Sep), and November (Nov); N s refers to the national scale; C s refers to the Caribbean slope; and P s refers to the Pacific slope; and p (no subscript) is the probability value where R N S , R C S , and R P S are the rainfall erosivity values (MJ mm ha 21 h 21 yr 21 ) at national scale, Caribbean slope, and Pacific slope, respectively. Figure 5 shows the relationship between the measured and calculated R-factors.

Discussion
Hastenrath (1968) noted the marked rainfall differences between the Caribbean and Pacific slopes of Central America. This differentiation is consistent with the pattern in Costa Rica, where the continental divide works as a topographic barrier that affects precipitation patterns (Calvo 1990). This study provided a strong argument to develop regional equations to estimate the R-factor. Additionally, a strong elevation effect is included in all equations estimating the R-factor across Costa Rica, with negative altitude coefficients in all equations. Goovaerts (1999) remarked on the potential use of elevation to calculate rainfall erosivity in Portugal and underlined the power of digital elevation models for this purpose. Even with high MAP, Costa Rica shows a constant reduction in the R-factor above 1000 masl on both slopes. This pattern is probably associated with the dominant effects of lowintensity orographic rainfall in the highlands, which contrast with the effects of the frequent high-intensity convective storms of the lowlands. Guswa et al (2007), using water isotopic signatures, showed the presence of orographic precipitation in Monteverde, Costa Rica (850 masl). This type of precipitation accumulates more than 2000 mm y 21 over 3000 masl (Horn 1989, cited by DeForest Safford 1999. High elevations in Costa Rica are characterized by unimodal rainfall seasonality, with 9 months of precipitation greater than 50 mm mo 21 (DeForest Safford FIGURE 4 Relations between mean annual precipitation and R-factor (A), mean annual precipitation and F index (B), and F index and R-factor (C) in Costa Rica, based on data from 106 stations.
FIGURE 5 Relationships between measured and calculated R-factor calculations based on Equations 1, 2, and 3. 1999). In these highlands, orographic precipitation is registered over long periods with low rainfall intensity, illustrating the prevailing effect of orographic precipitation on mountaintops.
The R-factor for the Pacific slope is strongly affected by September's rainfall due to the high water volume just after the short dry season that takes place between June and July. On the Caribbean slope, April and July have the most influence on the R-factor; this effect is related to the seasonal transition from dry to wet in this region. In April, the rains begin to increase in volume; they reach a peak in July. On a national scale, the R-factor shows a more complex interaction with monthly precipitation. The three months that are included in the national-scale equation mark the edges of the dry season (January and April) and the end of the rainy season (November).
Seasonal rainfall distribution has been shown to have an important effect on the R-factor in other regions around the world. Wilkes and Sawada (2005) proved the effect of this variable in the Great Lakes region of North America, while the effects of global phenomena such as El Niñ o and La Niñ a have been shown to modify the volume and intensity of rainfall in Peru (Romero et al 2007).
The R-factor has been positively correlated with MAP and the F index, but the manner of estimating varies among authors. Mannaerts and Gabriels (2000) suggested a linear relationship between MAP and the R-factor. Equation 6 depicts the procedure used by Ferro et al (1999) and Diodato (2004), employing the F index and MAP as predictor variables. Eltaif et al (2010) worked with MAP as an independent variable to estimate the R-factor on an annual basis in Jordan, as shown in Equation 7.
where x is the P or F index.
where x is the P, while a and b represent coefficients that describe the local conditions. In Costa Rica, P was not a significant variable in the prediction of the R-factor, probably because the R-factor is a combined effect of elevation, the passage of ITCZ, and the trade winds, which leads to a multitude of combinations of orographic, convective, and frontal storms during the year in each site. The F index, on the other hand, was included in the national-scale analysis as the best proxy to gauge the P seasonality.
R-factors have been frequently estimated for mountainous countries such as Switzerland (Meusburger et al 2012), Peru (Romero et al 2007), and Spain and Italy with their Mediterranean rainfall regime (Ferro et al 1991(Ferro et al , 1999Brath et al 2002;Diodato 2004Diodato , 2006Angulo-Martínez et al 2009;Diodato and Bellocchi 2009). Only a few R-factor estimates have been conducted for tropical countries such as Costa Rica (Vahrson 1990) or for arid and semiarid countries such as Jordan (Eltaif et al 2010) and Cape Verde (Mannaerts and Gabriels 2000). Extrapolation and prediction of the R-factor have also been performed through maps in South Africa (Smithen and Schulze 1982) and Costa Rica (Vahrson 1990).
The work performed in Costa Rica by Vahrson (1990) provides a base map of rainfall erosivity, but today's technology allows the spatial estimation of variables such as the R-factor, enabling the use of well-defined equations to interpolate this parameter. The estimation equations presented here can support the production of more accurate maps such as was done for the Great Lakes region by Wilkes and Sawada (2005). These equations can also facilitate the application of the USLE or RUSLE soil erosion equations by using digital elevation models and GIS.

Conclusions
Data availability is an important factor to consider in the analysis of hydrological processes. The lack of meteorological stations with long, reliable records in the tropics precludes the application of many hydrological models used worldwide. Empirical models for estimating soil erosion, such as USLE and RUSLE, require a good series of precipitation data to compute the R-factor. Because highlands are the areas most susceptible to soil erosion, an accurate R-factor estimation is highly necessary. The limited number of meteorological stations in Costa Rica's mountain regions underlines the importance of prediction equations for estimating R-factors.
This study provides an easy way to estimate the Rfactor using regression equations that require only simple and easy-to-collect geophysical data. The use of these equations in conjunction with soil maps, land-use and ground-cover maps, and digital elevation models will allow robust estimation and evaluation of soil erosion on a watershed scale. This will also improve the application of other hydrological models to estimation of suspended sediment yields, which requires soil erosion information as the main input variable.
The application of the equations developed in this study to other Central American countries that share a similar climate must be conducted with caution, due to microclimatic differences and the topographical characteristics of the places under evaluation. As a consequence, it is important to generate information on true R-factor for a few stations in order to test these equations and evaluate their precision and applicability.