1 Introduction

Over the period 1880 to 2012, almost the entire globe has experienced surface warming (Intergovernmental Panel on Climate Change (IPCC 2013)), with climate change expected to affect hydrological regimes at a global level (Arnell 1999). The changing patterns of extreme weather and climatic events are linked to an increase in the magnitude and frequency of natural hazards (Devkota and Bhattarai 2015), that will have a negative impact on global food security, since they may influence important crops that humanity depends on, such as wheat (Li et al. 2016) and rice (Huang et al. 2018).

The South Aegean Sea is a Greek sea that is located between the Greek mainland and Turkey and is dominated by 385 islands, which comprise the Region of South Aegean; at least 29 of these islands face water scarcity problems (Kaldellis and Kondili 2007). Additionally, since this region constitutes a top tourist destination worldwide, the latter problems are significantly intensified during summer months due to large tourist influx (HSA, 2020). To alleviate the problem, the Greek Government has been hiring tanker fleets all year round to transport potable water from the Greek mainland to the islands up to 2019, at an average annual cost of about 8.2 €/m3 (SGAIP, 2019a).

The purpose of this paper is to examine whether the impact of changing climate patterns, which may be attributed to the ongoing global warming and the tourist surge, are associated with the depletion of water resources in the Region of South Aegean. This phenomenon is evident due to the increasing need to transfer drinking water to the study area by utilising ships. The tendency of climatic extremes and the precipitation, temperature, drought, aridity and tourism time series for the Region of South Aegean were captured by using the Mann–Kendall trend detection test (M–K test) and the relationship between the aforementioned variables and the water transported by ships to the Region of South Aegean was investigated by using the Pearson (r) correlation coefficient. Moreover, Pettitt’s (1979) test was performed to detect the exact time that an abrupt change had occurred within the aforementioned time series.

To the best to the author's knowledge, the linkage between the ongoing changing climate patterns, tourism growth and the pressure on the regional capacities to cope with the increasing demand for water resources has garnered much global attention. For instance, Martínez-Ibarra (2015) attempted to explore the connection, but this research was based on local press releases, whereas the majority of the studies are focused on identifying the potential impact of the climate change on tourism in the future and the stress on water resources by using indexes (Sauri et al., 2013; Roson and Sartori, 2014) or hydrologic models (Klein et al., 2015). Furthermore, it should be highlighted that the applied methodology in this study is not only confined to the specific analysed dataset but it can also be utilised in any other water-stressed area of the world to capture the changes in climate patterns and tourism and their concomitant impact on any other significant water variables, such as the reservoir water levels, the river discharge magnitude, the groundwater table depth and other factors. Such studies would be extremely useful to enhance the awareness of the local communities and, subsequently, design an effective sustainable water management plan.

A thorough investigation of extreme weather and climatic events is provided through 27 climatic indices proposed by the Expert Team on Climate Change Detection and Indices (ETCCDI), which are computed based on daily precipitation, maximum and minimum temperature records, and are used to track trends of climate change. Τhe M–K trend detection test (Mann 1945; Kendall 1962) and the Theil–Sen estimator (Theil 1950; Sen 1968) are commonly applied to the aforementioned indices, so as to detect potentially significant trends and identify the magnitude of these trends respectively at a selected significance level (Kioutsioukis et al. 2010; Devkota and Bhattarai 2015; Li et al. 2016; Huang et al. 2018). The M–K test is considered the most popular trend detection test and has been widely applied so as to track such trends in hydro-meteorological time series (Gao et al. 2011; Lakhraj-Govender et al. 2017; Myronidis et al. 2018a,b).

The frequency and intensity of drought in the Mediterranean basin is likely to have increased and it is projected that this area will suffer from drought due to a global temperature increase (IPCC 2013) that could lead to higher aridity levels (Zarch et al. 2017). The tracking of drought severity can be performed through the Reconnaissance Drought Index (RDI) (Tsakiris and Vangelis 2005; Tsakiris et al. 2007) which is considered to be among the most credible and sensitive to climate variability (Khanmohammadi et al. 2018). It has low data requirements and is easy to compute (Mohammed and Scholz 2018). The latter advantages induced the authors to select the RDI index through its standardised form RDIst (Tigkas et al. 2012) in order to assess drought magnitude over the region (Myronidis et al. 2012). Additionally, as an extra validation, drought magnitude was also computed using the Standardized Precipitation Index (SPI) (McKee et al., 1993) at the same annual scale as the RDIst and the two drought indices outputs were then compared with each other. The SPI is the most popular meteorological drought index worldwide, which is promoted by the WMO organization to be used by the National Meteorological and Hydrological Services (NMHSs) worldwide so as to characterize drought (Hayes et al., 2011).

