Increase in surface runoff in the central mountains of Mexico: lessons from the past and predictive scenario for the next century

The hydrological response of a medium scale mountainous watershed (Mexico) is analysed over half a century. The hydrograph separation highlights an increasing surface runoff contribution since the early 1970’s. This increase is attributed to land use changes while the meteorological forcing (rains) remains statistically stable over the same period. As a consequence, the intensity of annual extreme floods has tripled up over the period of survey, increasing flood risks in the region. The paper ends with a climatic projection over the 21st century. The decrease of precipitation and the increase of temperature should accentuate the trend engaged since the 1970’s by reducing groundwater resources and increasing surface-runoff and associated risks.


Introduction
Over the last decades, Mexico has suffered from degradation of its surface water bodies which is imposing undeniable economical costs (Alcocer and Escobar, 1993). Nowadays Mexican water resources are commonly considered poor in quality and sparse in quantity (Vidal et al., 1985), a situation exemplified in Michoacán state. A recent study led by various local institutions (OrtizÁvila, 2009) stressed that water contamination, solid residuals management and drinkable water are the main environmental priorities of the Michoacán settlements. In terms of quantity, the state has been facing a decrease of about 70% in its surface water resources over the last century (Vargas Uribe, reported by Morales, 2007). This Correspondence to: N. Gratiot (nicolas.gratiot@ird.fr) evolution has been correlated with the high emigration rate in Michoacán (63% of total population), which contributed to the main soils and water resources use changes (López-Granados et al., 2006).
The present paper aims at investigating water cycle changes from 1956 to 2001 in the Cointzio watershed (Michoacán state), a medium scale mountainous basin representative of the Central Transvolcanic Mexican Belt.
Our approach considers the impact of climate and humaninduced environmental changes on water runoff in the watershed. To this purpose a hydro-meteorological database running from 1956 to 2001 was employed following the latter methodology: available rain and water discharge series were gathered (and digitalized when necessary) and criticized (missing data and quality control) to provide reliable series. A particular effort was paid to estimate the accuracy of the historical water discharge data in relation with sampling frequency (Sect. 3). In a second time, hydrograph separation technique was applied to water discharge records to define the baseflow/surface-runoff ratio and its evolution over decades.
The generated time series were tested with various statistical methods in order to provide accurate hydro-climatic trends and objective interpretations (Sect. 4). Finally, a climate model was applied to five meteorological stations to go towards scenarios of evolution for the coming century (Sect. 5).

