Modeling streamflow response to persistent drought in a coastal tropical mountainous watershed, Sierra Nevada de Santa Marta, Colombia

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.


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, 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.
ater 2018, 10, x FOR PEER REVIEW 3 of 2 and covers. This paper is structured as follows. First, we present a description of our study area ydrologic model setup and definition of drought scenarios. We then present the results of the ydrologic model calibration and validation, and the effects of meteorological drought on water yield or selected land covers in the watershed. We then discuss our findings and finish with some oncluding remarks.

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

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. 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).

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 km 2 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  show minimum mean monthly discharge values from January through April (~5.6 m 3 s −1 ), and maximum values from August through November (~21.4 m 3 s −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 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).

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 km 2 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 m 3 s −1 ), and maximum values from August through November (~21.4 m 3 s −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].

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.

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].

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).

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).

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). Table 2. Annual and monthly drought scenarios used to simulate streamflow response in SWAT.

Annual scenarios
Rainfall year 1

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].

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).

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 km 2 ; Figure 1), shrub (19% or 58 km 2 ), or shade coffee (17% or 51 km 2 ). Dry forest is limited to the lower basin, covering 6% of the total basin area (17 km 2 ), while pasture is distributed in patches along the upper margins of the Rio Frio (10% or 30 km 2 ). Paramo grasslands (2% or 6 km 2 ) are limited to the upper reaches of the basin, above 3100-3200 masl.

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

Discharge Calibration and Validation
Our final model had a satisfactory performance rating for both the calibration and validation periods ( Table 3 Table 3. Selected SWAT model performance statistics for the calibration and validation periods. Pand 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.

Discharge Calibration and Validation
Our final model had a satisfactory performance rating for both the calibration and validation periods ( Table 3 Table 3. Selected SWAT model performance statistics for the calibration and validation periods. Pand 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.

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 (Figures 6 and 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 oneyear 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.

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 (Figures 6 and 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.  Figure1); 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 2018, 10, x FOR PEER REVIEW 13 of 23 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.
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.

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).
Water 2018, 10, x FOR PEER REVIEW 15 of 23 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 fourmonth 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.

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.

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

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.

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 km 2 ), 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, Figures 1 and 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).

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 km 2 ) 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.

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. 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.