Impact of Hydroclimatic Changes on Water Security in the Cantareira Water Production System, Brazil

1. Introduction

The availability of freshwater is essential for the sustainability and well-being of societies, and its efficient management becomes even more critical in the face of the challenges posed by climate change. In this context, the Cantareira Water Producer System (CWPS), which supplies the Metropolitan Region of São Paulo (MRSP), stands out as an important subject of study to understand the impacts of these changes on water security.

The CWPS is responsible for supplying water to approximately 9 million people in the MRSP. Its strategic relevance is undeniable, considering the fundamental role that water plays in socioeconomic development and quality of life. However, in recent years, the CWPS has faced significant challenges due to climate change, jeopardizing the water security of the region.

Climate change has affected the hydrological cycle and precipitation patterns worldwide, bringing direct consequences to water availability and quality. In the case of CWPS, an increase in climatic variability is observed, with extreme droughts and heavy rains becoming more frequent. This climatic variability, combined with increasing water consumption, puts CWPS in a vulnerable position, susceptible to periods of scarcity and instability in supply [1].

The operating rules of CWPS established by the current water resource management system do not consider the impacts of climate change. The hypothesis of this study was based on the idea that even with control over demands and land use in the watersheds, the current infrastructure and operating rules will not be able to ensure water security in the future, as these challenges require more complex management measures.

Computational simulations are hybrid scientific tools that allow both the reproduction of a phenomenon and intervention in its behavior. Evaluating the impact of climate change on the hydrological regime has been presented by researchers as useful techniques for addressing the problem [2,3]. Reliable information on potential changes to future hydrological conditions require knowledge of the influence of heterogeneity in soil-vegetation-atmospheric relationships [4], and basic data related to climate, topography, land cover and soil characteristics are necessary for the hydrological simulation of basins [5].

The main objective of this paper is to analyze the impact of climate change on the water security of CWPS. To achieve this goal, multidisciplinary approaches will be used, including analysis of historical climate data, hydrological modeling, and evaluation of existing water resource management systems.

2. Study Area

The study area is in Brazil and comprises a portion of the Paraíba do Sul River Basin, which flows into the Atlantic Ocean, and another portion of the Tietê River Basin, which is a tributary of the Paraná River, ultimately forming the River Plate Basin in the context of South America (Figure 1).

The study area shown in Figure 2 is a delineation of four hydrographic basins representing the CWPS: the Jaguari/Jacareí River, with its outlet at the Buenópolis fluviometric station, and the Atibaia River, with its outlet at the Valinhos Hydrometric Station, the both rivers are source of the Piracicaba River, the Juqueri River, the both a tributary of the Tietê River, with its outlet at the Paiva Castro Dam; and the Jaguari River, a tributary of the Paraíba do Sul River, with its outlet at the Jaguari Hydroelectric Plant.

