Sediment loss and its cause in Puerto Rico watersheds

Abstract. A major environmental concern in the Commonwealth of Puerto Rico is increased sediment load to water reservoirs, to estuaries, and finally to coral reef areas outside the estuaries. Sediment deposition has significantly reduced the storage capacity of reservoirs, and sediments, with their associated contaminants and nutrients that are adsorbed, can stress corals and negatively impact reef health. To prevent and manage sediment loss it is therefore important to understand local soil erosion and sediment transport processes. The main objective of this study was to determine the influence of landscape characteristics on sediment loss. We analyzed available precipitation and sediment data collected in Puerto Rico during the past three decades, as well as information on land use, soil properties, and topography. Our partial least squares analysis was not very successful in identifying major factors associated with sediment loss due to the complexity of the study's watersheds; however, it was found that topography and rainfall factors do not play a leading role. Sediment loss from the ridge watersheds in Puerto Rico was mainly caused by interactions of development, heavy rainfall events (especially hurricanes), and steep mountainous slopes associated with the ridges. These results improve our understanding of sediment loss resulting from changes in land use/cover within a Puerto Rico watershed, and allow stakeholders to make more informed decisions about land use planning.


Introduction
Water bodies and coastal areas around the world are threatened by increases in upstream sediment and nutrient load, which influence drinking water sources, aquatic species, and other ecologic functions and services of streams, lakes, and coastal water bodies (Haycock and Muscutt, 1995;Verhoeven et al., 2006). Puerto Rico (PR) faces considerable challenges regarding sustainable land use and current land use effects on adjacent coastal ecosystems and their services. Previous studies in PR have shown that sediment contaminants have increased 5-to 10-fold since pre-colonial levels, with a 2-to 3-fold increase in the last 40-50 years (Sturm et al., 2012). The increased sediment contamination could originate from anthropogenic activities such as agriculture and urban development, or from natural erosion (Tong and Chen, 2002;Gellis, 1993Gellis, , 2013. The primary concern regarding increased sediment load to reservoirs is that sediment depo-sition significantly reduces storage capacity. This sedimentation can also reduce photosynthetic activity of aquatic plants and algae in reservoirs, and increase water treatment costs for domestic and industrial uses (Estades Hernández et al., 1997). Sediments and pollutants adsorbed to them can ultimately reach offshore reef areas and stress or kill the corals comprising the reef. Reducing sediments that reach reservoirs and offshore reefs is key to managing and conserving natural resources (Morgan, 1986).
Watershed-scale studies of the potential effect of land use changes on water quality are essential to minimizing water pollution. Various studies have linked stream pollutants to landscape variables using process-based hydrological models (Jha et al., 2010;Ullrich and Volk, 2009;Hu and Yuan, 2013) and/or statistical methods (Lenat and Crawford, 1994;Nie et al., 2011;Mbonimpa et al., 2014). Process-based hydrologic models have been used successfully to character-ize watershed processes and sources of stream pollutants, but they require detailed input data, which may not be available for some areas. For instance, Hu and Yuan (2013) showed the difficulty of calibrating a SWAT model for the Guánica Bay, PR, watershed due to limited data for numerous reservoirs and dams in the basin. Other studies, however, have demonstrated statistical relationships between landscape metrics and water quality. For example, when Lenat and Crawford (1994) analyzed water samples from three watersheds having different dominant land use (forest, urban, or agricultural) in the Piedmont ecoregion of North Carolina, they found urban land use was the greatest contributor to sediment loss. Mbonimpa et al. (2014) identified urban land use and agricultural land growing corn as the main factors that caused increases in total suspended sediment and total phosphorous in streams, using partial least squares (PLS) regression analysis. The objective of our study was to evaluate the influence of landscape characteristics on levels of suspended sediment and identify factors that impact sediment loss, using statistical methods. Identifying major factors causing sediment loss will help stakeholders make better planning decisions about future land use and sediment control.

Study area
PR is located in the western Atlantic and Caribbean region. The island was almost entirely covered with forests prior to European colonization in the 16th century, and it was not until the beginning of the 19th century that forests were cleared for planting sugarcane, coffee, cotton, and tobacco, with coffee as the primary agricultural crop. Many watersheds are located in the major coffee-growing zone.
PR receives an average of 1651 mm of precipitation annually, delivered by infrequent but high-intensity rainfall events. The top mountains receive a much higher amount than coastal areas. The island is located in the Caribbean hurricane belt, and hurricanes are the number one weather threat because of damaging high winds, waves, and large volumes of rain. Six to 10 hurricanes threaten this region annually, and several powerful ones (Hugo -1989;Hortense -1996;Georges -1998;Lenny -1999;and Irene -2011) caused catastrophic damage to the environment. Various parts of PR also experience severe yearly storms that cause floods, tree falls, and landslides (Nagle et al., 1999).
Rain intensity, duration, frequency, and areal extent are the most important factors contributing to erosion (Jiang et al., 2008); in addition, runoff generated from storm events of high intensity or long duration transports large quantities of suspended sediment (SS); annual SS loss in some PR watersheds can be as high as 130 t ha −1 year −1 . Furthermore, landslides triggered by heavy rainfall contribute an average of 3 t ha −1 year −1 into river channels (Zack and Larsen, 1994). Adding to the transport of SS is severe stream bank ero- sion, which often occurs in sun-grown coffee-growing areas (CWP, 2008).

Data acquisition and preliminary analysis
Monthly precipitation data were obtained from the National Oceanic and Atmospheric Administration's National Climatic Data Center (NOAA-NCDC; http://www.ncdc.noaa. gov/cdo-web/search). Thirty-three NOAA-NCDC weather stations are located on the main island of PR (Fig. 1). The annual total for each water year was calculated by monthly precipitation data and expressed in mm year −1 . There are 31 United States Geological Survey (USGS) stations monitoring stream flow and SS in PR watersheds ( Fig. 1). To meet the assumption of independence of watersheds (observations) for regression analysis, nested watersheds (i.e., those situated within larger watersheds) were not included in this analysis, leaving 20 independent watersheds with USGS monitoring. Further investigation of these 20 watersheds revealed that some stations have very short monitoring periods and that similarities of monitored SS concentrations and load exist between different watersheds. We therefore eliminated nine stations with short monitoring periods and/or similar SS concentrations and load (Table 1), leaving 11 watersheds in this study.
Annual (water year) statistics of discharge (cubic feet per second), suspended sediment concentration (milligrams per liter), and suspended sediment discharge (US tons per day) were downloaded for the 11 monitoring stations. For SS load, annual statistics values were first converted to annual t year −1 , then normalized by dividing by the drainage area and expressing as t ha −1 .
Thiessen polygons, created from available weather stations using ArcMap (ArcGIS10), were overlaid onto the watershed boundary. Areal average precipitation was calculated for each watershed based on where the polygons fall.
Distributions (percent of total watershed area) of land use type (Table 2) and soil type (Table 3) for each watershed were determined by overlaying the National Land Cover Database (NLCD) 2001 land use map and SSURGO soil map onto the watershed boundaries. The primary land covers in the studied watersheds were evergreen forest and grassland/herbaceous (ranging from 50 to 100 %) ( Table 2). Other land uses include developed/urban (ranging from 0 to 47 %), pasture/hay (0-31 %) and scrub/shrub (0-7 %). Shade coffee plantations and secondary forests regenerated from abandoned pastures and coffee plantations are both classified as shrub because of their low canopies.
Soil properties such as erodibility, texture, and hydraulic conductivity were retrieved from the SSURGO soil database (Table 4). In summary, soils in the studied watersheds varied from very deep to shallow, well drained to poorly drained, and slowly permeable to rapidly permeable soils; the majority are very to moderately deep, well drained, and very to moderately permeable.
We calculated the correlation between the calculated watershed precipitation and SS concentrations and load. Land use (Table 2) and soil distribution (SSURGO soil types, Table 3), slope (Table 5), and calculated precipitation subsequently formed predictors for loss of sediments from each watershed.
Annual average SS load and concentration were calculated for each monitoring station based on available data from 1983 to 2011. USGS stations do not have measured data in exactly the same time period, but they overlap. Annual average of areal rainfall for each watershed was also calculated using data from 1983 to 2011.

Correlation and partial least squares analysis
Correlations among the variables were calculated (R statistical software: https://www.R-project.org) to assess relationships between variables. Correlation analysis was also performed between precipitation and SS load/concentration to evaluate the impact of precipitation on SS load/concentration for individual watersheds. PLS was used to find the association between measured SS load, rainfall, and landscape characteristics (e.g., land use type, soil type, and topography). A small sample size, a large number of predictors, and collinearity between predictors prevented the use of standard multivariate regression (Yeniay and Goktas, 2002); however, Table 3. Soil type distribution (in percent of total watershed drainage area).
Soil type (% of total area) USGS s9538 s9537 s9542 s9532 s9539 s9533 s9536 s8369 s9531 s9535 s9553 s9530 s9541 s9543 s9518 s9540 s9552 s9529 s9534 s9544 PLS has been shown to work well under these constraints (Abdi, 2010;Nash and Chaloud, 2002). PLS extracts orthogonal factors (latent variables), requiring these latent variables to explain as much of the covariance between measured response (SS) and predictors (landscape variables) as possible. This is followed by regression of predictors on response (Helland, 1988;Höskuldsson, 1988). Predictor coefficients (magnitude and direction) from the PLS regression can define the role and influence of predictor variables on response.
Magnitude of the coefficient indicates the weight and degree to which the predictor influenced the response; direction of influence is "+" for an increase in the coefficient and "−" for a decrease.

The relationship of SS load/concentration with precipitation
Temporal variations in SS concentration and SS load were analyzed for each study monitoring station; the highest SS concentration occurred in 1989 (Fig. 2) and the highest SS load occurred in 1998 (Fig. 3) at monitoring station 3. Hurricanes Hugo (1989) and Georges (1998) generated serious soil erosion and sediment loss; hurricanes are PR's top weather threat because of the damage they cause (Nagle et al., 1999). Hurricane Georges may also have been responsible for higher SS loads that occurred at monitoring stations 1, 10, 11, and 9 (in order from high to low) in 1998, and also for higher SS concentrations at monitoring stations 10, 1, 9, and 11 (in order from high to low). No monitoring information was available for 1989 for monitoring stations 1, 10, 9, and 11. Although SS monitoring started in 1986 at monitoring station 11, monitoring information was not available for 1989 due to equipment malfunction caused by the hurricane.  As expected, precipitation was a strong predictor for SS concentration and SS load for most monitoring stations, and it explains 19-94 % variation in SS concentration and 19-83 % variation in SS load (Table 6).

Comparison of SS load among different USGS monitoring stations
The watershed monitored by station 3, where the highest SS load occurred, (Fig. 3) has clay soils with relatively low erodibility (Tables 3 and 4); it also has the lowest mean slope (Table 5) and an annual average rainfall of 4506 mm, which is substantially lower than the mean of the study watersheds (Table 7). All these factors suggest that watershed 3 should have a much lower SS load; however, it has the highest percentage of developed land (47 % - Table 2), which may out-weigh the other factors (Lenat and Crawford, 1994;Mbonimpa et al., 2014).
The driver for SS load in the watershed of monitoring station 1 (second highest load, Table 7) was not developed lands (only 8 % - Table 2), but likely was associated with the much steeper slopes (Table 5). This is further borne out by the SS load from the watershed of monitoring station 10 (Fig. 3), which has clay soils with low erodibility (Tables 3 and 4) and a low percentage of developed land (4 %) but also has the highest mean slope (41 % - Table 5).
The lowest SS load occurred in the watershed at monitoring station 6 ( Table 7). Although the primary soil was moderately erodible (Tables 3 and 4), with the highest annual average rainfall (Table 7) and moderately steep slopes (Table 5), the watershed had 100 % land cover by evergreen forest (Table 2). Land cover for the watershed of monitoring station 7 also had a high degree of natural cover (70 % evergreen forest and 26 % grassland/herbaceous), as well as a very low SS load (Table 7).

Results of partial least squares analysis
PLS regression coefficients indicated the degree and direction of association for land uses, soil types, and slope characteristics (Fig. 4). Of the land use categories, developed and barren lands appear to greatly increase SS load (Fig. 4), while land covered by evergreen forests, scrub, grass, and pasture decreases the load. The high regression coefficient associated with developed land use shows it is the strongest predictor for SS load. Watershed 3 had the highest developed land use and the highest SS load (Fig. 3), which is discussed in Sect. 3.2.
Slope can also have a strong and directional influence on SS loading (Fig. 4). For watersheds with relatively low developed land use (Table 2), the highest sediment loads came from watersheds 1 and 10 (Table 7). Both have very high mean percent slope (Table 5), which may be responsible for their higher SS loads ( Table 7). Lack of steep slopes may also reduce the amount of SS loading, as in watershed 5 (Table 7). Generally, soil loss is linearly related to the sine of slope angle for slopes ranging from 9 to 55 % (Liu et al., 1994). The higher the erodibility of a soil, the higher the potential for increased sediment loading. For example, the primary soil type of s9542 (Pellejas-Lirios) in watershed 1 has an erodibility coefficient of 0.17, which could increase SS load (Fig. 4). Soil erosion potential rises as soil particle size content larger than 0.05 mm increases, since the material is less cohesive (Morgan, 2001).
Watershed 8 had a high SS load (Table 7) but also received the highest rainfall (Table 7). On the other hand, watershed 6 received the same amount of rainfall but had the lowest SS load of any of the watersheds. Therefore, even with increased rainfall, other factors associated with land use and soil type may have reduced SS loading. Watershed 6 has 100 % evergreen forest land cover (Table 2) and a primary soil type of s9542, which has an erodibility coefficient of 0.17. It may be that large amounts of forested land cover and soil that is not high erodible can lead to decreased SS load in a watershed, even in high rainfall areas. Neither of these watersheds had a weather station, however, and the same weather stations were assigned to them. The weather station is closer to monitoring station 8 than monitoring station 6. We therefore cannot rule out that rainfall data collected at the weather station did not represent actual rainfall on either watershed.
Additional regression analysis was subsequently used to examine watersheds with high vegetation cover (< 10 % developed and barren land) and low vegetation cover (> 10 % developed and barren land) ( Table 2). Precipitation had less impact on SS load in the highly vegetated watersheds. Slope of the regression for highly vegetated land suggests that vegetative cover can retard soil loss, while watersheds with reduced vegetative cover appeared to have increased soil loss.
PLS analysis was not efficient in identifying key factors associated with sediment loading from the study watersheds. This was likely due to the complex attributes of these watersheds that often acted jointly, in opposition, or in both directions simultaneously to affect SS loading from the watersheds. Land use was identified as a strong predictor for sediment loss from the study watersheds, however, and land use is a factor that can be managed. Vegetation can be an important factor affecting soil erosion processes, and in areas with serious soil erosion and sediment loss, the natural vegetation has often been destroyed (Gyssels et al., 2005;Shi and Shao, 2000;Zhou et al., 2008). Natural vegetative cover may be able to compensate for the erosive potential of high precipitation and steep slopes (e.g., watershed 6, Table 7). Where development must occur in PR, areas with highly erodible soils, steep slopes, and high precipitation should be avoided.

Conclusions
Development was the one controllable factor associated with SS load and SS concentration increases in the PR study watersheds. SS loads were also influenced by precipitation, steep terrain, and soils with higher erodibility. In watersheds with high percentages of natural vegetative cover (e.g., watershed 6, Table 7, evergreen forest), SS loading was low, even in steep terrain with high precipitation. We found PLS to be unsuccessful in identifying the main factors causing high SS loading in such complex watersheds. Future studies are needed to examine the spatial distribution of different land cover types within a watershed that mitigate SS loads from complex watersheds.