Study area
The Cointzio watershed is located in the hydrological region of Lerma-Chapala, within the Central Transvolcanic Mexican Belt, in the state of Michoacán (Fig. 1). It drains a surface of about 650 km 2 , ended by the Cointzio reservoir (4 km 2 , 65 Mm 3 ). The latter provides about 22% of drinking water distributed in Morelia, capital of the state situated 13 km downstream. Water demand of the city has been growing over the last decades because of increasing individual water consumption coupled with a severe urban growth: Morelia experienced an augmentation of its population of 600% during the period 1975-2000(López-Granados et al., 2001 and counts now over 700 000 inhabitants (INEGI, 2006).
The climate of the region is temperate sub-humid, characterized by a rainy season from May to October and a dry season the rest of the year (Rubio, 1998). Mean annual rainfall is 770 mm in Morelia, ranging from 400 to 1.100 mm/y (Carlón Allende and Mendoza, 2007 and Fig. 5 in this paper). The seasonality implies that 77% of the drainage system consists of temporary watercourses active only during rainy season (Susperregui et al., 2009). The main river of the watershed is the Rio Grande de Morelia whose source lies about 25 km upstream of the Cointzio reservoir.
The Cointzio basin is underlain by igneous rocks (both lavas and pyroclasts) originating from Quaternary volcanic activity. Soils and landforms developed in most of the watershed have been derived from these volcanic materials (Carlón Allende et al., 2009). Such soils present susceptibility to erosion (Poulenard et al., 2001). They are mainly Andisol in the headwater areas, Acrisol in the hillsides and Luvisol in the plain (INEGI, 2002).
According to López-Granados et al. (2001), the basin can be classified as agro-sylvo-pastoral. Using GIS to represent data of land cover and land use, Mendoza and López-Granados (2007) were able to identify the major changes in the Cointzio watershed over the period 1975-2000. The main changes in surface were increase of scrublands (9.6%), recovery of forests (6.2%), deforestation (5.5%), degradation of forests (4.1%) and urbanization (1.3%). The major part of these changes occurred during the 1986-1996 period.

The hydrometeorological historical database and its limitation
Rain series were extracted from the Mexico Climatological Station Network Data (CLICOM). This database reports all of the meteorological stations surveyed by the Servicio Nacional de Meteorología de México. A pre-analysis was realized by Vinson (2008 1955-60 1961-70 1971-80 1981-90 1991- Following this pre-processing of data, we decided to focus the analysis to the stations of Contzio and Santiago Undameo. These latter are the only consistent time series within the watershed (see Vinson, 2008 for more details). The three other stations have been considered during the modelling exercise presented in Sect. 5 to offer a regional projection of potential climate variation.
The gauging station of Santiago Undameo (SU in Fig. 1) constitutes the ultimate control upstream the reservoir of Cointzio: it drains 628 km 2 . Its monitoring has been launched in 1939, year of construction of the reservoir of Cointzio. A Parshall flume was built, allowing a control of the hydraulic section, and a stage-discharge rating curve was established by the Comisión Nacional del Agua (CNA). The monitoring ended in 2001.
The 1939-1985 period was previously digitalized by the CNA (BANDAS database), while the 1985-2001 period was digitalized as part of our investigations. Although overall database covers the 1939-2001 period, we decided to focus on 1956-2001 series since the first flume was destroyed during a major event in 1953. The re-building and its calibration ended during 1955. The CNA archives at our disposal do not report any change of the flume section (by sedimen-tation, erosion or vegetation growth); therefore the stagedischarge rating curve was assumed to remain constant between 1956 and 2001. The water discharge database is presented in Fig. 2b. Some minor missing period events were identified in 1975, from the 17 to 20 October because of a flood event. Four days of missing data are reported in 1992 and one more in 1997. With a period of missing data of several weeks in September and October, the year 1998 could not be considered in our analysis.
From the beginning of 2008, the station has been updated to estimate instantaneous water and sediment discharges. Water level was surveyed at a five minute time-step with a Thalimede OTT water-level gauge. Water discharge timeseries were determined through the stage-discharge rating curve. Our aim was to take advantage of this high-frequency monitoring to evaluate the accuracy of daily discharge values provided by historical time-series, and thus to carry out a long-term analysis without misinterpretation.
Daily discharge data was historically deduced by the CNA from a minimum of three manual measures (respectively at 06:00 LT, 12:00 LT and 18:00 LT). Taking into account the occurrence of flash flood events in the basin, the historical sampling frequency is questionable. To validate the methodology, a sub-sampling of the real-time-series acquired in 2008 was generated. Resulting mean daily discharges were then compared to reference mean daily discharges derived from the high-frequency time-series (Fig. 3).
Reference discharge values and sub-sampled ones exhibit a strong linear relationship (Fig. 3a). Distribution of the relative error presented in Fig. 3b exhibits a Gaussian shape, centred on zero and characterized by a low standard deviation of about 7%. Such pattern demonstrates the reliability of the historical sampling. Loss of accuracy remains very limited if focusing on short-term dynamics as well as on seasonal water budgets calculation (maximum under-estimation of 2% of reference volume). The further processing of 1956-2001 historical time-series is thereby validated.

Hydrograph separation technique
We aimed to consider the effects that may have induced an alteration of physical characteristics of the watershed, such as land use change, on the hydrological cycle. We subsequently focused on the characterization of annual baseflow/surfacerunoff ratio and its evolution over the past fifty years at Santiago Undameo.
Stream flow hydrographs were separated into Annual BaseFlow (ABF) and Annual Surface-Runoff (ASR) components. The baseflow component has traditionally been associated with groundwater discharging into the stream and the surface-runoff component with precipitation that enters the stream as overland runoff (Sloto and Crouse, 1996;Chapman, 1999). The aim of this paper is not to work on the physically based concept of hydrograph separation; the technique was applied as a tool for detecting trends in water discharge behaviour. The reader can refer to Chapman (1999) and Nathan and McMahon (1990) to get an overview of techniques commonly used by engineers to quantify the baseflow contribution of a watershed.
A variant of a hydrograph separation method was programmed using Matlab© software. The algorithm used is based on the concept developed by Pettyjohn and Henning (1979) and is commonly referred as the "smoothed minima technique" in the literature (Brodie and Hostetler, 2005). This method can be described as connecting with straight lines the minima of fixed intervals τ of the hydrograph. The sequence of these connecting lines defines the baseflow hydrograph. In case of missing data, a mean value calculated from nearest neighbours was used. Calculated Base Flow delimitation (CBF) was thus obtained as follows: with d: day of the year τ : interval parameter (in days and necessarily pair) Q d : mean daily water discharge (in m 3 .s −1 ) Annual Base Flow (ABF) and Annual Surface Runoff (ASR) were then estimated by: In our analysis, we tested a broad range of τ values from four to fourteen days. From a general point of view, it is clear that CBF (and thus ABF) will decrease for increasing τ .
Such behaviour is illustrated in Fig. 4a and b for two contrasted hydrological years that occurred in 1993 and 1994.
Year 1993 was hydrologically active with many significant flooding periods lasting several days. CBF function decreases significantly for increasing τ values from six to fourteen days. Conversely, in the case of year 1994 (Fig. 4b), periods with high water discharge were scarce and never exceeded a couple of days. As expected, CBF function does not vary significantly with τ . The estimation of baseflow level requires a minimal period of one day without significant precipitation. By examining the precipitation database at our disposal, it appeared that the longest time-interval required to reach a day without significant rainfall event (P d <5 mm/day) is approximately of five days. Thus a τ interval of six days is best suitable in Fig. 4a, while coinciding with both rainfall pattern and hydrograph visual examination. However, to prevent any misinterpretation, the algorithm was also run with both a shorter (four days) and a longer (ten days) τ value.

Trends detection from statistical methods
Statistical significance of gradual trends was detected by applying the rank-based Mann-Kendall test (Mann, 1945;Kendall, 1975) and magnitude of trends was estimated from Sen's method (Sen, 1968).
The Mann Kendall test is a non parametric statistic that has been widely used to assess the significance of monotonic trends in hydro-meteorological time-series (e.g. Lettenmaier et al., 1994;Marengo and Tomasella, 1998;Jiang et al., 2007;Zhang et al., 2008;among others). The test assumes that there is no serial correlation in the data. Such assumption is reasonable for the rainfall and runoff records presented in this paper.
The null hypothesis H 0 is that the sample of data is independent and identically distributed. The alternative hypothesis H 1 is that a monotonic trend exists in the time-series. Mann-Kendall method was applied by considering the statistic S as: Where x i and x j are the sequential data values, n is the length of the time-series and sign(x i − x j ) is −1 for (x i − x j )<0; 0 for (x i −x j ) = 0 and 1 for (x i −x j )>0.
In the absence of ties, the variance Var [S] of the statistic S was obtained as (Kendall, 1975): The standardized statistical test Z was computed by: Positive value of zindicates an increasing trend while negative zindicates a decreasing trend. When testing monotonic trends at an α significance level, H 0 was rejected for absolute value of Z greater than Z (1−α/2) , where Z (1−α/2) is the value of the standard normal distribution with a probability of α/2. In this work, significance level of α=0.01 (99% confidence) was applied and null hypothesis was rejected if |Z| > Z 0.995 =2.575.
Sen's method is a non-parametric statistic used in determining the presence and magnitude of a trend slope. This test proceeds by calculating the slope as a change in measurement per change in time. Trend slopes magnitudes were obtained following the method of Hirsch et al. (1982): Where x j and x i are data points measured at times j and i, respectively.
Mann-Kendall test and Sen's method were applied on precipitation, surface runoff and water discharge time series. Results are presented in Sect. 4.3.

Pattern of the precipitation time-series from 1956 to 2001
Our analysis is depicted in Fig. 5. It is based on the examination of three statistical indicators, namely the total annual precipitation (in mm/y), the maximum daily precipitation (in mm/day, from year to year) and the sum of the 5% days undergoing maximum rainfall events (in mm/18 days, from year to year). The latter basically corresponds to a virtual "top eighteen days" each year (sum of the rain data exceeding 95th percentile). Annual precipitation series fluctuates significantly in the range 650-1200 mm/y; however there is no evidence of tendencies over the last fifty years. Cumulated top 5% precipitations correspond to about 400-500 mm of total annual precipitation, which means that 50% of annual volumes are precipitated in 5% of the time. Such ratio highlights the heavy rainfall pattern characterizing the region during wet season.
Maximal precipitation exhibits the same pattern than annual volumes and top 5% precipitations. Although maximum daily precipitation fluctuates from year to year, it basically remains stable in the range 35-60 mm and roughly corresponds to 5% of annual precipitation.

Pattern of the Annual Surface Runoff (ASR) from 1956 to 2001
Annual Surface Runoff (ASR) and Annual Base Flow (ABF) were calculated from 1956 to 2001 by applying Eqs.
(1) and (2) with τ intervals of four, six and ten days. As shown in Fig. 6a, ASR volume series remains globally constant over 1956-2001 period, in the range [1-4] 10 7 m 3 y −1 . ABF volume is always predominant. It fluctuates in the range [3-8] 10 7 m 3 y −1 . Unlike ASR, ABF volume exhibits clear inter-decades trends. It shows a continuous decrease from the beginning of the seventies to the end of the eighties that could be correlated to the dry cycle discussed by Carlón Allende and  and Metcalfe et al. (2007).
Another interesting pattern is depicted by %ASR and %ABF series presented in Fig. 6b. Unlike ASR volume, %ASR has been increasing substantially since the seventies. It remains almost constant (nearly 25%) until 1970 and then increases to reach about 40% by the end of the eighties to 2001. This significant water balance change does not depend on the hydrograph separation processing as it is observed for every τ value (grey shading in Fig. 6b). The latest values remain in a typical range of values encountered in North American watersheds (Neff et al., 2005).
Looking simultaneously at volumes and percentages it appears that surface runoff versus baseflow is not affected by the precipitation pattern. This is well illustrated by couples of years 1956-1957 and 1993-1994 1956 and 1957, the total water discharge tripled up, from 3.9 10 7 m 3 to 11.3 10 7 m 3 (Fig. 6a) without modifying the partition between runoff and baseflow. As shown in Fig. 6b, %ASR (or %ABF) remained remarkably stable with a value of 28.5% for the dry year (1956) and of 33% for the wet year 1957. A similar pattern is observed four decades latter for years 1993 and 1994. The total water discharge dropped from 15.3 10 7 m 3 to 4.4 10 7 m 3 while %ASR (or %ABF) did not change significantly (34% to 36%). Accordingly, a pronounced stability is observed in the partition of water for consecutive years whatever the annual amount of precipitation (rainy or dry year).
The last two paragraphs point out that the ratio ASR/ABF and its evolution over years does not depend on meteorological forcing. Hence, it is the partition of rainfall input between surface-runoff and baseflow that has changed gradually in Cointzio. Between 1956 and 2001 the watershed has progressively been generating more surface runoff and sustaining lower baseflow. In other words, water discharge series progressively showed narrower peaks over years. This change in hydrological response is very likely to be associated with physical alterations of the watershed, such as land use change and surface impermeabilization.
During the 80's, the state of Michoacán experienced very strong migration fluxes which contributed to the abandonment of many rain-fed agricultural lands in steep relief (López-Granados et al., 2006). Inhabitants have been leaving countryside while Morelia (capital of Michoacán) has been growing tremendously, from 100 000 inhabitants in 1950 to 700 000 in 2005. As a consequence, water bodies now suffer higher turbidity level than in the past, as a response of deforestation and growing urbanization (Rendón López et al., 2007).

Evolution of water discharge extreme events
The goal of the present section is to quantify hydro meteorological trends and evaluate the relative impacts of both meteorological and human-induced changes on floods. The Mann Kendal test was applied to annual rain (Cointzio weather station), %ASR series and water discharge. Results are presented in Table 2 and trends are depicted in Fig. 7.
As specified in Sect. 4.1, rain series does not exhibit any trend (Fig. 7a, Table 2). Hence, it can be considered that water input in the watershed did not change significantly during the study period. Conversely, %ASR increases significantly between 1956 and 2001. Applying Sen's method, the resulting trend slope has amplitude of 0.33, which corresponds to an increase of 69% in the contribution of surface-runoff to overall water discharge. Looking into more details, the series can be split into two distinct periods: from 1956 to 1970 values are slightly decreasing and from 1970 to 2001 there is a clear increase ( Table 2). The human-induced changes at the origin of the %ASR increment (see Sect. 4.2) have very probably contributed to the increase in extreme floods intensity depicted in Fig. 7c. This graph displays the annual flood peak Q d max. The series exhibits a clear increase with trend slope amplitude of 0.5. It corresponds to an increase rate of Q d max as high as 232% over half a century. The series also shows a clear increment of the variability since 1975 but there is no doubt that flood risk has increased over years.  (period 1961-1990, centred in 1975), and for a climate change scenario for decades centred in 2030, 2060 and 2090. (b) Estimated mean annual temperature for contemporary climate (period 1961-1990, centred in 1975), and for a climate change scenario for decades centred in 2030, 2060 and 2090. (c) Annual aridity index for contemporary climate (period 1961-1990, centred in 1975), and for a climate change scenario for decades centred in 2030, 2060 and 2090.

Potential evolutions over next decades
In terms of water management, it is crucial to evaluate whether the engaged trend will continue or not.
Future climate change prediction was thus conducted for the five weather stations at the vicinity of the watershed. Estimates of mean annual precipitation and temperature presented in Fig. 8a and b are deduced from a spline climatic model developed for a normalized period (years 1961 to 1990, nominated "contemporary climate") and updated with weighted outputs from a Global Circulation Model (Canadian Centre for Climate Modelling and Analysis), emission scenario A2, for the decades centered in the years 2030, 2060and 2090(Crookston, 2008Sáenz-Romero et al., 2009). Following Sáenz-Romero et al. (2009), precipitation and temperature series are combined to provide an aridity index (ratio of the square root of annual degree days >5 • C to annual precipitation) that represents the potential for moisture stresses to develop in the vegetation (Fig. 8c).
Average estimations among the five weather stations for the contemporary climate and compared to the global change scenario, indicate an expected decrement in precipitation of 15.4, 19.1 and 27.7%, and an increment in mean annual temperature of 1.6, 2.5 and 4.4 • C, for the decades centred in the years 2030, 2060 and 2090, respectively.
In terms of extreme events, it would be inappropriate to associate the annual precipitation decreasing pattern with a decrement in flood intensity: flood dynamics is governed by instantaneous rainfall intensity rather than by annual budgets.
The combination of increment of temperature and decrement of precipitation evidently causes an increment of aridity, measurable through the annual aridity index (Fig. 8c, where smaller values mean a climate colder and moister and larger values indicate a climate warmer and dryer). Since the aridity index is closely related to the type of vegetation (Rehfeldt, 2006;Rehfeldt et al., 2006), it is reasonable to expect that climate change will cause a decrement of the veg-etation coverage, and consequently, an increment of runoff contribution, as previously reported by Ranzi et al. (2002) and García-Ruiz et al. (2008).

Conclusion and perspectives
Over half a century, significant changes have occurred in the water balance of Cointzio, a medium scale watershed representative of the mountainous highlands of Central Mexico. Surface-runoff has increased of about 70% and has consequently led to a severe increase in extreme flood events (magnitude has tripled up over . During the same period, precipitation exhibited no statistical trend in the region, attesting the absence of significant modification of the annual budget. The main annual surface runoff increase occurred from the 1970's and is believed to be principally driven by human-induced changes. Predicting in which proportion water cycle disequilibrium engaged since the 1970's will extend in the coming decades is a key issue. First, substantial climate modifications are likely to occur: numerical simulations presented herein indicate a drastic reduction in rain budget (−28% in 2090) and an increase of the aridity. These changes are expected to alter vegetation coverage and consequently, accentuate the surface-runoff contribution. Secondly, the growth in water consumption related to urbanization will reinforce the pressure on local water resources.
At a global scale, Bates et al. (2008) recently reported that flash-floods and inundation should become more frequent worldwide. A regional analysis of precipitation intensity would be required to further improve our understanding of flood risk in the highlands of Central Mexico. especially M. Esteves and C. Prat who coordinate these projects. The authors wish to thank the Centro de Investigaciones en Ecosistemas, particularly K. A. Oyama. Our thanks are extended to the Comisión Nacional del Agua for providing us with the long-term hydro-climatologic data and for all people who contributed to the digitalization (P. Burgos, C. Hernandez Wences, and S. Van Hecke).
Edited by: F. Laio The publication of this article is financed by CNRS-INSU.