Objective assessment of ecosystem hydrological services in tropical areas : A Colombian experience in arid and semi-arid zones

This study presents a methodology to address the challenge of objectively demonstrating Ecosystem Hydrological Services (EHS). A case study is used in the region of the La Guajira Peninsula (Colombia), with a focus on the EHS of water flow regulation. The proposed methodology hypothesizes that EHS that have not been objectively demonstrated lead to failures in the implementation of guidelines of Payments for Ecosystem Services (PES). Following this idea, we have tried to understand and quantify the relationship between vegetation coverage and streamflow regulation. To prove this relationship and the existence of the mentioned EHS in the La Guajira Ecosystem, we determined land cover changes from 2000 to 2013 using the Normalized Difference Vegetation Index (NDVI), and we also quantified the streamflow regulation using hydrological and meteorological time series in the study area. The analysis methods used were insufficient to determine the influence of vegetation on hydric regulation EHS; nevertheless, a greater influence of morphometry was observed in medium-and large-sized basins. Another important finding shows the relevance of selecting an adequate spatial and time resolution when quantifying water flow regulation services and its relationship with land cover characteristics. In this way, this exercise shows the complexity of quantifying EHS. Furthermore, we highlight some aspects that must be taken into account to properly quantify streamflow regulation due to vegetation coverage.


INTRODUCTION
Ecosystems provide essential services for the well-being of humanity.However, the continuous degradation of ecosystems caused by activities such as agriculture, cattle raising and mining have threatened the supply of these services.One of the main works regarding this concern is the Millennium Ecosystem Assessment (MEA), a program created to "provide an integrated evaluation of the consequences of the changes in ecosystems to human well-being and to analyze available options to improve the preservation of ecosystems and their contribution to meet human needs" (MEA, 2005).According to the MEA (2005), ecosystem services are defined as the benefits that the population obtains from ecosystems and are classified into four categories: (i) provision services that provide goods for the direct use for humankind (i.e., food, drinking water, wood and fibers); (ii) regulation services that maintain a world in which it is biophysically possible to exist (i.e., crop pollination, water damage reduction, climate stabilization); (iii) cultural services that make the world a place where people would like to live (i.e., recreation and aesthetic, intellectual and spiritual inspiration); and (iv) support services, defined as the subjacent ecosystem processes that provide the direct services described above.In particular, the services related to water that are offered by ecosystems have been called Ecosystem Hydrological Services (EHS) and cover the following four categories: water supply, water damage reduction, hydrological and cultural services, and support services related to water (Brauman et al., 2007).
Several authors have defined that the integrated management of ecosystem services in a region requires one to first identify, quantify and evaluate the ensemble of services that can be provided by different ecosystems (Harrison et al., 2010;Holland et al., 2011;Lele, 2009;MEA, Objective assessment of ecosystem hydrological services in tropical areas: … Rev. Ambient.Água vol. 12 n. 3 Taubaté -May / Jun.2017 2005;Quintero, 2010).The results of this process will provide the information base for the development of instruments whose aims would be the conservation of the ecosystems that provide these services, thus assuring human well-being for present and future generations (Fletcher and Breitling, 2012;Kinzig et al., 2011).In the same way, it will allow a more effective integration of these mechanisms for the institutional and political framework when making decisions about land management.(Locatelli et al., 2010;Quintero, 2010;Ferraro, 2011;Lü et al., 2012).One of the most widespread economic instruments for the management of ecosystems services is the Payments for Ecosystem Services (PES), which corresponds to a "voluntary transaction of an ecosystem service where at least a buyer and a seller participate; in that, the ecosystem service is 'well defined' and the service purveyor assures its provision" (Wunder, 2005).Despite its boom, it has been reported that in general, the hydrological PES schemes have not been preceded by a hydrological analysis that evaluates the effects of the activities and decisions promoted by the proposed PES scheme that ensures the supply of the studied service (Roumasset and Wada, 2013).It has been said about the EHS that in several of the world's regions, the abuse and wrongful use of resources, specifically those related to deforestation and forest degradation, along with contamination, increasingly threaten the availability and quality of the EHS, in particular that of the water supply (Calder et al., 2007;MEA, 2005;Quintero, 2010).In this context, the relationship between the vegetation cover, land use and hydrological cycle is widely recognized (Bosch and Hewlett, 1982;Brown et al., 2005;Calder, 1992;2007;Calder et al., 1997;Farley et al., 2005;Huber, 2008;Iroumé and Huber, 2002;Lara et al., 2009;Little et al., 2009;Ward and Trimble, 2004); but despite recent advances, the understanding of this relationship is still controversial.The results vary between geographical latitudes and can be influenced by different factors (climate seasonality, bio-geographical basin characteristics, and the time and space resolution of developed studies, among others).In general, the results reported by different authors have agreed in asserting that the loss of forest coverage for the establishment of agricultural land or cattle pastures leads to an increment of the peak river discharges during rainy seasons and to a water retention capacity reduction and an increased production of sediments as well.However, these conclusions, especially the ones referring to the services of water supply and water flow regulation, still feed into an enormous discussion (FAO and CIFOR, 2005;Quintero, 2010).
Recently, the Colombian legal system has incorporated laws that favor the function and recognize the benefits of PES; in 2007, Law 1151 was enacted, a framework law for the National Development Plan 2006 -2010 of Colombia.The structure of this law supports the management of the environment to achieve a sustainable development and includes the commitment of municipalities and departments to dedicate no less than 1% of its budget to purchase areas of high water availability or to implement PES instruments to ensure their preservation as water providers (Ruiz-Agudelo, 2011).Nevertheless, to implement these regulations appropriately, the formulation of PES schemes should be based on an appropriated quantification of the demanded EHS, and failing in this negatively affects the decisions made for natural resources management (Daily et al., 2009).Accordingly, this study establishes conceptual hydrological elements for the quantification of the EHS of water flow regulation (river regulation), taking as a case study the basins of the San Salvador, Negro, Jerez, Lagarto and Rancheria Rivers and the upper basin of the Cesar River in the Guajira Department, Colombia.