The study area was further analyzed using the Digital Elevation Model (DEM) developed by [6], employing machine learning techniques to remove buildings and forests from the Copernicus DEM to produce a global elevation map with a grid spacing of 1 arc second (∼30 m). For the extraction and analysis of hydrographic information from the DEM, the TauDEM software was utilized (available at https://hydrology.usu.edu/taudem). This tool, employing flow direction methods, facilitated the delineation of contribution areas and determination of stream networks [7,8].

As described in Figure 3, the Jaguari/Jacareí basin, with its outlet being the reservoir of the same name (JAG), discharges into the downstream Buenópolis basin with its outlet at the control point of the same name (BUE). The Cachoeira and Atibainha basins, with outlets at their respective reservoirs (CAC and ATA, respectively), flow into the Atibaia basin with a control point of the same name (ATI). Subsequently, the waters proceed to the Valinhos control point (VAL), incorporating the contribution area of the same name.

As shown in Figure 2 and Figure 3, the reservoirs are interconnected by tunnels and conduits, with the main demand points being the control points, the Santa Inês Pumping Station (SIPS), and the Hydroelectric Power Plant (HPP). Table 1 presents the drainage area of the incremental basins within the study area.

3. Materials and Methods

The study’s method began with information on climate change scenarios, simulated the behavior of the CWPS for each scenario, and generated results in terms of water security indicators. Figure 4 shows the flowchart of the method.

Since the goal is to intervene in the behavior of the CWPS to achieve acceptable water security indicators, the methodology checks the sensitivity of operational variables. The chapter is divided into the presentation of the CWPS study area, description of the models and databases used for calibration and simulations.

2.1. Data

The observed meteorological data were extracted from [5], who generated a dataset of daily gridded meteorological data with a spatial resolution of 0.1° for the period from 1961 to 2020 for the entire Brazil. The information includes minimum and maximum temperatures, precipitation, solar radiation, wind speed, relative humidity, and potential evapotranspiration. Observed data from 11,473 rain gauges and 1,252 meteorological stations were used.

Simulated meteorological data were extracted from [9], who provided a dataset based on a set of 19 CMIP6 climate models with bias correction for projections over Brazilian territory, based on the SSP2-4.5 and SSP5-8.5 scenarios. The Quantile Delta Mapping approach was used to correct biases in daily time series of precipitation, maximum and minimum temperature, solar radiation, wind speed, and relative humidity. The bias-corrected dataset is available for both historical (1980-2013) and future (2015-2100) simulations at a spatial resolution of 0.25°.

The calculation of potential evapotranspiration from the climate models was processed using the PyEt library [10] employing the Hargreaves methodology [11] presented in equation 1.

$$\mathrm{E}\mathrm{P}=\frac{k}{\lambda }{R}_a\left(\frac{T_{max}-T_{min}}{2}-\mathrm{17,8}\right)\sqrt{T_{max}-T_{min}}$$

(1)

where k is an empirical constant, with a adopted value of 0.0135, Ra is the extraterrestrial solar radiation in MJ.m⁻²d⁻¹, λ is the latent heat of evaporation in MJ.kg⁻¹, Tmax is the maximum air temperature in °C, Tmin is the minimum air temperature in °C. The values of Ra and λ are calculated using the latitude grid values.

Discharge data from the stream gauges were obtained from the hydrological database of São Paulo State, and information about the inflow discharge of the reservoirs was obtained from the reservoir monitoring system of Brazilian National Water Agency.

The data were collected, extracted for the area of interest, and stored in a dimensional database [12]. This facilitated data aggregation for generating information and statistics. The database was built on the PostgreSQL platform, organizing information by basin, daily time basis, parameter, model, and scenario.

2.3. Evaluate climate projection models

Multicriteria decision analysis was used to evaluate climate projection models with bias correction based on the proposal by [13]. The analysis is based on comparing the precipitation from the models with observed data in the study area. The applied criteria aim to assess the overall ability of models, both with and without bias correction, to replicate the main statistics of observed data relevant to hydrological studies. The indicators include the time series (a) daily, (b) monthly, (c) annual, and (d) hydrological year (from October to September). Include average seasonality (e) day of the year, and (f) month of the year. Wet extreme about (g) day of the year, and (h) month of the year. Dry extreme of (i) day of the year, and (j) month of the year.

The performance of the ten indicators was evaluated using the modified Kling-Gupta Efficiency KGE’ coefficient initially proposed by [14] as an improvement over the Nash-Sutcliffe coefficient, and later modified by [15]. The coefficient was applied to evaluate projected precipitation in Brazil, showing satisfactory results [16]. Equation 2 shows KGE’.

$$\mathrm{K}\mathrm{G}{\mathrm{E}}^{\prime }=1-\sqrt{{\left(r-1\right)}^2+{\left(\frac{{\mu }_{sim}}{{\mu }_{obs}}-1\right)}^2+{\left(\frac{{\sigma }_{sim}.{\mu }_{obs}}{{\sigma }_{obs}.{\mu }_{sim}}-1\right)}^2}$$

(2)

Where: r (dimensionless) is the correlation coefficient, μ (mm) is the mean precipitation, σ (mm) is the standard deviation of the precipitation series, and the indices obs and sim indicate the observed and simulated series, respectively. The KGE’ varies from -∞ to 1, with 1 indicating the best fit. The KGE’ values corresponding to the selected criteria were scored according to the categories presented in Table 1 [17].

The sum of the scores for each indicator was used for ranking climate models, from highest to lowest. The maximum value is 20, where all 10 indicators are classified as high, and the minimum value is zero, where all indicators are low.

2.4. Hydrological model

The hydrological model used in this paper is based [18], a conceptual models based on some physical reasoning that represent catchment processes by several interconnected buckets, which mimic water storage and transfer within a homogeneous area. This model has been successfully applied in Brazil [19]. This deterministic model operates continuously over time and considers a daily time step. It employs a lumped approach at the watershed scale and adopts a conceptual framework incorporates linear transfer functions, considering the direct runoff, soil moisture depletion, unsaturated zone dynamics, groundwater recharge, and three linear reservoirs representing surface, soil, and aquifer retention characteristics (Figure 5).

As shown in Figure 5, a fraction of the precipitation (P in mm) is conveyed as direct runoff (DR in mm) using equation 3. This DR is routed through a linear reservoir reaching the basin outlet as routed direct runoff (RDR in mm) using equation 4.

$$\mathrm{D}\mathrm{R}=\frac{{(\mathrm{P}-\mathrm{I}\mathrm{A})}^2}{\mathrm{P}-\mathrm{I}\mathrm{A}+\mathrm{S}\mathrm{A}\mathrm{T}-\mathrm{R}\mathrm{S}\mathrm{O}\mathrm{I}\mathrm{L}}$$

(3)

$$\mathrm{R}\mathrm{D}\mathrm{R}=(1-\mathrm{K}\mathrm{D}\mathrm{R})\times \mathrm{R}\mathrm{S}\mathrm{U}\mathrm{R}\mathrm{F}$$

(4)

Where IA is initial abstraction (mm), SAT is the soil saturation capacity (mm), RSOIL is the level of the soil reservoir (mm), RSURF is the level of the surface reservoir (mm), and KDR is the direct runoff recession constant (d⁻¹).

The remaining water depth (P − DR) is subject to depletion at the potential evapotranspiration rate (PE in mm). The excess water (P – DR − PE) infiltrates into a linear reservoir that represents the upper soil horizon (unsaturated zone). Moisture is lost from this zone at an evapotranspiration rate (ER in mm) proportional to the moisture content (MC in %, shown in equation 5) and PE. The output from the unsaturated zone reservoir corresponds to the recharge (REC in mm) of the groundwater reservoir. If RSOIL exceeds the field capacity (FC in mm), REC occurs following the equation 6.

$$\mathrm{M}\mathrm{C}=\frac{\mathrm{R}\mathrm{S}\mathrm{O}\mathrm{I}\mathrm{L}}{\mathrm{S}\mathrm{A}\mathrm{T}}$$

(5)

$$\mathrm{R}\mathrm{E}\mathrm{C}=(\mathrm{R}\mathrm{S}\mathrm{O}\mathrm{I}\mathrm{L}-\mathrm{F}\mathrm{C})\times \mathrm{M}\mathrm{C}\times \mathrm{C}\mathrm{R}\mathrm{E}\mathrm{C}$$

(6)

Where CREC (%) is recharge coefficient, a parameter related to the movement of water in the unsaturated soil zone, the water transferred to the groundwater reservoir regulated by this parameter. The output of the groundwater reservoir is the base flow component (BF in mm) shown in equation 7.

$$\mathrm{B}\mathrm{F}=(1-\mathrm{K}\mathrm{B}\mathrm{F})\times \mathrm{R}\mathrm{G}\mathrm{R}$$

(7)

Where RGR is the level of the groundwater reservoir (mm), and KBF is the base flow recession constant (d⁻¹). Runoff (RO, in m³s⁻¹) is obtained by equation 8.

$$\mathrm{R}\mathrm{O}=(\mathrm{R}\mathrm{D}\mathrm{R}+\mathrm{B}\mathrm{F})\times \mathrm{A}\times 86.4$$

(8)

Where A is the drainage area in km². Following the above procedure, the levels in each reservoir are continuously updated for each simulation day. This calculation applies to the headwater basins. To calculate the discharge for the downstream basins, an additional method is employed to routed upstream discharge (RUD in m³s⁻¹), calculated by the routing model (McCarthy, 1938) describes in the equation 9.

$${\mathrm{R}\mathrm{U}\mathrm{D}}_{\mathrm{i}}=\frac{\left(1-2KX\right){\mathrm{U}\mathrm{D}}_{\mathrm{i}}+\left(1+2KX\right){\mathrm{U}\mathrm{D}}_{\mathrm{i}-1}+\left(2K\left(1-X\right)-1\right){\mathrm{R}\mathrm{U}\mathrm{D}}_{\mathrm{i}-1}}{2K\left(1-X\right)+1}$$

(9)

Where the RUD is routed upstream discharge (m³s⁻¹), K is the parameter to represent the time release is delayed (d) and X is the parameter referring to the reduction of the release peak in the stream. Index i is about the value on a time interval and i-1 is a previous interval value. The downstream discharge (UD in m³s⁻¹) is the sum of RO and RUD. To ensure calculation stability the channel is divided into sections to K fits within the limits shown in equation 10.

$$\frac{1}{2(1-X)}<K<\frac{1}{2X}$$

(10)

The hydrological model calibration was performed using data from October 2011 to September 2019 (8 hydrological years) for the sub-basins. The choice of this period aims to cross-reference the available observed daily flow data with rainfall and evapotranspiration data in the watersheds. Unlike the traditional approach, no specific period was selected for model validation, as suggested by [20], who evaluated the performance of two conceptual hydrological models in 463 watersheds using 50 different data splitting schemes. They demonstrated that it is more robust to use the complete available dataset for calibration and skip model validation.

The parameters were optimized using the minimize function from the scipy.optimize library [21], which provides a unified interface for finding local minima of non-linear optimization problems. Within a trust region, a local model of the objective function is built based on first and second derivative information. It approximates the best point until reaching a local minimum of the original objective function. The parameters varied SAT [100, 2000], FC [30, 50], IA [2.5, 5], CREC [0, 20], KDR [0.2, 5], KBF [30, 180], K [1, 10] and X [0.25, 0.35].

The objective function was based on KGE on basins headwater basins, for better representing the volume flowing into reservoirs and downstream basins are performed Log-NSE, for better representing the limit flows at control points. The indicators NSE, R2 and PBIAS were also analyzed.

4. Discussion

Daily flow is greatly influenced by variations in soil moisture, infiltration and precipitation. The physical water characteristics of the soil influence infiltration, water storage, and aquifer recharge. These variables were adjusted in the calibration of the hydrological model by parameters such as the direct flow recession constant (KDR), the base flow recession constant (KBF), and the recharge coefficient (CREC). The intra-annual distribution of precipitation is an important condition for the performance of the hydrological model.

We selected ten indicators for a representative climate model, four for the accumulated rainfall of the time series (day, month, year, hydrological year), two for the daily and monthly seasonality, and four for the daily and monthly extremes in wet and dry periods. Of the 19 models surveyed, the GFDL-CM4 presented the best performance for the study area in terms of the modified Kling-Gupta efficiency coefficient (KGE’). According to [23], the model appears to produce a reasonable multidecadal modulation of the El Niño-Southern Oscillation due to a slight weakening of the annual cycle, relatively small biases in seasonal spatial patterns of top-of-atmosphere fluxes, surface temperature, and precipitation.

GFDL-CM4 without bias correction represents better the time series of annual and hydrological years as well as the maximum monthly precipitation than the other models. In the CWPS basins, the model represents precipitation well during the wet period (October to March) and tends to underestimate it by 23.3% during the dry period (April to September). SPI and DM indicate a higher projected drought frequency, especially for the scenario SSP5-8.5. These results are consistent with the study by [25], which identified the GFDL model for the Piracicaba River basin using the tool developed by [26].

Performance indicators indicated a satisfactory adjustment of the hydrological model calibrated to the time series observed with the objective functions KGE and log-NSE. The projected flows in the basins without a reservoir were estimated using parameters calibrated with an objective function Log-NSE, which best adjusted the minimum flows. The basins with reservoirs were simulated with parameters calibrated with KGE.

The simulated flows were strongly correlated with forecast precipitation for both climate scenarios. The SSP5-8.5 scenario indicates the lowest water availability in all basins, and VAL is the one that shows the greatest reduction in minimum flows. Considering the frequency of the average of seven consecutive days of the P/Q ratio, JPS, CAC, PAI and ATA basins show faster reactions in flows to the P deficit in both climate scenarios. For a more pronounced deficit (P/Q = 0.2), the VAL basin showed greater sensitivity in the SSP2-4.5 scenario. According to [27], basins with low P/Q values are subject to longer streamflow droughts, possibly due to lower subsurface storage.

In systems with high water demand such as the CWPS, a deficit of seven consecutive days is enough to cause operational disruptions and ultimately supply conflicts. The expected moderate (P/Q = 0.5) and severe (P/Q = 0.2) deficits are high.

5. Conclusions

It is expected that the results of this study will provide a comprehensive and well-founded understanding of the challenges faced by CWPS in the context of climate change, as well as suggest strategies and adaptive measures to ensure water security for MRSP. Understanding these impacts and identifying sustainable solutions are of utmost importance to guide decision-making and formulate public policies aimed at efficient water resource management and adaptation to climate change.