WEAP21 based modelling under climate change considerations for a semi-arid region in southern-central Chile

The water balance of an irrigation project in the Chilean semi-arid zone was modelled using the WEAP21 code. A combined deterministic/stochastic protocol was applied for defining the hydrological prediction intervals of the model. Climate Change potential implications on the useful life of the irrigation project (year 2070) were predicted considering the variability in precipitation and temperature forecasts. Thereto, it was considered the scenario A1B of the 2007 Intergovernmental Panel on Climate Change report, generated by the Regional Climate Model PRECIS-ECHAM, specially developed for Chilean conditions. The study revealed that most of the inspected hydrological model parameters are insensitive to model predictions and the associated simulation limits may be categorized as acceptable. Nevertheless, the structure of WEAP21 had difficulties representing low flows because of the apparent inability to mimic the complex hydro-physical characteristics of the shrink-swell granitic soils which are predominant in the study basin. Even though the original Climate Change projections (CChP) of the RCM were refined, using observations of the historical period, it is important to underline that significant uncertainties may remain and as such the current results should be handled with care. With respect to historical records, mean annual climate forecasts suggest a maximum temperature increment of about +1.1oC and a maximum reduction in precipitation of 20.7%. The hydrological modelling suggests a maximum reduction in the mean annual streamflow of 49.7% and a reduction in the magnitude and frequency of streamflow peaks. Bearing in mind the potential uncertainties attached to CChP, the irrigation project will most probably be significantly affected in terms of water availability and crop water consumption since rainfall is expected to decrease and temperature, and as such evapotranspiration, to increase.


INTRODUCTION
Studying the effect of climate change on water resources is important because the change in climate leads already to devastating impacts across the planet (Sanderson, Jaitech, Levy, Redford, Wannebo et al., 2002;Barnett, Pierce, Hidalgo, Bonfils, Santer et al., 2008), which are expected to worsen due to the future population growth (Palmer, Bernhardt, Chornesky, Collins, Dobson, et al., 2004;Wagener et al., 2008).Easterling, Aggarwal, Batima, Brander, Erda et al. (2007) state that global warming affects crop and pasture yields with a positive effect in temperate regions and a negative impact in semi-arid and tropical regions, and this as a result of changes in the regime of both precipitation (Palmer et al., 2004;Wagener et al., 2008) and temperature.The expected changes in climate will make agricultural zones very vulnerable (Deressa, Hassan, & Ingler, 2008), and justify examining its impact on the management and planning of irrigation and drainage projects.
However, this is not a trivial task since there is a significant uncertainty (Hingray & Saïd, 2014;Willems, Mora, Vansteenkiste, Teferi Taye, & Van Steenbergen, 2014) in climate change projections (CChP) using Global Circulation Models (GCMs).Consequently, ensembles of GCMs are preferred, rather than the use of a single one, in conjunction with statistical downscaling or downscaling through multivariate stochastic and decision support techniques, for assessing the level of uncertainty on the different projections (Makropoulos & Butler, 2004;Hublart, Ruelland, Dezetter, & Jourde, 2013;Willems et al., 2014).An alternative is dynamic downscaling by means of Regional Climate Models (RCMs) nested within a GCM (Christensen, Carter, Rummukainen, & Amanatidis, 2007;Fowler, Ekström, Blenkinsop, & Smith, 2007), which has important advantages over GCM-based scenarios as the physical processes are better represented and at a higher resolution (Dankers & Feyen, 2008).In addition, an ensemble of RCMs enables estimation of the uncertainty level for the different selected projections.
This paper presents the results of a water balance model used for the planning of an irrigation project located in a Chilean river system with intermittent flow conditions.To simulate water availability Revista semestral de la DIUC 127 https://doi.org/10.18537/mskn.08.02.10 under historical and climate change conditions, a lumped hydrologic model was applied, on a monthly basis.The modelling protocol included a deterministic/stochastic procedure, based on Monte Carlo simulations, for determining the prediction limits.As top boundary a single CChP was used derived from a RCM that was developed for Chilean conditions.Although this projection has important advantages over a GCM-based one (i.e., the physical processes are better represented and at a higher resolution), some uncertainty is likely to remain, which was difficult to assess owing to the absence of other RCMs.Nevertheless, as an attempt to remove as much as possible the referred uncertainty it was decided to run additionally a bias-removal process.
The hypotheses that led to this research are that: (i) modelling of the monthly water balance associated to the planning of an irrigation project is possible applying a lumped model within the context of a deterministic/stochastic protocol; (ii) the hydrological model is capable of acceptable estimating streamflow time series at the location of the projected irrigation dam, where no observations are available, despite the fact that the predominant soils are shrink-swell clays, difficult to model; and (iii) climate change implications on water availability can be assessed fair by using the CChP of a RCM after bias noise has been removed through a downscaling approach.
In this context, (i) the application of a lumped model was encouraged since the vertical gradients of the meteorological variables in most of the study area are not that significant; and (ii) the assessment of the prediction limits of the hydrological model is an innovative element in the analysis in comparison to similar studies that assume uncertainty only be present in the CChP and not in the hydrological model.

The study site
The study site is in southern-central Chile in the Bío-Bío region, an agroecological zone known as the "Secano" (dry land), characterized by semi-arid Mediterranean climate and granitic soils.Both features made that the region evolved to a non-irrigated agricultural area characterized by significant fragmentation of the land, low agricultural output, and significant soil erosion as a consequence of poor soil management practices (Shinomi & Pérez, 2004).Because of these constraints, combined with the irregular topography and the low investment capacity of the farmers (Grau, 2007), missed the zone the benefits of the significant agricultural development that Chile exhibited in the last 30 years.To reactivate the agricultural economy, the Chilean government has been studying development projects, and herein, is the Lonquén irrigation project proposed for the communities of Ninhue, San Nicolás and part of San Carlos (Fig. 1).Considering that Ninhue is one of the poorest communities of Chile (Bengóa, 2003), the irrigation project will be a very important leverage for the improvement of the community's well-being.
An obstacle for the modelling study is that no hydrological gauging station is present near the spot where the irrigation dam is projected in the Lonquén river (Fig. 1).As a consequence, the design discharge time series had to be generated through the hydrological modelling of the upper contributing subcatchment, considering the streamflow information from the "Lonquén" and the "Trehuaco" gauging stations, both located downstream of the dam (Fig. 1).The contributing area of the projected irrigation reservoir, the "upper subcatchment" (LU 2) is 297 km 2 , whilst the corresponding contributing area of the gauging stations, the entire "catchment" has a surface of 1178 km 2 .The surface difference between both areas, being the surface of the "lower subcatchment" (LU1) is 881 km 2 .
The Lonquén River with an approximate length of 686 km (until Trehuaco), originates in the mountains of the coast and is one of the major tributaries of the Itata River.The catchment is covered by grassland (40.3%), shrubs (34.9%), forest (16.6%), and crop land (7.4%).The elevation of the study catchment, with an average slope of 16%, ranges from 22 to 902 m above sea level (a.s.l), and nearly 90% of the area is situated below 300 m a.s.l.Around 80% of the annual rainfall in the basin occurs between March and August (DGA, 1987).The isohyet method (Chow, Maidment, & Mays, 1988) suggests that the average annual rainfall is approximately 868 mm, with a dry season in the period November-March (final part of the spring, the whole summer and the first part of the autumn) and a wet season in the period April-October (final part of the autumn, the whole winter and the first part of the spring).The average temperature is about 13.9 o C. The areal crop reference potential evapotranspiration (ET0) upstream of the projected dam, calculated using the Penman method, is 1290 mm y -1 (CNR, 1997).According to the Chilean Water Directorate (DGA, 1987) is the average annual pan evaporation in the order of 1500 mm.
The soils of the catchment are characterized by a low infiltration rate, low permeability, and high water retention capacity due to the presence of shrink-swell clays and the high rate of degradation.The predominant soils are considered moderately deep; however, deep soils are also present on the granite and metamorphic sediments.According to its texture is the surface layer of the soils (up to a depth of 30 cm) classified as clay-loam, and the deeper layers as muddy-clay-loam (CIREN, 1998;Riquelme, Perez, & Shigehiko, 2004;SMI, 2012).The soils shrink in the dry season to the point that the structure of the top layer changes, and swells in the wet season; an annual cycle (Riquelme et al., 2004).The soil properties can cause high runoff rates in the wet season and close to nil rates in the dry season; runoff rates grow and decline quickly as a function of precipitation events (Selker, Rupp, Leñam, & Uribe, 2000;Uribe, Perez, & Okuda, 2004).These characteristics also hinder groundwater recharge whose values are estimated to range from 1 to 6% of the annual rainfall (Uribe, Arumí, González, & Salgado, 2003).Chemical and mechanical weathering, and the effects of the vegetation and climate, resulted in an altered top layer of the bedrock.The transition area between the altered layers and the underlying unaltered rocks contributed to the formation of several shallow wells, which by lack of surface water are considered the main water supply (Uribe et al., 2003).

Hydro-meteorological data
Hydrometeorological records were provided by the DGA.The Lonquén catchment upstream the Trehuaco outlet is equipped with one discharge and two rainfall stations (Fig. 1).The moderate conditions of the discharge gauging station for measuring low flows in the dry period suggest that the historical hydrological information is affected by significant uncertainties.Historical records of meteorological variables, such as temperature or wind, were collected at neighbouring stations by lack of a meteorological station within the study catchment.The main characteristics of the stations from where data were borrowed are listed in Table 1.Additionally, rainfall records were also taken from neighbouring gauging stations.The main characteristics of those stations are included in Table 1.
The average low flow in the Lonquén River during dry periods varies around 0.32 m 3 s -1 ; although there are historical dry periods where no discharge was observed.The average discharge during the wet season is approximately 20.6 m 3 s -1 , whilst the mean annual discharge is 12.3 m 3 s -1 .Consideration of these values suggests that the river does not receive important groundwater contributions, although further research would still be needed to confirm this assumption.Several studies in the basin revealed that the first precipitation events of the year, up to the first 100 mm of the annual rainfall, completely infiltrate and only thereafter is run-off produced (DGA, 1987;Selker et al., 2000).Figure 2 depicts the average seasonal discharge of the Lonquén River at the Trehuaco gauging station, for the period April 1986-March 2011.The figure clearly shows a significant difference in streamflow between the wet and the dry seasons.For most of the studied years, the mean discharge of the dry season is nearly zero; indeed, only in four hydrologic years a mean discharge significantly different from zero was recorded.Further, in six years of the available annual records the maximum discharges were in the order of 5 m 3 s -1 .

The digital elevation model of the study site
The Shuttle Radar Topography Mission 90-m resolution digital elevation model (DEM) was used to define the topography and assess the spatial distribution of the temperature within the catchment.The STRM 90 m DEMs have a resolution of 90 m at the equator and provide elevation data for over 80% of the world (CGIAR-CSI, 2008).

