Unravelling Precipitation Trends in Greece since 1950s Using ERA5 Climate Reanalysis Data

Abstract

Precipitation is one of the most variable climatic parameters, as it is determined by many physical processes. The spatiotemporal characteristics of precipitation have been significantly affected by climate change during the past decades. Analysis of precipitation trends is challenging, especially in regions such as Greece, which is characterized by complex topography and includes several ungauged areas. With this study, we aim to shed new light on the climatic characteristics and inter-annual trends of precipitation over Greece. For this purpose, we used ERA5 monthly precipitation data from 1950 to 2020 to estimate annual Theil–Sen trends and Mann–Kendall significance over Greece and surrounding areas. Additionally, in order to analyze and model the nonlinear relationships of monthly precipitation time series, we used generalized additive models (GAMs). The results indicated significant declining inter-annual trends of areal precipitation over the study area. Declining trends were more pronounced in winter over western and eastern Greece, but trends in spring, summer and autumn were mostly not significant. GAMs showcased that the trends were generally characterized by nonlinearity and precipitation over the study area presented high inter-decadal variability. Combining the results, we concluded that precipitation did not linearly change during the past 7 decades, but it first increased from the 1950s to the late 1960s, consequently decreased until the early 1990s and, afterwards, presented an increase until 2020 with a smaller rate than the 1950–1960s.

1. Introduction

Countries across the Mediterranean basin are susceptible to climate change, as water availability is expected to decrease in the future with a simultaneous expansion of aridity [1,2,3]. The climate change influences the alternation of seasons by increasing the disparity between wet and dry ones [4,5,6], while it changes the spatiotemporal distribution of precipitation, thus strongly affecting ecosystems and human activities [2,4,7]. This phenomenon is expected to be strong in countries such as Greece, which is characterized by large landscape variability additionally affecting the local climatic characteristics [6,8]. Although the typical climate of Greece is generally characterized by dry-hot summers and mild-wet winters, there are many Greek regions that present intense spatial and seasonal variabilities of precipitation [9,10]. The complex orography of Greece forces precipitation to follow an irregular pattern with a strong gradient between the western part and the other regions [11]. The topographic characteristics of Greece with steep mountains and sea–land transitions play a key role in the climatic pluralism holding, for example, more precipitation at northern–western mountainous areas in comparison to drier southern–eastern lowlands [6,12]. Therefore, the investigation of inter-annual precipitation trends is a difficult task in areas such as Greece, which includes many mountains, sea areas, islands, and several ungauged regions.

Previous studies have examined the temporal and/or spatial distribution of precipitation in Greece, applying statistical methods and using data retrieved from meteorological stations (e.g., [8,13,14,15,16,17]). It is noteworthy that previous studies found statistically significant declining trends for precipitation in Greek areas during the previous decades, mainly in the second half of the 20th century (e.g., [18,19,20]). For example, Feidas et al. [19] found decreasing trends in precipitation from 1955 to 2001 by statistically analyzing data from 24 Greek meteorological stations. In addition, Pnevmatikos and Katsoulis [21] implemented a statistical analysis using data from 36 Greek meteorological stations that resulted in clear evidence that most Greek territories have undergone shifts in their rainfall regimes. Apart from meteorological stations, another valuable source that can be used, especially at regions with poor observational data, is the spatially distributed precipitation data derived after reanalysis and reprocessing [22]. Furthermore, Philandras et al. [20] reported decrease of precipitation and rainy days in Greece from 1951 to 2009 by analyzing precipitation data from stations and spatially distributed datasets. In the study of Hatzianastassiou et al. [23], satellite-based monthly precipitation data obtained from the Global Precipitation Climatology Project (GPCP) were evaluated against measurements from 36 stations and then were used to examine the spatial and temporal variations of precipitation in Greece and surrounding regions from 1979 to 2004. This study presented decreasing and increasing trends of precipitation while it highlighted the value of using spatially distributed precipitation data in temporal analyses [23]. Therefore, the lack of dense, continuous, and long-term observational networks usually poses the necessity to employ long-term spatially distributed datasets for the investigation of inter-annual precipitation trends over continuous grids, including even ungauged areas. It is also important to consider the past 2 decades (i.e., 2000s and 2010s) in trend analyses due to their significance regarding climate change. However, the period since 2000 was not completely included in the previous relevant studies as also noted by Markonis et al. [8] and Katsafados et al. [24].

