Assessment of the Water, Environmental, Economic and Social Vulnerability of a Watershed to the Potential Effects of Climate Change and Land Use Change

: In semi-arid regions, where hydrological resources are very vulnerable and where there are water shortages in many regions of the world, it is of great importance to assess the vulnerability that a system is facing or will face to the potential impacts of climatic changes and changes on the use of land. For that reason, this research focuses on evaluating the global vulnerability of a hydrological basin, taking into consideration these changes. Being different from the existing methodologies that assess the vulnerability, our methodology interconnects through a new interface a distributed hydrological model, global climate models, climate change scenarios, land use change scenarios and the largest number of system variables calculated with information from ofﬁcial sources. Another important point of our methodology is that it quantiﬁes the global vulnerability of the system, taking into consideration hydrological, environmental, economic and social vulnerabilities. The results obtained show that the proposed methodology may provide a new approach to analyze vulnerability in semi-arid regions. Moreover, it made it possible to diagnose and establish that the greatest current and future vulnerabilities of the system are the result of activities in agricultural areas and urban centers.


Introduction
Vulnerability is a concept that emerges as a necessity for climate science and policies in the face of the potential impacts of climate change. During the last decade, the efforts to assess vulnerability to the potential impacts of climate change have contributed to the development of theories and assessment methodologies that can be found in the reports of the Intergovernmental Panel on Climate Change (IPCC) [1]. The IPCC defines vulnerability to climate change as a function of a system exposure and stimuli-response to climate and its ability to adapt to these effects [2]. According to [3], vulnerability is an indicator of possible future damages. Historically, the term "vulnerability" has been used in several contexts, such as: food security [4,5], livelihood [6], natural disasters [7] and general risk management [7,8]. In that same context, climate change has caused global concern due to adverse effects on the global ecosystem, the economy and society [9]. Therefore, in recent years, the tendency has been to assess the environmental [10], economic [11] and social vulnerability [10,11]. Other studies have focused on assessing the world's hydrological vulnerability to the possible effects of climate change, such as those carried out by [12,13]. The foregoing is due to considering that inland water resources have been identified as one of the most important global change problems [14]. According to [1], assessing vulnerability to the potential impacts of climate change has the objective of developing policies that reduce the risks associated with climate change. According to [15], climate change has become a threat by applying greater pressure on already stressed hydrological systems and resources. Moreover, in the last decades, the impacts of climate change have already been visible, mainly in the increase of temperatures, and the increase and/or decrease of rainfall in several regions of the world [16,17]. According to [14], these impacts are of extreme importance, considering that the hydrological resources are very vulnerable, and that at present, many regions of the world face water scarcity. Under this reality, it is transcendental to estimate the vulnerability that a system is currently facing or will face taking into account future changes [18]. It is for that reason that this article focuses on assessing hydrological, environmental, economic and social vulnerability in a spatially distributed global context. Global vulnerability combines and distributes at the cell scale the four types of vulnerability described above, in which hydrological, economic, social and environmental factors converge.
According to [19], when analyzing different studies, we can observe that the methodology mostly used to quantify the vulnerability of a system by is using a composite index that includes a set of indicators. The same authors comment that these indicators represent the vulnerability of a system and are mathematically combined into a single composite index. A composite index may be easier to interpret than the tendencies of single indicators [19]. In literature, different methodologies might be found to quantify the composite indicator, such as: multiple linear regression models [20], principal components analysis and factor analysis [21], neutralization of correlation effect [22], distance to targets [23], distance to objects [23], experts opinion [24], public opinion [25], spatial multi-criteria analysis [26], Analytic Hierarchy Process [27], Bayesian model [28], Artificial Neural Networks [29], and fuzzy logic techniques [29]. According to [15], the vulnerability calculated using the existing methodology does not consider all the variables involved in the system, thus, limiting its accuracy.
In the same context, the methodologies do not introduce land use change in the vulnerability assessment. According to [30], the impacts of land use change (Urbanization) on a hydrological system produce an increase in volume and surface runoff rate, causing a decrease in groundwater recharge and base flow. According to [31,32], the land use change and climate change have been identified as two of the main factors affecting the ecosystem services. Studies like that of [33] highlight the importance of using land use change in the environmental assessment. For that reason, the authors consider to include it in the quantification of the global vulnerability. In addition, we believe that land use change will have an important effect, because the study area is located in a region where availability studies accentuate a shortage of water resources and over-exploitation of groundwater. Finally, the hypothesis of this article is that the vulnerability of a basin to the potential impacts of climate change and land use change can be evaluated through the connection of mathematical tools that will contribute to decision making, with the aim of mitigating the impacts and knowing the resilience of the system.