The overall expression and significance of climatic change in bio-climatic terms is best expressed by the Aridity Index (AI) defined by the ratio of precipitation (P) to potential evapotranspiration (PET), whilst PET is preferably estimated, in terms of accuracy, by using the Penman (1948) formula (Kafle and Bruins 2009). However, since the necessary meteorological parameters for the calculation of the Penman formula are usually lacking for most parts of the world, the calculation of PET in the AI formula is usually carried out by using the simpler Thornthwaite approach (Nastos et al. 2013; Cheval et al. 2017; Şarlak and Agha 2018). This simplified AI calculation method has been adopted by both the United Nations Environment Programme (UNEP 1992) and the Food and Agriculture Organization (FAO 1993).

In a comparison with similar research in the study area, Kioutsioukis et al. (2010) investigated climate extremes over Greece for the period 1955–2002 and Nastos et al. (2013) captured the temporal variability of the AI index in Greece during the period 1951 to 2000. The present study focuses on the South Aegean region, accounts for climatic extreme indices and the AI index computation for the recent period 1991–2020 and considers four additional meteorological stations compared to the aforementioned studies. The water shortage problem in the Aegean archipelago islands has also been highlighted by Kaldellis and Kondili (2007) but their study did not examine the association between climatic variables and water transfers. Finally, although Tigkas (2008) captured drought conditions in the Cyclades region, he only used one representative station and his analysis was focused on a previous period (1955–2002).

2 Study area and description of datasets

The South Aegean region features more than 385 magnificent islands covering an area of 5286 km2 with a total coastal length of sandy beaches equal to 4223 km (Fig. 1). The islands are divided into the Dodecanese and Cyclades prefectures and, as the 2011 official population census revealed, are populated by 308,975 inhabitants; however, in summertime, their population may be five times higher than in winter (Kaldellis and Kondili 2007).

Fig. 1
figure 1

Location map of the selected observation stations within the South Aegean Sea (Station index numbers are based on WMO standards)

Full size image

The aforementioned islands face water resource shortage problems, which are intensifying over time, mainly due to the limited and non-uniform rainfall distribution, as well as the high temperatures they experience during summer, which result in high evapotranspiration rates (Fig. 2). Additionally, only four islands, namely Naxos, Andros, Rhodes and Kos seem to have sufficient local water resources (Kaldellis and Kondili 2007), while during the period 1997–2018, 17.8 Mm3 of drinking water were transported to the Region with ships from the mainland at a total cost of 124.2 million euros (SGAIP, 2019a). The drinking water quantities were transferred to the Region by private ships hired out by the SGAIP and each have about 2100 m3 water tank capacity.

Fig. 2
figure 2

Observation station’s climograph covering the following period: 1991 – 2020

Full size image

The water transfer to the islands at the Prefecture of Cyclades was done from the urban water facilities of the Athens Water Supply and Sewerage Company (EYDAP S.A.) at Lavrio port, whilst the water transfer to the islands at the Prefecture of Dodecanese was done from the water drillings located at Kalathos village (Fig. 1). The transported water was then conveyed from the vessels through closed pipelines to appropriate storage facilities in the ports through which it was pumped to the arid islands central water tanks. Furthermore, due both to the growth in tourism and the limited rainfall in recent years, the water demand now exceeds the water availability, as a result of which the Region has been suffering from intense water shortages.

Daily time series of precipitation totals, along with maximum and minimum temperatures, were obtained from six meteorological stations (M.S.) operated by the Hellenic National Meteorological Service (HNMS); their spatial locations in the region are displayed in Fig. 1. The stations are identified by their World Meteorological Organization (WMO) number and name and information is provided regarding their location, the altitude of each station and the mean annual temperature and precipitation. They are situated in lowland areas, their altitude ranging from 7 to 167 m (Table 1).

Table 1 List of Meteorological Stations analyzed in this study
Full size table

The stations’ datasets have different starting dates and a common continuous reference period spanning precisely over the past 30 years, from 01/01/1991 to 31/12/2020, with the exception of Milos M.S. and Naxos M.S. which have a dataset ending in 12/2010 and 12/2017, respectively. The missing dataset for Milos for 01/2011—06/2012 was provided through daily measurements from two neighbouring stations, Santorini M.S. and Naxos M.S., following the application of the normal ratio (NR) method. The NR method was selected so as to estimate the missing daily precipitation and daily temperature values since similar studies have found that it outperformed other daily rainfall infilling methods during the use of two to five index gauges (Mair and Fares, 2010). It is also a recommended infilling method for timeseries of daily temperature data (Shabalala et al. 2019) and is expressed as follows:

$$P_{x} = \frac{1}{n}\sum\limits_{i = 1}^{i = n} {\frac{{N_{x} }}{{N_{i} }}} P_{i}$$
(1)

where Px is the missing value of precipitation or temperature at meteorological station x, Pi is the observed precipitation or temperature value for the same period at the ith neighbouring station, Nx is the normal annual precipitation or temperature of the station x, Ni is the normal annual precipitation or temperature of the ith station and n is the number of neighbouring stations. Moreover, the Milos M.S. dataset for the period 07/2012—12/2020 was obtained from a new automated weather station located 9 km away from the old HNMS station and is operated by the National Observatory of Athens (NOA) (Lagouvardos et al. 2017).

The measurements for the period 01/2018—12/2020 for Naxos M.S. were obtained from another NOA automated weather station in the island located 12 km away to the HNMS station at about the same elevation. Additionally, the NR method was employed in order to fill some missing temperature records for Santorini M.S., Mykonos and Naxos M.S. whilst the percentage of incomplete records for these stations were 1.11%, 1.67% and 4.44%, respectively. Furthermore, the 30 years record length of the dataset is considered sufficient in order to ensure the statistical validity of the trend test (Burn and Elnur 2002) as well as to track drought (Horion et al. 2012).

A homogeneity assessment of the daily precipitation and minimum and maximum temperature datasets was performed using RHtest V4 software, which is freely available from http://etccdi.pacificclimate.org (Wang and Feng 2013). The latter software was used without a reference time series in order to locate step changes in the record that are significant at the 95% confidence level. Figure 3 illustrates the precipitation homogeneity test for Milos M.S., in which the continuous red line implies that the data is homogenous while the smallest shift in 2011 is not statistically significant at the 95 per cent uncertainty range.

Fig. 3
figure 3

Precipitation homogeneity test for Milos M.S

Full size image

The same homogeneity test was applied for all understudy variables in all six meteorological stations. The applied homogeneity testing revealed no sudden jumps in the entire data series, which denotes that the climatic series under study can be considered homogenous.

By employing the Thiessen polygon method (Myronidis et. al. 2012), the areal mean annual rainfall and temperature were found to be equal to 492.5 mm and 18.7 °C respectively whilst, based on the Koppen (1936) classification, the study area climate is characterized as Mediterranean with mild winters and long dry summers. The average summer precipitation is less than 2.5 mm and the rainiest month is January with 100.7 mm, while the mean monthly temperature ranges between 11.8 °C in January and 26.3 °C in August. Finally, the annual aggregate volumes of water transferred to the region and the related cost for the period 1997–2018 were collected from the SGAIP (2019a) and are presented in Fig. 4.

Fig. 4
figure 4

Water transported to the region, aggregate volume and total cost (SGAIP, 2019a)

Full size image

At the Region of South Aegean, the irrigation in agricultural lands and the domestic water use are based on the groundwater abstractions. The latter mainly cause the drain of groundwater level, the intrusion of seawater in aquifers and the important deviation of drinking water from good quality, which is illustrated as excess levels on water quality parameters beyond the maximum admissible concentrations standards (Karavoltsos et al. 2008), as defined by the drinking water directive (EU 1998).

Between 1991 and 2019, the number of tourist arrivals (NA) and the number of the nights that they spent in hotels as well as similar establishments (NS) at the South Aegean Region increased more than forth-fold and triple-fold, respectively, whilst the relevant annual registry is displayed on Fig. 5. However, it should be marked that the 2020 tourist records are not considered in this study, due to COVID-19, as well the number of tourists arriving by cruise ships in the Region, which touched 2.26 million in 2019 (SETE, 2020), since there are no track records prior to 2013.

Fig. 5
figure 5

Tourist arrivals and the nights that they spent in hotels and similar establishments in the Region between 1991 and 2019

Full size image

Besides the tourism sector, the agriculture sector also uses an essential segment of the available water resources. However, since these data are scarce, just seven farm structure surveys from 1993–2016, they could not be used to estimate trends with the M–K test or for correlation analysis. The latter datasets (HSA, 2020), as a simple indication, revealed a strong negative trend on the irrigated land area between 1993 and 2016 (Fig. 6).

Fig. 6
figure 6

Irrigated land area in the Region between 1993 and 2016

Full size image

3 Methodology

3.1 Estimation of climate extreme indices

The recommended 27 climate extreme indices by ETCCDI (Table 2) which are suitable for climate change monitoring can be separated into two categories: 16 temperature related indices and 11 precipitation related indices. They can be obtained from daily temperature and precipitation datasets respectively, using the RClimDex V1 software (Gao and Shi 2015; Zhang et al. 2005); the latter software is used prior to the indices’ computation in order to quality control the precipitation and temperature time series (Abbas et al. 2018). Additional details regarding the indices’ definitions can be found in (Sillmann et al. 2013; Wang et al. 2013), whilst Tmax and Tmin in Table 2 are the abbreviations of the daily maximum temperature and daily minimum temperature respectively. Furthermore, in Table 2, the applied rainstorm events’ classification in the form of R10mm, R20mm and R25mm indices refer to the annual count of days where daily precipitation is more than 10, 20 and 25 mm, respectively (Li et al., 2016; Zhang et al., 2014). Once the 27 climatic indexes were derived for each meteorological station, the regional average climate extreme indices time series were computed as an unweighted mean of the stations’ values (Wang et al. 2013; Zhang et al. 2014).

Table 2 List of the 27 ETCCDMI core climate indices
Full size table

3.2 Aridity and drought estimation

In order to calculate aridity and drought, monthly PET values were required and estimated using the Thornthwaite (1948) formula:

$${\text{PET}} = 16 \cdot F \cdot \left( {\frac{{10T_{m} }}{J}} \right)^{a}$$
(2)

where Tm is the average monthly temperature (°C), F is the adjustment factor that depends on latitude and the month of the year, J is the heat index which is calculated based on the Tm values and a is a coefficient that depends on J; more details regarding the calculation of the heat index and the a coefficient can be found in (Vangelis et al. 2013).

The temporal evolution of the aridity index on an annual scale was examined by using the AI index (UNEP 1992; FAO 1993), which is expressed as follows:

$$AI = \frac{P}{PET}$$
(3)

where P and PET are the annual total precipitation (mm) and the potential annual evapotranspiration (mm) values respectively.

The RDI index was used in its standardized form RDIst to track the drought phenomena over the region. Initially, the value a0 was computed for the i-th year on an annual basis as follows (Tsakiris et al. 2007):

$$a_{0}^{(i)} = \frac{{\sum\nolimits_{j = 1}^{12} {P_{ij} } }}{{\sum\limits_{j = 1}^{12} {PET_{ij} } }},\;{\text{i}} = {\text{1 - N and j }} = { 1 - 12}$$
(4)

where Pij and PETij are the precipitation and potential evapotranspiration of the j-th month of the i-th year, and N is the total number of years of available data.

The RDIst, in case that the lognormal distribution is applied, can be computed as follows:

$$RDI_{st}^{(i)} = \frac{{y^{(i)} - \overline{y}}}{{\hat{\sigma }y}}$$
(5)

where RDIst is the standardised RDI, \(y^{(i)}\) is the \(\ln (a_{0}^{(i)} )\), \(\overline{y}\) and \(\hat{\sigma }y\) are its arithmetic mean and standard deviation, respectively.

Although the RDI computation is based on P and PET values, no significant influence has been reported on RDI by any of the selected Hargreaves, Thornthwaite, Blaney-Criddle and FAO Penman–Monteith PET estimation methods (Vangelis et al. 2013). Moreover, in spite of the fact that the a0 value satisfactorily follows both the lognormal and the gamma distributions, the RDIst could be better performed in several cases by fitting the gamma probability density function (pdf) to the given frequency distribution of a0 (Tigkas et al. 2012). More details about this process can be found in Tsakiris et al. (2007) and Vangelis et al. (2013).

The SPI index estimation compared to the RDI index calculation, which is based on the ratio P/PET, is less data demanding since it is based solely on long term monthly precipitation records, whilst both indices could assess drought on multiple timescales. In order to calculate the SPI index, the gamma probability density function is initially fitted on the frequency distribution of the aggregate monthly precipitation. Then, through appropriate mathematical transformations, the gamma cumulative probability distribution is converted to the standardized normal distribution of Z with mean and standard deviation equal to 0 and 1 respectively, while the Z value represents the SPI index. For more details regarding the SPI derivation, please refer to Edwards and McKee (1997). Furthermore, the RDI and the SPI indices were computed for an accumulation period of 12 months (RDI-12, SPI-12) by using the DrinC Software developed by the National Technical University of Athens (http://drought-software.com/).

Once the annual values of P, PET, AI, RDIst and SPI were derived for each meteorological station, the relevant regional average time series was calculated by considering the region as one catchment area and using the Thiessen polygon methodology (Tsakiris et al. 2007; Pashiardis and Michaelides 2008). Additionally, the basic linear regression test was used (Feidas 2007) so as to identify the linear patterns of drought and aridity. Finally, the dryland zone classification, as defined by the AI index (UNEP 1992; FAO 1993), and the categorization of drought severity based on both the RDIst (Vangelis et al. 2013) and the SPI (McKee et al., 1993) indices, are illustrated in Table 3.

Table 3 Dryland zone categories based on the AI and drought categorization based on RDIst and SPI indices
Full size table

3.3 Trend detection analysis and abrupt change detection (ACD)

The M–K test was selected so as to investigate potential trends in the climatic variables under study and in order to test the statistical significance of the tracked trends at the 95% confidence level (Myronidis et al. 2018a,b). The Z statistical sign (positive or negative) indicates whether there is an increasing or decreasing trend respectively and is computed as follows:

$$Z = \left\{ \begin{gathered} \frac{S - 1}{{\sqrt {{\text{VAR}}(S)} }},\;if\;S > 0 \hfill \\ \;\;0,\;\;\;\quad \quad \;if\;S = 0 \hfill \\ \frac{(S + 1)}{{\sqrt {{\text{VAR}}(S)} }},\;if\;S < 0 \hfill \\ \end{gathered} \right.$$
(6)

in which,

$$S = \sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^{n} {{\text{sgn}} (x_{j} - x{}_{k}} )}$$
(7)
$${\text{sgn}} (x_{j} - x_{k} ) = \left\{ \begin{gathered} + 1,\;if\;x_{j} - x_{k} > 0 \hfill \\ \;\;0,\;if\;x_{j} - x_{k} = 0 \hfill \\ - 1,\;if\;x_{j} - x_{k} < 0 \hfill \\ \end{gathered} \right.$$
(8)
$${\text{VAR}}(S) = \frac{{n \cdot (n - 1) \cdot (2n + 5) - \sum\nolimits_{i = 1}^{m} {t_{i} \cdot (t_{i} - 1) \cdot (2t_{i} + 5)} }}{18}$$
(9)

where xj and xk are the data values for the years j and k respectively, with j > k; n is the total number of years, sgn() is the sign function, VAR(S) represents the variance of S which has a zero mean for n ≥ 8 and follows the normal distribution, ti represents the number of ties for the ith value and m represents the number of tied data (Amirataee et al. 2016).

According to this test, the null hypothesis H0 indicates no trend in the data, and the alternative hypothesis HΑ shows that the data follows a monotonic trend over time. The null hypothesis is either accepted or rejected if the statistical value of Z is lower or higher respectively than the value of the normal distribution function Za/2 for the selected significance level a = 0.05, which is obtained from the normal distribution tables and is equal to 1.95996. Furthermore, the trend’s magnitude was obtained by using the Theil–Sen approach (Theil 1950; Sen 1968). Compared to linear regression, the latter method is preferable since it limits the outliers’ influence on the computed slope (Hirsch et al. 1982). Furthermore, before applying the M–K test, all datasets have to be checked for the presence of autocorrelation using a two-tailed test and the lag-1 autocorrelation coefficient (r1) at a 5% significance level must be calculated (Myronidis et al. 2018a). The effect of serial correlation can then be eliminated by applying the modified Mann–Kendall test (Hamed and Rao 1998); more details regarding the application of the modified M–K test can be found in Amirataee and Montaseri (2017).

The Pettitt (1979) test provides the most probable change point when the exact time of the change is unknown (Caloiero et al., 2010). In this test, the null hypothesis H0 indicates no abrupt change in the time series, whilst the alternative HA shows that the time series has a change-point. The test statistic Ut, is given by

$$U_{t} = \sum\limits_{i = 1}^{t} {\sum\limits_{j = t + 1}^{n} {{\text{sgn}} (x_{i} - x_{j} )} }$$
(10)

where sgn (x) = 1 if x > 0, 0 if x = 0, -1 if x < 0.

The most significant change point is located at time t, when \(\left| {U_{t} } \right|\) is at its maximum, and is provided by Eq. 11, whilst the related confidence level (ρ) for the sample length (n) is computed using Eq. 12:

$$K_{n} = \mathop {\max }\limits_{1 \le t \le n} \left| {U_{t} } \right|$$
(11)
$$\rho = \exp \left( {\frac{{ - 6K^{2} }}{{n^{3} + n^{2} }}} \right)$$
(12)

when the ρ value is greater than the specific confidence level, the null hypothesis is rejected and the change point is considered significant; more details about the method can be found in Jaiswal et al. (2015).

3.4 Relationship between water transfers, climatic variables, tourism variables and indices

With a view to tracking the strength of the relationship between the aggregate annual volume (V) and cost (C) of the water transported to the Region with the climatic extreme indices, precipitation, temperature, aridity, drought and tourism variables for the period 1997–2008, the Pearson correlation coefficient (r) (Myronidis et al. 2018b) was utilized. The latter period was selected in order to capture the aforementioned relationship, since 32 desalination plants were constructed in the region after 2008, which provide about 48,000 m3/day (SGAIP, 2019b). Due to the desalination plants operation, water transfers into the region were drastically reduced after 2008. Therefore, correlating these water transfers to the region with the understudy variables after 2008 was not considered statistically valid. Moreover, during 2019 and 2020, some limited amounts of water were transported to only three remote islands of the Region of South Aegean by the Hellenic navy, where the construction of desalination units continues to be in progress.

Prior to the estimation of the r coefficient magnitude between the climatic variables, indices and the water transfer dataset under study, it was necessary to reduce the impact of potentially confounding non-climatic factors by applying the first differential technique (Yang et al. 2015). According to this technique, the differences in values between two consecutive years were initially calculated for the whole dataset and then the Pearson’s coefficient magnitude was computed for all possible relationships; more details about this method can be found in Huang et al. (2018) and Li et al. (2016).

4 Results and discussion

4.1 Trends of temperature extreme indices and ACD

Initially, the lag-1 autocorrelation magnitude was computed for the 16 temperature extreme indices. It was observed that only the TR, TN10p, TX10p, TN90p and DTR had an r1 value exceeding the confidence interval bounds (0.32 to −0.39). Thus, the latter variables are considered to be serially dependent and the modified Mann–Kendall test by Hamed and Rao (1998) was applied to them, whilst the M–K test standard form was utilized for all the remaining variables (Myronidis et al. 2018a,b). In all cases, the 95% significance level was selected, which is typical for trend analysis studies of climatic extreme indices (Sensoy et al., 2013; Dashkhuu et al. 2015). The results of the M–K and Pettitt tests are presented in Table 4.

Table 4 M–K test Z value, Sen’s slope β value and Pettitt’s test abrupt change timing
Full size table

Table 4 illustrated insignificant positive trends regarding the growing season length (GLS), the monthly maximum value of daily minimum temperature (TNx), and the monthly minimum value of daily minimum temperature (TNn). Meanwhile, the number of ice days (ID) in the area is equal to zero for the entire climatic archive. Insignificant negative trends were observed in relation to the indices for the number of frost days, summer days, the monthly maximum value of daily maximum temperature (TXx), and the monthly minimum value of daily maximum temp (TXn). Insignificant trends for the number of summer days were also captured by Kioutsioukis et al. (2010) for a small fraction of the study area. Insignificant decreasing trends for winter and insignificant increasing trends for the remaining seasons of the year were observed for the neighboring coastal areas of Turkey by Turkes et al. (1996) but they concern an older period (1930–1993) of their study.

Significant positive trends were noted concerning the number of tropical nights, warm nights, warm days and the warm spell duration index, which have increased by almost 9 days, 4 days, 2 days and 1.4 days per decade, whilst Pettitt’s test identified an abrupt change in the aforementioned variables between 2006/07. A statistically significant increase of the latter variables in neighboring Turkey by 3.7, 1.5, 1.4 days and 2.3 days per decade respectively was also reported by Sensoy et al. (2013), who highlighted that the greatest trends involved the coastal areas. The diurnal temperature range showed significant positive trends by 0.2 °C per decade with a step change at 2006, whilst long-term temperature data in neighboring Cyprus revealed a significant increase of the DTR by 0.1 °C per decade (Price et al. 1999). Significant negative trends were observed regarding the number of cool nights, cool days and the cold spell duration index (CSDI), which has decreased by 4.8 days, 3.8 days and 1.1 days per decade respectively, while the change-point year for these variables was 2006. A significant decrease of the cool night, cool day and CSDI figures was also observed by Sensoy et al. (2013) in neighboring Turkey’s coastal areas.

4.2 Trends of precipitation extreme indices and ACD

The M–K and Pettitt test processes and runs were repeated in order to identify any changes in the precipitation extreme indices, with the relevant results displayed in Table 5.

Table 5 M–K test Z value, Sen’s slope β value and Pettitt’s test abrupt change timing
Full size table

Table 5 illustrated that 10 of the 11 precipitation related indices exhibited a downward trend, thus highlighting that drier conditions have been prevailing over the region during the recent decades. The same negative trends in all precipitation extremes indices over the eastern Mediterranean region have also been marked by Kostopoulou and Jones (2005). Sensoy et al. (2013) also reported a decrease in the number of heavy precipitation days and the PRCPTOT for Turkey’s coastal areas, whilst Kioutsioukis et al. (2010) also tracked a decrease in the annual total wet-day precipitation. Moreover, significant positive trends were not noted in the region, while significant negative trends were observed for the number of consecutive dry days which were reduced by about −12 days per decade, with the timing of the abrupt change being the year 2001. Nastos and Zerefos (2009) captured a statistically significant downward trend for the CDD by -7.6 days per decade for the Rhodes M.S. whilst Kioutsioukis et al. (2010) found no significant trend for the R95pTOT index. Their data however concerns an older period (1956–2001) and only a small fraction of the region.

4.3 Precipitation, temperature, aridity, drought, tourism trends and ACD

The r1 value is within the confidence interval for the precipitation, aridity and drought variables. Therefore, the M–K test standard form was applied for the latter variables. Moreover, since the r1 value for the tourism variables and temperature was found to transcend the upper confidence bound, 0.33 and 0.32 respectively, at the significant level of 95%, the modified M–K test by Hamed and Rao (1998) was applied for these variables and the emerged results are illustrated on Table 6.

Table 6 M–K test Z value, Sen’s slope β value
Full size table

Table 6 presents an insignificant downward trend for precipitation and an insignificant upward trend for the temperature with a magnitude of −34.7 mm and 0.1 °C per decade respectively. Annual temperature warming trends of no statistical significance were also reported in the region by Mamara et al. (2016), while the annual temperature increase by 0.1 °C per decade is in good accordance with the IPCC (2013) calculation of 0.12 °C per decade for the period 1951–2012. Insignificant decreasing trends of precipitation for the central Aegean Sea were also documented by Feidas (2007).

The AI illustrates an insignificant negative trend, while the basic linear regression test (Feidas 2007) reveals that the region’s climate characterization is changing from dry sub-humid to semi-arid conditions (Fig. 5a). Türkeş (2003) by using the AI index in neighboring Turkey’s coastal area documented a significant change from humid conditions to dry sub-humid or semi-arid climatic conditions. Meteorological stations also documented a shift of the climate from humid conditions to dry sub-humid or semi-arid climatic conditions, whilst the same tendency over the Cyclades prefecture has also been reported by Nastos et al. (2013). The two drought indices revealed normal to mild dry conditions over the Region for most of the time. Insignificant negative trends that suggest that drought phenomena have slightly intensified over the Region, with the linear regression tests illustrating the same shift (Fig. 7b). These findings are in accordance with the global drought patterns which suggest that there have been insignificant and minor changes in drought over the recent decades. Moreover, Tigkas (2008), by using both linear regression and the Mann–Kendall test, noted that there is no statistically significant trend of drought over any one station in the Cyclades prefecture, while Giannikopoulou et al. (2014) observed that mild drought conditions were prevalent in three Cycladic islands.

Fig. 7
figure 7

AI index a and RDIst and SPI indices b variation and their linear trend

Full size image

The modified M–K trend test revealed significant positive trends at the 0.05 significance level for the NA and the NS variables, which have increased by 0.7 and 2.4 million per decade, respectively (Table 6). Moreover, the Pettitt's test captured the timing of this abrupt change in the NA and the NS variables in the years 2006 and 2009 correspondingly.

4.4 Relationship between water transfers, climatic variables, tourism variables and indices

Initially, in order to reduce the effect of potentially confounding non-climatic factors on the climatic variables, indices and the water transfer dataset under study, a first-difference time series was constructed for the entire dataset, whilst the magnitude of the association between the variables was captured by means of the Pearson correlation analysis method (Huang et al. 2018). Moreover, the IBM SPSS Statistics v.25 (IBM Corp., Armonk, NY, USA, 2017) software was used so as to create the r matrix between the water transfer dataset and the climatic variables; it is portrayed in Fig. 8. The r coefficient ranges between -1 and 1, and the association between the two parameters can be categorized as weak, moderate and strong, when the absolute value of the coefficient ranges between (0—0.3), (0.31—0.60) and (< 0.60) respectively (Gerber and Finn 2005).

Fig. 8
figure 8

Pearson (r) correlation matrix value between the water transfer dataset and climatic variables for the period 1997—2008

Full size image

The correlation analysis indicates that there is a strong positive correlation, between the water transfer dataset and the diurnal temperature range at a significance level of 0.05, which suggests that as the diurnal temperature variation over the region increases, so does the transported water volume and related costs. On the contrary, there is a moderately negative relationship between the water transfer dataset and the RDIst, SPI, CDW, CWD, Rx5day, PRCPTOT, R95pTOT, SDII, R10, aridity and precipitation. The latter suggests that a) as the monthly minimum value of the daily minimum temperature reduces, the water transfer levels increase b) as 7 of the 11 precipitation related indices and annual precipitation decreased, the transported water volume and costs rose and c) as drought phenomena and aridity magnified over the region, more water and money were needed in order to transfer potable water by ships to the region. The negative impact of precipitation extreme indices on explanatory variables such as the rice yield (Huang et al. 2018; Wang et al. 2016) or wheat (Li et al. 2016; Zarei and Mahmoudi 2020) has also been documented in China.