MATERIALS AND METHODS
According to the above, this document applies the following analysis stages: For stages (a) and (b), the concepts of Normalized Difference Vegetation Index and regulation coefficient are introduced, but basin morphometric characteristics are also presented, allowing for the further determination of the morphometric characteristics that control the surface runoff and water flow regulation processes.

Study area
The study area is approximately 7.000 km² and belongs to the "La Guajira" peninsula in Colombia.The range of geographical characteristics is quite varied and goes from the altitude of 5.575 meters above the sea level in the Sierra Nevada de Santa Marta to the coastal zone of the Caribbean Sea.The ocean influence and the fact that this territory is a part of a peninsular zone, in which the largest desert area of Colombia is located, makes it a region with a very high hydrological and ecosystem diversity.The variety of hydrological regimes within a relatively small area facilitates the identification of the ecosystem characteristics that influence the hydrological variables that are also the basis for EHS quantification.
To establish the dependency of water flow regulation variables on vegetation coverage, 23 watersheds were studied.These watersheds were obtained throughout supervised basin delineation using the ArcHydro tool incorporated in the extension HEC-GeoHMS tools in ARCGIS 10 (USACE and HEC 2011).As outlet points, hydrological stations were used (red points, Figure 1).Precipitation in the area is monitored by a meteorological network.Rainfall stations are shown in black hexagons.Objective assessment of ecosystem hydrological services in tropical areas: … Rev. Ambient.Água vol. 12 n. 3 Taubaté -May / Jun.2017

The Normalized Difference Vegetation Index (NDVI)
This index varies between values of -1 and +1, where dense vegetation presents values between 0.5 and 0.7 (Holben, 1986).In this sense, the high index values are directly related to the content of chlorophyll, the phenological dynamics, the amount of CO2, the amount of rainfall received and the potential evapotranspiration (Chuvieco, 2006).It is a reliable index that has extensive use in the description of biophysical parameters of vegetation cover.The NDVI is calculated by the Near-Infrared Band (BIRC) and the visible red band (BR) as Equation 1.