The main scope of this study was the identification of statistically significant (or simply significant hereafter) inter-annual trends of precipitation over Greece and neighboring areas for a 71-year period from 1950 to 2020, including even ungauged land regions as well as sea areas where precipitation has remarkable effects (e.g., on sea state and sea temperature [25,26]). To this aim, we used ERA5 reanalysis precipitation data, which are consistent in space and time covering 7 decades and which have recently been used by other studies, such as those by Mastrantonas et al. [27] and De Luca et al. [28], to identify spatiotemporal characteristics of precipitation in the Mediterranean. The innovation of our study lies in the processing of the abovementioned data with robust statistical techniques to explore statistically significant trends of precipitation over the study area and, subsequently, to quantify nonlinear variabilities that have occurred within the past approximately 7 decades. In this way, our study can reveal areas that have been mostly affected by climate change and, thus may be more vulnerable to climate change in the future.

2. Materials and Methods

In order to investigate statistically significant precipitation trends over Greece and neighboring areas we estimated Theil–Sen linear trends and Mann–Kendall significance. Then we used generalized additive models (GAMs) to analyze and model nonlinear variabilities of precipitation over the study area. Precipitation data, study area (Figure 1) and statistical methods used in this study are described in the following subsections. The main data processing and calculation stages are briefly illustrated in a flowchart (Figure 2) at the end of the section.

2.1. Climatological Precipitation Data and Study Area

The historical precipitation data used in this study were obtained from the ERA5 reanalysis dataset [29]. ERA5 climatological data were selected, as they are considered the most modern and reliable datasets that can be used in a wide range of long-term analyses and multidisciplinary studies (e.g., [27,30,31,32,33,34]). ERA5 data were produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) and they are available online by the Copernicus Climate Change Service (C3S, [35]). ERA5 data cover both land and sea areas at a horizontal resolution of 0.25° × 0.25°. ERA5 multiyear reanalysis facilitates inter-annual and climatic analyses as they cover a long period from 1950 to nearly present (2–3 months before today). Until recently, ERA5 has been covering a period from 1979 to nearly the present, but currently they also cover the period 1950–1978, after their recent back extension [36]. It is important to note that this back extension of ERA5 for the period 1950–1978 is considered as preliminary due to an issue related to the representation of tropical cyclones. This issue concerns tropical areas away from our study area, and therefore preliminary 1950–1978 ERA5 data can be normally used in this study as a back extension of the dataset for the period from 1979 to nearly the present. An advantage of ERA5 data is that they provide high-quality climatological information on hourly or monthly averaged basis. It is easy to retrieve ERA5 data in friendly GRIB and NetCDF formats, and they are also available on custom regular latitude–longitude grids anywhere on Earth, thus facilitating various analyses and applications in any country, including Greece [12,37].

First, we selected the study area covering the entire territory of Greece and neighboring areas enclosed in a region between 18.75° W–29.75° E and 34.5° N–42° N, as demonstrated in Figure 1. Figure 1 also illustrates the topography of the study area using Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 (ASTGTM) data at a spatial resolution of approximately 30 m [38]. The study area consists of 1395 ERA5 grid points. Consequently, we obtained ERA5 monthly precipitation data for each month in the 71-year period of the study, spanning from January 1950 to December 2020. The data were retrieved from the C3S service on a regular latitude–longitude grid covering the study area. They are known as ERA5 monthly averaged precipitation data and they are expressed in m day⁻¹, representing the average precipitation amount per day for each month. Then, the monthly averaged precipitation values were used to estimate the respective monthly accumulated precipitation (or simply monthly precipitation hereafter) values expressed in mm month⁻¹ or simply in mm, by multiplying the first ones by 1000 (conversion from m to mm) and, consequently, by the number of days of each month. The monthly precipitation values for each grid point were used as a basis for calculating annual values that were used to estimate Theil–Sen inter-annual trends. The monthly precipitation time series for each point were also employed for the estimations performed by GAMs.