Although the M–K trend test revealed significant positive trends for the NA and the NS variables, with the timing of the abrupt change for the study period 1991–2020 being 2006 and 2009, respectively (Table 6), the correlation analysis between water transfers and tourism variables for the period 1997—2008 (Fig. 8) indicated a weak negative relationship between the NA variable with the Volume (V) and Cost (V) of the transported water. In addition, a strong negative statistically significant relationship was found between NS with V and C. The latter suggests that—from 1998 to 2008, tourism declined in the Region whilst was an increasing quantity of water was transferred to the Region so as to match the water shortages At the same time, the irrigated land decreased progressively (Fig. 6), as the farmers were facing increasing difficulties while coping with the shrinking water resources due to the changeover of the Region of South Aegean's climate to warmer conditions with less precipitation. Since both tourism and irrigated land were contracted during the 1997–2008 period, the increase on water transfers could be primarily attributed to the changes on the Region’s climate patterns given that they are expressed by warmer conditions with reduced precipitation and the environment shift to semi-arid conditions.

Regarding the period 2009–2020, a clear strong surge in tourism was documented, due to the desalination plants’ operation, sufficient quantities of water were provided for drinking and domestic use in the Region. Thus, the transported water was reduced to 2020 up to some limited amounts of water to only three remote islands. Today, the desalination plants are operating at their limits, apparently due to changing climate conditions and increased tourism, which explains the urgent need to develop new solutions so as to cope with the water shortages in the Region, such as the dam construction, new desalination plants establishment or the assembly of an underwater pipeline in order to provide freshwater to the islands from the Greek mainland.

5 Conclusions

The IPCC’s Fifth Assessment Report (IPCC 2013) has highlighted the warming of the climate system, noting that any of the observed changes are unprecedented. Moreover, water resource availability in the Mediterranean has already been affected by environmental change (García-Ruiz et al. 2011). The reduced water resource availability over the Region of South Aegean in Greece led Greek governments to hire a tanker fleet so as to transport 17.8 Mm3 of drinking water to the region with ships from the mainland between 1997–2018 at a total cost of 124.2 million euros (SGAIP, 2019a).

The present study initially captured the climate extreme indices, aridity, and drought patterns over the region using the RClimDex software, aridity index, and the RDIst and SPI indices respectively, so as to investigate the changes in the latter variables. The aforementioned variables’ trend and abrupt change detection were tracked by means of the Mann–Kendall and Pettitt tests respectively and the research findings illustrate a tendency of the region's climate towards semi-arid conditions with an intensification of drought conditions, whilst warmer conditions with less precipitation have prevailed in the area, with the timing of the abrupt change being around 2006. The changes in climate patterns have a direct effect on drinking water transfers to the South Aegean Region by ships and leads to an overburden of the water resources regional capacities to cope with the growing water demand.

The volume and cost dataset of the water transferred to the region was connected to the above-mentioned variables through strong positive correlation trends regarding the diurnal temperature range and a moderately negative relationship regarding 7/11 of the precipitation related indices, annual precipitation, drought phenomena and aridity. The relevant research findings can be used by policy makers so as to design sustainable infrastructure that could withstand current and future anticipated warmer and drier conditions.