Case Study
The study area used in this research is the Turbio river sub-basin ( Figure 1). The sub-basin is located in the Lerma Santiago Hydrological Region, which includes the states of Guanajuato and Jalisco in Mexico. The sub-basin has an area of 2983 km 2 and is one of the most important for agricultural production of the country. The sub-basin is located between elevations 1726 and 2872 m. Accumulated average annual rainfall varies from 560 to 807 mm, with a very marked seasonality of the rainy season and an occurrence of 90% between the months of June and September. The maximum average temperature is 25.59 • C and the minimum average temperature is 9.56 • C. The variations in temperatures have a significant influence on evaporation rates in the study area with rates ranging from 0.25 mm day −1 to 12.45 mm day −1 .
In the last three decades, the Turbio river sub-basin has presented problems of over-exploitation of surface and groundwater. This is due to the population growth and increased agricultural activity, supplied with surface and groundwater. This over-exploitation has generated great uncertainty among decision-makers about the future of water availability in the sub-basin, and how it will be affected by the potential effects of climate change. This scenario makes our methodology more relevant, as it aims to respond to the concerns, such as the future of water availability from the point of view of hydrological, environmental, economic and social vulnerability.

Methodology
In contrast to the existing methodologies, in this article, a totally different new methodology is proposed, which combines a distributed hydrological model, global climate models, climate change scenarios, land use change scenarios and the largest number of system variables calculated on the basis of information provided from official sources. Another important point of our methodology is that it quantifies the global vulnerability of the system, taking into consideration the hydrological, environmental, economic and social vulnerability. It is important to emphasize that nowadays there are no studies that consider all the vulnerabilities (i.e., hydrological, environmental, economic and social) in the assessment of global vulnerability, connecting several mathematical models in a hydrological basin. This makes our methodology new with respect to the existing ones. In addition, the effects of land use change are included in the evaluation.
The methodology that we propose to carry out this research has included the development of a Predictive Vulnerability Interface that we have called MPDV1.0. It has been programmed in Visual Basic (Microsoft, Washington, DC, USA) language and allows quantifying the hydrological, environmental, economic, social and global vulnerability of a hydrological basin. MPDV1.0 integrates the distributed hydrological model TETIS [34,35], the Global Climate Models (GCM) of the Inter-Comparison of Coupled Models-Phase 5 (CMIP5) Project [36], and RCP4.5, RCP6.0 and RCP8.5 scenarios of the Intergovernmental Panel on Climate Change (IPCC) [37]. The foregoing is following the scheme shown in Figure 2. Unlike the methodologies proposed by [38][39][40][41], MPDV1.0 calculates vulnerability in a distributed way, including all the processes and storages of the hydrological cycle through its modeling in TETIS. In addition, the effects of climate change and land use change are included as well. In the hydrological modeling with TETIS, the future land use changes are introduced, estimated with the Terrset software (Clark Labs, Worcester, MA, USA) [42]. This makes the MPDV1.0 a different and innovative way of estimating vulnerability. Another difference from the existing methodologies is that it allows quantifying five types of vulnerabilities (i.e., hydrological, environmental, economic, social and global). The MPDV1.0 executes the connected tools and the input maps following the flowchart shown in Figure 3. The input maps have been processed in ArcGIS (ESRI, Sacramento, CA, USA), considering the available information, obtained from different governmental departments. In the research, the error on the model and calculation of vulnerability is not estimated because the input maps are official information in Mexico and are carefully validated. These maps are introduced in ascii format in the MPDV1.0. The MPDV1.0 conducts information cross-checks to obtain hydrological, environmental, economic, social and global vulnerability.
The MPDV1.0 calculates global vulnerability, based on the IPCC's conception of vulnerability [43], which generally explains vulnerability by applying the following equation: where GV is the global vulnerability, HV is the hydrological vulnerability, EV is the economic vulnerability, SV is the social vulnerability, AV is the environmental vulnerability and n is the number of vulnerabilities. The vulnerabilities of the Equation (1) are obtained through the indicators shown in Table 1. Each vulnerability map will be calculated from the normalized indicators involved and their corresponding weights with the following equation: where X i is the normalized indicators i, w i is the weight of the normalized indicator i and n is the number of indicators for each type of vulnerability. The indicators are normalized to eliminate the units and to vary the values from 0 to 1. The equations implemented to perform the standardization are the following: where x i , x min and x max are the minimum and maximum values of the data set x i . If the indicator has the capacity to adaptation, Equation (3) shall be used; and If the indicator has a degree of exposure or sensitivity, Equation (4) shall be used. For each normalized indicator, a weight will be obtained using the following equation: where σ 2 i is the standard deviation of the set of the indicator values i, and n is the number of selected indicators.
In this research, the MPDV1.0 has been applied to obtain present and future vulnerability maps (for the year 2035). The foregoing is with the purpose to assess the potential impacts of climate change and land use change.