NDVI
(1) For the calculation, we used a total of 298 images from the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor, downloaded using the GloVis (USGS Global Visualization) tool.The images corresponded to the product MOD13Q1 taken from February 18, 2000 through January 17, 2013, every 16 days, with a spatial resolution of 250 meters.
To understand the spatial variability and dynamics of the NDVI vegetation coverage fields, temporal series were constructed in each one of the watersheds.The geoprocessing tasks performed to build such series included a) transformation into the spatial reference system MAGNA-SIRGAS and b) obtaining the proportion of four main categories of land coverage in each watershed according to the NDVI thresholds shown in Table 1 (Holben, 1986;Chuvieco, 2006).These tasks allowed us to obtain 144 time series of the NVDI, 4 for each of the 36 afferent areas and each with a length of 298 registries.This analysis was performed using scripts that were developed in Python 2.6 and its scientific programming modules Numpy, Scipy, OGR and GDAL.Finally, we considered linear trend tests, the least squares test and the Mann-Kendall test (Salmi et al., 2002).

Regulation coefficient
The regulation coefficient is a quantitative measure of the streamflow regulation using data from a hydrological station, which is defined as Equation 2.
(2) where: Ap and Aт are obtained from the flow duration curve (Figure 2), Efrain Antonio Domínguez Calle et al.Ap represents the partial area contained below the line of the mean flow, and Aт is the total area under the whole curve.The regulation coefficient is always less than or equal to one, and according to Samokhin and Salaviov (1980), its values of 0.6 to 0.8 are characteristic of the surface currents regulated by lakes or reservoirs (artificial or natural) and by the currents located in geographical zones with a high level of rainfall saturation.The coefficient values of regulation from 0.2 to 0.1 are characteristic of currents with a very low regulation capacity.For this article, the authors propose a classification of the different regulation values as follows: Very high (0.81 -1.00), High (0.71 -0.80), Medium (0.51 -0.70),Low (0.31 -0.50) and Very low (0.00 -0.30).Taking as a base point the previous definition, the regulation coefficient was calculated for each one of the hydrological series available using data at a daily resolution.However, there are different ways to calculate the regulation coefficient that depend mainly on how the flow duration curve is built.One way to calculate the regulation coefficient is to build only one flow duration curve for each station, using the full record of the time series of daily flow, and the multi-annual regulation coefficient is obtained as a result.In a second method, one flow duration curve is built for each year and each series, and the annual dynamic of the regulation coefficient is obtained as a result.The third method is similar to the second, but for every hydrologic station case, one regulation coefficient is calculated for each month based on the construction of a flow duration curve for every month.In this way, the monthly dynamic of the regulation coefficient was obtained.
Following these three methods of calculation, we obtained the following: a) For the first calculation method, 23 multi-annual regulation coefficients, one for each hydrological station.b) For the second method, 23 series that show the annual temporal evolution of the regulation coefficient on each hydrological station.c) For the third method, 23 series that show the monthly temporal evolution of the regulation coefficient on each hydrological station.
Objective assessment of ecosystem hydrological services in tropical areas: … Rev. Ambient.Água vol. 12 n. 3 Taubaté -May / Jun.2017 To compare the regulation coefficients from the different methodologies, three versions for the multi-annual regulation coefficients were used.To obtain the multi-annual values from the second and third methods, the respective series were aggregated.

Possible control factors of the regulation coefficient
Having a quantitative characterization of the streamflow regulation for the watershed calculated on the outlet points (hydrological stations), we proceeded to look for the factors that may affect this magnitude.A correlation analysis was performed between the regulation coefficient values, the NDVI index and the morphometric characteristics of each watershed.The morphometric parameters for each of the 23 watersheds were determined using the SRTM digital elevation model of a 90-meter cell size (CGIAR-CSI, 2012; Sanders, 2007).Scripts in Python and R were developed to calculate the morphometric parameters.In Table 2, the Pearson correlation coefficient between the multi-annual regulation coefficient and the morphometric parameters is shown (p-value ≤ 0,05).A detailed explanation of the calculation of these parameters is given in Samokhin and Saloviov (1980).Subsequently, the collinearity analysis was conducted between the 9 morphometric parameters with statistically significant linear correlations.Then, the hydrographic network density (D), the mean height of the basin (H ) and the mean flow (Qmed) parameters were selected to build a multiple regression to determine the regulation coefficient.

