Assessing climate change impacts on water availability in southern Portugal using a daily sequential soil-water budget model

The Mediterranean, a typically low precipitation region with high interannual variability, has been recurrently pointed out as highly vulnerable to extreme events resulting from climate change. Hosting more than 500 million people, the region has been continuously exposed to water stress and some of those extreme events are already being experienced today.

In the southern region of Portugal, a decrease in rainfall ranging between 5 and 20 mm/year per decade is expected (Portela et al. 2020a; IPCC 2021), resulting in a likely increase in drought intensity (Espinosa & Portela 2022). The occurrence of warmer and drier summers (Giorgi & Lionello 2008; Dias et al. 2020) will result in lower soil moisture and increased evapotranspiration (Pan et al. 2015; Riedel & Weber 2020). Although the mean annual rainfall is expected to decrease in the Mediterranean region, IPCC (2021) reports a likely intensification of spatiotemporal rainfall variability with a probable increase in the intensity of heavy rainfall events.

At the soil level, the decrease in rainfall and the increase in evapotranspiration are likely to result in the decrease of available water that flows through the soil, thus generating a lower amount of diffuse groundwater recharge (Oliveira 2006). In the long term, this can lead to declining groundwater levels as the groundwater storage (reserve) may suffer from a lack of replenishment volumes.

The decrease in the natural water availability coupled with an expected increase in the water irrigation demand and, in some areas, tourism intensification, will lead to more intense water-stress situations that may ultimately drive land-use and land-cover transformations (García-Ruiz et al. 2011). This intensification in water demand may have consequences not only in surface water sources but also in groundwater bodies and groundwater-dependent ecosystems (Oliveira et al. 2007; Novo et al. 2018; Erostate et al. 2020). Also, the transition to a drier climate with more time-concentrated extreme rainfall episodes may intensify flash-flooding events (Bisselink et al. 2018).

For management planning, it is relevant, using climate-change projections, to predict the expected impacts until the end of the 21st century. Although with some uncertainty, these predictions may support the characterization of the effects of climate change on the aquifer recharge and discharge patterns (Stigter et al. 2014) and project adaptation measures. Also, understanding the trends in potential evapotranspiration (PET) changes is critical for the agricultural sector, which relies heavily on water (Rocha et al. 2020).

Hydrological modelling is an important tool to assess climate-change impacts and has been extensively used to study water availability, with the results commonly agreeing on the future trend of decreased availability (Milano et al. 2013), combined with a strong contrast between the wet and dry seasons. These models should be able to provide clear results to be translated into effective policies (Dias et al. 2020).

The uncertainty underlying future scenarios' climate models results in different predictions for the impact of climate change on groundwater recharge (Smerdon 2017; Atawneh et al. 2021). Several studies quantify the effects of climate change on hydrological processes paying particular attention to the aquifer recharge effects. For example, in the US Southern High Plains, using the SWAP model, the mean annual recharge ranged from −75% to +35% resulting from the decline in mean annual rainfall and an increase in PET (Ng et al. 2010). A systematic literature review pointed out that most studies, using different models like MODFLOW, SWAT, SWAP, or HYDRUS-1D, predict declines in groundwater recharge (Atawneh et al. 2021). In the European region, under different scenarios and climate models, groundwater recharge is likely to decrease from 30% in northern Italy to 83% in southern Spain (Treidel et al. 2011). Specifically, for mainland Spain, under a pessimistic scenario, a 12% reduction in net aquifer recharge was obtained, with aquifer recharge reduction detected in 99.8% of the territory (Pulido-Velazquez et al. 2018). For Portugal, changes in runoff will likely impact economic activities (Teotónio et al. 2017; Canuto et al. 2019; APA 2022a, 2022b).

This paper presents an estimation, at a regional level, of the effects that the changes in rainfall and PET can have on hydrological processes based on a daily sequential water budget at the soil level.

The main objective is to understand the impacts of those changes in the watersheds within the study area and characterize the extent and degree of that impact. A set of meteorological data series – rainfall and PET – were assembled as inputs. These series included (1) for the reference period, measured gap-filled time series and (2) for the future scenarios, bias-corrected EURO-CORDEX data under different representative concentration pathways (RCPs). The daily outputs of the implemented model – actual evapotranspiration (AET), direct runoff and deep infiltration, defined as the water that vertically percolates below the soil zone subject to evapotranspiration – are presented, at an annual time-scale. A comparison between the reference period and future scenarios was conducted to evaluate the degree of increase or decrease for each computed hydrological process.

