Evaluating the BFAST method to detect and characterise changing trends in water time series: A case study on the impact of droughts on the Mediterranean climate

Evaluating the BFAST method to detect and characterise changing trends in water time series: A case study on the impact of droughts on the Mediterranean climate Maria Paula Mendes a , Victor Rodriguez-Galiano b, ⁎,David Aragones b,c a CERIS, Civil Engineering Research and Innovation for Sustainability, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisbon, Portugal b Departamento de Geografía Física y Análisis Geográfico Regional, Universidad de Sevilla, 41004 Seville, Spain c Remote Sensing and Geographic Information Systems Lab (LAST-EBD), Estación Biológica de Doñana, C.S.I.C., 41092 Seville, Spain HIGHLIGHTS GRAPHICAL ABSTRACT •BFAST allowed detection of abrupt and gradual changes in the freshwater resources. •BFAST enabled the detection and quantification of recovery and drawdown periods. •BFAST identified trends in groundwater and surface water due to weather fluctuations. •BFAST assessed impact of water management on both surface and groundwater resources. ABSTRACTARTICLE INFO Editor: Martin Drews Keywords: Water resources LOESS, BFAST Time series analysis Climate variability Climate resilience Mediterranean climate regions are facing increased aridity conditions and water scarcity, thus needing integrated management of water resources. Detecting and characterising changes in water resources over time is the natural first step towards identifying the drivers of these changes and understanding the mechanism of change. The aim of this study is to evaluate the potential of Breaks For Additive Seasonal and Trend (BFAST) method to identify gradual (trend) and abrupt (stepchange) changes in the freshwater resources time series over a long-term period. This research shows an alternative to the Pettitt's test, LOESS (locally estimated scatterplot smoothing) filter, Mann– Kendall trend test among other common methods for change detection in hydrological data, and paves the way for further scientific investigation related to climate variability and its influence on water resources. We used the monthly accumulated stored water in three reservoirs, the monthly groundwater levels of three hydrological settings and a standardized precipitation index to show BFAST performance. BFAST was successfully applied, enabling: (1) assessment of the suitability of past management decisions when tackling drought events; (2) detection of recovery and drawdown periods (duration and magnitude values) of accumulated stored water in reservoirs and groundwater Science of the Total Environment 846 (2022) 157428 Abbreviations: BFAST, Breaks for Additive Season and Trend; LOESS, locally estimated scatterplot smoothing; (CF), Campina de Faro; (MS), Meridional Rim of Sotavento; (MOSUM), Moving sum; (OLS), Ordinary least squares; (QS), Querença-Silves; (SPI), Standartized Precipitation Index; (STL), Seasonal Trend decomposition using Loess. ⁎Corresponding author. E-mail address: [email protected] (V. Rodriguez-Galiano). http://dx.doi.org/10.1016/j.scitotenv.2022.157428 Received 24 February 2022; Received in revised form 11 July 2022; Accepted 12 July 2022 Available online 19 July 2022 0048-9697/© 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv bodies after wet and dry periods; 3) measurement of resilienceto drought conditions;(4) establishment of similarities/ differences in trends between different reservoirs and groundwater bodies with regard to drought events. 1. Introduction The Mediterranean region and Portugal mainland are being faced with increasing aridity conditions (Sousa et al., 2011;Vicente-Serrano et al., 2014) that fit well into predicted future climate trends, with higher frequency and intensity of droughts and floods in regions already prone to these extreme weather events (IPCC, 2018, 2019). A resulting decrease in water resources is foreseen that will seriously affect the availability of freshwater resources for drinking and irrigation purposes (Hiscock et al., 2008). Climate variability, in general, and precipitation variability, in particular, can have substantial effects on the rate of aquifer recharge (Green et al., 2011) and the availability of surface water resources stored in dammed reservoirs (Deng et al., 2020). Sustainable groundwater management requires the evaluation of groundwater resources to ensure a balance between abstraction and groundwater recharge. Water resource availability is dependent on interannual short-term climate fluctuations (e.g. precipitation and evapotranspiration), but also on long-term changes such as climate change (Pinto et al., 2021). Furthermore, abrupt changes, such as those resulting from floods or overexploitation, or gradual changes resulting from modifications in land cover, management, as well as socio-economic drivers (Antunes Aghsaei et al., 2020;Lian et al., 2020;Mendes and Ribeiro, 2014)caninfluence fluctuations in groundwater levels and surface water reserves. Traditionally, different problems related to hydrology have been addressed applying machine learning for predictive modelling (Adnan et al., 2021; Cardenas-Martinez et al., 2021;Mendes et al., 2021;Pandey et al., 2020), numerical models (Casillas-Trasvina et al., 2019;dos Santos et al., 2019; Sondermann and de Oliveira, 2021), and statistical methods (Młyński et al., 2021;Peña-Angulo et al., 2020;Vicente-Serrano et al., 2019; Zeleňáková et al., 2017) among others. The establishment of gradual or abrupt changes of temporal hydrological dynamics can help to establish and identify any probable causes, e.g. to distinguish natural climate variability from anthropogenic changes, especially if supported by descriptive characteristics, such as drought indices or climate teleconnections (Anyah et al., 2018;Secci et al., 2021;Wang et al., 2021;Xu et al., 2021). Parametric approaches require features to follow specific statistical distributions, such as a normal distribution, and they are very sensitive to outliers. Non-parametric approaches do not require any prior assumptionabout the distributional characteristics of the samples (Teegavarapu, 2019). The non-parametric Mann–Kendall trend test (Gilbert, 1987)iswidely used in hydrology, since it is simple, outlier-resistant and can handle missing values or values that are below detection limits. However, this method canonly beused in the case of monotonic trends and it assumes that sample data are serially independent. Block bootstrapping (Carlstein, 1986a; Kunsch, 1989a) will be used to overcome the serial dependence enabling the checking the robustness of the trend assessment. On the other hand, LOESS (locally estimated scatterplot smoothing) is less outlier-resistant, but allows examination of both monotonic and nonmonotonic trends (Cleveland and Devlin, 1988). LOESS is recommended for assessing statistically significant trends at the groundwater bodies level, within the scope of the European Union River Basin Management Plans (Grath et al., 2001). Although, abrupt and gradual changes can coexist (Sharma et al., 2016), these two abovementioned tests cannot detect abruptchanges, since it assumes that thetrend component variessmoothly. The Pettitt's Test (Pettitt, 1979) is widely used in hydroclimatological studies for addressing abrupt change or a step change within the time series, allowing the detection of a single shift at an unknown time t(Mallakpour and Villarini, 2016). The BFAST (Breaks for Additive Season and Trend) method was developed to identify gradual and abrupt changes, allowing the detection of multiple breakpoints, while explicitly considering seasonal variations (Verbesselt et al., 2010a). BFAST permits characterisation of the time, magnitude and direction of change, flagging disturbance events in groundwater and surface water resources. The B-FAST method can identify gradual and multiple abrupt changes in a time series, but it has rarely been used in hydrological applications as few exceptions can be found in the literature. For instance, Liu et al. (2021) used BFAST to evaluate the trends of runoff change in the upper and middle reaches of the Yellow River, Coladello et al. (2020) to assess the over-abundance of aquatic macrophytes in tropical reservoirs and, Deng et al. (2020) to identify and supplement newly impounded reservoirs during 2000–2018 in the Yellow River. The aim of this work is to evaluate the potential of BFAST to quantify the impacts of disturbance events, such as droughts periods, on groundwater and surface water resources. The specific objectives are to: (1) Detect abrupt changes such as hydrological droughts in the seasonal and trend components of piezometric levels and accumulated stored water in dammed reservoirs in a Mediterranean climate region; (2) Identify trends in groundwater and surface water due to weather fluctuations, discerning these from effects induced by human activities such as overexploitation of groundwater resources or the filling of a reservoir; and (3) Compare the performance of the BFAST method with that of the Pettitt test and Mann- Kendall. 2. Methods and data 2.1. Water resources The southernmost region of mainland Portugal, named Algarve (with 16 municipalities and around a resident population of 467,475; CENSUS, 2021)(Fig. 1), has high urban pressures along the coastline, which began with the tourism boom in the 1960s. Along with tourism, the expansion of irrigated agriculture also started during the 60's, leading to an exponential rise in water demand. Since 1990, Algarve has expanded its provision and range of golf courses (nowadays around 39), and facilities extensively (Barros et al., 2010). Public water supply, with an volume of 69.5 hm 3 (in 2018) comes from surface-water sources (76 %) and groundwater (24 %). Public water supply is characterized by a strong seasonality, with water consumption in the summer months being more than one and a half times the average annual consumption. Groundwater was the dominant source for public water supply, irrigated farming and golf courses until the end of the 20th century (Stigter et al., 2009), which led to a severe decrease in the piezometric levels in the majority of the aquifers located in the coastal zones (Maia and Silva, 2009;Mendes and Ribeiro, 2014). Nowadays, groundwater remains the primary water source for irrigated farming (≈68 %) and golf courses (≈65 %). Surface water represents around 40 %of the total water consumption in this region (APA, 2016). The Bravura and Beliche dams were built in the late 1950s and in the 1980s (Fig. 1), respectively, to cope with water shortages in months with the greatest tourist influx (June to September) (Antunes do Carmo, 2019; Carriço et al., 2014). In 2002, >80 % of the total public supply was sustained by surface water after the construction of the Odeleite dam in late 1990's and the rehabilitation of the Beliche and Bravura dams (Table 1). The Odeleite and Beliche reservoirs have a storage capacity of 117 hm 3 and 47.6 hm 3 , respectively. These two dams were connected through a tunnel, constituting a unique water supply system. The Odeleite-Beliche system is located in the Guadiana river basin, and serves the irrigation and water consumption needs of the eastern part of the study region. Odeleite and Beliche have, respectively, an average annual precipitation of 722 mm and 644 mm. M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 2 The Bravura reservoir (storage capacity of 32.33 hm 3 ) and two other reservoirs (Arade and Funcho) supplied the western part alone until construction of the Odelouca dam, with the capacity of these three storage reservoirs being lower than the Odeleite-Beliche system (Maia and Silva, 2009). The Bravura reservoir is located in Odeáxere river basin and has an average annual precipitation of 821 mm. The severe drought that occurred in 2004–2005 caused a major water shortage in these storage reservoirs, with water use restriction measures being imposed for the agriculture sector. This issue was addressed by reactivating some municipal boreholes that had previously been abandoned. Additionally, new boreholes located in the western part of the Querença-Silves (QS) aquifer system (Fig. 1) were used (Maia and Silva, 2009). The QS aquifer system is the largest (324 km 2 ) and the most important groundwater resource (mean yield of 11.1 l/ s) in the region (Stigter et al., 2009). This aquifer system is comprised of karstified carbonate rock and has a significant recharge (long-term mean annual recharge of approximately 111 hm 3 /year) (ARH Algarve, 2012a;Hugman et al., 2013; Nicolau, 2002). Its western sector has a small (approximately 1.5 %) hydraulic gradient, a predominant E–Wflow direction and discharge occurs at Estombar springs located near the mouth of the Arade river (Almeida et al., 2000;Neves et al., 2016;Stigter et al., 2009). The Campina de Faro (CF) multi-aquifer system (Fig. 1) is composed of three aquifer systems: the deeper Cretaceous limestone/marls, the overlying Miocene sandy limestone aquifer system, and the sand and gravel Plio-Quaternary upper aquifer. This latter aquifer system is located in the centre and south, with an average thickness of 50 m, and can be considered phreatic (Stigter et al., 2006). Recharge occurs directly via precipitation and is approximately 9.95 hm 3 /year, which is the same as the estimated available resources (ARH Algarve, 2012). The Meridional Rim of Sotavento (MS) has a low productivity and is composed of biocalcarenites and sands. The mean annual recharge is 11.08 hm 3 /year with mean available resources of 9.97 hm 3 /year. These three groundwater bodies (where two of the three are aquifer systems) mirror distinct hydrogeological settings with quick response to precipitation events. Three piezometers located at the Querença-Silves aquifer system (QS), Campina de Faro aquifer system (CF) and Meridional Rim of the Sotavento groundwater body (MS) were selected (Table 1), each with good data quality (i.e. good length and integrity of water level records), as well as three dammed reservoirs (Bravura, Beliche and Odeleite; Table 1), in order to illustrate the advantages of using BFAST for detecting both gradual and abrupt changes in the trend component. 2.2. Climate The Algarve region has a warm-summer Mediterranean climate (Csa) under the Köppen-Geiger classification, with a well marked dry summer with nearly five-month-long. 64 % of the precipitation falls between November and February, with a low mean annual precipitation of around 500 mm/year (1981–2010 climate normal). The highest monthly accumulated precipitation occurrs inDecember (115.6 mm) followed by November (83.5 mm). The rivers show typically seasonality, with high winter discharges and low summer discharges, generally resulting of the seasonality of the rainfall exacerbated by the high tempertures of the summer months (Guerreiro et al., 2017). The daily average temperature ranges between 12 °C (January) and 24 °C (July and August). The average number of tropical nights (i.e., days with low temperatures over 20 °C) for the summer months are 18 days (June), 8.6 days (July), 8.7 days (August) and 4.8 days (September) (IPMA, 2021). Due to the orography, coastal areas receive on average <70 % of precipitation than northern areas, making the south coast one of the aridest regions on the Portuguese mainland (Neves et al., 2019). Only the Querença-Silves aquifer is located at the foothills of the Serra do Caldeirão, one of the rainiest regions in the Algarve (mean annual precipitation around 1000 mm), where episodes of heavy rainfall are most frequent and torrential nature (Fragoso and Gomes, 2008). 2.3. Standardized Precipitation Index (SPI) The standartized Precipitation Index (SPI; Mckee et al., 1993)isoneof the drought indices that is extensively used worldwide since is easy to Fig. 1. a) Location of Portugal in Europe; b) the four major hydrogeological units and; location of the six water wells (the three reservoirs, Bravura, Beliche and Sotavento) and the three wells (Sotavento, Silves e Faro) and the two aquifer systems. Table 1 Location of water points. ID Database Name Coordinates WGS 84 Altitude (m) Surface (S)/Groundwater (G) body Data Source NW 600/134 Sotavento 37.17 −7.56 33.51 G- Meridional Rim of Sotavento (MS) https://snirh.apambiente.pt 610/167 Faro 37.06 −8.00 4.13 G- Campina de Faro (CF) https://snirh.apambiente.pt 595/215 Silves 37.18 −8.44 63.76 G- Querença-Silves (QS) https://snirh.apambiente.pt Odeleite Odeleite 37.33 −7.52 47.00 S- Guadiana river basin https://snirh.apambiente.pt Beliche Beliche 37.28 −7.51 39.00 S- Guadiana river basin https://snirh.apambiente.pt Bravura Bravura 37.20 −8.70 96.00 S- Ribeiras do Algarve https://snirh.apambiente.pt M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 3 calculate, is effective in analysing wet-dry periods/cycles and, only requires precipitation as an input parameter (WMO, 2012). Calculation of the SPI is based on the long-term precipitation data for a certain period. SPI values are units of standard deviation from the long-term mean, allowing the SPI to be used to compare any precipitation anomalies for any location and for any number of time scales. When calculated for longer periods, the SPI is more related to hydrological droughts with low frequency. The 14-month SPI was used because SPIs of these longuer timescales are usually tied to streamflows, reservoir levels, and groundwater levels (WMO, 2012). In this study, it serves as an indicator for reduced reservoirs and groundwater recharge. The 14-month SPI values were classified according to seven categories: extremely dry (Min ≤SPI ≤-2.0), very dry (−2.0 < SPI ≤-1.5), moderately dry (−1.5 <SPI ≤-1.0), normal precipitation (−1.0 <SPI ≤1.0), moderately wet (1.0 <SPI ≤1.5), very wet (1.5 < SPI ≤2) and extremely wet (2.0 <SPI ≤Max) (EDO, 2020). 2.4. Methods to evaluate the temporal behaviour of time series The Mann-Kendall test (Kendall, 1955;Mann, 1945) is a nonparametric trend test, widely used in hydrology. The Mann-Kendall test is based on the calculation of Kendall's tau measure of association between two samples, which is itself based on the ranks with the samples. This test requires the datato be independent,meaning that a positive/negative serial correlation might affect the power of the test, rejecting or accepting the null hypothesis (H 0 ) of monotonic trend absence even if true (Fatichi et al., 2009). For instance, a positive serial correlation among the observations can be perceived by the test as the data having a trend, even in the absence of a trend (error of type I; Cox and Stuart, 1955). Non-parametric bootstrap strategies such as block bootstrapping (Carlstein, 1986b;Kunsch, 1989b) can be used to cope with the serial dependence, breaking the data set into smaller chunks for sampling purposes that might wreck any correlations that exist in the larger data set (Arteche, 2021;WCAP, 2004). By separating the observations in blocks far enough, in time, they are nearly uncorrelated and can be treated as exchangeable, if the time series has length nand can be factored as n=b∗l,whereband lare integers (Chernick and LaBudde, 2014). Using block bootstrap (Canty and Ripley, 2021), time-series data will be first deseasonalized, and then block bootstrap will run with a fixed block length (l) equal to 12 and a (b) number of blocks for the assessment of a Mean Mann Kendall Tau, i.e., trend component. Seasonal Trend decomposition using Loess (STL) was proposed by Cleveland et al. (1990), and it is a filtering procedure for decomposing a seasonal time series at a period vinto three components: trend (T v ), seasonal (S v ) and remainder (R v ) (or irregular or error), for v=1to N: Yv¼TvþSvþRvð1Þ STL consists of an inner loop nested inside an outer loop (in the case of STL robustness estimation). If Y v indicates no non-Gaussian behaviour (i.e. aberrant behaviour on S v and T v ), then only the inner loop is used. The STL algorithm (Verbesselt et al., 2015) is tuned using six parameters: n (p) − number of observations in each cycle of the seasonal component related to periodicity; n (i) −number of passes through the inner loop; n (0) : number of robust interactions of the outer loop, e.g. if robustness is not neededthen: n (0) =0;n (l) −smoothing parameter for the low-pass filter, e.g. if robustness is not needed then: n (i) =2;n (t) −smoothing parameter for the trend component; and n (s) −smoothing parameter for the seasonal component. Selecting the first four parameters is straightforward since they are based on empirical rules/default values and can be selected in an automated way. The n (s) is always odd and its choice can be ambiguous as it depends on the expert knowledge of the characteristics of the data being analysed. Nevertheless, STL considers subseries of values at each position of the seasonal cycle (called cycle-subseries), revealing changes in timing, amplitude and variance that occur in the seasonal cycle (Sellinger et al., 2008). However, STL assumes that the T v and S v components change slowly and smoothly, and does not permit the detection of abrupt changes within the time series (Verbesselt et al., 2010a). Pettitt's test (Pettitt, 1979) is a rank based nonparametric statistical test method for determining a shift (i.e., change point) in the mean of a time series. Pettitt's test can be used to detect the statistical significance and position of an abrupt change or a step change in the time series. Traditionally, this test has been assumed to be less sensitive to outliers and skewed distribution. The null hypothesis (H 0 ) of the test is that, when arbitrarily splitting the sample into two segments, there is no change in the mean value of each segment (Xie et al., 2014). Pettitt's test normally only allows detection of a unique changing point in a time series. The Pettitt test uses the Mann-Whitney Utest to find a single change point in the median of a series and returns a p-value for the change point (Pettitt, 1979). All the aforementioned statistical methods are most frequently used in hydrological studies and can either assess the trend or an abrupt change. This work presents the BFAST method, which is commonly used in remote sensing studies (Verbesselt et al., 2010a), but not in hydrological field, and has some advantages, as it allows the smooth and abrupt changes detection of time series. BFAST is an additive decomposition model like STL (Eq. (1)) that iteratively fits a seasonal model (S v ) and a piecewise linear trend (T v ). T v is determined with mbreakpoints υ 1 ∗ ,…,υ m ∗ , with defined υ 0 ∗ =0and υ m+1 ∗ =N, and corresponding segment-specific intercept α j and slope β j , such as: Tv¼αjþβjv υ j−1<v≤υ j;where j ¼1;…;m:ð2Þ This way, the trend component has every time series segment fixed between breakpoints, but can differ across breakpoints. The magnitude of an abrupt change at a breakpoint is derived by the diffrence between T v at υ j−1 ∗ and υ j ∗ ,suchas(Verbesselt et al., 2010a): Magnitude ¼αj−1−αj  þβj−1−βj  tð3Þ A harmonic seasonal pattern (S v ) is calculated with pbreakpoints υ 1 # ,…, υ p # with defined υ 0 # =0and υ p+1 # =N, and corresponding segment-specific harmonic a i,k and phase δ i,k ,forKharmonic terms (k=1,2,3)and frequency (f;e.g.f=12 annual observations for monthly time series) (Verbesselt et al., 2010b): Sv¼X K k¼1 αi;ksin 2πkv fþδi;k  υ# i−1<v≤υ# i;where i ¼1;…;p ð4Þ Hence, the seasonal component can be different for every time series segment comprised between two different breakpoints. The segment decomposition model is not straightforward, as the number of required segments in the trend (m+ 1) and seasonal (p+ 1) components have to be determined. First, the ordinary least squares (OLS) residuals, based moving sum (MOSUM) test (Chu et al., 1995), is used to test whether one or more breakpoints occur. If the test is significant (p< 0.05), the optimal position of breakpoints and the optimal number of breaks can be estimated by minimising the residual sum of squares, and minimising the Bayesian Information Criteria (BIC), respectively. The method starts the iteration decomposition by estimating a seasonal component b Svfrom a standard season-trend decomposition, followed by the iterative estimation of the parameters until the number and position of the breakpoints are stable: Step 1- the OLS-MOSUMis used to test the breakpoints ofthe trend component. If they occur, the number and position of these breakpoints (υ 1 ∗ ,…,υ m ∗ ) are estimated using least squares from the seasonally adjusted data; Step 2- a robust regression based on M-estimation is performed to estimate the α j and β j parameters, i.e. trend changes are estimated based on Eq. (2); M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 4 Step 3- the OLS-MOSUM is used to test the breakpoints of the season component. If they occur, the number and position of these breakpoints (υ 1 # ,…,υ p # ) are estimated using least square from the seasonally adjusted data; Step 4- a robust regression based on M-estimation is performed for the estimation of the a i,k and δ i,k parameters, i.e. the seasonal changes based on Eq. (4). Some limitations of this method have been reported in previous studies (Masiliūnas et al., 2021;Saxena et al., 2018): i)the time seriesmust be completed without gaps; ii) the user has to choose a priori the number of breakpoints based on the duration and pattern of the time series; only the strongest breaks will be identified if the number of breaks in the time series is greater than the number of breaks specified by the user and; iii) the processing time can be high in the case of big data, as it needs to converge in both the seasonal and trend components. All statistical analyses were performed in the R statistical environment (R Core Team, 2018). 3. Results and discussion The widelly used Mann-kendall test was performed to test monotonic trend. The monthly accumulated stored water in the three reservoirs showed nonstatistical evidence of monotonic trend before and after the use of block bootstrap (Table 2). If the Mann-Kendall test was performed without considering the serial correlation of the time series, a false significant monotonic trend (α= 5 %) would be detected in all three piezometers. Block bootstrap Mann-Kendall test was performed detecting nonmonotonic behaviour of the trend for all-time series (Mean Mann Kendall Tau ≅0), illustrating the need to examine within the hydrological time series both monotonic and non-monotonic trends. The trend (i.e., smoth change) decomposed by BFAST and LOESS showed similar behaviour in each freshwater source(Fig. 2). The difference between these two methods lies in the fact that BFAST enables identifying several distinct time periods, quantifying the slope and its significance in each segment. Moreover 14 months SPI data were closely related to the behaviour of the trend of both freshwater sources. Piezometers Sotavento (MS) and Faro (CF) presented a similar trend behaviour (LOESS, Fig. 1), following the dry/wet conditions (14-month SPI; Fig. 2), with valeys during the drought years and peaks during the wet years. Both piezometers (Sotavento and Faro) had two similar significant breakdates (p<0.05) for groundwater levels in October/November 1983 and December 1995. The years of 1983 and 1995 were classified as drought years, ending these periods in October 1983 and December 1995, respectively. BFAST defined the biggest change of the trend component (recovery) from November to December 1995, where magnitudes of 4.71 m and 1.67 m, were respectively obtained for the piezometers Sotavento and Faro. From July 1995 to June 1996, piezometers Sotavento and Faro were highly correlated (R 2 ≅1), presenting an upward trend (slope = 32 cm and slope = 10 cm, respectively). With regard to the piezometer Sotavento, a significant change point (i.e., at 95 % confidence level) in the mean value of piezometric series was detected by the Pettitt's test in December 1995 (Table 3), where the mean value increased from 25.8 m to 26.6 m. This change point also corresponds to the highest magnitude breakpoint detected by BFAST. Concerning the piezometer Faro, the Pettitt's test detected the change point (i.e., a level shift) in May 1999 (Table 3), five months later than the breakpoint detected by BFAST in December 1998. According to Pettitt's test results, the mean went down from 1.32 m to 0.96 m. Silves groundwater levels also closely reflected the SPI 14 months. During the period from January 1990 to December 1992, there was a drawdown of 1.69 m (1992 was a drought year with six months classified as severe drought, four months as moderate drought, and one month as extreme drought), followed by a negative slope of −1.21 m for the period of July 2003 to June 2006 (2005 was a drought year, with three months classified as a moderate drought, eight months as severe and one month as extreme). After the drought event of 2004/05, no recovery periods occurred in this piezometer (Fig. 2). The significant breakdates (p<0.05) of November 1989, December 1992 and 1995, and January 2010 detected by BFAST in piezometer Sotavento were similar to those observed in piezometer Silves, although with a delay in 1989 (1 month), 1995 (2 months) and 2010 (2 months). Moreover, the piezometer Silves had also similar datebreaks of February 1996 and 1999, June 2006 and March 2010 to the piezometer Faro, albeit with a delay in 1995 and 1998 (2 months) and 2010 (3 months), and advancement in 2006 (4 months). Mallakpour and Villarini (2016) showed that the Pettitt test detects more easily a break if the variability of the data decreases. Observing Table 3, Sotavento was the piezometer with the lowest coefficient of variation (4.64 %). For this piezometer, the datebreak determined by the BFAST coincided with the breakpoint of the most significant magnitude detected by the BFAST.The majority of these breakpoints showed a connection between drought episodes (14 months SPI) and shifts in piezometric levels. As a result, the main drivers for the trend behaviour during these common periods could be related to climate variability. The most negative slopes of the three piezometers (Fig. 2) occurred during the period when groundwater assumed a relevant role in the Algarve region (i.e., until 1998; ARH Algarve, 2012), showing that the uptake of water intensified during the periods of drought. Some significant breakpoints were detected that could not be associated with climate variability, such as: December 2000 in piezometer Sotavento (slope of −0.69 m), December 1998 and 2001 (respectively, slopes of −0.38 m and 0.12 m) in piezometer Faro, and February 1999 and June 2006 (respectively, slopes of 0.19 m and 0.22 m) in piezometer Silves (Fig. 2). These breakpoints can be convincingly related to human activities. For instance, during the period of March 1996 and June 2003, the piezometer Silves showed two positive slopes (0.19 m and 0.22 m, respectively) related to a recovery period since, after the 2000s, surface water resources began to serve as the primary freshwater source (Fig. 2). However, due to the drought from 2003 to 2005 (Fig. 3), exploration of the Querença- Silves aquifer system resumed in October 2004 with an extraction rate of 500 l/s, which was halved in July 2005 due to a risk of saltwater intrusion (Do Ó and Monteiro, 2006). Concerning the water accumulated at the end of the month, Beliche and Bravura reservoirs showed significant breakpoints in the trend in December of 1995, after the drought of 1995 (Fig. 3), with the highest turn of the trend component (recovery) from November to December Table 2 Mann-kendall results of trend evaluation. Tau Kendall Two-sided p-value Block Bootstraping Result H 0 (α=5%) Mann Kendall Tau estimator standard error Mean Mann Kendall Tau Sotavento 0.125 6.425E-05 0.0960 −0.0042 H 0 = not reject: Nonstatistical evidence of monotonic trend Faro −0.132 4.07e-05 0.0863 −0.0026 H 0 = not reject: Nonstatistical evidence of monotonic trend Silves 0.130 0.00011 0.1096 0.0016 H 0 = not reject: Nonstatistical evidence of monotonic trend Odeleite −0.065 0.10874 0.1092 −0.0084 H 0 = not reject: Nonstatistical evidence of monotonic trend Beliche −0.077 0.03014 0.0983 −0.0063 H 0 = not reject: Nonstatistical evidence of monotonic trend Bravura −0.023 0.52052 0.1059 −0.0052 H 0 = not reject: Nonstatistical evidence of monotonic trend M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 5 1995 (i.e., magnitudes of 35.3 hm 3 and 24.7 hm 3 , respectively). Regarding the water stored in the Beliche reservoir, Pettitt's testdetected a lowering of the mean of the trend in March 1993, being detected a breakpoint by the BFAST three months earlier, i.e., December 1993 (Fig. 2). The Bravura reservoir exhibited negative slopes during a moderate drought between 1990 and 1995 (−4.21 hm 3 ) and a severe-to-extreme drought between 2003 and 2005 (−6.74 hm 3 ). All positive slopes can be related to normal precipitation periods, i.e. January 1999–December 2002 (slope = 2.98 hm 3 ) where 77 % of the months were normal and the early months were moderate drought (16.7 %), January 2006–January 2012 (slope = 1.54 hm 3 ) where approximately 85 % are normal years and, February 2012–April 2015 (slope = 2.43 hm 3 ), all classified as normal years (Fig. 2 and Fig. 3). Odeleite dam was built after the two other dams in 1997, and BFAST signalized four breakdates (November 1999 and 2009, December 2005 and 2016, February 2013), that were, all related to drought episodes. During the study period, a significant positive slope (26 hm 3 ) was only found between December 1996 and November 1999 (during moderate drought; Fig. 3). This positive trend corresponded to the filling of the reservoir. The highest magnitude (in this case a lowering of storage levels) was in December of 2005 at the end of a severe drought period (≅18 hm 3 ) (Fig. 3). The Pettitt's test detected a significant point of change in the mean in August 2014 (Table 3), that could not be related to the breakpoints signalized by BFAST, and to drought periods Beliche showed a similar pattern of storage water levels, with breakdates for October 2005, February 2010 and 2013, and November 2016. These similarities are related to the fact that these two reservoirs constitute a unique system. After the 2003–2006 drought, the first recovery period of the Bravura reservoir lasted approximately six years at 0.62 hm 3 /month (January 2006 to January 2012), and Odeleite reservoir took four years and four Fig. 2. Results of the BFAST application to hydrological series. Rectangle on the top:Trend resulting from BFAST (black line) and LOESS (blue line) for the piezometric levels (Silves, Sotavento and, Faro) and monthly accumulated stored water (Bravura, Beliche, and Odeleite); Rectangle at the middle (black line): Seasonal component resulting from BFAST; Rectangle at the middle (Red line): Remainder; Rectangle at the bottom: Trend resulting from BFAST with respective segments slopes (β)significant (p<0.05). Table 3 Results of the Pettitt's test. Name Two-sided p-value Probable change point at time Coefficient of variation (%) Pettitt's test Sotavento 4.628E-11 December 1995 4.64 Faro 2.20E-16 May 1999 36.00 Silves 1.71E-11 August 1992 24.46 Odeleite 3.01E-07 August 2014 19.81 Bravura 2.95E-07 March 1996 22.80 Beliche 3.29E-05 March 1993 26.51 M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 6 months at 2.22 hm 3 /month (November 2005 to February 2010). The BFAST method enabled the detection and quantification of recovery and drawdown periods, something that is not possible using common methods (LOESS, Pettitt's Test and Mann-KendallTest).A correspondence was found between the mean change point detected by Pettitt's test and breakpoints detected by BFAST in all data, except for Beliche reservoir. This might be a consequence of data variability (Table 3) and the extension of the record which is the shortest of the three reservoirs (Fig. 2)(Mallakpour and Villarini, 2016). Storage capacity, independent of the medium, had an influence on the resilience of water resources to drought conditions. For instance, both phreatic aquifers reacted quickly to weather conditions (14-month SPI), as the recovery periods matched to the end of drought periods and the upwards were related to the starting of the normal/wet period of rain. The detection of abrupt changes in the trend components of accumulated stored water in dammed reservoirs allowed us to establish relationships between similar trends. As expected, the Bravura reservoir followed the standard 14 month SPI more closely than the Odeleite-Beliche system as the latter is a unique system, located in two different locations. Moreover, Beliche reservoir, with its lower storage capacity, sometimes showed a delay or advance in the breakpoints of the trend when compared to the Odeleite reservoir. Nevertheless, some differences found between the trends of both reservoirs only occurred within one month. None of the data (groundwater and surface water) showed seasonal breakpoints. 4. Conclusions The performance of BFAST in detecting changing trends was evaluated in water resources over long-term period, (between 26 and 37 years). BFAST enabled robust detection of abrupt and gradual changes. The BFAST method outperformed the widely used Mann-Kendall and Pettitt tests in hydroclimatological studies, showing that monotonic trends and a single abrupt change are not adjustable to this type of data. The 14- Fig. 3. Z-score values of the trends resulting from BFAST (dotted line) and LOESS (solid line), and 14 months SPI (Standardized Precipitation Index) (yellow area) for the zscore values of piezometric levels (Faro, Sotavento, and Silves) and z-score values of monthly accumulated stored water (Odeleite, Beliche, and Bravura). M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 7 month SPI could be tied to reservoir storage levels and groundwater levels, showing BFAST some common datebreaks of surface water and groundwater resources (e.g., a common decline of the levels of the piezometers Silves and Sotavento and Beliche reservoir levels during the drought episode of 1992). BFAST enabled assessment of the suitability of past management decisions when tackling drought events, showing, for instance, that during the period of 2003–2005, the three reservoirs demonstrated negative slopes (Beliche = −6.83 hm 3 ;Odeleite=−18.15 hm 3 ,and Bravura = −6.74 hm 3 ), where Odeleite (storage capacity of 117 hm 3 ) started showing this negative trend earlier in January 2003, followed by Bravura (storage capacity of 32.33 hm 3 ) in November 2003. This negative trend observed in the Odeleite storage volumes can be attributed to the fact that the severity of the hydrological drought was greater in that basin (i.e. the 14-month SPI classified six months as extreme drought) than that observed in Bravura (i.e. only one month was classified as extreme drought). Furthermore, Odeleite is located in a basin with average annual precipitation (722 mm/year) lower than observed in Bravura (821 mm/year). Piezometer Silves only showed a negative trend in July 2003, with this downward period ending in June 2006. This gap shows that the conjunctive use of groundwater and surface water resources was a good decision to reduce the effect of this drought, as there was a shortage of surface water resources. In the Mediterranean region, drought events were found to increase substantially in summer. This fact can have a greater effect on society than the long-term rise in the average temperature, since adaption to weather extremes has to be faster. Quantification of the recovery and decline periods of water resources and their significance can help to establish mitigation measures for drought events. In the European Union, although it is not mandatory, the Drought Management Plan (EC, 2007) is a planning document that would benefit from the use of BFAST, since it can contribute to the establishment of a reliable early warning system based on hydrological indicators including stored surface reservoir volumes and aquifer water levels. CRediT authorship contribution statement Maria Paula Mendes: Conceptualization, Data curation, Research, Writing-original draft, Software Victor Rodriguez-Galiano.: Conceptualization, Methodology, Writing, Visualization, Reviewing and editing. David Aragones-Borrego: Software. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Victor Rodriguez-Galiano reports financial support was provided by Ministerio de Ciencia e Innovación. Victor Rodriguez-Galiano reports article publishing charges was provided by University of Seville. Acknowledgements Maria Paula Mendes would like to thank Fundación Carolina for her two-month mobility grant for researchers at the University of Seville, Spain. The authors are grateful for the financial support provided by the projects RTI2018-096561-A-I00 and US-1262552, funded by the “Agencia Estatal de Investigación”of the “Ministerio de Ciencia e Innovación”and “Consejería de Economía y Conocimiento”of the “Junta de Andalucía”, respectively. References Adnan, R.M., Mostafa, R.R., Islam, A.R.M.T., Gorgij, A.D., Kuriqi, A., Kisi, O., 2021. Improving drought modeling using hybrid random vector functional link methods. Water 13 (23), 3379. https://doi.org/10.3390/w13233379. Almeida, C., Mendonça, J.J.L., Jesus, M.R., Gomes, A.J., 2000. SISTEMA AQUÍFERO : QUERENÇA –SILVES ( M5 ). Antunes Aghsaei, H., Mobarghaee Dinan, N., Moridi, A., Asadolahi, Z., Delavar, M., Fohrer, N., Wagner, P.D., 2020. Effects of dynamic land use/land cover change on water resources and sediment yield in the Anzali wetland catchment, Gilan, Iran. Sci. Total Environ. 712, 136449. https://doi.org/10.1016/j.scitotenv.2019.136449. Antunes do Carmo, J.S., 2019. The changing paradigm of coastal management: the Portuguese case. Sci. Total Environ. 695, 133807. https://doi.org/10.1016/j. scitotenv.2019.133807. Anyah, R.O., Forootan, E., Awange, J.L., Khaki, M., 2018. Understanding linkages between global climate indices and terrestrial water storage changes over Africa using GRACE products. Sci. Total Environ. 635, 1405–1416. https://doi.org/10.1016/j.scitotenv. 2018.04.159. APA, 2016. Plano de Gestão de Região Hidrográfica das Ribeiras do Algarve (RH8). Parte 2: caracterização e diagnóstico. ARH Algarve, 2012. Plano de Gestão das Bacias Hidrográficas que Integram a Região Hidrográfica das Ribeiras do Algarve (RH8). Volume I - Relatório. Parte 2- Caracterização Territorial e Fisiográfica: Tomo 1A- peças escritas. Arteche, J., 2021. Bootstrapping long memory time series: application in low frequency estimators. Econ. Stat. https://doi.org/10.1016/j.ecosta.2021.06.002. Barros, C.P., Butler, R., Correia, A., 2010. The length of stay of golf tourism: a survival analysis. Tour. Manag. 31, 13–21. https://doi.org/10.1016/j.tourman.2009.02.010. Canty, A., Ripley, B., 2021. Bootstrap Functions (Originally by Angelo Canty for S). Packag. R. Cardenas-Martinez, A., Rodriguez-Galiano, V., Luque-Espinar, J.A., Mendes, M.P., 2021. Predictive modelling benchmark of nitrate vulnerable zones at a regional scale based on machine learning and remote sensing. J. Hydrol. 603 (Part C), 127092. https://doi. org/10.1016/j.jhydrol.2021.127092. Carlstein, E., 1986a. The use of subseries values for estimating the variance of a general statistic from a stationary sequence. Ann. Stat. 14, 1171–1179. https://doi.org/10.1214/ aos/1176350057. Carlstein, E., 1986b. The use of subseries values for estimating the variance of a general statistic from a stationary sequence. Ann. Stat. 14, 1171–1179. https://doi.org/10.1214/ aos/1176350057. Carriço,N.J.G.,Gonçalves,F.V.,Covas,D.I.C.,doCéuAlmeida,M.,Alegre,H.,2014. Multi-criteria analysis for the selection of the best energy efficient option in urban water systems. Proc. Eng. 70, 292–301. https://doi.org/10.1016/j.proeng.2014. 02.033. Casillas-Trasvina, A., Zhou, Y., Stigter, T.Y., Mussáa, F.E.F., Juízo, D., 2019. Application of numerical models to assess multi-source saltwater intrusion under natural and pumping conditions in the Great Maputo aquifer, Mozambique. Hydrogeol. J. 27 (8), 2973–2992. https://doi.org/10.1007/s10040-019-02053-5. CENSUS, 2021. Resultados Provisórios. https://www.ine.pt/scripts/db_censos_2021.html. (Accessed 1 May 2022). Chernick, M.R., LaBudde, R.A., 2014. An Introduction to Bootstrap Methods with Applications to R. Wiley Publishing, Inc. Chu, C.-S.J., Hornik, K., Kuan, C.-M., 1995. MOSUM Tests for Parameter Constancy. 82. Biometrika, pp. 603–617. https://doi.org/10.2307/2337537. Cleveland, W.S., Devlin, S.J., 1988. Locally weighted regression: an approach to regression analysis by local fitting. J. Am. Stat. Assoc. 83, 596–610. https://doi.org/10.2307/ 2289282. Cleveland, R.S., Cleveland, W.S., Mcrae, J.E., Terpenning, I., 1990. STL: a seasonal-trend decomposition procedure based on loess. J. Off. Stat. 6, 3–73. Coladello, L.F., Galo, M.D.L.B.T., Shimabukuro, M.H., Ivánová, I., Awange, J., 2020. Macrophytes' abundance changes in eutrophicated tropical reservoirs exemplified by Salto Grande (Brazil): trends and temporal analysis exploiting Landsat remotely sensed data. Appl. Geogr. 121, 102242. https://doi.org/10.1016/j.apgeog.2020.102242. Cox, D.R., Stuart, A., 1955. Some quick sign tests for trend in location and dispersion. Biometrika 42 (80), 95. https://doi.org/10.1093/biomet/42.1-2.80. Deng, X., Song, C., Liu, K., Ke, L., Zhang, W., Ma, R., Zhu, J., Wu, Q., 2020. Remote sensing estimation of catchment-scale reservoir water impoundment in the upper Yellow River and implications for river discharge alteration. J. Hydrol. 585, 124791. https://doi.org/ 10.1016/j.jhydrol.2020.124791. Do .A., Ó, Monteiro, J.P., 2006. Estimação da procura real de água no Algarve por sectores, in: Proc. V Congresso Ibérico Sobre Gestão e Planeamento Da Água. dos Santos, F.M., de Oliveira, R.P., Mauad, F.F., 2019. Evaluating a parsimonious watershed model versus SWAT to estimate streamflow, soil loss and river contamination in two case studies in Tietê river basin. J. Hydrol. Reg. Stud. 29, 100685. https://doi.org/10. 1016/j.ejrh.2020.100685. EC, 2007. Drought Management Plan Report, Including Agricultural, Drought Indicators and Climate Change Aspects. Technical Report 2008–023. Luxembourg. EDO, 2020. EDO INDICATOR FACTSHEET: Standardized Precipitation Index (SPI). Fatichi, S., Barbosa, S.M., Caporali, E., Silva, M.E., 2009. Deterministic versus stochastic trends: detection and challenges. J. Geophys. Res. Atmos. 114 (D18). https://doi.org/ 10.1029/2009JD011960. Fragoso, M., Gomes, P.T., 2008. Classification of Daily Abundant Rainfall Patterns and Associated Large-Scale Atmospheric Circulation Types in. 544, pp. 537–544. https://doi.org/ 10.1002/joc. Gilbert, R.O., 1987. Statistical Methods for Environmental Pollution Monitoring. United States. Grath, J., Scheidleder, A., Uhlig, S., Weber, K., Kralik, M., Keimel, T., Gruber, D., 2001. The EU Water Framework Directive: Statistical Aspects of the Identification of Groundwater Pollution Trends, and Aggregation of Monitoring Results. Vienna. Green, T.R., Taniguchi, M., Kooi, H., Gurdak, J.J., Allen, D.M., Hiscock, K.M., Treidel, H., Aureli, A., 2011. Beneath the surface of global change: impacts of climate change on groundwater. J. Hydrol. 405, 532–560. https://doi.org/10.1016/j.jhydrol.2011.05. 002. Guerreiro, S.B., Kilsby, C., Fowler, H.J., 2017. Assessing the threat of future megadrought in Iberia. Int. J. Climatol. 37, 5024–5034. https://doi.org/10.1002/joc.5140. M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 8 Hiscock, H., Sparkes, R., Hodgson, A., Martin, J.L., 2008. Evaluation of future climate change impacts in Europe on potential groundwater recharge. Geophysical Research Abstracts. EGU General Assembly. Hugman, R., Stigter, T.Y., Monteiro, J.P., 2013. The importance of temporal scale when optimising abstraction volumes for sustainable aquifer exploitation: a case study in semi-arid South Portugal. J. Hydrol. 490, 1–10. https://doi.org/10.1016/j.jhydrol. 2013.02.053. IPCC, 2018. Summary for Policymakers Global Warming of 1.5°C. An IPCC Special Report on the Impacts of Global Warming of 1.5°C Above Pre-industrial Levels and Related Global Greenhouse Gas Emission Pathways, in the Context of Strengthening the Global Response to the. IPCC, 2019. Technical Summary, 2019 Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems. IPMA, 2021. Ficha Climatológica Faro/aeroporto. 554 pp. 1971–2000. Kendall, M.G., 1955. Rank Correlation Methods. Griffin, London. Kunsch, H.R., 1989a. The Jackknife and the bootstrap for general stationary observations. Ann. Stat. 17, 1217–1241. Kunsch, H.R., 1989b. The Jackknife and the bootstrap for general stationary observations. Ann. Stat. 17, 1217–1241. Lian, J., Chen, H., Wang, F., Nie, Y., Wang, K., 2020. Separating the relative contributions of climate change and ecological restoration to runoff change in a mesoscale karst basin. CATENA 194, 104705. https://doi.org/10.1016/j.catena.2020.104705. Liu, W., Shi, C., Zhou, Y., 2021. Trends and attribution of runoff changes in the upper and middle reaches of the Yellow River in China. J. Hydro Environ. Res. 37, 57–66. https://doi.org/10.1016/j.jher.2021.05.002. Maia, R., Silva, C., 2009. DSS application at a river basin scale, taking into account water resources exploitation risks and associated costs: the Algarve Region. Desalination 237, 81–91. https://doi.org/10.1016/j.desal.2007.12.024. Mallakpour, I., Villarini, G., 2016. A simulation study to examine the sensitivity of the Pettitt test to detect abrupt changes in mean. Hydrol. Sci. J. 61, 245–254. https://doi.org/10. 1080/02626667.2015.1008482. Mann, H.B., 1945. Nonparametric tests against trend. Econometrika 13, 245–259. Masiliūnas, D., Tsendbazar, N.E., Herold, M., Verbesselt, J., 2021. BFAST lite: a lightweight break detection method for time series analysis. Remote Sens. 13 (16), 3308. https:// doi.org/10.3390/rs13163308. Mckee, T.B., Doesken, N.J., Kleist, J., 1993. The relationship of drought frequency and duration to time scales. Eighth Conference on Applied Climatology,17-22 January 1993. Anaheim, California. Mendes, M.P., Matias, M., Gomes, R.C., Falcão, A.P., 2021. Delimitation of low topsoil moisture content areas in a vineyard using remote sensing imagery (Sentinel-1 and Sentinel-2) in a Mediterranean-climate region. Soil Water Res. 16, 85–94. https://doi.org/10.17221/ 101/2019-SWR. Mendes, M.P., Ribeiro, L., 2014. The importance of groundwater for the delimitation of Portuguese National Ecological Reserve. Environ. Earth Sci. 72. https://doi.org/10.1007/ s12665-013-3039-y. Młyński, D., Wałęga, A., Kuriqi, A., 2021. Influence of meteorological drought on environmental flows in mountain catchments. Ecol. Indic. 133, 108460. https://doi.org/10.1016/j. ecolind.2021.108460. Neves, M.C., Costa, L., Monteiro, J.P., 2016. Climatic and geologic controls on the piezometry of the Querença-Silves karst aquifer, Algarve (Portugal). Hydrogeol. J. 24, 1015–1028. https://doi.org/10.1007/s10040-015-1359-6. Neves, M.C., Jerez, S., Trigo, R.M., 2019. The response of piezometric levels in Portugal to NAO, EA, and SCAND climate patterns. J. Hydrol. 568, 1105–1117. https://doi.org/10. 1016/j.jhydrol.2018.11.054. Nicolau, R., 2002. Modelação e Mapeamento da Distribuição Espacial da Precipitação-Uma Aplicação a Portugal Continental. Universidade Nova de Lisboa, Lisbon. Pandey, K., Kumar, S., Malik, A., Kuriqi, A., 2020. Artificial neural network optimized with a genetic algorithm for seasonal groundwater table depth prediction in Uttar Pradesh, India. Sustainability 12 (21), 8932. https://doi.org/10.3390/su12218932. Peña-Angulo, D., Vicente-Serrano, S.M., Domínguez-Castro, F., Murphy, C., Reig, F., Tramblay, Y., Trigo, R.M., Luna, M.Y., Turco, M., Noguera, I., Aznárez-Balta, M., García-Herrera, R., Tomas-Burguera, M., Kenawy, A. El, 2020. Long-term precipitation in Southwestern Europe reveals no clear trend attributable to anthropogenic forcing. Environ. Res. Lett. 15 (9), 94070. https://doi.org/10.1088/1748-9326/ab9c4f. Pettitt, A.N., 1979. A non-parametric approach to the change-point problem. J. R. Stat. Soc.: Ser. C: Appl. Stat. 28, 126–135. https://doi.org/10.2307/2346729. Pinto, F.S., de Carvalho, B., Marques, R.C., 2021. Adapting water tariffs to climate change: Linking resource availability, costs, demand, and tariff design flexibility. J. Clean. Prod. 290, 125803. https://doi.org/10.1016/j.jclepro.2021.125803. R Core Team, 2018. R: A Language and Environment for Statistical Computing. Saxena, R., Watson, L.T., Wynne, R.H., Brooks, E.B., Thomas, V.A., Zhiqiang, Y., Kennedy, R. E., 2018. Towards a polyalgorithm for land use change detection. ISPRS J. Photogramm. Remote Sens. 2014, 217–234. https://doi.org/10.1016/j.isprsjprs.2018.07.002. Secci, D., Tanda, M.G., D’Oria, M., Todaro, V., Fagandini, C., 2021. Impacts of climate change on groundwater droughts by means of standardized indices and regional climate models. J. Hydrol. 603, 127154. https://doi.org/10.1016/j.jhydrol.2021.127154. Sellinger, C.E., Stow, C.A., Lamon, E.C., Qian, S.S., 2008. Recent Water Level Declines in the Lake Michigan−Huron System. Environ. Sci. Technol. 42, 367–373. https://doi.org/10. 1021/es070664+. Sharma, S., Swayne, D.A., Obimbo, C., Ess, L.O.W., 2016. Trend analysis and change point techniques : a survey. Energy, Ecol. Environ. 1, 123–130. https://doi.org/10.1007/ s40974-016-0011-1. Sondermann, M.N., de Oliveira, R.P., 2021. A shared vision on the transboundary water management challenges of the Tagus River Basin. Water Resour. Manag. 35 (13), 4647–4664. https://doi.org/10.1007/s11269-021-02973-6. Sousa, P.M., Trigo, R.M., Aizpurua, P., Nieto, R., Gimeno, L., Garcia-Herrera, R., 2011. Trends and extremes of drought indices throughout the 20th century in the Mediterranean. Nat. Hazards Earth Syst. Sci. 11, 33–51. https://doi.org/10.5194/nhess-11-33-2011. Stigter, T.Y., Ribeiro, L., Dill, A.M.M.C., 2006. Evaluation of an intrinsic and a specificvulnerability assessment method in comparison with groundwater salinisation and nitrate contamination levels in two agricultural regions in the south of Portugal. Hydrogeol. J. 79–99. https://doi.org/10.1007/s10040-004-0396-3. Stigter, T.Y., Monteiro, J.P., Nunes, L.M., Vieira, J., Cunha, M.C., Ribeiro, L., Nascimento, J., Lucas, H., 2009. Screening of sustainable groundwater sources for integration into a regional drought-prone water supply system. Hydrol. Earth Syst. Sci. 1185–1199. Teegavarapu, R.S.V., 2019. Chapter 1 - methods for analysis of trends and changes in hydroclimatological time-series. In: Teegavarapu, R. (Ed.), Trends and Changes in Hydroclimatic Variables. Elsevier, pp. 1–89 https://doi.org/10.1016/B978-0-12-810985-4. 00001-3. Verbesselt, J., Hyndman, R., Newnham, G., Culvenor, D., 2010a. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 114, 106–115. https://doi. org/10.1016/j.rse.2009.08.014. Verbesselt, J., Hyndman, R., Zeileis, A., Culvenor, D., 2010b. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 114, 2970–2980. https://doi.org/10.1016/j.rse.2010.08.003. Verbesselt, J., Zeileis, A., Hyndman, R., 2015. Breaks For Additive Season and Trend (BFAST). Vicente-Serrano, S.M., Lopez-Moreno, J.-I., Beguería, S., Lorenzo-Lacruz, J., Sanchez-Lorenzo, A., García-Ruiz, J.M., Azorin-Molina, C., Morán-Tejeda, E., Revuelto, J., Trigo, R., Coelho, F., Espejo, F., 2014. Evidence of increasing drought severity caused by temperature risein southern Europe. Environ. Res. Lett. 9, 44001. https://doi.org/10.1088/1748-9326/9/4/ 044001. Vicente-Serrano, S.M., Peña-Gallardo, M., Hannaford, J., Murphy, C., Lorenzo-Lacruz, J., Dominguez-Castro, F., López-Moreno, J.I., Beguería, S., Noguera, I., Harrigan, S., Vidal, J.-P., 2019. Climate, irrigation, and land cover change explain streamflow trends in countries bordering the Northeast Atlantic. Geophys. Res. Lett. 46 (19), 10821–10833. https://doi.org/10.1029/2019GL084084. Wang, M., Jiang, S., Ren, L., Xu, C.Y., Menzel, L., Yuan, F., Xu, Q., Liu, Y., Yang, X., 2021. Separating the effects of climate change and human activities on drought propagation via a natural and human-impacted catchment comparison method. J. Hydrol. 603, 126913. https://doi.org/10.1016/j.jhydrol.2021.126913. WCAP, 2004. Development, Use and Application of the Hydrospect Data Analysis System for the Detection of Changes in Hydrological Time Series for Use in WCP-WATER and National Hydrological Services WCASP - 65. WMO/TD-No. 1240. WMO, 2012. Standardized Precipitation Index: User Guide. Xie, H., Li, D., Xiong, L., 2014. Exploring the ability of the Pettitt method for detecting change point by Monte Carlo simulation. Stoch. Env. Res. Risk A. 28, 1643–1655. https://doi. org/10.1007/s00477-013-0814-y. Xu, H., Wang, X., Zhao, C., Shan, S., Guo, J., 2021. Seasonal and aridity influences on the relationships between drought indices and hydrological variables over China. Weather Clim. Extrem. 34, 100393. https://doi.org/10.1016/j.wace.2021.100393. Zeleňáková, M., Vido, J., Portela, M.M., Purcz, P., Blištán, P., Hlavatá, H., Hluštík, P., 2017. Precipitation trends over Slosvakia in the period 1981–2013. Water 9 (12), 922. https://doi.org/10.3390/w9120922. M.P. Mendes et al. Science of the Total Environment 846 (2022) 157428 9