Variability of precipitation series
Simultaneously, there are local and global processes influencing rainfall regimes in different world regions (Huntington, 2006).All these processes of global change seem to influence not only the average magnitude of rainfall but also its variance (Haan Thomas, 1977;Kovalenko, 1993;Rozhdenstvenskiy and Chevotariov, 1974).The relationship between rainfall and surface runoff has been widely studied; in fact, this is the keystone that supports rainfall-runoff models (Chow et al., 1994;Clarke, 1973;Domínguez and Rivera, 2010).Changes in rainfall will likely affect streamflow.However, it is also well known that a basin works as a filter of the precipitation signal, and just a fraction of the precipitation variability reaches stream flows.
As a metric of how much a basin regulates the precipitation signal before it becomes a streamflow, we compared the precipitation regulation coefficients vs. the streamflow regulation coefficients.To achieve this type of analysis, the stations with rainfall registry were identified inside each analyzed watershed, and in cases without measurements, we associated the closest rainfall station.Then, to understand the role of basin regulation, the following procedure was developed: a) A linear correlation and regression analysis between the regulation coefficients of the stream time series and the regulation coefficients of the rainfall series was performed.
b) The correlation coefficient between each of the series of the proportion of NDVI and the streamflow series for each of the 23 areas of surface runoff was determined.This comparison was performed using a monthly and annual resolution of time series.c) From each of the regulation coefficients of the streams, we subtracted the obtained regulation coefficients from the corresponding rainfall series.In this way, a high difference means high basin regulation; in contrast, low differences mean low basin regulation.Then, a correlation analysis was performed using the series of regulation coefficient differences vs. NDVI series.

The Normalized Difference Vegetation Index (NDVI)
The NDVI values between 0.7 and 1.0 (Areas with very dense vegetation) appear to be dominant for the 2000-to 2012-time period, sharing 62% of the study area.In the range between 0.5 and 0.7 (areas with dense vegetation), 23% of the area is found on average; with regard to the values between 0.0 and 0.5 (areas without dense vegetation cover), the share reaches 15%.The NDVI values lower than 0.0 (Barren areas) occupy a low share, on the order of 0.2%.
The barren areas, despite registering a small increasing trend, lacked statistical significance, at least to reach a 95% confidence level.For areas with vegetation cover, the Mann Kendall test as well as the linear regression trend analysis show negative trends (α = 0.05); in this case, the decreasing rate was -0.5% annually.For areas with dense vegetation cover, even though the Mann Kendall test showed a significant negative trend (α = 0.05), the slope coefficient in this trend defined by linear regression showed a p-value of 0.06.As for the multi-annual monthly behavior of the vegetation cover, there is a clear pattern differentiating the vegetation cover in the months of February, March and April from the coverage pattern for the rest of the year.In these months, the proportion of the area with NDVI between 0.7 and 1.0 Objective assessment of ecosystem hydrological services in tropical areas: … Rev. Ambient.Água vol. 12 n. 3 Taubaté -May / Jun.2017 decreases on average to 33% of the basin's area.In contrast, areas with NDVI between 0.0 and 0.5 increase to 33%.This behavior seems highly conditioned by the rainfall regime, which has its first peak of the year between the months of March through April.Although the second high-water season occurs between the months of October and December, with higher rainfall values, there is no evidence of the same behavior in these months (Figure 3).

Regulation coefficient
The regulation coefficient magnitudes vary in accordance with the temporal resolution to which the flow duration curve was built.In general, the regulation coefficient magnitude increases with the flow duration curve's temporal resolution.In this way, the lesser regulation coefficient is obtained for the flow duration curve that is built for the multi-annual interval, and the biggest is obtained through the monthly regulation coefficient series, leaving an intermediate point of the regulation coefficient that is obtained as the average of the series of annual regulation coefficients.It must be noted that the regulation coefficients of these three methodologies scale linearly (Figure 4).Differences in the regulation coefficients obtained by the three methodologies seem to obey the fact that a monthly grouping of daily stream data has less variance than does one from annual and multi-annual groupings.It is widely known that the regulation coefficient depends in a significant way on the variance of the data used to calculate it.