Similar studies were conducted by other authors at a single-watershed scale (Oliveira et al. 2007; Mourato et al. 2014; Rocha et al. 2020) but, in this application, the daily sequential water budget model was applied sequentially covering all the southern Portuguese region watersheds making use not only of ground-based observational data but also corrected future series under different RCP scenarios.

The results may be valuable for defining adequate climate-change adaptation strategies. This is the case in the implementation of methodologies such as managed aquifer recharge (MAR), which, coupled with other water -resource management measures, may prove useful in supporting scarcity adaptation. One outcome of this research – the geographical distribution of the expected change in aquifer recharge – may serve as a water availability criterion in an analysis for selecting suitable areas for MAR implementation.

Economically, the Alentejo region underwent a significant transformation resulting from the construction and progressive expansion of the Alqueva dam multipurpose system, which increased water availability for agriculture and urban supply (Canuto et al. 2019; APA 2022a; EDIA 2022). In the Algarve, irrigation and urban supply rely heavily on groundwater (APA 2022b). Furthermore, in 2020, the Algarve showed the highest volume of water distributed per inhabitant in Portugal (119.6 m³/inhabitant), followed by the Lisbon metropolitan area (77.0 m³/inhabitant) (INE 2020).

According to CORINE Land Cover 2018 (Copernicus 2022a), 62% of the study area is agricultural areas followed by 36% as forests and semi-natural areas. The expansion of the irrigated areas in both Algarve and Alentejo with the goal of optimizing productivity is increasing the pressure on both surface and groundwater resources.

The simulations were conducted with the daily sequential water budget model at soil level BALSEQ (Lobo Ferreira 1981) adapted to account for the antecedent moisture conditions (AMCs) (Oliveira 2006). It estimates direct runoff, deep infiltration, and AET using daily rainfall (P) and PET, the US Soil Conservation Service Curve Number (CN), and the maximum amount of water available in the soil for evapotranspiration (AGUT). The last two parameters are originally derived from soil and land use maps based on the methods defined by Vermeulen et al. (1993). The assumptions underlying the conceptual model for the soil control volume of the BALSEQ model are the following: (1) the only water input is infiltration, given by the difference between rainfall and direct runoff; (2) soil flow is purely vertical downwards.

The model requires a small amount of input data, which is advantageous, namely when aiming at modelling large areas (Oliveira 2006; Martins et al. 2021b). Daily data is required, in opposition to decadal or monthly data, as the water budget at soil level is sensitive to the time step.

BALSEQ has been extensively used in several case studies for determining the potential groundwater recharge (e.g., Lobo Ferreira 1981; Paralta & Oliveira 2005; Oliveira et al. 2007, 2022). In comparison with other models, the BALSEQ structure does not have a user interface. However, it is easily scripted for sequential simulations in multiple watersheds. Python3 was used in this research.

Alternative models such as SWAT or USGS PRMS-IV require more complex parametrization with a large number of inputs and the process of conducting numerous simulations for different watersheds is not easily automated. Another technique that may allow groundwater recharge only to be computed is by conducting surface-flow hydrograph separation (Oliveira 2005, 2006). Although surface-flow information is available for the region, there were three main factors that resulted in the exclusion of this last method: (1) the poor and irregular regional coverage of the gauging stations within the study area; (2) the scarcity of surface-flow historical data (short series with many gaps) that would not allow the generation of adjusted surface-flow series for future scenarios; (3) the difficulty – if not almost impossibility – to conceive an accurate-enough model for hydrograph separation able to be generalized to numerous watersheds in a large region without accounting for their specific constraints.

Martins et al. (2021b) compared BALSEQ outputs with observed surface-flow data in different watersheds in Alentejo and Algarve with acceptable results, thus supporting the use of this model for the evaluation of water availability.

The study area was subdivided into the 21,335 second-order-stream-segment small watersheds, as defined by the Portuguese Environmental Agency (APA 2015). Each watershed, or hydrologic unit, was assigned a set of four variables (P daily series, PET daily series, CN, and AGUT). The outputs were aggregated for the seven main river basins in the study area.