Data Sets
The geographic information used in estimating the indicators (Table 1 has been obtained from the National Institute of Statistics and Geography (INEGI) available on https://www.inegi.org.mx/ and from the National Water Information System of the National Water Commission (SINA-CONAGUA) available on http://sina.conagua.gob.mx/sina/. They are governmental departments that allow free downloading of information at scales of 1:50,000 and 1:250,000. Table 1. Equations and sources of information used in the calculation of indicators [44].  Hydrological modeling calibration has been conducted using daily flow data measured at Las Adjuntas gauge, which is located at the mouth of the sub-basin (Figure 1). These hydrometric data are also freely downloaded from the website of CONAGUA National Data Bank of Surface Waters (https://www.imta.gob.mx/bandas). In the case of precipitation data, it was initially tested using daily data from 20 weather stations. However, the first results showed some inconsistencies in the interpolation of vulnerability maps, because spatial interpolation of daily rainfall gauge data using Inverse Distance Weighting in the TETIS model has shown errors of overestimation in rainfall interpolated in the cells adjacent to the border of the basin. Moreover, the results showed the formation of bull's-eye patterns produced by the spatial location of the weather stations. That is, the isohyets are circular in the stations without nearby neighbors, denoting a poor spatial distribution of the weather stations. To support the above, we decided to use the Multi Source Weighted-Ensable Precipitation (MSWEP) product from Princeton Climate Analytics, available on https://platform.princetonclimate. com/PCA_Platform/index.html. MSWEP is a precipitation mesh with a spatial resolution of 0.25 • and a temporal resolution of three hours [45][46][47][48]. MSWEP data have been correlated with daily rainfall, observed at weather stations for the period from October 2010 to September 2014. The foregoing is with the purpose of validating the product data and justifying its use. Figure 4 shows some examples of the obtained results; in general terms, the correlations from 0.5 to 0.8 were obtained. In the case of satellite rainfall products, correlation coefficients greater than 0.5 are considered acceptable [9,49]. Despite the overestimates observed in Figure 4, we decided to use the MSWEP, because the interpolated rainfall from the weather stations presents 13.5% overestimates in cells far from the stations compared to the MSWEP. We also evaluated in the analysis period the differences in the annual average accumulated precipitations and we obtained that the observed precipitations are greater than those of the MSWEP with an average of 175 mm year −1 . There are other products with precipitation information, such as: NOAA CPC Morphing Technique (CMORPH), Integrated Multi-satellite Retrievals for GPM (IMERG) and Climate Hazards Group Infrared Precipitation with Stations (CHIRPS). Studies such as those carried out by [46,[50][51][52] consider these high-precision products and have been implemented in various regions, obtaining correlations of 0.65 to 0.8. However, the information period for these products is limited. Therefore and considering the MSWEP has a longer reporting period and that similar correlations are obtained, we have decided to use it in the research.

Global
In the hydrological modeling, 36 cells of the MSWEP product were selected; their centroids have been named "virtual stations". By interpolating the precipitations of these stations in the period from October 2000 to September 2017, we obtained the distribution of accumulated average annual precipitations shown in Figure 5a. Figure 5b,c show the distribution of elevations in the sub-basin, the variation of precipitation in the high areas is due to an effect of latitude.

Automatic Calibration of the TETIS Model
The TETIS model uses nine initial parameters (Table 2), Figure 6 shows the three most important parameters, calibrated automatically through Correcting Factors (FC). For more information on the initial parameters of the TETIS model, see for example [34,35,53].
The optimization algorithm used is the Shuffled Complex Evolution (SCE-UA), developed at the University of Arizona [13,54,55]. The SCE-UA obtains better performances, given that it requires less iterations to reach optimum and it is more suitable for non-expert users, while Simplex, which is a more efficient algorithm, needs a step-by-step process supervised by an expert modeler. We decided to use as an objective function, the Nash-Sutcliffe Efficiency Index (NSE), because it is the most used in hydrological model calibration [56]. The NSE is calculated using the following equation: where q obs is the observed flow, q sim the simulated flow andq obs is the mean flow observed.

Downscaling Climate Model
In this research, we decided to use the CMIP5 GCMs because they are the most used models from the end of the eighties [57,58]. Moreover, according to [59], the CMIP5 models have a better performance to represent numerous phenomena include extreme precipitation, diurnal variations in precipitation and tropical cyclone intensity. Studies such as the one carried out by [60] have shown that climate models in CMIP5 are better in future projections, because the efforts of CMIP5 are enormous, with a larger number of more complex models run at higher resolution, with more complete representations of external forcings, more types of scenario and more diagnostics stored. The projections with the GCMs have been obtained using the Climatic Atlas platform developed by the Center for Atmospheric Sciences of the National Autonomous University of Mexico and are available on https://atlasclimatico.unam.mx/AECC/servmapas. This platform was developed exclusively for Mexico considering the models shown in Table 2. It allows calculating projections using each model individually and using the 15 models together. We decided to use them together through the Reliability Ensemble Averaging (REA) conducted by [61]. According to [62], the REA reduces the uncertainty in the projections by giving greater weight to the GCMs that present the smallest errors and biases with respect to the variables observed in a given mesh point. According to [61], the weight assigned to each model is based on two criteria, (1) trend (difference with observations) and (2) convergence (difference between simulations).
One of the problems of the GCMs is that they have a very large resolution. To solve this problem, the authors [61] divided Mexico into four regions: Northwest, Northeast, South and Southeast to compare the results of the REA with the data from the Climate Research Unit (CRU) [63]. The CRU database is available on http://www.cru.uea.ac.uk/data. In addition to the downscaling performed by the previous authors, we decided to apply a downscaling to the REA projections for the South region (specifically in the state of Guanajuato, which is the region where our study sub-basin is located). Downscaling has been made by combining statistically observed information for the 1971-2000 period, data from the CRU, REA and RCP4.5, RCP6.0 and RCP8.5 scenarios [64]. The authors [61] obtained correlations between REA and CRU higher than 0.8 for the South region. Moreover, significant reductions in the Mean Absolute Error (MAE ≤ 0.9 mm day −1 ) and Root Mean Square Error (RMSE ≤ 1.3 mm day −1 ). We obtained for the state of Guanajuato a correlation coefficient of 0.78, MAE ≤ 0.99 mm day −1 and RMSE ≤ 1.5 mm day −1 . Figure 7 presents the examples of maps with monthly anomalies obtained with downscaling for the entire state of Guanajuato. These anomalies are introduced as inputs to the TETIS hydrological model to simulate the effects of climate change.

Forecast of Land Use Changes
For the prediction of future changes in land use and vegetation cover, it has been decided to use the Terrset software [42]. TerrSet is an integrated geospatial software system for monitoring and modeling the earth system for sustainable development. We decided to use Terrset because it uses Markov Chain analysis to project the expected quantity of change and a competitive land allocation model to determine scenarios for a specified future date [65,66]. In addition, the class assignment is based on the contiguity of elements through multicriteria-multiobjective. Its main advantage is that it models with a stronger approach than the other models, since it models the suitability for transition rather than suitability for the ultimate cover class. The technique produced projections maps with an overall accuracy of more than 90% [67]. More information can be found at the link https://clarklabs.org/terrset/land-change-modeler/.
The main input information used in the modeling was the land cover and land use maps at a scale of 1: 250,000, obtained from INEGI. The land use map year 1997 was used in the construction of the Terrset model (Figure 8a), and the VI series (year 2016) was used for the validation of the Terrset model results (Figure 8b).

Hydrology Model Performance
The TETIS model adequately reproduces the flows observed at Las Adjuntas gauge for the periods used in calibration (from October 2012 to September 2014) and validation (from October 2005 to September 2008). We decided to use only two hydrological years in the calibration, because there is little information and since according to [68] there is no direct relationship between the number of years and the results of the models. Moreover, that computing times are significantly reduced. Consequently, in the selection of the calibration event, it was necessary to have a wet and a dry year to obtain robust parameters. In terms of efficiency, the model has achieved an average NSE of 0.67, thus, simulating a wet year and a dry year (Figure 9). Calibration is acceptable considering that according to [69], the value of NSE are considered acceptable if it is greater than 0.6 and excellent if it is greater than 0.8. Table 2 shows the average effective parameters obtained by automatic calibration with the SCE-UA algorithm. SMITH2013300.  Table 3 shows the average effective parameters obtained by automatic calibration with the SCE-UA algorithm, calibrated by correction factors (FC), which make a global correction of the initial parameters quantified with the information available [53]. Table 3. Average effective parameters obtained by automatic calibration with the Shuffled Complex Evolution (SCE-UA) algorithm.

Parameter Correction Factor Parameter Equation Effective Parameter
Static storage Vegetation cover index FC 2 = 0.500 Overland flow velocity Percolation capacity Interflow velocity FC 6 = 0.010 k * ss(i) = FC 6 · k ss 0.003 (mm h −1 ) Deep aquifer permeability Connected aquifer permeability River channel velocity The performance of the model in temporary validation, once again, shows an acceptable behavior of the TETIS model, with an NSE of 0.63 ( Figure 10). The NSE obtained in validation shows an acceptable capacity of the model to reproduce the flows in another period, which is different from the one used in the calibration. According to [69], validation is usually considered acceptable if NSE > 0.5 and very good if NSE > 0.7.

Climate Projections and Downscaling
When downscaling the results obtained by the GCM of CMIP5, it was possible to observe a significant reduction in precipitations of the year 2035 for the three scenarios of IPCC ( Figure 11). The foregoing depends on the control period used in the comparison. This control period has been generated with historical precipitation data from the MSWEP (period: 1982 to 2014). We have obtained that the REA-CMIP5 models predict an average decrease in rainfall of 1 mm day −1 for the IPCC scenarios. The greatest differences between future projections are presented in the months of the rainy season. The high variability observed in Figure 11 corresponds to the uncertainty calculated for the three scenarios, which is ±2 mm day −1 . We find the same trend in evaluating annual projections ( Table 4). The same pattern was observed with annual anomalies, which predict an average decrease in rainfall of 34.7 mm year −1 . The results in Table 4 show that in the near future, the RCP4.5 scenario represents the most critical conditions with a greater decrease in rainfall. However, the projections of REA-CMIP5 in the year 2100 show that RCP8.5 is the most critical with an average decrease in rainfall of 80 mm year −1 (period: 2075-2100).  The projections show a general increase in maximum temperatures of up to 2 • C ( Table 5). This will have a direct effect on the process of evaporation, evapotranspiration, and consequently, on hydrological availability. The estimated uncertainty for these projections is around ±1 • C for the three IPCC scenarios. The above uncertainty coincides with the results obtained by [61] for the South region of the country. As expected, the most critical scenario for the near future has been RCP8.5. Finally, the REA-CMIP5 project for the year 2100 and the same scenario an increase in temperatures of 4.1 • C.

Future Land Use Change Scenarios Modeling
The validation of land use change projections up to the year 2035 presents 85% of successes, when comparing the vegetation cover of the VI series ( Figure 7) and the Terrset results. Once the model was validated, the land use change scenario was generated for the year 2035 ( Figure 12). In the projection, it can be observed that the greater changes are seen in the crop land class ( Figure 12). Likewise, a significant increase is observed in the urban nuclei. The authors decided not to include the distant future (year 2095) in the research, because we observed inconsistencies in the Terrset prognosis when we distanced from the validation period. In the year 2095, Terrset predicts a decrease in urban areas of 30% compared to 2035 and the occupation of this area by irrigation agriculture, which is not logical considering patterns of large cities growth trend over the last years.

Assessment of Vulnerability
The results obtained with the proposed methodology and the developed MPDV1.0 show that the Turbio river sub-basin presents high vulnerability for the current conditions ( Figure 13). The greatest vulnerabilities have been obtained in urban centers and agricultural areas. This coincides with what has been reported in some studies, such as: [70,71]. We decided to analyze more exhaustively the results of hydrological modeling, and observed that infiltration and percolation processes were reduced to 34.5% and 50%, respectively, in urban areas. In addition, 83% of groundwater was extracted in urban centers and urban areas. All this enables us to conclude that the planning and management of the current hydrological availability is not as equitable and sustainable as it was expected. The previous results may change in the cells located in the valley of the Turbio river sub-basin, since in these cells the MSWEP overestimates rainfall by 0.88 mm day −1 . That is to say, the overestimations of rainfall are diminishing the effect of groundwater withdrawals in the formulations of the indicators involved, possibly reducing the vulnerabilities in the cells of the valley of the sub-basin.
The MPDV1.0 forecasts for the near future a generalized increase in vulnerabilities to the potential impacts of climate change and future land use changes (Figures 14 and 15). We expected this increase due to the decrease in rainfall (from 4.97% to 12.32%), predicted by the CMIP5 GCM. There is a potential effect of climate change on evapotranspiration processes, infiltration and percolation. Other authors have also detected these effects for other areas of the world (e.g., [72,73]).
By analyzing the vulnerability for each IPCC scenario, we can observe patterns influenced by the characteristics of the scenarios. In general, the greatest vulnerabilities are presented in the RCP6.0 and RCP8.5 scenarios, which surpasses the rest of the scenarios with 0.1% in the sub-basin area and with vulnerabilities of more than 0.5. We believe that the reason for this small increase is because by definition the most drastic scenarios in 2035 have a gradual rise in effects. Likewise, the increase in the use of urban centers and agricultural zones directly influences on the increase of these vulnerabilities. Another interesting conclusion is that according to the hydrological vulnerability rate, the worst conditions for the system are expected in water availability. In other words, practically 75% of the system area will be very vulnerable, with the possibility of not being able to meet the demands of different consumptive uses. The above results are giving indications that there is an urgent need for a restructuring of the management policies that are currently being implemented.

Conclusions
In this research, we assess the global vulnerability through hydrological, environmental, economic and social vulnerabilities, taking into consideration the potential impacts of climate change and land use change. In order to achieve this, a new methodology was designed and a Predictive Vulnerability Interface was developed. The results obtained in the study sub-basin show that under current conditions the system is already vulnerable. In the sub-basin, environmental and economic vulnerability is the most important. The areas with the greatest environmental vulnerability were urban centers and agricultural areas, as these areas are where the greatest water vulnerability occurs. This makes evident a restructuring of the management policies that are currently being implemented in the country. In the near future, the situation shows small changes influenced more by the change in land use than by the effects of climate change. The foregoing may be due to the fact that the same degree of exploitation of the sub-basin and the groundwater of the current conditions is considered. This is a non-conclusive result, which will need improvement and further testing. Furthermore, it is necessary to evaluate the uncertainty to detect and improve the variables that introduce the errors in the model. However, if an effect is observed in the sub-basin for each IPCC scenario. This is due to the fact that the Global Climate Models predict a decrease in rainfall and a significant increase in temperature. By integrating these effects of climate and land use changes into the Vulnerability Predictive Interface, there are the combined effects in the flow regimes, producing negative impacts on water resources; thus, a decrease in water availability in the aquifer is predicted. This gives evidence of the great importance of having a methodology such as the one proposed for planning and managing current and future water availability for the region.