Possible factors controlling the regulation coefficient
The analysis between the proportion series of the area occupied by NDVI ranges and the regulation coefficient series did not find statistically significant linear correlations.In the case of the correlation analysis of NDVI with the regulation coefficient through a monthly series, the result of a low correlation is not unexpected.To begin with, because of the positive bias of the regulation coefficient magnitude when it is calculated as a monthly resolution, it is correlated with a low variability magnitude (the monthly regulation coefficient) and with very high variability (the monthly NDVI).From this contrast, we can only expect a low linear correlation.On the other hand, in the case of the linear correlation analysis for the annual series of the regulation coefficient and the NDVI index, it might be explained by the annual NDVI being a bad indicator of the vegetation cycle.This result is a consequence of most of the vegetation dynamics occurring in the hydrological year, and using the mean NDVI value for the whole year could hide this property.
Additionally, even though the set of data of the NDVI has gone through a process of quality control, there are still problems that must be considered when using this type of study, such as the presence of mixed pixels (pixels that have more than one coverage type affecting the values of NDVI), registry deficiency (appearing as a result of mistakes when re-calculating the satellite position at the time of the image being recorded), comparison of values of the NDVI amongst pixels (factors such as plant architecture arrangements, canopy interactions, height, species composition, vegetation strength, leaf properties and vegetation stress can significantly affect the information of remote sensors, making it so that equal values of NDVI may represent different conditions for different vegetation communities, recommending previous knowledge in these elements of the study area) and quality information about spatial location (bigger mistakes may occur close to the equator because of the variations in the angle of the solar zenith in the majority of the satellite registry) (Pettorelli et al., 2005).
It is also important to take into consideration the fact that even if the image availability of MODIS-NDVI, which provides a high-quality series of time data, represents a great advance in the monitoring of the annual changes in the earth coverage and vegetation condition in large geographical regions, it also brings with it disadvantages and spatial resolution limitations (250 m) because minor event changes to approximately 1.5 hectare would have a low probability of being detected.This issue is problematic for monitoring the changes in zones of riparian dampening and other events of small scale conversion that could be associated with high value ecological resources such as the evaluated EHS in this study.
The good relationship between the drainage density parameters, basin mean height and mean flow with the regulation coefficient can be explained more easily in some cases than in others.As for the drainage density, the negative correlation seems very clear; it is to be expected that in the measure in which a basin shows a higher length of currents per unit of area, the superficial surface runoff will flow easily, concentrating in the stream quickly, which will create more severe/sharp peaks of flooding.In contrast, a low drainage density indicates that the water should drain for a longer time over the terrain surface, in which the roughness coefficients are usually bigger than in the main stream, creating in this way higher friction to the water movement, slowing its concentration in the main stream and thus conducting a smoothed peak raise (Chow, 2009).In the first instance, it is not expected that the basin's mean height would have a positive influence on the regulation coefficient; on the contrary, an elevated mean height Objective assessment of ecosystem hydrological services in tropical areas: … Rev. Ambient.Água vol. 12 n. 3 Taubaté -May / Jun.2017 may indicate high slopes in which the superficial surface runoff is developed with high speeds.However, on the analyzed study area, it was found that the weight of this problem is not very significant (7.2659E 05 vs. 0.77492865 corresponding to drainage density), and when analyzing the spatial distribution of the hydrological stations for which the regulation coefficient was calculated, it was found that the sloped and steep relief areas of the study area conserved vegetation coverage in such a way so that the small influence of the median height of the basin could be attributed indirectly to the EHS offered by the vegetation.
Finally, the multiple regression between the three selected parameters (D, H and Q ) proved to be statistically significant.According to this result, morphometrics can explain close to 77% of the regulation coefficient of the streams time series, and it is proposed that the remaining 23% should be explained by other factors such as vegetation and ground characteristics (Table 3).