Figure 1. Study area covering the entire territory of Greece and neighboring areas. The map also demonstrates topography (m) based on Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 (ASTGTM) data at a spatial resolution of approximately 30 m [38].

Figure 1. Study area covering the entire territory of Greece and neighboring areas. The map also demonstrates topography (m) based on Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 (ASTGTM) data at a spatial resolution of approximately 30 m [38].

Figure 2. Flowchart demonstrating the main data processing and calculation stages.

Figure 2. Flowchart demonstrating the main data processing and calculation stages.

2.2. Data Formatting, Spatial Analysis and Estimation of Inter-Annual Trends

We edited the 71-year ERA5 dataset to come in a format serving the purposes of this study. First, we estimated annual accumulated precipitation (or simply annual precipitation hereafter) values for the 1395 ERA5 grid points, summing up the 12 monthly precipitation values for each year. Consequently, we estimated seasonal accumulated precipitation (or simply seasonal precipitation) values for each point considering the four seasons. We estimated seasonal precipitation values for winter (i.e., December, January, and February; DJF), spring (i.e., March, April, and May; MAM), summer (i.e., June, July, and August; JJA) and autumn (i.e., September, October and November; SON). Having calculated annual and seasonal precipitation fields for each year, then, we estimated the 71-year average of annual and seasonal precipitation for each grid point. We also created maps with these 71-year averaged values to present the spatial distribution of precipitation across the study area, on annual and seasonal basis. To better demonstrate the values for each reanalysis point, average annual and seasonal precipitation were depicted on raster-based maps avoiding spatial interpolation.

Afterwards, we estimated Theil–Sen linear trends [39,40] for each ERA5 grid point in the study area, considering the 71-year period from 1950 to 2020. The Theil–Sen trend estimation method is one of the most popular nonparametric techniques for estimating a linear trend. It is insensitive to outliers, and it can be better than simple linear regression for skewed and heteroskedastic datasets, competing also well against nonrobust least squares even for normally distributed data in terms of statistical power. To estimate Theil–Sen trends, we used annual precipitation values for each point considering the 71 years. We used annual precipitation values instead of monthly ones to make deseasonalized time series. In this way, we avoided probable errors in our estimations related with seasonal variability of precipitation within each year. We also applied Mann–Kendall (M–K) nonparametric tests for monotonic trends [41,42,43] to check the significance of the 71-year trends for each point. The M–K test is distribution-free and is used to detect the presence of both linear and nonlinear trends in time series assessing if they are steadily increasing/decreasing or are unchanging. Consequently, we estimated Theil–Sen linear trends and M–K significance considering time series of annual precipitation spatially averaged over the study area (annual areal precipitation hereafter). The annual areal precipitation is translated to the average water amount expressed in mm that drops across the entire study area each year and was calculated as the arithmetic mean of the 1395 annual precipitation amounts from the respective grid points covering the study area. The estimation of Theil–Sen trends and M–K test was based on the respective functions of NCL (National Center for Atmospheric Research-NCAR Command Language [44,45]).

2.3. Generalized Additive Models for Linear and Nonlinear Trends

After the estimation of the inter-annual Theil–Sen trends and the Mann–Kendall significance, we used GAMs to analyze and model the monthly precipitation time series data aiming beyond a fitting based on linear trends. We chose GAMs because they are capable of modelling nonlinear relationships between the response variable and the predictors [46] and therefore they have many applications in environmental sciences [47,48,49]. They are also widely used in time-series analysis since they can capture nonlinear trends in time series and deal with the irregular spacing of samples in time [37]. The components of a time series are represented as smooth functions, which are nonlinear representations of the covariates, composed by the sum of K simpler basis functions [48]. A general form of generalized additive model is:

$$g\left(Y\right)=\beta +f_1\left(x_1\right)+f_2\left(x_2\right)+\dots +f_n\left(x_n\right)$$

(1)

where Y is the expected response value, β is the model intercept and f₁, f₂, f_n are smooth functions of the predictors x₁, x₂, x_n [50]. In our models the response variable Y is the monthly precipitation, and the predictors are the trend (named “Time”) and the intra-annual variation (named “Month”). The smoothness selection for each term was based on a restricted maximum likelihood approach (REML). Thus, we fitted 1395 models, one for each grid point, using for “Time” smooth term a cubic regression smoothing spline, while for “Month” we used the cyclic cubic spline with k = 12. We also fitted a model using the monthly areal precipitation values which were calculated as the arithmetic means of the monthly precipitation values from the 1395 grid points of the study area for each month to assess the inter-annual variability for the whole region of study. All models were fitted with the gam function of the “mgcv” package in R environment [51].

3. Results

3.1. Distribution of Average Annual and Seasonal Precipitation

The spatial distribution of 71-year averaged annual and seasonal precipitation over the study area is presented in Figure 3 and Figure 4, respectively. The highest average annual precipitation values are depicted at western Greece with peaks over 1350 mm detected at Epirus and the Ionian Sea (Figure 3). There are some additional areas that are characterized by high precipitation values. These include the eastern Aegean Sea, western Peloponnese, northern Thessaly, eastern Macedonia, western Thrace, and western Crete. On the other hand, the lowest average annual precipitation values are demonstrated at eastern–southern Crete, and at the southeastern–southern–central Aegean Sea as well as at Athens and eastern Peloponnese (Figure 3). Precipitation values about 350–450 mm are illustrated at these areas, while values even lower than 350 mm can be detected at the southeastern Aegean Sea and southeast of Athens. The analysis of the 71-year ERA5 dataset conducted in this study generally corroborates previous studies. In particular, the spatial distribution of average annual precipitation over land areas of Greece agrees with the respective one presented in a map included in the classic book of Mariolopoulos and Karapiperis [52] for precipitation in Greece. It is also in agreement with the respective distribution demonstrated in the climatic atlas of Greece of the Hellenic National Meteorological Service [53,54]. The spatial distribution of annual precipitation is mostly determined by the mountains and the passage of barometric depressions accompanied by fronts, mainly occurring during the cold period. In general, the barometric systems move in a west-to-east direction bringing larger water amounts over the windward slopes of mountains in comparison to the leeward ones (Figure 1 and Figure 3), thus sufficiently explaining local variations of precipitation across Greece [9]. It is noteworthy that Pindus Mountain range (Figure 1) is a typical example as it plays a key role on the formation of a wetter climate over the western continental Greece in comparison with the drier climatic conditions that prevail over the eastern continental Greece.

The seasonal variability of precipitation in Greece follows a typical pattern with alternations of wet, transitional, and dry periods. The peak of precipitation occurs in winter (DJF) (Figure 4a) illustrating similar pattern as the annual distribution (Figure 3). Smaller precipitation amounts are generally shown for spring (MAM) with peaks located over mountainous areas (Figure 4b). Summer (JJA) is too dry especially for the islands and over the sea (Figure 4c). Especially for the Aegean Sea, Etesian winds which are dry north winds prevail for long periods in the summer [9]. However, some precipitation amounts are presented over the mainland areas in summer (Figure 4c), which are associated with thermal thunderstorms, usually occurring from late spring to early autumn. Significant amounts of precipitation are observed in autumn (SON), mainly over the windward slopes of mountains in western Greece (Figure 4d). Overall, the passage of barometric systems and fronts explains most of precipitation amounts in winter, autumn (especially late) and spring (especially early), while rainless Etesian winds over the Aegean Sea and local convective storms over continental areas strongly determine precipitation distribution in summer.

3.2. Inter-Annual Precipitation Trends

The results of trend analysis of annual precipitation conducted for each ERA5 grid point are presented in Figure 5. The map of Figure 5 spatially combines the trends estimated for the 1395 grid points to highlight areas with increasing, unchanging, and declining trends, also adding information for the Mann–Kendall significance. Figure 5 demonstrates a general declining trend of annual precipitation across the study area during the 71-year period from 1950 to 2020. The trends are significant, especially in western Greece and the eastern Aegean Sea, with some additional limited areas at the northern Aegean Sea and the sea area west of Crete. It is interesting to note that there are some areas such as Athens, northeastern Peloponnese, central–eastern Crete, and the central–western Aegean Sea where slight increasing trends are pronounced (Figure 5). Even though these increasing trends are not significant, they at least do not imply declining trends. This is an optimistic finding for those regions that are characterized by relatively low precipitation (Figure 3) and might be susceptible to desertification in the future.

It is obvious that annual precipitation trends are determined by winter (DJF) that demonstrates a similar pattern at the eastern Aegean Sea and western Greece (Figure 6a). Additionally, more grid points with significant trends at the Ionian Sea are presented for winter (Figure 6a). In spring (MAM) and summer (JJA), there are few points characterized by significant trends either unchanging or declining (Figure 6b,c). In autumn (SON), there are more points with significant trends compared with spring and summer but fewer than winter and limited at the southeastern Aegean Sea (Figure 6d). Therefore, winter seems to mostly determine the inter-annual trends of precipitation. The analysis resulted in declining trends of precipitation during the past 7 decades that are more pronounced over the wettest areas of Greece (i.e., western Greece and the eastern Aegean Sea as illustrated in Figure 3). This finding could be considered encouraging because these areas are not as vulnerable to drought as some drier Greek areas such as the region around Athens, the eastern Peloponnese, and eastern Crete. Nevertheless, the reduction of precipitation even at the wettest areas might have significant implications for ecosystems, agriculture, the economy, and society.

It is also important to investigate how the areal precipitation over the entire study area changes over time. In this context, Figure 7a demonstrates time series of annual areal precipitation over the study area. The most striking result to emerge from the 71-year time series is that the areal precipitation is not characterized by a linear inter-annual trend. In more detail, we can identify an increasing trend from 1950 to the late 1960s, followed by a declining trend until the late 1980s to the early 1990s and a relatively small increasing trend afterwards that gradually turns to an unchanging trend during the past decade. Overall, a 71-year declining Theil–Sen trend of −1.02 mm/yr was estimated for the 71-year period with Mann–Kendall significance equal to 93%. The overall inter-annual trend is not significant at the 95% level. This may be attributed to the increased nonlinearity of time series. Although the estimated trend did not reach 95% significance level, it is noteworthy that our study provides considerable insight into the impact of climate change on precipitation, which is not straightforward. The inter-annual variability of areal precipitation is also captured by the time series of anomalies. Precipitation anomalies were constructed by subtracting the 71-year average precipitation from annual precipitation values for each year (Figure 7b). In general, the 1950s and especially the 1960s and 1970s were characterized by high precipitation amounts across the study area; however, the 1980s and 1990s had lower ones. During the 2 most recent decades, i.e., the 2000s and 2010s, alternations between dry and wet years occurred. Weak increasing and unchanging phases can be also detected during the past 20 years.

3.3. GAM Results

Consequently, in order to identify significant linear and nonlinear inter-annual trends of precipitation, we fitted a total of 1395 models, one for each grid point. More than half of the grid points (732 from 1395) showed significant trends (p ≤ 0.05). R² values for these models ranged from 0.06 to 0.66 with an average R² of 0.46. Most grid points showed variable trends with precipitation following an increasing trend from the 1950s to the middle 1960s, to a decline until the middle 1980s and then an increase again until today (e.g., Figure 8). For some other areas (e.g., Chania and Rhodes, Figure 8) a declining nonlinear trend from the 1950s to the early–middle 1980s is apparent before a slight increase until 2020. Non-significant trends were mostly linear (Figure 11), as it was the trend for Athens that showed an increasing linear curve (Figure 8).

Regarding intra-annual variations, the majority of selected locations are characterized by a typical seasonal distribution, following alternations illustrated in Figure 9. Most of the areas have a dry summer (adding also late spring and early autumn), wet winter (adding also late autumn and early spring) and two transitional seasons, including sub-periods of autumn and spring. This distribution strongly depends on cyclonic and frontal activity at this area [9]. It is noteworthy that Drama, Thessaloniki and Alexandroupoli, which are cities in northern Greece, additionally present high precipitation amounts in late spring and early summer. These are explained by the frequent occurrence of convective storms there during this period.

The “whole region” model had an R² of 0.65 and a significant trend (p < 0.001). The trend for the whole region (Figure 10) shows a highly variable trend with an increasing phase from the 1950s to the late 1960s, a decreasing period until the early 1990s and then again an increasing trend until 2020. However, the increasing precipitation during the past 3 decades has not managed to reach the levels of the period from the late 1950s to the early 1970s. This remains as a noteworthy issue for the coming years. Figure 10 also demonstrates the intra-annual variations of precipitation as represented by the “whole region” model, confirming the seasonal distribution shown in Figure 4a–d.

We also quantified the degree of nonlinear trends by using the effective degrees of freedom (edf), which reflects the nonlinearity of the fitted GAM curve [51]. Edf values close to 1 indicate a linear trend, whereas higher values show more complex and nonlinear trends. Figure 11 depicts the spatial distribution of significant (p ≤ 0.05) and nonlinear trends showing that areas of Western Greece, and the Ionian and Eastern Aegean are more likely to show significant and nonlinear trends (edf > 2.5). These findings agree with the results from the Theil–Sen trend analysis (Figure 5), although the analysis with GAMs revealed a notably higher number of cells with significant trends. This is attributed to the advantage of GAMs to better fit the data and deal with nonlinear trends.

4. Discussion

4.1. Variability and Trends of Precipitation in the Study Area

Our results corroborate previous studies that reported inter-annual trends of precipitation over the study area (e.g., [8,18,19,20,23]). For example, our study agrees with some findings from Feidas et al. [19], who found declining precipitation trends analyzing data from 24 Greek meteorological stations for the period from 1955 to 2001. This previous study also reported significant inter-annual trends of precipitation only for winter, which seems to mostly determine trends of annual precipitation in Greece over time. It is noteworthy that our study also provides considerable insight into the inter-decadal variabilities of precipitation, again substantiating the study of Feidas et al. [19]. They also stated that the decrease of precipitation during 1955–2001 is mostly attributed to the relatively drier period from the middle 1980s until 2001, which generally agrees with our findings. Our results are also in line with the study conducted by Philandras et al. [20]. They estimated declining precipitation trends across our study area for the period from 1951 to 2009, employing data from meteorological stations as well as lower-resolution gridded data. They reported significant declining trends for the rainy period (i.e., from October to March), that includes winter, more pronounced in western Greece, which is consistent with our findings. Our study is also in agreement with the study of Markonis et al. [8] highlighting that most of studied Greek regions (records from 136 stations) showed a decline in precipitation since 1950 and an increase since 1980, while they remained stable from the late 1990s to 2012. In addition, Hatzianastassiou et al. [23] reported an inter-annual variability of precipitation from 1979 to 2004, which largely fitted the one we presented in this study for the abovementioned period. They analyzed satellite-based data and found a decreasing trend from 1979 through the 1980s and an increase during the 1990s up to the early 2000s, which was followed by a decline up to 2004.

One of the innovative characteristics of our study is the use of ERA5 reanalysis precipitation data, which are consistent in space and time. Given the lack of a dense, continuous, and long-term network of meteorological stations in Greece, ERA5 reanalysis data can be a valuable basis to explore inter-annual precipitation variabilities even in ungauged areas. An additional advantage of our work is the long-term study period spanning 1950 to 2020. It is noteworthy that the study period includes the past 2 decades from 2000 to 2020. This is critical due to the growing impact of climate change on the spatiotemporal characteristics of precipitation during the past decades. It is worth mentioning that our study showed that precipitation was not characterized by a continuous declining trend from the 1950s up to 2020. It presented a small increasing trend during the past 2 decades, thus adding unexpected information in the research of inter-annual precipitation trends over Greece. It is also likely that a future continuation of the current increasing phase could potentially lead to precipitation levels of 1950 to the 1970s. In addition, we explored inter-annual trends of precipitation over Greece and neighboring areas, including sea areas, which are often excluded from precipitation trend analyses. Finally, we demonstrated the advantageous use of GAMs in analyzing inter-annual trends of precipitation data. In particular, we showed that GAMs can be employed to reveal hidden patterns in precipitation time series by using flexible and easy-to-apply functions with relatively high interpretability.

4.2. Possible Reasons of the Nonlinear Precipitation Variabilty over Time

The nonlinear variability of precipitation in Greece and surrounding areas is mostly determined by the general atmospheric circulation and large-scale climatic patterns that are affected by climate change over time. More specifically, the North Atlantic Oscillation (NAO), which is one of the main large-scale climatic patterns in the Northern Hemisphere, plays an important role in the long-term periodic variabilities of precipitation in Greece as highlighted by previous studies (e.g., [18,19,20,23,55]). NAO implies alternations of atmospheric pressure between the Azores high and the Icelandic low having either positive or negative phase (or index), while it is characterized by quasi-biennial to quasi-decadal and/or decadal time scales [23]. NAO impacts weather and climate conditions in various areas such as the Atlantic Ocean and Europe, causing strong precipitation variations even in the Mediterranean basin, as occurred, for example, during the dry period of the 1980s [56]. The declining precipitation trend in the 1980s, which was also evidenced by our analysis, was related to an increase in the NAO index [23]. On the other hand, the increasing precipitation trend in the 1990s, also identified by our study, was related to a decrease in the NAO index [23]. In general, precipitation trends in Greece during various periods seem to be anticorrelated with the NAO index as discussed by Feidas et al. [19] and Philandras et al. [20]. Except for NAO, studies have been conducted to investigate the effects of additional large-scale climatic patterns in parts of our study area (e.g., Greek Ionian Sea region [55]). These climatic patterns include the El Niño Southern Oscillation or ENSO, the East Atlantic or EA, the Scandinavian or SCAND, the East Atlantic–West Russian or EA–WR, and others (Kalimeris et al. [55] and references therein).

4.3. Implications of Precipitation Variability on Inland Freshwater Resources

The present findings have important implications for current and future freshwater availability, with particular focus set on the inland freshwater ecosystems and their delivered services. Several studies have reported declining river flows and changes in the water quantity during the past decades [57,58], but it remains unclear how reducing precipitation is linked with hydrological alterations. This is more prominent in the Mediterranean since river damming and water abstractions are major environmental problems that are responsible for changing flow patterns and reduced water quantity [59,60,61]. For instance, Mentzafou et al. [57] showed significant relationships between precipitation and discharge with similar fluctuation patterns for major Greek rivers but concluded that a more in-depth analysis that considers the water uses will reveal the actual causes of flow regime changes. Nevertheless, since our results showed small but significant declines of total precipitation for a considerably large share of the studied area, it becomes obvious that a decrease in water availability will impact both the anthropogenic water uses and the ecological integrity.

4.4. Necessity for Inter-Annual Analysis of Short-Term Extreme Precipitation

This study is a first step towards enhancing our understanding and interpreting of inter-annual precipitation trends over Greece focusing on annual and seasonal accumulated amounts. One of the major issues in the context of climate change in the Mediterranean area seems to be the change in frequency and intensity of heavy precipitation events [62]. Therefore, long-term analyses of short-scale precipitation are important to better explain how the large-scale climatic variabilities impact the spatiotemporal characteristics and the magnitude of precipitation extremes [15,17,27,63,64,65,66]. Furthermore, in order to quantify the potential impact of precipitation variability on the ecosystems, a more detailed and elaborate analytical approach of rainfall events is required that is based on the concept of individual rainfall episodes [67]. Particularly in the Mediterranean environments where precipitation shows high and irregular variability, the use of annual or monthly averages of precipitation might not be representative of the actual hydro-geomorphological processes [67,68,69]. Instead, precipitation intensity for short time scales (e.g., daily and hourly) might be more appropriate for assessing the impact on the water resources and their biotic components. Precipitation data at higher spatial resolution might also be more appropriate, especially in a mountainous country such as Greece.

In this context, the hourly ERA5-Land [70] reanalysis precipitation dataset masked over the land at a native horizontal resolution of 9 km × 9 km might be very useful. Just as the ERA5 data, ERA5-Land data were recently produced by the ECMWF and are freely available by the Copernicus-C3S platform [35]. A main advantage of ERA5-Land, compared to ERA5, is the higher horizontal resolution. Thus, considering the better representation of orography, ERA5-Land could facilitate a more accurate analysis of inter-annual trends of precipitation in a future study, calculating Theil–Sen trends and Mann–Kendall significance as well as executing GAMs. Moreover, an improved accuracy of hourly or daily ERA5-Land precipitation data can help to investigate inter-annual changes in the intensity and frequency of extreme precipitation events that often cause flash floods, having critical implications for freshwater ecosystems, the economy, and society [71,72,73,74,75,76,77,78]. It is also interesting that such future work using ERA5-Land data could facilitate applications of the Indicators of Hydrologic Alteration (IHA) methodology in river ecosystems to explore their ecological integrity and functioning [60].

5. Conclusions

This work presents a detailed overview of long-term trends of precipitation for the entire territory of Greece and neighboring areas, spanning 7 decades from 1950 to 2020. To this aim, we used ERA5 monthly averaged precipitation reanalysis data from January 1950 to December 2020. At first, we estimated Theil–Sen linear trends and Mann–Kendall significance considering annual precipitation time series for each ERA5 grid point and the areal precipitation over the study area. Consequently, we used GAMs to explore and analyze nonlinear inter-annual variabilities of precipitation.

Our results show a spatially variable pattern of changes over time. Significant Theil–Sen declining linear trends were limited mostly at the regions of western Greece and the eastern Aegean Sea. Declining linear trends were also estimated for winter more or less at the same areas, however, trends for spring, summer and autumn were not significant in the majority of areas. Overall, the results indicate large inter-decadal fluctuations of areal precipitation over the study area. The analysis of areal precipitation resulted to a declining linear trend, but significant at a 93% significance level. GAMs highlight that the trends are generally nonlinear. We identified significant nonlinear changes for more than half of the studied area. Both nonparametric and modelling analysis showed a complex pattern of precipitation trends for most of Greece. To summarize, overall precipitation over Greece was not characterized by a linear trend during the past 7 decades. Overall precipitation first increased until the late 1960s, consequently decreased until the early 1990s and then increased until 2020 with a relatively small rate.

We believe that such information can be exploited by other scientists for climate and interdisciplinary studies as well as by decision-makers and civil protection who will have to take into consideration precipitation variabilities when they make management plans. For example, our study sets a basis for long-term hydrological studies such as that of Mentzafou et al. [78]. Nevertheless, our study has identified some research issues in need of further examination. One of the most critical questions is whether there are inter-decadal variabilities of short-term extreme precipitation amounts, except for total annual and seasonal ones. In this context, the findings of this study would support in the future a trend analysis of the magnitude and frequency of extreme precipitation events during the past 7 decades. Such a work might unveil significant impacts of climate change on the characteristics of extreme precipitation having crucial implications on society and ecosystems.