Next Article in Journal
Development of an Ice Jam Flood Forecasting System for the Lower Oder River—Requirements for Real-Time Predictions of Water, Ice and Sediment Transport
Next Article in Special Issue
Voluntary Management of Residential Water Demand in Low and Middle-Low Income Homes: A Pilot Study of Soacha (Colombia)
Previous Article in Journal
The Use of Non-Conventional Water Resources as a Means of Adaptation to Drought and Climate Change in Semi-Arid Regions: South-Eastern Spain
Previous Article in Special Issue
Landscape Drivers and Social Dynamics Shaping Microbial Contamination Risk in Three Maya Communities in Southern Belize, Central America
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Streamflow Response to Persistent Drought in a Coastal Tropical Mountainous Watershed, Sierra Nevada De Santa Marta, Colombia

1
Department of History and Social Sciences, Universidad del Norte, Barranquilla 080001, Colombia
2
Instituto de Geología, Universidad Nacional Autónoma de México, Coyoacán 04510, Mexico
3
Smithsonian Tropical Research Institute, Panama City 32401, Panama
4
Global Water Center, University of Nevada, Reno, NV 89557, USA
5
Department of Geography, University of Vermont, Burlington, VT 05405, USA
6
Integrated Data Repository, Clinical and Translational Science Institute and UF Health, University of Florida, Gainesville, FL 32610, USA
7
Department of Civil Engineering, Universidad del Norte, Barranquilla 080001, Colombia
8
Department of Physics and Geosciences, Universidad del Norte, Barranquilla 080001, Colombia
9
Department of Geology, University of Regina, Regina, SK S4S 4P5, Canada
*
Author to whom correspondence should be addressed.
Water 2019, 11(1), 94; https://doi.org/10.3390/w11010094
Submission received: 26 November 2018 / Revised: 20 December 2018 / Accepted: 24 December 2018 / Published: 8 January 2019
(This article belongs to the Special Issue Current and Emerging Issues Surrounding Water in the Americas )

Abstract

:
Droughts constitute natural hazards that affect water supply for ecosystems and human livelihoods. In 2013–2016, the Caribbean experienced the worst drought since the 1950s, and climate projections for the southern Caribbean predict less rainfall by the end of the 21st century. We assessed streamflow response to drought for a watershed in the Colombian Caribbean by analyzing the effects of drought length and land cover on streamflow recovery. We generated a calibrated SWAT model and created annual and monthly drought scenarios from rainfall records. We used our model to predict water yield for selected land covers (wet forest, shade coffee, shrub, and dry forest) under drought conditions. Annual scenarios resulted in water yield reductions of ~15 mm month−1 (wet forest, coffee, and shrub) and 5 mm month−1 (dry forest) for the first month after a two-year drought. Maximum water yield reductions for monthly scenarios occurred after a 10-month drought and were ~100 mm month−1 (wet forest, coffee, and shrub) and 20 mm month−1 (dry forest). Streamflow recovered within nine months (annual scenarios), and two to eight months (monthly scenarios) after drought termination. Drought response seems to be conditioned by climatic factors (rainfall seasonality and spatial variability) and catchment properties.

1. Introduction

In the past centuries, human activities have been the principal factor driving environmental and climate change [1]. For coastal regions such as the Caribbean, these changes include land use change, extreme weather conditions, rising sea levels, increasing temperatures, and changes in precipitation patterns [2] which can have a detrimental effect on water supply. Within this context, coastal watersheds are critical in many regions of tropical America and the Caribbean [3]. Although the region as a whole has abundant freshwater resources, there are large disparities in water availability per person [3]. In terms of climate projections, models for the Caribbean predict changes in seasonal and annual precipitation, with a decrease in annual and summer (JJA) precipitation, and a slight increase in winter (DJF) precipitation by the end of the century [2,4,5]. The northern Caribbean is predicted to have more intense rainfall, whereas the southern Caribbean will experience less rainfall by the end of the century [6]. Temperature projections include a 1–5 °C increase by the end of the century [6,7,8]. In 2013–2016, the Caribbean experienced the worst drought since the 1950s related to El Niño [9]. Although the analysis of Herrera et al. (2018) [9] was limited to the Caribbean islands, it suggested that the projected drying trend for the region is already underway, with anthropogenic warming contributing to both its severity and spatial extent.
Drought, defined as “sustained period of below-normal water availability”, is classified as meteorological drought (i.e., rainfall deficit), soil moisture drought, hydrological drought (surface and subsurface water deficits) and socioeconomic drought (impacts of the above) [10]. Meteorological droughts may propagate to the soil and surface/subsurface components of the hydrologic system, affecting water supply for the agriculture, energy, industry, domestic, navigation, recreation sectors, as well as aquatic and terrestrial ecosystems [10]. Droughts tend to affect a large number of people, as they develop over large geographical areas and during extended periods (i.e., months to years). In addition, drought occurrence and development is complex due to the interaction of multiple natural and anthropogenic factors, including inter- and intra-annual climate variability, catchment geology and soil characteristics, vegetation and land use, and water extraction practices, all of which can either aggravate or alleviate their effects [11,12]. With regard to climate, regions in warm seasonal climates are particularly susceptible to droughts as rainfall deficits during the wet season may translate into a hydrological multi-year drought due to the low probability of recovery during the dry season [13].
The Colombian Caribbean region, located between the Caribbean Sea and the Andes mountains, is characterized by a warm seasonal climate that transitions into dry conditions to the north, and to wet conditions to the south, near the Panama border. This region supports around 20% of the country´s total population (~50 million in 2018) and exhibits a profound socioeconomic disparity when compared to the central Andean region [14]. Within the Colombian Caribbean, there is the Sierra Nevada de Santa Marta Massif (SNSM), a steep coastal triangular mountain range that rises to 5775 meters above sea level (masl) within 42 km from the coast (Figure 1). Due to its elevation gradient and proximity to the coast, it possesses nine types of vegetation biomes, at least 600 botanical genera, and no fewer than 3000 species of higher plants with a high level of endemism in biomes above 1000 masl [15]. The massif is drained by rivers that are critical for urban water supply (~1.5 million people), large and small-scale agriculture, fishing, tourism, and mining. Downscaled climate projections for the region show contradictory results for precipitation, with predictions for the 2011–2040 period ranging from 10–40% increase to 20–40% decrease relative to historic data [16,17]. Historical records show mixed results for trends in discharge and precipitation, and an overall increase in temperature on the order of 0.3 °C to 0.5 °C per decade [18,19,20,21,22,23]. Droughts in this region are exacerbated by strong rainfall seasonality and affect domestic water supply, agricultural and cattle production, fisheries (e.g., through increased salinity in coastal lagoons), fire occurrence, and the prevalence of certain diseases [24]. For instance, during the 2014–2016 drought, the northern Colombian Caribbean experienced the largest decrease in rainfall (44–70%) and widespread water availability issues, with > 80% of municipalities affected [24].
In this work, we analyzed streamflow response to persistent meteorological drought for a ~300 km2 watershed in the Colombian Caribbean region. We use the term ‘streamflow recovery’ to refer to streamflow levels returning to average conditions. Our specific objectives were to: (1) evaluate how meteorological drought length affects streamflow recovery, and (2) assess the effect of land cover on streamflow recovery. To do this, we used a hydrologic model (SWAT) combined with synthetic rainfall series representing droughts with increasing lengths, to predict water yield from different land covers. This paper is structured as follows. First, we present a description of our study area, hydrologic model setup and definition of drought scenarios. We then present the results of the hydrologic model calibration and validation, and the effects of meteorological drought on water yield for selected land covers in the watershed. We then discuss our findings and finish with some concluding remarks.

2. Methods

Our methods involved (1) setup, calibration and validation of the SWAT hydrologic model, (2) calculation of synthetic rainfall series to assess streamflow response to persistent dry conditions, (3) execution of the calibrated model with multiple synthetic rainfall series, and (4) analysis of model outputs (Figure 2). The following sections present a description of our study area and detailed methods.

2.1. Study Site