A set of 153 meteorological stations with at least five years of daily rainfall records were selected from the National Water Resources Information System database (SNIRH – https://snirh.apambiente.pt/) covering the study area. The meteorological stations' information and locations are presented in Supplementary Figure S1 and Table S1.

From the original set of meteorological stations with daily rainfall data existing in this database, stations were selected that have at least five years of data. Although other daily rainfall databases are available, the ground-based SNIRH daily rainfall database was found to be more adequate. The use of ground-based databases continues to be fundamental in hydrology studies, as verified by several authors (e.g., Sun et al. 2018; Mudelsee 2019; Dong et al. 2020; Djaman et al. 2020). In fact, the use of other meteorological databases, such as satellite or reanalysis databases, would still require the use of the ground-based type of data for their validation, since they are often associated with meshes with a spatial discretization that is very coarse, so they do not constitute a more precise alternative to those databases.

An alternative ground-based database for the region would be the Portuguese Meteorological Institute (IPMA) network. However, this network contains only 33 stations within the study area, the time-series length and frequency of missing data are unknown, and it is not an open database.

If, for a specific day, there is no rainfall value at station ⁠, a new station should be selected by following steps (2)–(6) again. An original dataset missingness map was produced and is presented in Supplementary Figure S2. The map shows that meteorological stations with more than 50% of records are evenly distributed in the study area, thereby enhancing the likelihood of using local information to fill the gaps; 80% of the gaps were filled using records from a station distanced up to 17.8 km, and 95% up to 29.4 km (see Supplementary Figure S3).

Additionally, to have a global overview of the appropriateness of the gap-filling process, a comparison was made based on dimensionless empirical non-exceedance probability curves of rainfall. To obtain this curve for each dataset (original non-filled and gap-filled) and time scale (year, month, and day), all the corresponding rainfall values were previously normalized (by subtracting the average and dividing the result by the standard deviation). The results, presented in Supplementary Figure S4, confirm a similar global pattern between the two datasets, supporting the suitability of the gap-filled dataset to perform the analysis.

Spatial weighted rainfall series were computed for each hydrological unit. These were produced by discretizing each hydrological unit in a 0.5 km cell-sized grid and computing the precipitation in each grid cell by applying the inverse distance squared weighting method. A subset of at least two rain-gauge stations selected on a predefined 5 km search distance from the centroid of the hydrological unit was used. When the minimum number of stations was not met, the search distance was expanded with 5 km increments until at least two rain-gauge stations were found. The precipitation series of each hydrological unit was computed by averaging the precipitation values from each grid cell.

Reference evapotranspiration data (mm) was collected from 33 monitoring stations (Supplementary Figure S1 and Table S1) from different sources and periods, namely: (1) EDIA meteorological database with monthly data from 2010 to 2019 (https://regante.edia.pt/suporteaatividade/meteorologia/SitePages/Home.aspx), (2) DRAP Algarve database with daily data from 2006 to 2019 (https://www.drapalgarve.gov.pt/pt/servicos-e-produtos/servicos/fitossanidade/avisos-agricolas), and (3) INAG (1998) with monthly data for the period of 1941 to 1991. The previous data was then used to assemble monthly data series for the period of 1941 to 2019. PET was assumed equal to reference evapotranspiration. Monthly data was converted to daily data by dividing it by the number of days of the month, under the assumption that PET daily variation is relatively small, and that the inaccuracy produced by this process has reduced impact on the water budget. Each missing value was assigned a mean PET, computed by using the other existing years' same-day records. A spatial weighted PET series was generated for each hydrological unit using the Thiessen polygons method.

For future scenario evaluation, meteorological data was extracted from the EURO-CORDEX database (Jacob et al. 2014; Benestadt et al. 2021) which contains European regional models for (1) historical observations, (2) RCP 4.5 scenario (characterized by ambitious emissions reductions, likely to cause a 1.1–2.6 °C global mean temperature (GMT) rise in the 2081–2100 period in relation to the average of the 1986–2005 period), and (3) RCP 8.5 scenario (business as usual, 2.6–4.8 °C GMT rise) (IPCC 2021). The grids of 0.11° of resolution represent downscaled regional climate models (RCMs) from general circulation models (GCMs – 0.44° resolution) that were developed by different research groups and spanned different time-periods. No reference evapotranspiration projections were available; hence, the Penman–Monteith method (Allen et al. 1998) was applied based on monthly projections of maximum and minimum temperatures, sunshine duration, relative humidity, surface wind, and pressure at sea level. The models and variables considered are represented in Table 1.

An ensemble mean of the three regional models was composed for rainfall and also for PET and a bias correction of the systematic error in climate projections was determined by adjusting the modelled dataset with observational data from a defined historical reference period (Benestadt et al. 2021). A complete series was extracted from the GCM/RCM grids for each rainfall and PET monitoring station, using non-gap-filled data for bias correction. Correction adjustment of bias was conducted based on the quantile–quantile calibration method to correct mean, variability, and shape errors of the RCM cumulative distribution functions (Boé et al. 2007; Amengual et al. 2012; Gudmundsson et al. 2012). The analysis procedure is based on the methodology presented in Duijnisveld (2018), making use of the qmap R package (Gudmundsson et al. 2012).

The CN and AGUT parameters were mapped using the method presented by Martins et al. (2021a), derived from soil (DGADR 2020) and CORINE Land Cover 2018 (CLC) (COPERNICUS 2020a) maps.

The water budget effects of karstic and urban areas were considered in the calculations. Areas of occurrence of soils associated with karsts were identified. In those areas, it is assumed that, after briefly flowing over the surface, all direct runoff will intersect a karstic feature and infiltrate through the soil directly below the depth subject to evapotranspiration. In these circumstances, all direct runoff is considered null and added to deep infiltration ⁠. These karstic areas represent 6.52% of the total study area. For the impact of urban areas (namely artificial surfaces), which represent 1.37% of the study area, the permeability ratio map (COPERNICUS 2020b) was used to compute the mean impervious area ratio for each anthropogenic area identified in CLC maps. These parameters are presented in Supplementary Table S2.

Given a mean impervious area ratio ⁠, direct runoff is assumed to be given by and deep infiltration by ⁠.

The BALSEQ model was run sequentially for each grouped set of variables and parameters (P daily series, PET daily series, CN, AGUT, ⁠, and existence of karstic features) in each small-watershed for the control period of 1941–2019 based on gap-filled observational data. Future-scenario simulations were conducted for the time-windows of 2020–2060 and 2061–2100 both for RCP 4.5 and RCP 8.5 conditions. No changes were considered in land use/land cover between 1941 and 2019 and future scenarios.

Rainfall is expected to decrease towards the end of the century both for RCP 4.5 and 8.5, with a more notable decrease in the latter. RCP 8.5 also shows a considerable decrease in the [2020–2060] period, which is larger than the RCP 4.5 decrease in the [2061–2100] period.

The RCP 8.5 scenario is particularly critical, and a significant rainfall reduction is observed in comparison with the reference time-window. In both future scenarios – RCP 4.5 and 8.5 – the general seasonal distribution of rainfall follows the historical behaviour with dry summers with almost no rainfall and wet autumn and winter months. This represents the typical rainfall annual distribution for the Mediterranean region.

This rise in PET is primarily due to the expected increase in air temperatures (Riedel & Weber 2020; IPCC 2021), which are the main variables that control the changes of this parameter. Similarly to the rainfall process, the variations in PET under climate-change scenarios are not expected to result in a major variation in this process seasonal pattern – PET is higher in warmer summer months and lower in winter months. As expected, the more perceptible PET increases occur in the summer months.

For the mean annual Ed in the control period (Figure 6(a)), 89% of the study area shows Ed below 50 mm/year (19,134 of the 21,335 small watersheds), 2% above 100 mm/year (454 small watersheds) and 0.1% above 200 mm/year (20 small watersheds). The computed higher Ed values are observed where the higher rainfall values were obtained, in the regions between Mira, Barlavento and eastern Arade, the southern part of Guadiana and the northern part of Sotavento and northern Guadiana (Figure 2).

Concerning mean annual Ip (Figure 7(a)), 70% of the region shows Ip below 50 mm/year (14,652 small watersheds), 6% between 100 and 200 mm/year (1,403 small watersheds) and 0.4% above 200 mm/year (88 small watersheds). The geographical location of higher values of Ip is observed to be in the Sado, Ribeiras do Alentejo, and a great part of Sotavento. These areas are associated with the location of the major sedimentary and karstic regions (Figure 1), whose geological characteristics are prone to higher infiltration.

Table 2 compares, in percentage of the area, the differences between the mean annual Ed and Ip, for each RCP scenario time-window, and the reference period. More than 93% of the total study area in all scenarios presents a decrease in Ed between 0 and 50 mm/year, and less than 0.1% of the total area shows an increase. The major decreases in mean annual Ed are observed in the areas where the mean annual rainfall is above 500 mm/year and is expected to show a higher decline in future.

Concerning Ip, Table 2 shows that most of the area presents differences between −50 and 0 mm/year, pointing to a general decrease in the region. However, there are areas where Ip is expected to increase: for the case of RCP 4.5, in 31% (period [2020–2060]) to 30% (period [2061–2100]) of the area and, for the case of RCP 8.5, in 18% (period [2020–2060]) and 5% (period [2061–2100]) of the area.

Water availability is regarded as the water that will flow naturally on the surface or in the underground saturated zone and may be exploited for human use. It may be regarded as the sum of direct runoff and deep infiltration if one can assume that all water that vertically flows down below the soil zone subject to evapotranspiration will result in groundwater recharge. In this case, is equal to ⁠.

Regarding the ratio between Etotal and rainfall (Figure 8(c)), Etotal is less than 20% of the rainfall in 93% of the total area (20,723 small watersheds). Etotal represents more than 50% of the rainfall in less than 1% of the area. When comparing the mean annual values obtained for all future scenarios with the values of the control period, more than 92% of the total area shows a decrease in Etotal (Table 2). These areas are predominantly located in the lower Guadiana, Sotavento, southern Mira and northern Barlavento and Arade (Figure 8(b1), 8(b2), 8(d1) and 8(d2)).

The decrease in rainfall, particularly pronounced in RCP 8.5 [2061–2100], results in a general decrease in almost all river basins of Ed and Ip. As mentioned further on, exceptions are observed for the case of Ip in Sado and Arade. Considering specifically the Ip in comparison with the reference period, for RCP 4.5 and respective time-windows of [2020–2060] and [2061–2100], the observed decreases were: Sotavento – 42% and 39%; Mira – 30% and 33%; Barlavento – 15% and 16%; Ribeiras do Alentejo – 15% and 17%; and Guadiana – 13% and 12%. For Sado and Arade, in the [2020–2060] time-window, an increase of 4% and 2%, respectively, was observed followed by a decrease in the [2060–2100] of 3% and 4%, respectively. For RCP 8.5, the observed decreases for the equivalent time-windows were: Sotavento – 51% and 60%; Mira – 50% and 58%; Barlavento – 32% and 49%; Ribeiras do Alentejo – 28% and 41%; Guadiana – 24% and 40%; Arade – 15% and 41%; and Sado – 12% and 29%.

Concerning Ed, it always decreases in RCP 4.5. Comparatively with its reference period and respectively for the time-windows of [2020–2060] and [2061–2100], the decreases were: Guadiana – 82% and 80%; Mira – 80% and 74%; Ribeiras do Alentejo – 70% and 64%; Sado – 80% and 76%; Sotavento – 71% and 69%; Barlavento – 70% and 64%; and Arade – 68% and 65%. The equivalent results for RCP 8.5 are: Guadiana – 87% and 83%; Mira – 86% and 83%; Ribeiras do Alentejo – 72% and 70%; Sado – 83% and 80%; Sotavento – 77% and 79%; Arade – 77% and 78%; and Barlavento – 76% and 78%.

It is possible to observe in the later future period that Guadiana, Mira, Ribeiras do Alentejo, and Sado show an increase in Ed in comparison with the first future period, for both scenarios. In RCP 8.5, Arade, Barlavento, and Sotavento show a decrease in Ed in the later future period when compared with the first half.

The evolution of the two hydrological processes, Ed and Ip, between the first and second future time-windows can be explained as follows: (1) for the Ip process, the decrease in rainfall will result in a decrease in the amount of water stored in the soil, which, combined with increased water used in evapotranspiration, will reduce the amount of water that contributes to deep infiltration; (2) the increase of mean Ed between the first time-window and the second one, even with the expected decrease in rainfall, can be a consequence of the change in rainfall patterns. Similar results were presented by Chang et al. (2017) and Gusev et al. (2019).

In fact, more frequent intense rainfall events are expected in the future, even under climate-change rainfall-decrease scenarios (Soares et al. 2017; Yin et al. 2018; Santos et al. 2019; Espinosa et al. 2024). Intense rainfall occurring in shorter periods and exceeding the infiltration capacity will result in higher Ed. Cumulatively to what was said about deep infiltration, they will also contribute to further diminish Ip, since they reduce the opportunity for infiltration.

The increase in PET leading to an increase in AET in certain river basins under the business-as-usual scenario RCP 4.5 highlights the sensitivity of these regions to rising air-temperatures. The results for RCP 8.5 simulations present an even more severe scenario, as expected. This increase in PET may result in a higher demand for water resources, which could impact various sectors, including agriculture, industry, and urban water supply.

The changes in the AET, Ed, and Ip processes within the main river basins of the study area may have important implications in water resources management, namely regarding groundwater depletion. These likely alterations are indicative of the potential challenges related to water scarcity that may arise in the region due to climate change.

This study allowed the identification of the areas where groundwater recharge is likely to be reduced as a result of PET increase and rainfall (P) decrease. Areas where deep infiltration (Ip) and direct runoff (Ed) are expected to decrease for all RCP scenarios coincide with those where rainfall is expected to suffer higher decreases. These decreases mean lower water availability. Mira, Barlavento, southern Guadiana, and Sotavento are the river basins where more differences are found between the control period and the RCP scenarios. Comparatively with the control period, AET is more likely to increase in almost all river basins for RCP 4.5 mainly due to the temperature rise. Although not expected, the results showed two basins (Arade and Sado) where Ip is likely to increase in the first time-window of RCP 4.5 as a consequence of the decrease in Ed.

Future methodology developments may include parameter calibration based on field data or benchmarking the adapted BALSEQ results with those from other models. Such developments would also allow a better understanding of the amount of Ip that returns as baseflow to the stream flow. Also, this study only considered Ip as a hydrological process, comparable to a potential diffuse recharge, but did not consider other types of recharge that occur naturally, e.g., fluxes between surface water and groundwater bodies, or inter-watershed flows that are not easily quantified.

Furthermore, it is important to emphasize that no changes were considered in the soil use/land cover throughout time. It is expected that, under future increased water stress, some adaptation measures will be required, particularly in agriculture, the most water-dependent sector, eventually resulting in land use changes. On the other hand, higher PET may also ultimately result in higher water demands, putting additional pressure on surface and groundwater resources, already overstressed under the expected effects of more intense and recurrent droughts. The higher water stress may force changes in agricultural practices, for example, by reducing the more intensive water-use crops, thus decreasing the water yields. Such changes in crop patterns will impact the general water budget and should be considered in future studies.

Although not addressed in the study, it is relevant to stress that there is uncertainty associated with the input meteorological and spatial data, the parameters' values and the downscaled climate projections for the adopted scenarios. This unquantified uncertainty may impair the use of the results presented herein, particularly from the perspective of the main stakeholders' (farmers, water utilities, etc.) decision-making processes. However, it may be useful as part of other applications aiming to support policy-making or management strategies, such as the identification of climate-change-impacted areas suitable for deploying adaptation measures, like MAR.

From the water-management point of view, the results of this study, particularly the identified areas where climate-change impacts are expected to be higher, will provide input to the decision-making process for the effective implementation of such adaptation measures. The results may also contribute towards the definition of adequate changes in land use, the outlining and fine-tuning of protection policies for maximum infiltration zones, and the implementation of alternative and integrated water-resources management solutions, all together, ensuring better use of water surpluses and improving groundwater storage efficiency.

T.N.M. thanks the Fundação para a Ciência e Tecnologia (FCT), Portugal, for the PhD Grant PD/BD/135590/2018 within the H2Doc program. This paper was developed under the framework of LNEC's 2013–2020 Research and Innovation Plan, Risk Management and Safety in Hydraulics and Environment (Process Nr. 0605/112/20383).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.