The WEAP21 model
The Water Evaluation And Planning System (WEAP21), widely used in Chile for climate change impact studies (Vicuña, Garreaud, & MacPhee, 2010;Poblete, Vicuña, Meza, & Bustos, 2012;Cortés, Schaller, Rojas, Garcia, Descalzi et al., 2012;Vargas, 2012;Ragettli et al., 2014), was selected for the modelling of streamflow and the water availability at the outlet of the upper subcatchment (LU 2) using respectively historical records and future weather estimates.An additional objective of the study was to evaluate whether the fixed structure of this model can perform well the modelling.The WEAP21 model, A Demand-, Priority-, and Preference-Driven Water Planning Model, simulates consecutively the water demand and supply, runoff, infiltration, crop water requirements, the water flow and storage, among other aspects (Yates, Purkey, Sieber, Huber-Lee, & Galbraith, 2005;Ingol-Blanco & McKinney, 2009;CCGUC-SEI, 2009;Sieber & Purkey, 2015).The hydrological module is based on a conceptual soil water content method (Yates et al., 2005) that represents the rainfall-runoff relationship by means of two linear reservoirs accounting explicitly for surface flow, interflow and baseflow, as illustrated in Fig. K1 and K2 are the hydraulic conductivities [mm month -1 ] in the upper and lower reservoirs.K1 and K2 represent the maximum possible water fluxes at full storage, when Z 1 and Z2 equal 100%.However, they should not be considered the saturated hydraulic conductivities in the strict sense (Rawls, Ahuja, Brakensiek, & Shirmohammadi, 1993;Yates et al., 2005).Further, Df is the preferential flow direction [-], which splits groundwater flow from the surface flow; and FR is the runoff resistance factor of the land cover [-].
A catchment can be subdivided into several fractional areas representing different land uses/soil types.The balance equations in every one of the reservoirs for every fractional area can be estimated using the Eqs.(1 & 2): where Z1max is the upper reservoir capacity [mm]; Pe is the rainfall and snow melt [mm]; and Z2max is the lower reservoir capacity [mm].ET 0 [mm d -1 ] is estimated through an implementation (Sieber & Purkey, 2015) of the Penman-Monteith equation (Maidment, 1993;Allen, Pereira, Raes, & Smith, 1998).Although this equation is not reported herein to keep the extent of the current manuscript within reasonable limits, it is important to mention that it is a function of, among other aspects, the net solar radiation (R n), the daily mean temperature (T), the wind velocity measured at 2 m above the ground surface (u2), and the saturated (es) and current (ea) vapour pressures.In Figure 3 the FAO-56 method (Allen et al., 1998) is used to estimate the crop potential evapotranspiration (ET c = Kc x ET0), which to define ET r is then affected by a factor, that in turn is depending on the current storage in the upper reservoir.

Initial parameterisation of the hydrological model
On the basis of the requirements of the irrigation project, a monthly hydrological model was established for the study site with the outlet at Lonquén, nearby the Trehuaco station (Fig. 1).The study catchment was represented by means of two (lumped) units, the lower subcatchment (LU 1 ) and the upper subcatchment (LU 2), and this with the intention of obtaining reservoir operation discharge time series WEAP21.Table 2 provides a description of the tuned parameters with their respective lower and upper limit.

Model calibration and validation
The calibration and validation periods were defined based on the hydrological year, respectively the period April 1986-March 2005 and April 2005-March 2011.Furthermore, a warming-up period was considered between April 1979 and March 1986.Given the extent of the warming up period, it is expected that the effects of the assumed initial conditions were reduced to a minimum, not affecting at all the modelling during the calibration and validation periods.
The Generalised Likelihood Uncertainty Estimator (GLUE) methodology, a deterministicstochastic approach, was used for the calibration and evaluation of the model.The GLUE method is based on the "equifinality" concept, which states that there are several model structures with different parameter sets that may be acceptable in reproducing the observed behaviour of the modelled system.The approach has been widely used by Binley, Beven, Calver, & Watts (1991), Beven & Binley (1992), Beven (1993), Beven (2006), Beven (2011), among other authors.
GLUE calculates prediction limits through Monte Carlo simulations that are defined using multiple parameter sets that are sampled out from the feasible parameter space on the basis of prior parameter density distributions.Model performance (i.e., likelihood measure) is then evaluated for every tried parameter set comparing predictions and observations.The parameter sets that produce predictions that exceed a specified threshold of performance are considered as behavioural and, as such, are retained for further consideration.Finally, the distribution functions of both the model parameters and the predicted variables can be calculated for a given level of confidence (Binley et al., 1991;Beven, 1993).
GLUE enables the integration of additional information into the error assessment by means of a Bayesian type approach to update the likelihood weights and estimated prediction limits, as used in this study: ( where Lo(Ω i) is the prior likelihood distribution of the i-th parameter set; LO(Ωi| O ) is the likelihood measure of the i-th parameter set, provided the new observations (O) and computed in the newer period of observations; Lp(Ωi | O ) is the posterior likelihood distribution of the i-th parameter set; and cGL is a scaling constant to ensure that the summation of the posterior likelihood measure of the behavioural simulations is equal to one.In the current research, it was decided to use in GLUE the Bayesian type updating of likelihoods.
Reflecting the lack of prior knowledge on the distribution associated to model parameters, independent uniform distributions were assumed for all of the studied parameters.For each run of the model a new set of parameter values was sampled.A total of 100,000 parameter sets were sampled, which can be considered as an appropriate number of Monte Carlo simulations according to the number of model parameters in the current analysis.
The likelihood measure that was computed to characterise the model performance is proportional to the Nash-Sutcliffe efficiency coefficient (Nash & Sutcliffe, 1970) and evaluated using a square root transformation of streamflow (EF 2sqrt).The latter is supposed to be less sensitive to the simulation of peak flows than the standard EF2 Nash-Sutcliffe coefficient (Legates & McCabe, 1999;Vázquez, 2003) and therefore more suitable to assess the performance of the general representation of the flows.
Thus: Q are respectively the observed and the simulated flow at step i, obs Q is the mean observed discharges and N is the total number of observations.EF2sqrt varies between -∞ and 1.This path was followed in the current research given it offers an acceptable measure of the combined systematic and random error.
The behavioural sets were defined by considering an efficiency threshold or cut off criteria (EF 2sqrt>0.5)for the calibration period.A likelihood weight was then determined for each behavioural parameter set on the basis of the EF2sqrt index.Upon these prior likelihood distributions, a first set of prediction limits (or prediction band) was defined.The WEAP21 structure was run again to obtain simulation results for a second (i.e., validation) period considering only the behavioural sets defined in the prior calibration period.The value of this additional information in the validation period, on refining the prediction limits, was also inspected by means of the Bayes type equation ( 3) that was used to update the prior likelihood distributions.

Historical meteorological data
Spatial interpolation methods were used to reconstruct the meteorological features in each of the two study subcatchments, given the meteorological stations of which data were used are located outside the limits of the study zone (Fig. 4a).Temperature was considered to be a function of elevation, although the vertical gradient is not that significant in most of the catchment area (Fig. 4b).The latter enabled estimating the multiyear mean monthly temperature for a given elevation spot and a given month of the year, for which the catchment extent was subdivided into 5 elevation bands/strips with similar area (Fig 4c ).
Revista semestral de la DIUC 134 https://doi.org/10.18537/mskn.08.02.10 The spatial discretization was achieved using Eq. ( 5 where, for month i and year j, Tastrip i,k and Tabi are the multiyear mean monthly temperature for the kth elevation strip and the base gauging station (i.e., a station that is used as the reference for the calculations), respectively; Tstrip i,j,k and Tbi,j are the mean monthly temperature for the k-th strip and the base gauging station, respectively.It should be noticed that Eq. ( 5) assumes for the base gauging station a direct proportionality between the multiyear ratio (Tastrip i,k /Tabi) and the current month ratio (Tstrip i,j,k /Tbi,j).The Parral station was chosen as the base gauging station due to its proximity to the study catchment.Finally, the mean monthly temperature in every subcatchment was calculated by means of an area-weighted average of the elevation strips present in every subcatchment.
As for the distribution of rainfall, the isohyetal method was used to derive mean annual values, since no vertical gradient was detected in the observed records.The precipitation isolines (Fig. 4d) were calculated with the Ordinary Kriging interpolation method considering the whole records of every station.Then, the monthly values in each subcatchment, for the different years of interest, were calculated by means of a reciprocal expression, that is similar to the one depicted in Eq. ( 5), among (i) the ratio between the mean annual rainfall over a given subcatchment and the respective annual rainfall recorded at a base (gauging) station in that subcatchment; and (ii) the ratio between the monthly rainfall (for a given year) in the subcatchment and the respective monthly rainfall recorded at the base station.The San Agutín de Puñal and the Mangarral stations (Table 1) were respectively the base (gauging) stations for the LU 1 and LU2 subcatchments.
The monthly wind values for the whole basin were calculated as a simple average of monthly values of the Digua Embalse and Coeihueco Embalse gauging stations.As for the distribution of relative humidity, mean monthly values were derived from the agroclimatic cartography of the Natural Resource Information Center of Chile (CIREN).Solar radiation (R s) was calculated using the the Angstrom formula (Allen et al., 1998), function of the relative cloudiness (n/N) and the extra-terrestrial radiation (R a).The relative cloudiness was obtained and validated by Petersen (2007) for the vicinity of the study catchment.Ra was derived on the basis of the latitude of the study site (Allen et al., 1998).Finally, it is important to state that neither basin transboundary contributions, nor snow processes are present in the study area.

Climate change data
For the generation of future climate projections monthly data of precipitation and temperature (forcing variables) were used as input in the PRECIS-ECHAM (http://mirasol.dgf.uchile.cl/PRECIS-ECHAM.html)model, a Regional Climate Model (RCM) developed for Chilean conditions by the Department of Geophysics of the Universidad de Chile.This RCM was the first dynamic downscaling exercise done in Chile, has a resolution of 25 x 25 km 2 , covers the whole extent of the country and integrates PRECIS (Providing Regional Climates for Impact Studies) whose boundaries were forced by the Global Circulation Model (GCM) ECHAM5.Further, it covers the period from 1950 to 2100, with a "historical period" from 1950 to 2000 and a "future period" from 2001 to 2100.As from the year 2001 on, the model simulates greenhouse gas emission according to the A1B scenario of the 2007 (4 th ) Intergovernmental Panel on Climate Change report (IPCC, 2007).This model and specific scenario have been used in the current study since it was developed particularly and thoroughly to support studies in the Bio-Bio region where the current study site is located (Rojas, 2013).
Even though RCMs have better spatial resolutions (Dankes & Feyen 2008) than GCMs, they are likely to contain biases that must be removed or at least acknowledged (Graham & Bergström, 2001).This bias can even be bigger when the assessment is at a local scale, such as in the current study.Therefore, a bias correction process for local scale, based on the Bias Correction and Spatial Downscaling (BCSD) method (Wood, Leung, Sridhar, & Lettenmaier, 2003;Wood & Edwin, 2002;Maurer, 2007) was applied.
The BCSD method was originally designed to process monthly data of distributed climatic models on three datasets: (i) observed climatic data; (ii) GCM (or RCM) forecasts for the historical period; and Revista semestral de la DIUC 135 https://doi.org/10.18537/mskn.08.02.10 (iii) GCM (or RCM) forecasts for the future period.In this method, the Bias Correction (BC) is applied on the first and second datasets, for which, both observed data and GCM (or RCM) forecasts are rescaled to a common spatial resolution.Cumulative Distribution Function (CDF) curves are generated at each grid cell for both datasets (since both datasets are spatially distributed), by plotting the sorted values versus the rank probabilities.Then, the CDFs for both observed and modelled meteorological data are related through probability threshold to define a quantile map.This quantile map is finally used to remove the bias from the third dataset (future period).It is worth noticing that in the BSCD method there are implicit assumptions in the sense that (i) the topographic and the climatic features determining the fine-scale distribution of large-scale climate will remain unchanged in future; as well as (ii) the climate model biases in both present and future periods follow the same pattern.
The original BC method is normally applied on distributed climatic models.However, in the current study the spatial distribution of forcing variables is based on a meteorological station.Thus, the BC method was applied only at the base (gauging) station for each forcing variable.Here the closest PRECIS-ECHAM forecasts to the locations of the base stations of interest were used.The base stations used in this study are (Table 1; Fig. 4a): Parral (temperature, ID: 15); Mangarral (rainfall, ID: 8); San Agustin de Puñal (rainfall, ID: 4).The considered historical period was from 1979 to 2000, whilst the Near Future (NF) period spanned from 2011 to 2040 and the Mid Century (MC) period spanned from 2040 to 2070.
The CDF curves of precipitation were modified to cope with the condition that the modelling deals with an intermittent stream, whereby in the summer season there are several months with no precipitation, so that if P o, being the probability threshold of having in the historical records a month with zero precipitation, the rainfall values of PRECIS-ECHAM having a probability threshold lower than P o were set to zero (Ahmed, Wang, Silander, Wilson, Allen et al., 2013).Doing so resulted to the fact that the CDFs of both the climatological modelled and observed data were more consistent.This approach is a direct consequence of the implicit assumptions of the BSCD method that the climatic properties remain the same in both present and future periods.
To adjust the monthly CDFs between the observed and modelled meteorological values, bias correction indexes were used.For which time series data are ranked from the lowest to the highest values.Then, for precipitation, a scale factor (PBC) was defined as the ratio between the observed value and the corresponding modelled meteorological value in the historical period.For temperature, a scale term (TBC) was defined as the difference between the observed value and the corresponding modelled meteorological value in the historical period.The historical and future modelled CDFs were adjusted by multiplying each modelled meteorological value of the monthly time series by the corresponding precipitation factors.For temperature, the adjustment operation implies adding rather than multiplying For values that lay outside the range of historical values two criteria were used.For temperature first the future trend is removed before applying BC, then, after applying BC, the future trend is added back in (Wood et al., 2003).At the other hand, the high precipitation value used the same PBC of the highest simulated historical value; likewise, the low precipitation value received the same PBC associated with the lowest simulated historical value (MetED, 2012).

Hydrological modelling under climate change conditions
With respect to the assessment of the potential implications of the climate change A1B scenario on the water cycle of the modelled catchment, it should be underlined that in addition to the causative factor of climate change, a change in climate might cause additional changes in the patterns and cycles of vegetation that directly affect the soil and vegetation parameters in the hydrologic model (e.g., K 1, FR, K c, Z1max, Df).It is likely that the intrinsic parameters linked to the deeper soil layers (e.g., K2 and Z2max) are less affected.Furthermore, it is expected that future potential climate scenarios will enhance the frequency and extent of water use conflicts.Defining corresponding changes in the value of these dynamic model parameters is rather difficult and uncertain.Therefore, all of the model parameters of interest were kept constant throughout the simulation of future conditions, assuming a steady behaviour of the physical properties of the modelled system regardless of the previously referred expected changes.As such, the study was limited to the assessment of the average hydrological effects related to the A1B Revista semestral de la DIUC 136 https://doi.org/10.18537/mskn.08.02.10 climate scenario, applying as a result of the GLUE analysis the best fitted parameterised structure of the WEAP21 model.

Magnitude and frequency of peak events
For the assessment of the frequency and magnitude of the hydrograph peak events, which are important for reservoir planning and operation, a partial duration time series (PDS) analysis (Maidment, 1993) was carried out at the outlet of the catchment where there are discharge observations, with the aim of selecting (monthly) independent peak events from the observed records as well as from the simulated time series.The analysis was based on the geometric investigation of the hydrographs to identify null derivatives and negative curvatures, as well as, on the screening of the catchment and flow properties such as the baseflow recession constant and discharge threshold value (defined as 50 m 3 s -1 ) above which a peak event is selected for further consideration.These enabled comparing the number and magnitude of simulated peak events with the ones that were observed.Foregoing might shed some insight on the implications of potential climate change on reservoir management.

RESULTS AND DISCUSSION
The GLUE analysis yielded during the calibration period 13,673 behavioural parameter sets, out of a sample of 100,000.In the validation period, upon the availability of newer evidence (i.e., observations), the updating of the parameter distribution by means of the Bayes type Eq. ( 3) reduced the number of behavioural parameter sets to 13,267.
The scatter (i.e., "dotty") plot for every studied parameter is depicted in Fig. 5.Each dot in each of the subfigures represents a behavioural parameter set in the Monte Carlo simulation analyses.The ranges in the horizontal axes coincide with the sampled parameter ranges.The scatter plots show that, in general, high likelihood values were obtained.Further, the subplots in Fig. 5 suggest that most of the parameters, except K 1 and Z1max, are insensitive to the model predictions because their scatter plots are flat-topped across the parameter subspace.The K c parameter also exhibits a marginal degree of sensitivity to model predictions suggesting that, at least in the current application, the WEAP21 model only has a marginal sensitivity to the ET p data (or to the ETo data).Thus, these results seem to confirm the premise of equifinality as it was not feasible to identify a single optimal parameter set that produced by far the best WEAP21 performance in reproducing the observed behaviour of the modelled system.
It is worth noticing that the dotty plots presented in Fig. 5 are simple projections of the complex parameter space and, as such, the assessment on the individual parameter sensitivity of model predictions requires a more suitable sensitivity test such as a sort of global sensitivity analysis (Spear & Hornberger 1980;Freer, Beven, & Ambroise, 1996;Vázquez, Beven, & Feyen, 2009).Figure 5 shows minimum values for K 1 that, according to the conceptual model of WEAP21 (Fig. 3), result in minimum rates of interflow and percolation and high rates of runoff.This hydrologic response of the basin has been observed in the wet season when soils reach saturated conditions (Uribe et al., 2004).In this context, Stewart, Majdi, Najm, Rupp, Lane et al. (2015) conducted research at plot scale and found that after a cumulative rainfall amount of 333 mm, the runoff ratio (runoff divided by rainfall) approached unity.They also observed that, under this condition, the surface crack networks, which are typical in dry shrink-swell clay soils (the predominant topsoil in the Lonquén basin), are mostly sealed, thus generating low values of hydraulic conductivity and, therefore, minimum rates of infiltration, interflow, and percolation.Figure 5 also suggests that the best model performances tend to occur for Z 1max values (maximum water holding capacity) between 250 and 450 mm, corresponding most of the time to Z 1 = 100% (Fig. 3), resulting in minimum values for K1 and high rates of runoff.In other words, according to the current hydrologic model, the average depth of the rooting zone in the upper reservoir of the study basin needs to reach a moisture value between 250 and 450 mm to produce direct runoff.This assumption needs to be confirmed by further research.Based on the above analysis, one may conclude that WEAP21 is able to represent reasonably well the soil hydro-physical properties and the hydrologic processes of the Lonquén basin in the wet season.However, the model characteristics for the wet season are not representative for the dry season when the soils shrink to the point of changing the structure of the surface layer and its hydro-physical properties.Stewart et al. (2015) found that, during dry and semi-wet conditions, the near-surface cracks and macropores strongly affect the soil hydraulic conductivity, resulting in an important rate of infiltration and interflow in the Lonquén basin soils'.WEAP21 is apparently unable to consider these processes, and overestimates the low flows in the dry and semi-wet periods, as can be seen in Fig. 6.This figure depicts the calculated 90% prediction limits using the posterior likelihood distribution and efficiency threshold or cut off criteria of EF 2sqrt ≥0.5.Under these conditions, and based on Stewart et al. (2015), WEAP21 should be capable of interpreting that the runoff process should not occur after a cumulative rainfall threshold between 120 and 170 mm and rainfall values less than 333 mm, the rates of infiltration and interflow are important.
Revista semestral de la DIUC 138 https://doi.org/10.18537/mskn.08.02.10 Figure 6.90 % streamflow prediction limits in the calibration and validation periods calculated using posterior likelihood distribution after conditioning based on observed streamflow available in the validation period.
The behaviour of the WEAP21 model in the scope of the GLUE framework responds to the condition of static soil parameters for the whole year.With other words, the values of the soil parameters are associated with the best model performances, which is during the wet season.Given the characteristics of the soils of the Lonquén catchment, a leaking bucket model where the cracks become smaller and possibly sealed with time seems to be more appropriate to represent the basin's hydrology (Stewart et al., 2015).
Despite the obvious limitations of the WEAP21 structure to simulate appropriately the lower flows in the hydrological year, the predictions are generally acceptable since the thickness of the band is relatively narrow (Fig. 6), even though in some cases the observed hydrograph crosses over the prediction limits, particularly for lower discharges.These simulation pitfalls are likely also due to the uncertainties attached to the observed discharges.
Notwithstanding the fact that in the current study forecasts of a RCM, developed specially for local conditions, were used, one should keep in mind that all of the results depicted in the following, might be affected by uncertainties that are intrinsic to future climate modelling.Unfortunately, owing to the lack of suitable information, no other CChP could be considered to assess this uncertainty.Therefore, the results of this study should be used carefully, as referential, before a decision entirely could be based on them.
Hence, the analysis of climate change revealed the presence of bias between the values of the selected gauging stations and the climate forecasts in the historical period, although the historical seasonality of the climate data was well captured.For temperature (Fig. 7a), there is a maximum difference of approximately 5 o C in terms of multiyear mean monthly values; similarly, considering mean annual values, the regional model overestimates the observations roughly with 1.6 o C (11.2%).On the contrary, for rainfall, there is not such a significant difference in terms of the multiyear mean monthly values (Fig. 7b and Fig. 7c); however, there are significant differences in terms of mean monthly values, so that the observed mean annual rainfall for the Mangarral and the San Agustin de Puñal stations are 845 and 953 mm, whilst their forecasted counterparts are respectively 876 and 927 mm.After the application of BC, the climate forecasts suggest an increment of the mean annual temperature (Fig. 8a) and a decrement of the mean annual rainfall (Fig. 8b and Fig. 8c).Herein, particularly the forecasted increment of temperature is congruent with the hypothesis of future climate change owing to gas emissions.Table 3 presents some results for the upper subcatchment.The historical values are the observed records in the case of temperature and precipitation, as well as the hydrological simulations in the case of streamflow and ET 0. It should be mentioned that in the historical period the hydrological model was run with the observations; never with the RCM forecast.In the future periods, the hydrological model was run with the RCM projection only after application of BC.Henceforth, the historical mean annual temperature of the Parral station (14.2 o C) would increase about 0.4 o C (2.8%) for the NF and +1.1 o C (7.7%) for the MC (Fig. 8a and Table 3).For rainfall, the mean annual value of (i) the San Agustin Puñal station (953 mm) would decrease with 135 mm (-14%) for the NF and 207 mm (-21.7%) for the MC (Fig. 8b and Table 3); and (ii) in the Mangarral station (845 mm) mean annual rainfall would decrease with 131 mm (-15.5%) for the NF and 169 mm (-20%) for the MC (Fig. 8c and Table 3).Considering the spatially interpolated climatic data (using base gauging stations and interpolation approaches), for the upper subcatchment, the mean annual temperature (13.5 o C) would increase about 0.4 o C (2.9%) for the NF and 1.1 o C (8.1%) for the MC.As for the mean annual rainfall (849 mm), it would decrease roughly with 128 mm (-15.1%) for the NF and 176 mm (-20.7%) for the MC (Table 3).
Further, WEAP21 estimated a mean ET 0 annual value of 1366 mm for the subcatchment LU2 which is bigger than the ET 0 calculated by CNR (1997), namely, 1290 mm.The difference of 76 mm is likely to be the direct consequence of using different methodologies for calculating ET 0 in CNR (1997) and this study.Therefore, the values obtained by WEAP21 should be used as reference, since WEAP21 uses a more modern and accepted methodology to estimate ET 0. With regard to the mean annual ET0, an increment of around 1.7% is projected for both subcatchments (LU 1 and LU2) for the NF period, whereas for the MC period the mean annual ET 0 would increase 5.3% for the LU1 subcatchment (Fig. 8d) and 5.1 % for the LU2 subcatchment (Table 3).
The flow duration curves for the three periods of interest using the GLUE based best model parameter set are depicted in Fig. 9. Figure 9a shows the observed and WEAP21 simulated time series of flow duration curves at the catchment outlet for the historical period (the validation period).The figure shows that the quality of the WEAP21 predictions are very acceptable, providing sufficient confidence to use the numerical model for predicting the hydrological conditions at the location of the projected dam (Fig. 1).In the historical period, the hydrological model was run with the observations.In the future periods the model was run with the RCM projection (after BC was applied).
Figure 9b shows the flow duration curves for the historical and future periods (NF and MC) at the location of the projected irrigation dam.According to this figure coincide the NF and MC curves but are significant different with the historical curve.The figure suggests that the current mean annual water availability (approximately 3 m 3 s -1 ) will be reduced around 40.6% (1.2 m 3 s -1 ) by the end of the NF period and 49.7% (1.5 m 3 s -1 ) by the end of the MC period (Table 3).
The PDS analysis reveals that the WEAP21 predictions demonstrate a considerable reduction in terms of the frequency and the magnitude of peak events.Whereas 15 peaks (maximum value 153.4 m 3 s -1 ) were observed in the historical period, considering the RCM forecasts as input to the hydrological model only 8 peaks (maximum value 106.6 m 3 s -1 ) were simulated for the FC period and 4 peaks (maximum value 116.0 m 3 s -1 ) for the MC period.Despite the significant uncertainty that might be associated to the only CChP used, the obtained results are likely to provide useful insight on future reservoir operation.

CONCLUSIONS AND RECOMMENDATIONS
The monthly water balance modelling for the planning of an irrigation project located in a semi-arid zone in southern-central Chile was carried out using historical data and predictions of climate change conditions.WEAP21 was used as modelling tool, and results revealed that the model faced difficulties in correctly representing low flows, despite the very good values of the EF 2sqrt model performance index, emphasising the oversensitivity of this index to the simulation of peak values.Indeed, although these results are not shown in this manuscript, the current study revealed that the EF 2sqrt index does not outperform the intrinsic limitation (i.e., being oversensitive to the simulation of peaks) of the Nash-Sutcliffe EF 2 index.The GLUE analysis proved the presence of equifinality, as there were several parameterisations of the WEAP21 code that could produce acceptable predictions, which questions the concept of attaining a single optimal set of parameters for model calibration.The present analysis showed that most of the inspected parameters are insensitive to model predictions, except the maximum water holding capacity (Z 1max) and the hydraulic conductivity of the upper (soil) reservoir (K1), whilst the Kc factor, affecting directly the value of the simulated evapotranspiration, exhibits a certain degree of sensitivity.The distribution of the parameters showed that the best model performances are associated to Z 1max values between 250 and 450 mm and near zero K1-values.These values mimic very well the hydrologic Revista semestral de la DIUC 142 https://doi.org/10.18537/mskn.08.02.10 processes of the Lonquén basin in the wet season, when the clayey soil in the surface layer swells.These parameter sets, however, failed to represent the hydrologic processes in the dry season, mainly due to shrinkage of the top soil layer impacting the physical properties of the surface horizon.These, along with the uncertainties attached to the observed discharges, resulted in a poor representation of the low flows by the fixed structure of the WEAP21 model.
In general terms provided the prediction limits a measure of the main uncertainty in the model output.Marked departures of the observed discharge from the estimated prediction band were especially observed for the lowest discharges, varying around zero, and which could not be well simulated by the WEAP21 structure.Moreover, the analysis depicted wider prediction bands during some peak events.Nevertheless, the current WEAP21 predictions may be regarded as being acceptable since, in general terms, the prediction limits are relatively narrow and the simulation of main peaks are very much acceptable, providing relevant information for the design of the irrigation infrastructure.
The BC approach was successful in bringing the original forecast to match the observations in the historical period at the base gauging stations, one for temperature and two for rainfall, which encouraged the use of the adjusted forecast for the future climate analyses.Nevertheless, in spite of the fact that the forecast is the product of a RCM developed for local conditions, one should keep in mind that all of the obtained results might be affected by uncertainties that are intrinsic to future climate modelling.Unfortunately, owing to the lack of suitable information, no other CChP could be considered to assess this uncertainty.
As expected, the mean annual temperature forecast shows an incremental evolution (positive tendency) with a maximum increment of about +1.1 o C for the Mid Century (MC) period with respect to the historical period.In this context, WEAP21 predicts for the upper subcatchment a maximum increment for the mean annual ET 0 of about 5.1% for the MC period.The RCM mean annual rainfall forecasts for the upper subcatchment show a maximum reduction of about 20.7% for the MC period.
For the analysis of the potential consequences of the A1B scenario on the hydrology of the modelled catchment, using the RCM forecast of the forcing variables (temperature and rainfall), the hydrologic predictions suggest a general reduction of the available streamflow for the projected irrigation system up to 49.7% for the MC (2041-2070) period.Furthermore, the PDS analysis showed that there might be a considerable reduction of both the frequency and the magnitude of the peak events, which might be at the end beneficial for reservoir operational rules since the reservoir input would be apparently more regular, with more attenuated peaks.However, one should keep in mind that all the results depicted herein, might be affected by uncertainties that are intrinsic to future climate modelling.Thus, these results should be used carefully, as referential, before a decision could be based entirely on them.Consequently, future work should aim at assessing the uncertainty attached to the climate change projections, and the development of a model that allows the hydrological simulation of the Lonquén catchment during the dry and wet season.

Figure 1 .
Figure 1.Location of the Lonquén river catchment and the projected irrigation district.

Figure 3 .
Figure 3. Schematic representation of the structure of the hydrological module of the Water Evaluation and Planning System (WEAP21) model for a hydrologic response unit.

Figure 4 .
Figure 4. Illustration of the methodology used to derive the spatial distribution of the meteorological data.(a) Location of the used hydro-meteorological stations; (b) example (for March) of a linear regression analysis on the elevation-temperature (multiyear monthly average) relationship; (c) elevation strips; and (d) average yearly rainfall isohyets.

Figure 5 .
Figure 5. Dotty plots of the likelihood measure (EF2sqrt) as a function of the studied WEAP21 parameters for the calibration period.

Figure 7 .
Figure 7. Assessment of the multiyear mean monthly data for the historical period 1979-2000, respectively for the observed and raw PRECIS-ECHAM forecasts: (a) temperature at the Parral gauging station; (b) rainfall at the San Agustin de Puñal gauging station; and (c) rainfall at the Mangarral gauging station.

Figure 8 .
Figure 8.Time series plots of historical data (temperature, rainfall and reference crop potential evapotranspiration-ET 0 simulated with the WEAP21 code) and future projection (PRECIS ECHAM forecast after the BC method was applied): (a) temperature at the Parral gauging station; (b) rainfall at the San Agustin de Puñal gauging station; (c) rainfall at the Mangarral gauging station; and (d) ET 0 for the lumped unit LU 1.

Figure 9 .
Figure 9. (a) Observed and simulated flow duration curves for the historical period at the outlet of the catchment; and (b) simulated flow duration curves for the three periods of interest at the projected dam location (considering the best hydrological model parameter set, characterised by EF 2sqrt = 0.97).In the historical period, the hydrological model was run with the observations.In the future periods the model was run with the RCM projection (after BC was applied).

Table 1 .
Location and main characteristics of the hydro-meteorological stations used in this study.

Table 2 .
Range of variation of the hydrological model parameters.
at the location of the projected irrigation dam (Fig1).The same value for the model parameters of both lumped units were applied, given the spatial homogeneity observed with regard to the main catchment properties, such as slope, soils, geology, etc.The variation ranges of the model parameters considered in this analysis were defined based on the irrigation project technical studies, the main catchment characteristics, and previous applications of Revista semestral de la DIUC 132 https://doi.org/10.18537/mskn.08.02.10

Table 3 .
Range of variation of the hydrological model parameters.