The Río Frío watershed is located on the southwestern slopes of the SNSM (Figure 1). The watershed covers an area of ~300 km2 with elevation ranging between 30 to 3400 masl on steep terrain with mean slopes of 60%. Downstream from the massif, the Río Frío runs on recent fluvial and marine lowland deposits, until it reaches the Santa Marta Lagoon system, the largest coastal lagoon in the country and part of the RAMSAR wetland convention due to its ecological importance (Figure 1). Currently, the river runs dry before reaching the lagoon as most of its water is used for irrigation in the lowlands.
Regional climate is characterized by latitudinal and elevation gradients, as well as by a strong intra- and inter-annual variability. Annual precipitation on the southwestern flank of the massif increases from north to south, and from lower to higher elevations. For instance, annual precipitation totals ~450 mm near the city of Santa Marta (Figure 1), compared to ~1400 mm further south (Prado Sevilla station, Figure 1). Intra-annual seasonality is controlled by the meridional movement of the Intertropical Convergence Zone (ITCZ) and related northeasterly trade winds. There is a marked dry season from December through March which transitions into a wet season running from May through November with a short drier season or ‘veranillo’ from mid-June through mid-August. Streamflow varies accordingly, i.e., streamflow records for the Río Frío (1967–2014) show minimum mean monthly discharge values from January through April (~5.6 m3s−1), and maximum values from August through November (~21.4 m3s−1). The contribution of baseflow to total streamflow also varies throughout the year, with maximum values during the dry season (e.g., mean of 82% in January) and minimum values during the wet season (e.g., mean of 62% in September). Interannual rainfall and discharge variability is related to several ocean-atmosphere oscillations, including El Niño/Southern Oscillation (ENSO), the Tropical North Atlantic index (TNA), the Atlantic Meridional Oscillation (AMO) and Pacific Decadal Oscillation (PDO) [25]. Annual mean temperature at sea level is 28 °C and decreases to 14 °C at 2200 masl.
Watershed relief is characterized by long and steep slopes developed on metamorphic and igneous rocks, with 70% of the watershed on slopes > 50%. Soils in the watershed are mostly shallow (<0.5 m) to moderately deep (~1.5 m) Inceptisols and Entisols with medium textures (i.e., loam to sandy loam) and some fine textured horizons. Soils are moderately to well drained, except where the bedrock is close to the surface [26]. Native vegetation varied according to precipitation amount and seasonality, and elevation [27]. At low elevations (0 to ~1200 masl), vegetation cover included lowland humid tropical forests, tropical dry forests, and xerophytic and sub-xerophytic scrub. Subandean forests covered the mid-elevations (1300 to 2000–2300 masl), while andean forests covered elevations up to 3300–3500 masl. Higher elevations were dominated by paramo grasslands up to ~4800 masl. Most of the native vegetation, however, has been transformed. Current land use/land cover in the basin is the result of dynamic social and economic processes that have taken place since colonial times [27]. With regards to coffee, the first accounts of coffee plantations in the region date back to the 1850s, but coffee expansion did not take place until the 1940s and 1950s, with the arrival of peasants fleeing from political violence in the central Andes. Forest cover in the SNSM has also been dynamic as a result of multiple processes, including forest clearing by colonos (1950s onwards) and cultivation of illegal crops. Specifically, the international trade of marijuana (Cannabis sativa) during the 1970s brought large-scale deforestation to the SNSM, as it became a prime spot for cultivation, with an estimated 70% of its forests (~150,000 ha) felled between 1975 and 1980 [27]. Forest recovery within the basin seems to have taken place after the mid-1980s due to a lower demand for marijuana in the international markets, and the arrival of illegal groups [28]. The extent of deforestation from that period within the Río Frío basin is uncertain but local inhabitants confirm that large patches of forest in the upper basin were cleared (Meza, pers. comm.). Current land use/land cover includes a mosaic of forests at different successional stages, secondary growth, pasture for cattle, and crops (mostly coffee and subsistence annual crops). Paramo grasslands are found above ~3100 masl [29], while rapidly receding glaciers are found above 4900–5100 masl [30].

2.2. Hydrologic Model

We used the SWAT (soil and water assessment tool) hydrologic model (rev 627), to represent the main hydrologic processes within the watershed. SWAT is a continuous-time, process-based, semi-distributed model that simulates the major water balance components for each subwatershed within a watershed. Further model details can be found in [31]. Since its release in the early 1990s, SWAT has been extensively used to assess multiple aspects of water resource management in river basins, such as non-point source pollution, land cover/land use changes, soil erosion, alternative management practices and climate change ([32] and references therein). However, examples from tropical America are limited (e.g., [33,34]). Our model selection was based on considerations such as its ability to represent the physical processes associated with water movement, support documentation and additional software for model calibration and scenario simulation.

2.2.1. Data Sources and Processing

Input data for SWAT were gathered from a variety of sources (Table 1) and were processed as follows. Digital topographic maps with a scale of 1:25,000 from the National Geographic Institute (IGAC) were used to generate a digital elevation model (DEM). Contour interval was 25 m up to 600 masl, and 50 m above that elevation. The DEM was generated using the GRASS module within QGIS v2.18 [35], with an output spatial resolution of 30 m.
Digital soil mapping units were obtained from the national ecosystems map [29]. Soil profile properties required by SWAT were gathered from the state´s soil survey [26]. Properties not included in the soil survey were calculated as follows. Saturated soil hydraulic conductivity was estimated from soil texture and organic matter content [36]. Bulk density and available water capacity were missing for certain soil units, and were also calculated using the relationships from [36]. Soil erodibility was estimated with the software KUERY v1.4 [37] based on the soil´s texture, organic matter content, and regional climate [38]. For each soil unit, we compared the distribution´s mean and modal values, selecting the higher erodibility of the two values. Soil albedo was estimated from the soil´s color [39]. Most soil mapping units within the watershed were associations (i.e., two or more components at the subgroup level), while some were consociations. For the former, we used the properties of the most extensive unit, while for the latter we used the properties of the dominant component.
Land cover data were generated from a Landsat 8 image from January 2015. The land cover classification was performed using a random forest classification [40]. Prior to classification, the Landsat image was atmospherically corrected using the Cos(t) model [41] available in IDRISI v17.02 (Clark Labs). Subsequently, bands 2 through 7 from the atmospherically corrected image were subsetted to the area of interest. The following ancillary raster data were generated to aid in the classification: (1) normalized difference vegetation index (NDVI), (2) ecosystems from the national ecosystem map [29], (3) topographic data including elevation, slope, aspect, and hillshade derived from the DEM, (4) coffee elevation mask where pixels with elevations between 800 and 1800 masl were assigned a value of 1 and others a value of zero, (5) distance data including distance to the ocean and coastal lagoons, distance to major cities and distance to towns. All rasters were subsetted and aligned to the Landsat subsetted bands, and were generated to match the Landsat spatial resolution. A set of training points was generated for each land cover class, with a mininum of 30 points per class except for aquatic vegetation which was limited to small areas and had only 24 training sites. For each training point, values of all Landsat subsetted bands and ancillary data were extracted, as well as its X and Y coordinates. Spatial data processing was performed in IDRISI and ArcGIS (ESRI). The final point dataset was used as input for the random forest R package [42]. This tool generates a large number of decision trees (i.e., forest) based on random subsets of the training samples. Subsequently, each pixel in need of classification is put down each tree and assigned a class. The final class is selected as the most common, or the one having the most votes. A 10-fold cross validation was performed internally during the generation of the decision trees [43]. Final land cover classes were assigned to the equivalent SWAT land cover code as follows: Subandean and Andean forest (referred to as wet forest from hereafter) to FRST, tropical dry forest to FRSD, coffee to COFF, shrub or secondary growth to RNGB, pasture to PAST, paramo grassland to RNGE, banana crops to BANA, and bare to BARR (codes used from hereafter).
Daily precipitation and temperature data from nearby stations (<15 km), as well as daily discharge data at the watershed outlet, were obtained from the National Environmental Institute (IDEAM, Figure 1, Supplementary Materials S1). Finally, leaf area index (LAI from hereafter) data from MODIS were used during the calibration stage to modify SWAT´s default vegetation growth cycles (see Section 2.2.3 below) [44].

2.2.2. Model Setup

We constructed the SWAT model using the ArcSWAT extension (v2012.10_0.15) within ArcGIS. The watershed was divided into 15 subbasins (Figure 1) spanning areas of 26–4497 ha, using a minimum flow accumulation area of 1500 ha. Subbasins were further divided into hydrologic response units (HRUs), each defined as a particular combination of land cover class, soil type, and slope, using the soil data and Landsat-derived land covers described in Section 2.2.1, and two slope classes separated by a threshold of 36%. After further specifying a minimum threshold area of 1 ha for land cover, soil, and slope, the subbasins were finally divided into a total of 271 HRUs constituting 99.9% of the total watershed area.
We forced SWAT using daily precipitation from 4 climate stations within ~12 km of the watershed, and daily maximum and minimum air temperature from 3 stations within ~15 km of the watershed (Figure 1, Supplementary Materials S1). Statistics from these stations were used in SWAT’s internal WGENX weather generator to estimate daily solar radiation, and relative humidity, needed for potential ET (Priestly–Taylor method selected), and for gap-filling of missing precipitation and temperature data. Additional climate stations were used to assign subbasin-values of precipitation lapse rate (Supplementary Materials S1 and S2) as follows: (1) 500 mm km−1 from 0–600 masl, and (2) 262 mm km−1 from 600–2000 masl. Temperature lapse rate was found to be −6.3 °C km−1 (Supplementary Materials S1 and S2).

2.2.3. Model Calibration and Validation