Variability of precipitation series
According to the trend and the frequency distribution of the quotients between stream and rainfall, the regulation coefficients of the stream series are two to five times bigger than the regulation coefficients measured on the rainfall time series (Figure 5).This result, together with the influence of the morphometries on the regulation coefficient variability of the streams time series, indicates that in the study area, the basin as a system executes an attenuating filtering action that diminishes the rainfall variability that comes into the system.According to the proposed hypothesis, it was expected that as a result of steps (b) and (c), part of the unexplained variability could be assigned to the basin's vegetation coverage.In this sense, the applied correlation analysis did not show any significant linear correlation between the NDVI indexes and the regulation coefficient differences determined in step (c).Despite this result, it cannot be sustained that vegetation coverage does not have any effect on streamflow regulation, and it can only be said that there is no overwhelming evidence found to prove this ecosystem service.However, given the relationship between morphometric and streamflow regulation and, in particular, the positive linear relationship of the mean height of the basin with the regulation coefficients, it can be indirectly attributed to vegetation.This assertion takes into account the fact that mountain and hill landscapes have a mean height in the study area that is higher than flat landforms, and they have been higher in the last transformation of vegetation coverage.This result indicates that vegetation coverage also has an influence on the regulation coefficient magnitude of the stream series.On the other hand, this study did not contemplate the influence of the soil characteristics over the magnitudes of the regulation coefficient.The soil-vegetation system establishes a unit that surely produces synergy into the streamflow regulation of the basins; however, for the studied basins, the existing studies do not provide adequate information on the soil and its humidity dynamic to allow one to tie together the soil properties.Recent publications suggest that the satellite information on land humidity is beginning to have precision levels that convert it into relevant information for the analysis of the influence on the soil-vegetation system over the regulation coefficient magnitudes.

CONCLUSIONS
This study quantitatively determines the streamflow regulation and shows the preponderance of influences that have morphometric basin properties.This study also shows the difficulties of establishing a relationship between vegetation coverage and streamflow regulation as an ecosystem service.This same analysis is not conclusive in the determination of additional factors to the morphometric that control the regulation coefficient magnitudes obtained for the analyzed hydrologic stations.The analysis methods used were insufficient to determine the influence of vegetation on the hydric regulation EHS.This result can be attributed to congruence problems in the temporal resolution in the regulation coefficient of the series and the NDVI and to the presence of mistakes in the values of the latter due to mixed pixels and the other elements mentioned above, together with its limitations in spatial resolution.Finally, the findings and conclusions elaborated during this study allow us to suggest the existence of a vacuum at the time of the formulation of Payments for Ecosystem Services schemes.It is, therefore, necessary to ponder the hypothesis that is usually applied to sustain these schemes.According to the experience shown here, it is suggested that the decision makers in this matter include a project stage that looks for a clear and overwhelming demonstration of the ecosystem hydrological services in a permanent manner.
a) Characterization of the state and changes in vegetation land cover through the NDVI; b) Hydrologic quantification of streamflow regulation; Efrain Antonio Domínguez Calle et al. c) Search of streamflow regulation-NDVI relationships; and d) Ranking of ecosystem factors that influence the analyzed streamflow regulation.

Figure 1 .
Figure 1.Spatial distribution of hydrological stations and watersheds studied.

Table 1 .
NDVI ranges used to classify the type of coverage according to the proportion and type of vegetation.

Figure 2 .
Figure 2. Determination of the regulation coefficient from the flow duration curve.

Figure 3 .
Figure 3. Monthly multi-annual distribution of area proportions by ranges of NDVI.

Figure 4 .
Figure 4. Regulation coefficients obtained from the three-stream aggregation methodologies of the series.

Figure 5 .
Figure 5. Frequency distribution of the quotients C R-Streams/ C R rainfall .

Table 2 .
Correlation between the morphometric parameters on each of the areas of surface runoff and the regulation coefficient.

Table 3 .
Multiple regression parameters for the calculation of C R .