We performed model calibration/validation and sensitivity/uncertainty analysis using the sequential uncertainty fitting (SUFI-2) algorithm [45,46] available in the SWAT-CUP software package [47]. A key aspect of this algorithm for this study is the prediction of a 95% uncertainty (95PPU) band found through an iterative calibration process. Our process of calibrating and validating models followed these steps (details in Supplementary Materials S2):
  • Step 1, Data quality assurance and control. Daily observed values of stream discharge were screened for quality. Values that were both quality-flagged and considered outliers were removed from the record, and remaining values used to construct monthly discharge for use in simulation and analysis.
  • Step 2, Simulation period definition. We selected calibration and validation periods that included wet, average and dry years, and had a relatively low number of quality flags. Accordingly, we defined 2002–2008 as the calibration period, with a three-year warm-up from 1999–2001, and 1985–1991 as the validation period, with a three-year warm-up from 1982–1984.
  • Step 3, Determination of flow-path goals. To identify aspects of the modeled flow paths that needed improvement (hereafter, “flow-path goals”), we performed an initial comparison of literature-based observations to the following outputs from SWAT in its default configuration: fraction of total runoff as baseflow, calculated from daily discharge using a baseflow filter (https://engineering.purdue.edu/mapserve/WHAT/), ratios of total ET/total precipitation and surface runoff/total streamflow as measured in the central Colombian Andes [48,49,50]. Based on this comparison, we determined that the calibration process should, relative to SWAT defaults, decrease total runoff and increase total baseflow and ET.
  • Step 4, SWAT calibration of LAI. We found that the default configuration of SWAT provided a poor representation of LAI dynamics. We therefore calibrated SWAT-simulated LAI to MODIS-derived, eight-day LAI composites using the SWAT-T module for tropical vegetation growth [51,52]. This SWAT module has been shown to improve the representation of shifts between vegetation dormancy and growth in tropical regions by using soil moisture, rather than day-length, to represent important phenological thresholds. For more information, see Supplementary Materials S3.
  • Step 5, Sensitivity analysis. To identify important model parameters to estimate via calibration, we performed a global sensitivity analysis of monthly simulated streamflow to 16 selected parameters. This was done by ranking trend p-values for each parameter in question, with other parameters varying simultaneously, using 500 simulations in SWAT-CUP. We selected the 10 most influential parameters (listed in Supplementary Materials S4) to estimate via calibration. During the calibration process described next, we also performed local, one-at-time sensitivity analyses to select initial parameter ranges for SUFI-2 calibration that were compatible with our flow-path goals (Step 3, Supplementary Materials S5).
  • Step 6, Calibration. We calibrated monthly streamflow with SUFI-2 using 1000 simulations per iteration (latin hypercube sampling) and the Nash–Sutcliffe (NS) coefficient as the “goal function” for suggested parameter range-centering. We manually limited parameter ranges input to SUFI-2 by considering our flow-path goals (Step 3), uncertainty statistics (p-value and r-value) [45,46], and performance criteria [53]. This was done by repeating Steps 5 and 6 on a trial-and-error basis until acceptable values were attained.
  • Step 7, Validation. Once we obtained satisfactory calibration results, we tested the model against streamflow using the parameter ranges obtained from calibration (steps 4 and 6).

2.3. Scenarios

2.3.1. Synthetic Rainfall Series

We created synthetic rainfall series to represent drought conditions at two different temporal scales (i.e., annual and monthly) in order to study streamflow response to consecutive dry years (months) (Figure 2). We generated the annual drought series for each rainfall climate station by (1) selecting years where annual rainfall was below the 10th percentile, (2) averaging daily rainfall values across those selected years to generate a ‘dry’ rainfall series of daily data for an entire year. In addition, an annual reference series was generated for each station by (1) selecting years with annual rainfall between the 45th to 55th percentiles, (2) averaging daily rainfall values across all years to generate a ‘reference’ series of daily data for an entire year. We then created multiple annual scenarios each of seven-year standard length and containing between one and four consecutive dry years (Table 2). We created the monthly synthetic series for each station by (1) selecting the driest month on record for each month (i.e., driest January, February, etc.), (2) creating a series with increasing number of dry months, up to 12 consecutive dry months (Table 2).

2.3.2. HRU Selection

Prior to running SWAT for each scenario, we selected four HRUs, one for each of the land cover types calibrated with MODIS data. In order to minimize sources of variability other than the land cover type, we selected HRUs that had the same rainfall station, soil type, and slope class. We were able to accomplish this for wet forest, shrub and coffee (HRUs in subbasins 8 and 14, rainfall station Vista Nieves, Figure 1). All dry forest HRUs, however, were located within subbasin 15 (i.e., lower basin), and were associated with a different rainfall station (El Enano, Figure 1). For each scenario, we ran one iteration in SUFI-2 with the same number of simulations (n = 1000) and parameter ranges as we did for the model calibration in order to represent parameter variability. Simulation output included water yield for the selected HRUs, defined as the total amount of water leaving the HRU and entering the main channel at each time step [54].

2.3.3. Statistical Analyses

For each scenario, water yield difference was calculated as the difference between predicted water yield and water yield under the reference scenario. As predictions were based on 1000 simulations, 1000 values of water yield difference were obtained per month, which were summarized through probability density functions (PDFs). For each scenario, monthly PDFs were used for assessing the following aspects:
  • Water yield recovery time after drought termination: For each selected HRU, we plotted PDFs starting in month 1 after drought termination.
  • Water yield decrease magnitude as indicated by the PDF median value for each month/scenario (value of water yield decrease that divides the area under the PDF in half).
  • Probability of water yield decrease measured as the area under the PDF and to the left of zero. It is worth mentioning that high probabilities of water yield reduction do not necessarily imply a severe reduction (i.e., high probability can be associated with low magnitude).

3. Results

3.1. Land Cover Classification

The overall land cover classification accuracy was 96%, and ranged from minimum values of 85% (shrub) and 88% (bare) to > 90% for all other categories. According to our results, the majority of the basin´s area (~80%) is under wet forest (45% of the basin, or 136 km2; Figure 1), shrub (19% or 58 km2), or shade coffee (17% or 51 km2). Dry forest is limited to the lower basin, covering 6% of the total basin area (17 km2), while pasture is distributed in patches along the upper margins of the Rio Frio (10% or 30 km2). Paramo grasslands (2% or 6 km2) are limited to the upper reaches of the basin, above 3100–3200 masl.

3.2. SWAT-T leaf Area Index (LAI) Calibration

Seasonality of MODIS LAI series had a bimodal pattern for the major land covers analyzed (wet forest, coffee, shrub, and dry forest; Figure 3). Wet forest and coffee had peak LAI values (~6.0) in January and June, and minimum values (3.5 and 5.0) in March and October-November, respectively. Shrub had maximum LAI values (6.5) in January and June, and minimum values (5.5) in March and September. Dry forest LAI values peaked at ~4.5 in June and November, and reached minimum values (1.0) in February. There was no significant long-term trend in LAI for coffee and shrub. Wet forest LAI showed a decreasing trend until 2004, while dry forest LAI had a decreasing trend starting in 2004. Calibrated SWAT-T mean LAI values were consistent with LAI ranges for each land cover, but failed to represent LAI´s bimodal seasonality as SWAT-T is capable of modeling only one wet/dry transition within the year. On the other hand, default SWAT LAI mean values failed to represent both LAI seasonality and range, for all major land covers (Figure 3).

3.3. Discharge Calibration and Validation

Our final model had a satisfactory performance rating for both the calibration and validation periods (Table 3, Figure 4), while the major water flux components were consistent with our flow-path goals and reference values (Supplementary Materials S5). Final parameter ranges and best-fit parameter values are included in the Supplementary Materials (Supplementary Materials S6).

3.4. Scenarios

3.4.1. Annual Drought Scenarios

Results for all annual scenarios were analyzed starting at the time of meteorological drought termination (i.e., end of “D” periods in Table 2), regardless of drought length. Annual rainfall for the selected HRUs decreased by ~25% (coffee, shrub, and wet forest) and 40% (dry forest) during dry years compared to the reference scenario (Figure 5a). Rainfall deficits affected water yield for all land covers analyzed, and these effects extended for ~9 months after drought termination as indicated by probabilities higher than 95% of water yield decrease (Figure 6 and Figure 7). Water yield reductions were largest during the first months after the drought had ceased with median values of ~12 mm month−1 for coffee, shrub and wet forest, and 5 mm month−1 for dry forest during the first month after a one-year drought (Figure 6). Water yield reduction was largest during the first month after the drought, reaching median values of ~12 mm month−1 for coffee, shrub, and wet forest, and 5 mm month−1 for dry forest (Figure 6). The effect of drought length was evidenced by longer droughts causing larger water yield reductions (i.e., two-year drought compared to one-year drought). However, there were no significant differences between two-year, three-year and four-year droughts (Supplementary Materials S7).
The probability of water yield reduction was very high (>0.95) for all land covers during the first year after the drought (Figure 7 and Supplementary Materials S7). Although probabilities remained high (~0.8) beyond two years, the magnitude of water yield decrease was very low after month ~12. With regards to differences between land cover types, water yield decrease had a similar behavior in wet forest, coffee, and shrub. For instance, median water yield decrease values in these land covers ranged from ~12 mm month−1 to 4 mm month−1 in the first 12 months after a 1-year drought (Figure 7). By comparison, water yield decrease in dry forest ranged from ~5 mm month−1 to 2 mm month−1 for the same period.

3.4.2. Monthly Drought Scenarios

Results for monthly scenarios were analyzed starting in month 1 after meteorological drought termination (i.e., effects of a one-month drought were assessed starting in month 2, and so on). Rainfall for the selected HRUs decreased by ~69% (coffee, shrub, and wet forest) and 95% (dry forest) during the drought year compared to the reference conditions (Figure 5b). Meteorological droughts that included the initial three months did not have an effect on water yield (Figure 8). Meteorological drought affected water yield for ~2–8 months after drought termination in coffee, wet forest, and shrub, as indicated by probabilities > 95% of water yield decrease (Figure 9 and Supplementary Materials S8). This range was related to land cover type, drought length, and timing. For instance, a 4-month drought affected wet forest water yield for 3 months after the drought had ceased, while a 10-month drought affected wet forest water yield for ~8 months after drought termination (Figure 9a,b). Water yields from dry forest were little affected by drought, except after 5-, 9-, and 10-month droughts. In addition, the effect was short lived, with water yield recovering within one month after the end of the meteorological drought.
The magnitude of water yield reduction was affected by both drought length and timing. For instance, when droughts occurred during months that were usually dry (i.e., months 1, 2, 3, 12, Figure 5b), their effect was either non-significant (e.g., months 1 through 3) to reduced (e.g., month 12) (Supplementary Materials S8). On the other hand, droughts that affected wet months (e.g., month 10, Figure 5b) had the largest effect on water yields. For instance, a 10-month drought resulted in median water yield decrease values between 90, 108, and 127 mm month−1 for shrub, coffee, and wet forest respectively, and 20 mm month−1 for dry forest during the first month after drought termination (Figure 9b). By comparison, a 12-month drought resulted in median water yield decrease values below 50 mm month−1 for all land covers during the first month after the drought had ceased (Figure 9c).
Water yield decrease probability was very high (>0.9) for five-month and longer droughts in all land covers (Figure 9). Dry forest had slightly lower probabilities (>0.7) of water yield reductions for a four-month drought (Figure 9a).

4. Discussion

In this section, we discuss our results within the context of research on streamflow response to drought conditions and the role of different land cover types in moderating such response. Our annual scenarios are illustrative of longer, but less severe meteorological drought conditions, while the monthly scenarios represent shorter, but extremely severe rainfall deficits. We recognize that results from our longer annual scenarios (three- to four-year droughts) may be limited as they did not incorporate long-term drought related vegetation changes such as tree mortality and changes in species composition. We therefore focus our annual scenario discussion on one- and two-year droughts.

4.1. Drought Propagation

Studies dealing with the propagation of meteorological drought show that the following effects take place as it moves through the hydrologic system [10,55]: (1) pooling, which is the combination of meteorological droughts into a longer term hydrological drought; (2) attenuation, or tempering of the drought signal by water storage components; (3) lag, or time between the occurrence of a meteorological drought and soil moisture/hydrological drought; and (4) lengthening, meaning that droughts become longer as they move through the hydrologic system. Catchment characteristics control lag and attenuation, while climate and catchment properties control pooling and lengthening. In our case, we observed lengthening of the drought signal for both annual and monthly scenarios. For the annual scenarios, the hydrological drought extended for ~9 months after the meteorological drought had ended, while for the monthly scenarios, it extended for up to ~8 months. In addition, results from the monthly scenarios indicated a short lag (i.e., ~1 month) between rainfall deficits and reductions in water yield. The magnitude of this response was conditioned by the severity of rainfall deficits (i.e., greater in monthly droughts vs. annual droughts), and timing of the meteorological drought (i.e., larger water yield reductions in droughts that included regularly wet months). Attenuation of the meteorological drought signal was observed for the annual and monthly scenarios. For instance, annual scenarios had maximum monthly rainfall deficits of 80%, while maximum monthly water yield deficits of ~50% for wet forest, coffee, and shrub. Dry forest had more severe rainfall (100%) and water yield deficits (80%). The same effect was observed for monthly scenarios, although the deficits were larger.
Based on the above, our basin´s response to meteorological drought was controlled by climate, specifically by strong rainfall seasonality, and by catchment characteristics that led to fast water yield response. With respect to climate, the long-term climatology of our study region provides an example of the “wet-to-dry season” drought type [55]. These droughts occur in warm regions with strong seasonal variation in precipitation, where most recharge takes place during the wet season. Therefore, rainfall deficits during the wet season affect storage and extend drought conditions into the following dry season, as the probability of recharge during the latter is very low [13]. Slow water pathways (i.e., baseflow, groundwater contribution to discharge) are critical in sustaining streamflow during the dry season under such climates [10]. This is evident in our case, where baseflow contribution to discharge reaches > ~80% during the dry season under normal conditions. Baseflow is also the main contributor to discharge under drought conditions, as faster pathways (i.e., surface runoff and lateral-flow/interflow) are first depleted during the propagation of a meteorological drought [10]. Our monthly scenarios showed that droughts during the first four months did not affect streamflow, as recharge during the previous wet season maintained baseflow contribution to streamflow. Longer droughts did affect streamflow, particularly when they extended to 10 consecutive months. We believe that reduced streamflows during the dry season immediately following a 12-month drought were related to decreased groundwater recharge and resulted from a combination of (1) a delayed start of the wet season and (2) decreased rainfall during the wettest months of the year.
The main catchment characteristics that affect streamflow response to drought are related to its storage capacity, and include soil properties, geology, groundwater system, presence of lakes, glaciers and swamps, vegetation, and topography [10]. In particular, groundwater response time seems to play a critical role on hydrological drought duration and severity [11]. Fast groundwater responding systems tend to exhibit shorter and more intense hydrological droughts, while slow responding systems tend to exhibit longer and more attenuated hydrological droughts [11,13]. Catchment response, however, seems to be controlled by a combination of multiple factors [56]. Based on our catchment characteristics and results, we consider the following aspects to play a significant role in drought response: (1) catchment area and topography, (2) vegetation, and (3) soils and geology.
Catchment size is considered an important aspect of drought response, as large catchments may have areas that respond to drought in different ways, modulating the overall response [10]. On the other hand, steep topography in tropical catchments leads to high rainfall spatial variability [57]. In our case, although the Río Frío catchment is fairly small (~300 km2), it exhibits a heterogeneous response to drought in terms of water yield decrease magnitude, which we attribute to the steep elevation gradient (20 masl to ~4200 masl) and its effect on rainfall. Rainfall records report mean annual totals of 1000 mm at 25 masl to more than 2500 mm at ~2000 masl. The upper basin therefore receives significantly more rainfall than the lower, as evidenced by differences in precipitation between wet forest (higher elevation), shrub and coffee (intermediate), and dry forest (lower, Figure 1 and Figure 5). This spatial rainfall pattern controls the magnitude of drought throughout the basin. For instance, our results showed that water yield decrease for dry forest was minimal compared to other land covers. Since dry forest is located in the lower and drier part of the watershed, reductions in water yield are small as they are limited by low precipitation inputs (Figure 5).

4.2. Land Cover Effects

The distribution of vegetation within our catchment is controlled by spatial rainfall patterns as well as altitudinal temperature changes. Vegetation, in turn, plays a role on meteorological drought intensification and hydrological drought recovery through its effects on land-atmosphere feedbacks [10] and hydrologic regulation. For example, in warm seasonal climates, potential evapotranspiration (PET) may increase as drought develops leading to higher actual ET which further depletes soil moisture and prevents groundwater recharge [10]. In the specific case of forests, their role in water regulation is critical for sustaining streamflow during dry periods. Data from a three-year paired catchment study in the Panama Canal basin showed that a forested catchment sustained higher baseflow during the dry season compared to pasture and mosaic (agriculture, secondary forest, and pasture) catchments, supporting the “sponge-effect hypothesis” [58]. Data from the central Colombian Andes showed that forests at different successional stages provided increased hydrological regulation through lower variability in soil moisture storage, compared to pasture and croplands [50]. Similar findings are reported for other tropical forest sites and are related to increased soil infiltration capacity during early successional stages ([59] and references therein). In addition to water regulation, forests also display higher resilience to drought due to a greater diversity of plant hydraulic strategies to cope with dry conditions [60]. In our case, wet forests cover a significant portion of the watershed (45% or 136 km2) between 800 and 4200 masl, with an average elevation of 2300 masl. Based on the recent land use history of this region, it is likely that wet forests have varying degrees of intervention and successional stages. Data from the central Colombian Andes show that hydrologic regulation recovery (e.g., greater soil moisture storage and decreased overland flow) starts early in the successional process [50], while data from central Mexico indicates that the hydrologic behavior of lower montane forests was restored within 20 years after disturbance [61,62]. Recovery may be hindered by topsoil erosion due to severe soil degradation [63]. Shrub in the Río Frío catchment displayed a drought response similar to wet forest, which is consistent with the mentioned studies. With regards to coffee, all coffee in the SNSM is grown as shade coffee, as opposed to other regions of Colombia, due to the marked dry season of December-March. Research comparing lower montane rainforest forest and coffee agroforestry systems in eastern Mexico showed that forest canopies retained more water due to higher leaf area, canopy projection, presence of epiphytes, vertical stratification, and density of individuals [64]. On the other hand, coffee systems had increased throughfall and stemflow (i.e., more water for infiltration). Studies in the central coffee region of Colombia found slightly higher interception in shade-grown coffee systems compared to secondary subandean forest, while no significant differences in infiltration and overall low surface runoff under both systems [65]. Our results are consistent with the latter, as coffee´s response to drought was similar to that of wet forest, suggesting that they both have comparable hydrologic behavior.
Finally, although we did not explicitly explore the effects of soil type and underlying geology on drought response, we consider both of them important as they partially control the amount of water that infiltrates and eventually reaches the groundwater system. Specifically, soil survey data [26] indicates that soils in the upper basin display greater infiltration rates (i.e., hydrologic group B) than soils in the mid and lower basin (i.e., hydrologic group C).
Based on the above, we consider that shrub, coffee and wet forest in the mid and upper basin are critical for water regulation and provisioning during droughts. In such areas, the interaction between climate (e.g., higher rainfall and lower PET), land cover and soils leads to greater hydrological regulation. Additional analyses on baseflow in selected land covers (not shown) indicates that the effects of a 12-month drought extend for ~8 months after drought termination, suggesting that an average wet season is sufficient to replenish all components of the system (soil and groundwater). We also identified several lines for future model improvement, mainly: (1) better representation of LAI patterns through the incorporation of a bimodal cycle and calibration of satellite data with field measurements, (2) incorporation of long-term effects of drought on vegetation, (3) better depiction of rainfall´s spatial variability considering the complex topography of our catchment, and (4) representation of preferential infiltration and percolation pathways, very likely to play an important role considering the local soil and rock properties.

5. Conclusions

Our study watershed exhibited a fast streamflow response to meteorological drought propagation and recovery. Although drought scenarios resulted in a significant decrease in discharge (e.g., 68% during the first month following a 12-month drought), streamflow recovered within ~9 months after drought termination, with recovery time varying according to the severity and timing of rainfall deficits. We hypothesize that water yield response was controlled by a combination of climatic factors and catchment properties. The influential climatic factors included large rainfall seasonality, with a marked four-month dry season that concentrates only 5% of the total annual rainfall under average conditions, and high spatial rainfall variability related to a steep elevation gradient. The influential catchment properties included catchment size, slope steepness, soil infiltration rates, and density of land covers. For this watershed, the hydrological regulation function of wet forests, shade coffee, and shrub in the mid and upper basin, with high rainfall amounts, is considered critical for sustaining streamflow during and after droughts.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/1/94/s1, Supplementary Materials S1. Relevant characteristics of climate stations used for the SWAT model. Stations listed from low to high elevation above sea level. Supplementary Materials S2. SWAT calibration and validation methods. Supplementary Materials S3. SWAT parameters used for LAI calibration. Parameters LAI_INIT, BIO_INIT and PHU_PLT are required when the land cover status code (IGRO) is set to 1, indicating that vegetation is already growing at the beginning of the simulation. Supplementary Materials S4. Global sensitivity analysis. Sensitive parameters are indicated by a high t-statistic value (in absolute terms) and a low p-value. Parameters are listed from high to low sensitivity. CANMX values for forest and coffee were set to reference values prior to the sensitivity analysis. Supplementary Materials S5. Comparison of SWAT average water flux components prior and after calibration with reference values. Calibration period from 2002 to 2008. Supplementary Materials S6. Parameter ranges and best-fit parameter value for the selected SWAT model. Scaling type: v (absolute) indicates that the parameter is replaced by the given value, r (relative) indicates that the parameter is multiplied by [1 + (given value)]. The latter preserves the parameter´s spatial variability. Supplementary Materials S7. Effect of (a) one-year, (b) two-year, (c) three-year, and (d) four-year meteorological drought on water yield (WY) (mm month-1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (minimum value, dashed line) of water yield decrease from month 1 through month 36 after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from higher (red) to lower (blue) probability. Supplementary Materials S8. Effect of monthly meteorological droughts on water yield (WY) (mm month-1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (dashed line) water yield decrease for subsequent months after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from high (red) to low (blue) probability.

Author Contributions

Conceptualization, N.H., A.C.-M., B.W., J.E., and J.C.R.; Formal analysis, N.H., A.C.-M., S.M.J., B.W., M.M., and R.D.; Funding acquisition, N.H., B.W., J.E., and M.I.V.; Methodology, N.H., A.C.-M., S.M.J., B.W., S.V., M.M., and R.D.; Project administration, N.H., J.E., and M.I.V.; Writing–original draft, N.H., S.M.J., and J.E.; Writing–review and editing, N.H., A.C.-M., S.M.J., B.W., M.M., J.E., and J.C.R.

Funding

This research was partially funded by the Inter American Institute for Global Change Research (IAI) CRN 798 3038 to establish the SAFER network and U.S. National Science Foundation Award #1336839. NH was partially supported by the Fulbright Visiting Scholar Program from the Fulbright Commission in Colombia while on research leave from Universidad del Norte. On-site resources for NH while a visiting scholar were provided by the Department of Geography at the University of Vermont. J.E. and N.H. were partially funded by The Canadian Queen Elizabeth II Diamond Jubilee Scholarships (QES), a partnership among Universities in Canada, the Rideau Hall Foundation (RHF), Community Foundations of Canada (CFC). The QES-AS is made possible with financial support from IDRC and SSHRC.

Acknowledgments

We gratefully acknowledge the support provided in the field by Hector Meza, Doña Mirosalba Meza and Wilington Meza. The MODIS LAI data were retrieved from the NASA Earthdata Search Tool, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://search.earthdata.nasa.gov/).We also thank the editor and anonymous reviewers for their comments which helped to improve our manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Steffen, W.; Persson, A.; Deutsch, L.; Zalasiewicz, J.; Williams, M.; Richardson, K.; Crumley, C.; Crutzen, P.; Folke, C.; Gordon, L.; et al. The anthropocene: From global change to planetary stewardship. Ambio 2011, 40, 739–761. [Google Scholar] [CrossRef] [PubMed]
  2. IPCC. Climate Change 2007: The Physical Science Basis; Cambridge University Press: Cambridge, UK, 2007; p. 987. [Google Scholar]
  3. FAO. AQUASTAT. Available online: http://www.fao.org/nr/water/aquastat/countries_regions/americas/indexesp3.stm (accessed on 3 October 2018).
  4. Nurse, L.A.; Sem, G.; Hay, J.E.; Suarez, A.G.; Wong, P.P.; Briguglio, L.; Ragoonaden, S. Small island states. In Climate Change 2001: Impacts, Adaptation, and Vulnerability; McCarthy, J.J., Canziani, O.F., Leary, N.A., Dokken, D.J., White, K.S., Eds.; Cambridge University Press: Cambridge, UK, 2001; pp. 843–875. [Google Scholar]
  5. Neelin, J.D.; Münnich, M.; Su, H.; Meyerson, J.E.; Holloway, C.E. Tropical drying trends in global warming models and observations. Proc. Natl. Acad. Sci. USA 2006, 103, 6110–6115. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Campbell, J.D.; Taylor, M.A.; Stephenson, T.S.; Watson, R.A.; Whyte, F.S. Future climate of the Caribbean from a regional climate model. Int. J. Climatol. 2011, 31, 1866–1878. [Google Scholar] [CrossRef]
  7. Singh, B. Climate changes in the greater and southern Caribbean. Int. J. Climatol. 1997, 17, 1093–1114. [Google Scholar] [CrossRef] [Green Version]
  8. Angeles, M.E.; Gonzalez, J.E.; Erickson, D.J.I.; Hernandez, J.L. Predictions of future climate change in the Caribbean region using global general circulation models. Int. J. Climatol. 2006, 27, 555–569. [Google Scholar] [CrossRef]
  9. Herrera, D.A.; Ault, T.R.; Fasullo, J.T.; Coats, S.J.; Carrillo, C.M.; Cook, B.I.; Williams, A.P. Exacerbation of the 2013–2016 Pan-Caribbean Drought by Anthropogenic Warming. Geophys. Res. Lett. 2018, 45, 10619–10626. [Google Scholar] [CrossRef]
  10. Van Loon, A.F. Hydrological drought explained. Wiley Interdiscip. Rev. Water 2015, 2, 359–392. [Google Scholar] [CrossRef] [Green Version]
  11. Van Lanen, H.A.J.; Wanders, N.; Tallaksen, L.M.; Van Loon, A.F. Hydrological drought across the world: Impact of climate and physical catchment structure. Hydrol. Earth Syst. Sci. 2013, 17, 1715–1732. [Google Scholar] [CrossRef]
  12. Van Loon, A.F.; Gleeson, T.; Clark, J.; Van Dijk, A.I.J.M.; Stahl, K.; Hannaford, J.; Di Baldassarre, G.; Teuling, A.J.; Tallaksen, L.M.; Uijlenhoet, R.; et al. Drought in the Anthropocene. Nat. Geosci. 2016, 9, 89–91. [Google Scholar] [CrossRef] [Green Version]
  13. Van Loon, A.F.; Tijdeman, E.; Wanders, N.; Van Lanen, H.A.J.; Teuling, A.J.; Uijlenhoet, R. How climate seasonality modifies drought duration and deficit. J. Geophys. Res. Atmos. 2014, 119, 4640–4656. [Google Scholar] [CrossRef] [Green Version]
  14. Aldana-Domínguez, J.; Montes, C.; Martínez, M.; Medina, N.; Hahn, J.; Duque, M. Biodiversity and ccosystem services knowledge in the Colombian Caribbean: Progress and challenges. Trop. Conserv. Sci. 2017, 10, 1940082917714229. [Google Scholar] [CrossRef]
  15. Tribin, M.C.D.G.; Rodriguez, N.G.E.; Valderrama, M. Sierra Nevada de Santa Marta: A Pioneer Experience of a Shared and Coordinated Management of a Bioregion; UNESCO: Paris, France, 1999; p. 40. [Google Scholar]
  16. Ruiz, F.; Rodriguez, A.; Armenta, G.; Grajales, F. Informe de Escenarios de Cambio Climático para Temperatura y Precipitación en Colombia; IDEAM: Bogotá, Colombia, 2013; p. 81.
  17. IDEAM; PNUD; MADS; DNP; CANCILLERÍA. Escenarios de Cambio Climático para Precipitación y Temperatura para Colombia 2011–2100: Herramientas Científicas para la Toma de Decisions—Estudio Técnico Completo: Tercera Comunicación Nacional de Cambio Climático; IDEAM: Bogotá, Colombia, 2015; p. 278.
  18. Ochoa, A.; Poveda, G. Distribución espacial de señales de cambio climático en Colombia. In Proceedings of the XXIII Congreso Latinoamericano de Hidráulica, Cartagena, Colombia, 2–6 September 2008. [Google Scholar]
  19. Cantor, D.C. Evaluación y Análisis Espacio Temporal de Tendencias de Largo plazo en la Hidroclimatología Colombiana; Universidad Nacional de Colombia: Medellín, Colombia, 2011. [Google Scholar]
  20. Pabón, J.D. Cambio climático en Colombia: Tendencias en la segunda mitad del siglo XX y escenarios posibles para el siglo XXI. Rev. Acad. Colomb. Cienc. Exactas Fís. Nat. 2012, 36, 261–278. [Google Scholar]
  21. Carmona, A.M.; Poveda, G. Detection of long-term trends in monthly hydro-climatic series of Colombia through Empirical Mode Decomposition. Clim. Chang. 2014, 123, 301–313. [Google Scholar] [CrossRef]
  22. Hurtado, A.F.; Mesa, O.J. Cambio climático y variabilidad espacio-temporal de la precipitación en Colombia. Rev. EIA 2015, 12, 131–150. [Google Scholar] [CrossRef]
  23. Pierini, J.O.; Restrepo, J.C.; Aguirre, J.; Bustamante, A.M.; Velásquez, G.J. Changes in seasonal streamflow extremes experienced in rivers of Northwestern South America (Colombia). Acta Geophys. 2017, 65, 377–394. [Google Scholar] [CrossRef] [Green Version]
  24. UNGRD. Fenómeno El Niño. Análisis Comparativo 1997–1998//2014–2016; Unidad Nacional para la Gestión del Riesgo de Desastres: Bogotá, Colombia, 2016; p. 143.
  25. Restrepo, J.C.; Ortiz, J.C.; Pierini, J.; Schrottke, K.; Maza, M.; Otero, L.; Aguirre, J. Freshwater discharge into the Caribbean Sea from the rivers of Northwestern South America (Colombia): Magnitude, variability and recent changes. J. Hydrol. 2014, 509, 266–281. [Google Scholar] [CrossRef]
  26. International Global Atmospheric Chemistry (IGAC). Estudio General de suelos y Zonificación del Tierras: Departamento del Magdalena, Escala 1:100000; IGAC: Bogotá, Colombia, 2009.
  27. Fundación Pro-Sierra Nevada de Santa Marta. Plan de Desarrollo Sostenible de la Sierra Nevada de Santa Marta; Fundación Pro-Sierra Nevada de Santa Marta: Santa Marta, Colombia, 1997; p. 228. [Google Scholar]
  28. Uribe, E. Natural Resource Conservation and Management in the Sierra Nevada of Santa Marta: Case Study; Universidad de los Andes: Bogotá, Colombia, 2005. [Google Scholar]
  29. MADS; IDEAM; IAvH; INVEMAR; IIAP; SINCHI; PNN; IGAC. Mapa de Ecosistemas Continentales, Costeros y Marinos de Colombia Version 1.0; IDEAM: Bogotá, Colombia, 2015.
  30. IDEAM (Instituto de Hidrología, Meteorología y Estudios Ambientales). Informe del Estado de los Glaciares Colombianos; Instituto de Hidrología, Meteorología y Estudios Ambientales: Bogotá, Colombia, 2018.
  31. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. Soil and Water Assessment Tool Theoretical Documentation, Version 2009; Texas Water Resources Institute: College Station, TX, USA, 2011; p. 618. [Google Scholar]
  32. Arnold, J.G.; Moriasi, D.N.; Gassman, P.W.; Abbaspour, K.C.; White, M.J.; Srinivasan, R.; Santhi, C.; Harmel, R.D.; van Griensven, A.; Van Liew, M.W.; et al. SWAT: Model use, calibration, and validation. Trans. ASABE 2012, 55, 1491–1508. [Google Scholar] [CrossRef]
  33. Abe, C.A.; Lobo, F.L.; Dibike, Y.B.; Costa, M.P.F.; Dos Santos, V.; Novo, E.M.L.M. Modelling the effects of historical and future land cover changes on the hydrology of an Amazonian basin. Water 2018, 10, 932. [Google Scholar] [CrossRef]
  34. Montecelos-Zamora, Y.; Cavazos, T.; Kretzschmar, T.; Vivoni, E.R.; Corzo, G.; Molina-Navarro, E. Hydrological modeling of climate change impacts in a tropical river basin: A case study of the Cauto river, Cuba. Water 2018, 10, 1135. [Google Scholar] [CrossRef]
  35. QGIS Development Team. QGIS Geographic Information System; Open Source Geospatial Foundation Project, 2.18.6; QGIS Development Team, 2017. [Google Scholar]
  36. Saxton, K.E.; Rawls, W.J. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil Sci. Soc. Am. J. 2006, 70, 1569–1578. [Google Scholar] [CrossRef]
  37. Borselli, L. KUERY. Global Erodibility Database Query, 1.4; San Luis Potosí, Mexico, 2013. [Google Scholar]
  38. Borselli, L.; Torri, D.; Poesen, J.; Iaquinta, P. A robust algorithm for estimating soil erodibility in different climates. Catena 2012, 97, 85–94. [Google Scholar] [CrossRef]
  39. Post, D.F.; Fimbres, A.; Matthias, A.D.; Sano, E.E.; Accioly, L.; Batchily, A.K.; Ferreira, L.G. Predicting soil albedo from soil color and spectral reflectance data. Soil Sci. Soc. Am. J. 2000, 64, 1027–1034. [Google Scholar] [CrossRef]
  40. Breiman, L. Random forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef]
  41. Chavez, P.S. Image-based atmospheric corrections—Revisited and improved. Photogramm. Eng. Remote Sens. 1996, 62, 1025–1036. [Google Scholar]
  42. R Core Team R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017.
  43. Breiman, L.; Cutler, A. Random Forests. Available online: https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm (accessed on 10 February 2016).
  44. Myneni, R.; Knyazikhin, Y.; Park, T. MOD15A2H MODIS/Terra Leaf Area Index/FPAR 8-Day L4 Global 500m SIN Grid V006; NASA EOSDIS Land Processes DAAC: Sioux Falls, SD, USA, 2015.
  45. Abbaspour, K.C.; Johnson, C.A.; van Genuchten, M.T. Estimating Uncertain Flow and Transport Parameters Using a Sequential Uncertainty Fitting Procedure. Vadose Zone J. 2004, 3, 1340–1352. [Google Scholar] [CrossRef]
  46. Abbaspour, K.C.; Yang, J.; Maximov, I.; Siber, R.; Bogner, K.; Mieleitner, J.; Zobrist, J.; Srinivasan, R. Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT. J. Hydrol. 2007, 333, 413–430. [Google Scholar] [CrossRef]
  47. Abbaspour, K.C.; Srinivasan, R. SWAT-CUP, 5.1.6.2; EAWAG: Dubendorf, Switzerland, 2013. [Google Scholar]
  48. Suárez de Castro, F.; Rodríguez Grandas, A. Investigaciones Sobre la Erosión y la Conservación de los Suelos en Colobia; Federación Nacional de Cafeteros: Bogotá, Colombia, 1962; p. 471. [Google Scholar]
  49. Jaramillo-Robledo, A. El balance hídrico. In Clima andino y café en Colombia; Cenicafé: Chinchiná, Colombia, 2005; pp. 107–123. [Google Scholar]
  50. García-Leoz, V.; Villegas, J.C.; Suescún, D.; Flórez, C.P.; Merino-Martín, L.; Betancur, T.; León, J.D. Land cover effects on water balance partitioning in the Colombian Andes: Improved water availability in early stages of natural vegetation recovery. Reg. Environ. Chang. 2018, 18, 1117–1129. [Google Scholar] [CrossRef]
  51. Strauch, M.; Volk, M. SWAT plant growth modification for improved modeling of perennial vegetation in the tropics. Ecol. Model. 2013, 269, 98–112. [Google Scholar] [CrossRef]
  52. Alemayehu, T.; van Griensven, A.; Woldegiorgis, B.T.; Bauwens, W. An improved SWAT vegetation growth module and its evaluation for four tropical ecosystems. Hydrol. Earth Syst. Sci. 2017, 21, 4449–4467. [Google Scholar] [CrossRef] [Green Version]
  53. Moriasi, D.N.; Arnold, J.G.; Van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  54. Arnold, J.G.; Kiniry, J.R.; Srinivasan, R.; Williams, J.R.; Haney, E.B.; Neitsch, S.L. SWAT Input/Output Documentation, Version 2012; Texas Water Resources Institute: College Station, TX, USA, 2012; p. 654. [Google Scholar]
  55. Van Loon, A.F.; Van Lanen, H.A.J. A process-based typology of hydrological drought. Hydrol. Earth Syst. Sci. 2012, 16, 1915–1946. [Google Scholar] [CrossRef] [Green Version]
  56. Van Loon, A.F.; Laaha, G. Hydrological drought severity explained by climate and catchment characteristics. J. Hydrol. 2015, 526, 3–14. [Google Scholar] [CrossRef] [Green Version]
  57. Wohl, E.; Barros, A.; Brunsell, N.; Chappell, N.A.; Coe, M.; Giambelluca, T.; Goldsmith, S.; Harmon, R.; Hendrickx, J.M.H.; Juvik, J.; et al. The hydrology of the humid tropics. Nat. Clim. Chang. 2012, 2, 655. [Google Scholar] [CrossRef]
  58. Ogden, F.L.; Crouch, T.D.; Stallard, R.F.; Hall, J.S. Effect of land cover and use on dry season river runoff, runoff efficiency, and peak storm runoff in the seasonal tropics of Central Panama. Water Resour. Res. 2013, 49, 8443–8462. [Google Scholar] [CrossRef] [Green Version]
  59. Bruijnzeel, L.A. Hydrological functions of tropical forests: Not seeing the soil for the trees? Agric. Ecosyst. Environ. 2004, 104, 185–228. [Google Scholar] [CrossRef]
  60. Anderegg, W.R.L.; Konings, A.G.; Trugman, A.T.; Yu, K.; Bowling, D.R.; Gabbitas, R.; Karp, D.S.; Pacala, S.; Sperry, J.S.; Sulman, B.N.; et al. Hydraulic diversity of forests regulates ecosystem resilience during drought. Nature 2018, 561, 538–541. [Google Scholar] [CrossRef] [PubMed]
  61. Muñoz-Villers, L.E.; Holwerda, F.; Gómez-Cárdenas, M.; Equihua, M.; Asbjornsen, H.; Bruijnzeel, L.A.; Marín-Castro, B.E.; Tobón, C. Water balances of old-growth and regenerating montane cloud forests in central Veracruz, Mexico. J. Hydrol. 2012, 462–463, 53–66. [Google Scholar] [CrossRef]
  62. Muñoz-Villers, L.E.; McDonnell, J.J. Land use change effects on runoff generation in a humid tropical montane cloud forest region. Hydrol. Earth Syst. Sci. 2013, 17, 3543–3560. [Google Scholar] [CrossRef] [Green Version]
  63. Cavelier, J.; Aide, T.M.; Santos, C.; Eusse, A.M.; Dupuy, J.M. The Savannization of Moist Forests in the Sierra Nevada de Santa Marta, Colombia. J. Biogeogr. 1998, 25, 901–912. [Google Scholar] [CrossRef]
  64. Ponette-González, A.G.; Weathers, K.C.; Curran, L.M. Water inputs across a tropical montane landscape in Veracruz, Mexico: Synergistic effects of land cover, rain and fog seasonality, and interannual precipitation variability. Glob. Chang. Biol. 2010, 16, 946–963. [Google Scholar] [CrossRef]
  65. Jaramillo-Robledo, A.; Cháves-Córdoba, B. Aspectos hidrológicos en un bosque y en plantaciones de café (Coffea arabica L.) al sol y bajo sombra. Cenicafé 1999, 50, 97–105. [Google Scholar]
Figure 1. (a) Regional location of the study area showing the Santa Marta Massif (SNSM) in northern Colombia; (b) Río Frío watershed on the southwestern flank of the SNSM, climate stations used for daily temperature, precipitation and discharge records, and major nearby city (Santa Marta). (c) Río Frío watershed land cover and subbasins. Land cover codes are: FRST (wet forest), FRSD (tropical dry forest), RNGB (shrub or secondary growth), COFF (coffee), RNGE (paramo grassland), PAST (pasture), BANA (banana crops), and BARR (bare).
Figure 1. (a) Regional location of the study area showing the Santa Marta Massif (SNSM) in northern Colombia; (b) Río Frío watershed on the southwestern flank of the SNSM, climate stations used for daily temperature, precipitation and discharge records, and major nearby city (Santa Marta). (c) Río Frío watershed land cover and subbasins. Land cover codes are: FRST (wet forest), FRSD (tropical dry forest), RNGB (shrub or secondary growth), COFF (coffee), RNGE (paramo grassland), PAST (pasture), BANA (banana crops), and BARR (bare).
Water 11 00094 g001
Figure 2. Summary of steps followed, including model setup, calibration and validation, generation of synthetic rainfall series, and simulation of drought conditions. Drought series were generated at two time scales, annual (i.e., dry years selected as those with annual precipitation below the 10th percentile for each station), and monthly (i.e., dry months selected as the driest respective month in the historic record for each station). Reference series refer to years with average rainfall conditions (45th to 55th percentile).
Figure 2. Summary of steps followed, including model setup, calibration and validation, generation of synthetic rainfall series, and simulation of drought conditions. Drought series were generated at two time scales, annual (i.e., dry years selected as those with annual precipitation below the 10th percentile for each station), and monthly (i.e., dry months selected as the driest respective month in the historic record for each station). Reference series refer to years with average rainfall conditions (45th to 55th percentile).
Water 11 00094 g002
Figure 3. Calibrated LAI (area weighted HRU mean) compared to bfast-filtered MODIS values and default SWAT values for (a) wet forest, (b) coffee, (c) shrub, and (d) dry forest. Calibrated SWAT results were obtained using the Alemayehu et al. (2017) modified vegetation growth module, while default SWAT results were obtained using the default vegetation growth algorithm based on day length. Bfast-filtered MODIS values include the seasonal and trend components of the LAI time series, and exclude the random component.
Figure 3. Calibrated LAI (area weighted HRU mean) compared to bfast-filtered MODIS values and default SWAT values for (a) wet forest, (b) coffee, (c) shrub, and (d) dry forest. Calibrated SWAT results were obtained using the Alemayehu et al. (2017) modified vegetation growth module, while default SWAT results were obtained using the default vegetation growth algorithm based on day length. Bfast-filtered MODIS values include the seasonal and trend components of the LAI time series, and exclude the random component.
Water 11 00094 g003
Figure 4. Total monthly rainfall and average monthly observed and simulated stream discharge for the calibration (2002–2008) and validation periods (1985–1991). Mean monthly discharge for October 2004 is missing. Total monthly rainfall from SWAT output, calculated from records from four climate stations (Supplementary Materials S1).
Figure 4. Total monthly rainfall and average monthly observed and simulated stream discharge for the calibration (2002–2008) and validation periods (1985–1991). Mean monthly discharge for October 2004 is missing. Total monthly rainfall from SWAT output, calculated from records from four climate stations (Supplementary Materials S1).
Water 11 00094 g004
Figure 5. Monthly rainfall for selected HRUs under reference and dry conditions for (a) annual and (b) monthly scenarios. Numbers within each figure are annual rainfall totals for reference (top) and dry (bottom) conditions. Selected coffee, shrub, and wet forest HRUs are associated with the same rainfall station (Vista Nieves, Figure 1); differences in their rainfall totals are due to differences in elevation between their respective subbasins. Dry forest is associated with a different rainfall station (El Enano, Figure 1).
Figure 5. Monthly rainfall for selected HRUs under reference and dry conditions for (a) annual and (b) monthly scenarios. Numbers within each figure are annual rainfall totals for reference (top) and dry (bottom) conditions. Selected coffee, shrub, and wet forest HRUs are associated with the same rainfall station (Vista Nieves, Figure 1); differences in their rainfall totals are due to differences in elevation between their respective subbasins. Dry forest is associated with a different rainfall station (El Enano, Figure 1).
Water 11 00094 g005
Figure 6. Effect of annual drought on water yield (WY) (mm month−1) for the major land cover types analyzed. Each shaded curve shows a probability density function (PDF) representing the 95% prediction uncertainty (95PPU) of water yield difference values (drought minus reference scenario) for a particular annual scenario (i.e., one-year, two-year, three-year, and four-year drought). The vertical dashed line at zero represents no change relative to the reference scenario, while negative values in the x-axis represent a decrease in water yield for the drought scenario relative to the reference scenario. The y-axis represents the number of months after the termination of the meteorological drought.
Figure 6. Effect of annual drought on water yield (WY) (mm month−1) for the major land cover types analyzed. Each shaded curve shows a probability density function (PDF) representing the 95% prediction uncertainty (95PPU) of water yield difference values (drought minus reference scenario) for a particular annual scenario (i.e., one-year, two-year, three-year, and four-year drought). The vertical dashed line at zero represents no change relative to the reference scenario, while negative values in the x-axis represent a decrease in water yield for the drought scenario relative to the reference scenario. The y-axis represents the number of months after the termination of the meteorological drought.
Water 11 00094 g006
Figure 7. Effect of (a) one-year and (b) two-year meteorological drought on water yield (WY) (mm month−1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (minimum value, dashed line) of water yield decrease from month 1 through month 36 after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from higher (red) to lower (blue) probability. Figure for all annual scenarios in Supplementary Materials S7.
Figure 7. Effect of (a) one-year and (b) two-year meteorological drought on water yield (WY) (mm month−1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (minimum value, dashed line) of water yield decrease from month 1 through month 36 after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from higher (red) to lower (blue) probability. Figure for all annual scenarios in Supplementary Materials S7.
Water 11 00094 g007
Figure 8. Effect of monthly drought on water yield (WY) (mm month−1) for the major land cover types analyzed. Each shaded curve shows a probability density function (PDF) representing the 95% prediction uncertainty (95PPU) of water yield difference values (drought minus reference scenario) for a particular monthly scenario (i.e., one-month drought, two-month drought, etc.). The vertical dashed line at zero represents no change relative to the reference scenario, while negative values in the x-axis represent a decrease in water yield for the drought scenario relative to the reference scenario. The stepladder represents increasing number of dry months, up to 12 consecutive dry months. For example, at step 7 there are six PDFs that represent the effects of a one-month drought, two-month drought, etc. up to a six-month drought. The effect of a specific drought can be assessed by following its PDF in subsequent months. For example, water yield decrease from a four-month drought in coffee can be assessed by following the blue arrow in the left panel (i.e., the effect of a four-month drought is first seen in month 5). Similarly, the effect of a six-month drought in coffee can be assessed by following the orange arrow in the left panel.
Figure 8. Effect of monthly drought on water yield (WY) (mm month−1) for the major land cover types analyzed. Each shaded curve shows a probability density function (PDF) representing the 95% prediction uncertainty (95PPU) of water yield difference values (drought minus reference scenario) for a particular monthly scenario (i.e., one-month drought, two-month drought, etc.). The vertical dashed line at zero represents no change relative to the reference scenario, while negative values in the x-axis represent a decrease in water yield for the drought scenario relative to the reference scenario. The stepladder represents increasing number of dry months, up to 12 consecutive dry months. For example, at step 7 there are six PDFs that represent the effects of a one-month drought, two-month drought, etc. up to a six-month drought. The effect of a specific drought can be assessed by following its PDF in subsequent months. For example, water yield decrease from a four-month drought in coffee can be assessed by following the blue arrow in the left panel (i.e., the effect of a four-month drought is first seen in month 5). Similarly, the effect of a six-month drought in coffee can be assessed by following the orange arrow in the left panel.
Water 11 00094 g008
Figure 9. Effect of a (a) 4-month, (b) 10-month, and (c) 12-month meteorological drought on water yield (WY) (mm month−1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (dashed line) water yield decrease for subsequent months after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from high (red) to low (blue) probability. Figure for all monthly scenarios in Supplementary Materials S8.
Figure 9. Effect of a (a) 4-month, (b) 10-month, and (c) 12-month meteorological drought on water yield (WY) (mm month−1) for selected HRUs representative of the study area. For each figure, the left panel shows the median (continuous line) and 95% probability (dashed line) water yield decrease for subsequent months after drought termination. The vertical line at zero represents no change relative to the reference scenario. Water yield decrease values to the right of the 95% probability line are unlikely. The right panel shows probabilities of water yield decrease with colors scaled from high (red) to low (blue) probability. Figure for all monthly scenarios in Supplementary Materials S8.
Water 11 00094 g009
Table 1. SWAT input data sources and characteristics. For details on rainfall and temperature stations refer to Supplementary Materials S1.
Table 1. SWAT input data sources and characteristics. For details on rainfall and temperature stations refer to Supplementary Materials S1.
DataSource 1Relevant Characteristics
ElevationIGAC1:25,000 scale
Soil mapping units[29]1:100,000 scale
Soil properties[26]Soil profiles from state´s soil survey
Land coverLandsat 8 image from 11 January 201530 m spatial resolution
Leaf area indexMODIS LAI products
MOD15A2Hv006 (January 2002–June 2002)
MCD15A2Hv006 (July 2002–December 2008) [44]
500 m spatial resolution, 8-day composites
PrecipitationIDEAMTotal daily precipitation
Daily maximum and minimum temperatureIDEAMMaximum and minimum daily temperature
DischargeIDEAMAverage daily discharge
Period 1968–2015
1 IGAC: National Geographic Institute (Instituto Geográfico Agustín Codazzi), IDEAM: National Environmental Institute.
Table 2. Annual and monthly drought scenarios used to simulate streamflow response in SWAT.
Table 2. Annual and monthly drought scenarios used to simulate streamflow response in SWAT.
Annual scenariosRainfall year 1
1234567
ReferenceNNNNNNN
Scen1DNNNNNN
Scen2DDNNNNN
Scen3DDDNNNN
Scen4DDDDNNN
Monthly scenariosRainfall month 1
12345678910111213 to 84
ReferenceNNNNNNNNNNNNN
Scen1DNNNNNNNNNNNN
Scen2DDNNNNNNNNNNN
Scen3…DDDNNNNNNNNNN
Scen12DDDDDDDDDDDDN
1 N: ‘normal’, daily series for each station calculated as the daily average for years with annual rainfall between the 45th and 55th percentiles. D: ‘dry’, for annual scenarios, daily series were calculated as the daily average from years with annual rainfall below the 10th percentile. For monthly scenarios, daily series were obtained from the driest month (January, February, etc.) in each station´s record.
Table 3. Selected SWAT model performance statistics for the calibration and validation periods. P- and r-factors evaluate model uncertainty while the Nash–Sutcliffe coefficient (NS), percent bias (PBIAS) and “root mean square error-observations standard deviation ratio” (RSR) assess the performance of the best-fit simulation.
Table 3. Selected SWAT model performance statistics for the calibration and validation periods. P- and r-factors evaluate model uncertainty while the Nash–Sutcliffe coefficient (NS), percent bias (PBIAS) and “root mean square error-observations standard deviation ratio” (RSR) assess the performance of the best-fit simulation.
Modeling PeriodEvaluation Statistics for Model UncertaintyEvaluation Statistics for Best-Fit Simulation 1Performance Rating 2
p-factorr-factorNS 1PBIASRSRR2
Calibration (2002–2008)0.700.570.790.20.460.79Very good
Validation (1985–1991)0.710.600.72−3.40.530.73Good
1 NS: Nash–Sutcliffe efficiency (selected objective function), PBIAS: Percent bias, RSR: RMSE-observations standard deviation ratio, R2: Coefficient of determination. 2 Based on NS, PBIAS and RSR, after [53].

Share and Cite

MDPI and ACS Style

Hoyos, N.; Correa-Metrio, A.; Jepsen, S.M.; Wemple, B.; Valencia, S.; Marsik, M.; Doria, R.; Escobar, J.; Restrepo, J.C.; Velez, M.I. Modeling Streamflow Response to Persistent Drought in a Coastal Tropical Mountainous Watershed, Sierra Nevada De Santa Marta, Colombia. Water 2019, 11, 94. https://doi.org/10.3390/w11010094

AMA Style

Hoyos N, Correa-Metrio A, Jepsen SM, Wemple B, Valencia S, Marsik M, Doria R, Escobar J, Restrepo JC, Velez MI. Modeling Streamflow Response to Persistent Drought in a Coastal Tropical Mountainous Watershed, Sierra Nevada De Santa Marta, Colombia. Water. 2019; 11(1):94. https://doi.org/10.3390/w11010094

Chicago/Turabian Style

Hoyos, Natalia, Alexander Correa-Metrio, Steven M. Jepsen, Beverley Wemple, Santiago Valencia, Matthew Marsik, Rubén Doria, Jaime Escobar, Juan C. Restrepo, and Maria I. Velez. 2019. "Modeling Streamflow Response to Persistent Drought in a Coastal Tropical Mountainous Watershed, Sierra Nevada De Santa Marta, Colombia" Water 11, no. 1: 94. https://doi.org/10.3390/w11010094

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop