Modeling Return Levels of Non-Stationary Rainfall Extremes Due to Climate Change

Abstract

Global warming increases evaporation and atmospheric water vapor, leading to more extreme events in both spatial and temporal domains. This study conducts a non-stationary extreme value analysis of the annual daily maximum at 36 meteorological stations over Iran from 1960 to 2021. We applied stationary and non-stationary Generalized Extreme Value (GEV) models within a Bayesian framework to estimate return levels for rainfall extremes, along with 90% confidence intervals. Our findings indicate that non-stationary models are not prominently evident based on AIC at most stations; however, non-stationary Generalized Extreme Value (GEV) models outperform stationary models based on RMSE and NSE evaluation criteria that sufficiently capture variations in extremes. Furthermore, most observed changes in extreme events exhibit a non-stationary pattern. Non-stationary analysis indicates that the frequency and severity of rainfall extremes have shown both increasing and decreasing trends, characterized by inconsistent spatial patterns.

1. Introduction

Climate change is one of the most significant factors in altering rainfall patterns. Increased greenhouse gas concentrations due to climate change can intensify hydrological cycles, leading to changes in water resources and seasonal variations in rainfall worldwide. These climatic disturbances have resulted in more extreme rainfall in the latter half of the 21st century. However, it is important to note that global warming will not uniformly increase rainfall everywhere. Changes in pressure patterns mean some regions may also experience decreased rainfall.

Extreme rainfall poses serious risks, particularly in the context of natural disasters like floods. Understanding how this phenomenon behaves across different climates and over time is essential for effective planning and infrastructure design, including drainage systems, dams, and sewage networks. Accurate knowledge of extreme rainfall patterns is critical for these efforts. Climate change is likely to change extreme events, especially extreme rainfall. Extreme rainfall causes floods, leading to severe infrastructure damage and loss of lives. Therefore, it is vital to know and model the spatial and temporal variability of rainfall extremes to control floods and prevent damage [1].

Understanding extreme rainfall patterns is valuable for effective water resource management and the design of flood-resistant infrastructure [2,3]. Understanding the relationship between patterns and the occurrence of droughts and floods is crucial, particularly in the context of climate change. It is important to note that climate change is expected to intensify the frequency and magnitude of extreme events in numerous regions globally [4]. Global warming and rising temperatures enhance evaporation, leading to increased atmospheric water vapor. Consequently, this alteration can intensify extreme rainfall events, potentially resulting in more occurrences in tropical and high-latitude regions while decreasing in hot and arid areas adjacent to the tropics. The IPCC’s sixth report predicts that maximum daily rainfall intensity may rise by approximately 7% for every 1 degree Celsius [5,6].

In recent years, the occurrence of non-stationary extreme rainfall has presented challenges due to the unreliability of assumptions in analyzing the frequency and intensity of extreme rainfall. Changes in climate as a result of natural variability or human intervention make the stationary assumption for the analysis of the frequency and intensity of extremes unreliable [7,8,9,10].

The current design of civil infrastructure systems is based on the assumption that rainfall characteristics will remain relatively constant over time. However, with climate change, the intensity and frequency of extreme rainfall events will likely become non-stationary. As a result, the intensity–duration–frequency (IDF) curves must be based on a non-stationary, time-varying approach [2,9,11]. Katz proposed the approach of modeling the temporal variability of GEV distribution parameters as a function of time [12]. Statistical models are used to derive the variability and frequency of extremes based on historical data [12]. Statistical models for study extremes are categorized as stationary and non-stationary [13,14]. In the stationary model, the statistical parameters of the probability distribution function do not change over time or are related to another covariate, while, in the non-stationary model, the parameters of the probability distribution function are not constant but can be changeable regard to time or a covariate [15,16]. This research is the first comprehensive investigation into the non-stationary trends in rainfall extremes across Iran. We analyzed the changes in the magnitude and frequency of rainfall extremes over time using time as a covariate. A Generalized Extreme Value (GEV) distribution was employed to model rainfall extreme events at each station in each year. Our objectives include (1) analyzing the trend of rainfall extremes over time, (2) understanding the non-stationary behavior of extreme rainfall events, (3) estimating return levels under stationary and non-stationary conditions, (4) determining the frequency and return period changes of rainfall extremes events over time, and (5) analyzing the non-stationary predictive probability density distributions to assess the tail behavior of rainfall extremes. The result of this study will be crucial for improved flood risk assessment, enhanced water resource management, and understanding of the behavior of rainfall extremes in a changing climate. This research enhances flood risk assessments by accurately characterizing changes in rainfall extreme trends. This information is vital for developing effective flood mitigation strategies and designing infrastructure that is resilient to more frequent and intense rainfall. Additionally, understanding the non-stationary behavior of rainfall extremes is crucial for sustainable water resource management. It enables precise predictions of water availability and fundamentals for planning water allocation, irrigation, and drought mitigation.

2. Materials and Methods

2.1. Data

We have utilized daily rainfall data collected from 36 weather stations across Iran, from 1960 to 2021. The geographical distribution of the stations is given in Figure 1. Two primary criteria guided the selection of these stations:

• Representation of diverse climatic regions: The chosen stations encompass a variety of climatic zones throughout Iran, ensuring a comprehensive overview of the country’s meteorological patterns.
• Data availability: Preference was given to stations with the longest continuous records, allowing us to analyze trends and variations in rainfall over an extended period.

We used MATLAB R2015a (version 8.5) (MathWorks, Natick, MA, USA) for data processing and modeling and ArcGIS 10.8.1. (Esri, Redlands, CA, USA) for spatial visualization.

Figure 1. Geographical distribution of the 36 meteorological stations in Iran.

Figure 1. Geographical distribution of the 36 meteorological stations in Iran.

2.2. Trend Test

Statistical tests frequently assume that the data conform to a normal distribution. However, in cases where the data deviate from normality, non-parametric tests are generally more appropriate and effective. This study used the non-parametric Mann–Kendall test to analyze trends in extreme rainfall. The Mann–Kendall test is a non-parametric test for identifying trends in time series data [17,18]. The null hypothesis in the M-K test represents no monotonic trend in the time series, while the alternative hypothesis indicates that the trend exists. This trend can be positive (increasing), negative (decreasing), or non-null. The Kendall Tau, or Kendall rank correlation coefficient, estimates the monotony of the slope. Kendall’s Tau varies between −1 and 1. When the trend is increasing, the Tau value is positive. On the contrary, negative values in the Tau test state a decreasing trend in the series.

2.3. Stationary and Non-Stationary Analysis of Extreme Values

Extreme value theory offers a robust framework for analyzing climate extremes and estimating their return levels [19]. Generalized Extreme Value (GEV) distribution is flexible for modeling different behaviors of extremes with three distribution parameters: θ = {$\xi$;$\mu$;$\sigma$}. The GEV distribution function is used to model time series of block maxima. The cumulative distribution function (CDF) of the GEV is given as

$$\begin{gathered} f\left(x∣\mu ;\sigma ;\xi \right)={\displaystyle \frac{1}{\sigma }}{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]}^{-\frac{1}{\xi }-1}\exp \left\{-{\left[1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]}^{-\frac{1}{\xi }}\right\};\xi \ne 0 \\ 1+\xi \left({\displaystyle \frac{x-\mu }{\sigma }}\right)>0 \end{gathered}$$

(1)

In this distribution, the random parameters of the random variable $f\left(x\right)$ are controlled by three parameters: location ($\mu$), scale ($\sigma$), and shape ($\xi$).

By letting the parameters of the cumulative distribution function of the stationary GEV be a function of a given covariate,$y$, the non-stationary GEV model can be described as

$$f\left(x∣y\right)=\exp \left\{-{\left[1+\xi \left(y\right)\left({\displaystyle \frac{x-\mu \left(y\right)}{\sigma \left(y\right)}}\right)\right]}^{-\frac{1}{\xi }}\left(y\right)\right\}$$

(2)

In the time-varying GEV distribution, the covariate explains changes in the distribution’s location, scale, and shape parameters.

In this study, location and scale parameters vary as a function of the time covariate, while the shape parameter remains constant for the non-stationary GEV model. Two models were used to analyze the time-varying patterns of extreme rainfall.

In stationary model M0, three GEV parameters are kept constant over time:

$${\mu }_t=\mu ; {\sigma }_t=\sigma ;{\sigma }_t= \xi$$

In non-stationary model M1, location and scale parameters vary as a function of the time covariate $\mu ={\mu }_0+{\mu }_{1 }{\times y}_{t }+\sigma ={\sigma }_0+{\sigma }_1 \times y_{t }+{\xi }_t= {\xi }_0.$

2.4. Parameter Estimate Using Bayesian Inference

Non-stationarity in hydrology is highly unpredictable, making deterministic methods unsuitable for calculation. Instead, Bayesian Inference provides an effective framework for modeling predictive uncertainty. In this study, the parameters of stationary and non-stationary GEV models were estimated using the Bayesian approach. We estimated the parameters of the posterior probability distribution by incorporating the prior probability. This process combined the prior distribution with the likelihood function based on observed data. The posterior distribution was estimated using the MCMC approach, which incorporated both the prior distribution and the likelihood function. We utilized the MLE method to estimate the GEV parameters. In the Markov chain calculation, five chains were used, each with 10,000 iterations, resulting in a total of 9000 samples after discarding the burn-in samples.

2.5. Model Evaluation

To evaluate the goodness of fit of the Generalized Extreme Value (GEV) model to the extreme rainfall time series, we used graphical diagnostic methods such as Probability–Probability (P-P) and Quantile–Quantile (Q-Q) plots, as well as the two-sample Kolmogorov–Smirnov (K-S) test. Additionally, to select the best-fitting GEV model, we utilized various metrics, including the Akaike Information Criterion (AIC) [20], root mean square error (RMSE), and Nash–Sutcliffe efficiency (NSE) [21] coefficient. Therefore, an optimal model is determined not solely based on goodness of fit but also considering model complexity and minimum residuals. This methodology aims to identify the model that best characterizes and estimates the temporal variation in rainfall extremes. Models with lower AIC and RMSE values and higher NSE generally exhibit a better fit to the data. In this study, the optimal model was identified by the AIC value, with a lower value indicating a better fit. If the AIC value of the non-stationary model is higher than that of the stationary model, then the model with the lower RMSE and higher NSE will be selected as the best-fit model.

2.5.1. Model Selection Based on Model Complexity

The Akaike Information Criterion (AIC) is formulated as follows:

$$\mathrm{A}\mathrm{I}\mathrm{C}=2\times (\mathrm{D}-\widehat{\mathrm{L}})$$

(3)

where D is the number of parameters of the statistical model; $\widehat{\mathrm{L}}$ is the log-likelihood function evaluated. The model with lower AIC is considered to have a better fit.

2.5.2. Model Selection Based on Minimum Residual

Root mean square error (RMSE) and the Nash–Sutcliff efficiency (NSE) coefficient are valuable and widely used tools in hydrological and climatological modeling. Both methods aim to minimize the residuals.

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{R}\mathrm{S}\mathrm{E}}_{\mathrm{i}}^2}{\mathrm{n}}}}$$

(4)

$$\mathrm{N}\mathrm{S}\mathrm{E}=1-{\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{R}\mathrm{S}\mathrm{E}}_{\mathrm{i}}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{z}}_{\left(\mathrm{i}\right)}-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}({\mathrm{Z}\left)\right)}^2}}$$

(5)

where RSE_i is the residual for the i-th data point, z_(i) is the observed value for the i-th data point, and mean(Z) is the mean of the observed values. A perfect fit is associated with RMSE = 0 and NSE = 1, given RMSE ∈ [0,$\mathrm{\infty }$) and NSE ∈ [$-\mathrm{\infty }$, 1].

2.6. Stationary and Non-Stationary Return Level Estimates

Extending the stationary return level to non-stationary conditions, where it changes over time and with a physical covariate, causes uncertainties in both the return period and the return level [22,23]. In a stationary model, the parameters of the probability distribution function remain constant over time and are not influenced by another covariate. On the other hand, an m-year return level corresponds to an m-year return period [23]. The return period T_i is the inverse of the exceedance probability q_i for quantile Q_i (i.e., T_i = 1/q_i). Due to a one-to-one relationship between return level and return period, in stationary statistical models, return level curves are defined as follows:

(T_i; Q_i,) T_i > 1 yr, i = 1, …

where T_i is the return period and Q_i is the return level (extreme rainfall intensity) [24,25].

If return levels vary over time and with physical variables, it causes significant uncertainties in both the return period and return level [23,24]. To analyze non-stationary return levels, two concepts are suggested: expected waiting time [25] and effective return level [12]. This study utilized effective return levels to estimate the return periods of extreme rainfall and examine changes in the intensity of these rainfalls over time. The effective return level is defined as q where Q changes as a function of a covariate (time). Therefore, effective return level curves are calculated for a fixed value RP = 1/q, where q is the annual exceedance probability. The effective return level can be obtained using the following equation [9,12]:

$$\left(\right({\mathrm{x}}_{\mathrm{C}}; {\mathrm{Q}}_{\mathrm{q}} \left({\mathrm{x}}_{\mathrm{c}}\right)); \mathrm{q}\in [0;1\left]\right)$$

(6)

where x_C is the covariate, and Q_q (x_c) is the q-quantile [12,13,14,15,16,17,18,19,20,21,22,23,24].

The magnitude of return period changes in rainfall extremes was estimated using expected waiting time. The expected waiting time until an event of magnitude Q_i is exceeded, calculated based on the probability of exceedance over time. Olsen et al. [26] (for the first time) and, later, Salas and Obeysekera [25] extended it, providing a comprehensive mathematical description of the suggested concept by Wigley [27]. The probability of the event to exceed Q_q0 at time m is given by Salas and Obeysekera [25]:

$$\left(\mathrm{f}\mathrm{m}\right)={\mathrm{q}}_{\mathrm{m}}·{\prod }_{\mathrm{t}=1}^{\mathrm{m}-1}\left(1-{\mathrm{q}}_{\mathrm{t}}\right)$$

(7)

where f(1) = q1 and f(m) = 1.

The event Q_q0 is defined as the event with the exceedance probability at time t = 0 equal to $\mathrm{q}0$. Under non-stationary conditions, at time $\mathrm{t}=1,$ the probability of exceedance of ${\mathrm{Q}}_{\mathrm{q}0}$ will be $\mathrm{q}1$; at time $\mathrm{t}=2$, it will be $\mathrm{q}2$, and so on. The cumulative distribution function (cdf) of a geometrical distribution based on the above equation gives the return period under non-stationary conditions.

$$\mathrm{T}=\mathrm{E}\left(\mathrm{X}\right)=1+{\sum }_{\mathrm{x}=1}^{\mathrm{x}\mathrm{m}\mathrm{a}\mathrm{x}}{\prod }_{\mathrm{t}=1}^{\mathrm{x}}\left(1-\mathrm{p}\mathrm{t}\right)$$

(8)

3. Results

3.1. Trend Analysis

To investigate if climate change has affected the patterns of extreme rainfall intensity in the study area, we applied the MK trend analysis to historical data. According to the Mann–Kendall trend test, extreme rainfall shows a significant trend in two stations. Other stations show an insignificant trend, with 55.6% and 44.4% showing positive and negative trends of extreme rainfall, respectively (Figure 2). This demonstrates that extreme rainfall is influenced by climate change in the study areas. According to Porporato and Ridolfi [28], even a weak trend (a small but consistent change over time) can significantly affect the probability analysis. Ganguli and Coulibaly [29] also emphasize that weak trends can influence the outcomes of the statistical models used in analysis. Thus, we evaluated the performance of both non-stationary and stationary models at all stations to determine the best GEV model for frequency analysis of extreme rainfall, even if the rainfall extremes at those stations did not exhibit significant trends or non-stationarity.

3.2. Model Selection

The result indicated no significant difference between stationary and non-stationary models. However, based on the goodness-of-fit test, the non-stationary model shows better performance in fitting the GEV distribution. We observed better performance in fitting the tail part of the non-stationary GEV distributions based on quantile plots. Figure 3 presents probability and quantile plots for three representative stations: Anzali (north), Zabol (east), and Kerman (southeast). As shown in Figure 3, the p-p and q-q plot non-stationary model provides a better fit to the extreme rainfall data than the stationary model. Considering changes in model parameters over time enhances the model’s accuracy in representing data.

Figure 4 illustrates the spatial distribution results of the best-selected GEV model based on the value of RMSE, NSE, and AIC for each station, indicating that the behavior of extreme rainfall intensity at most stations was modeled using a non-stationary, time-covariate approach. According to Figure 4 and

Table 1. Statistical results include three criteria metrics: RMSE, NSE, and AIC.

Station	Model	RMSE	NSE	AIC	Station	Model	RMSE	NSE	AIC
Abadan	M0	3.09	0.92	496.3	Khoramabad	M0	1.7	0.97	513.4
	M1	2.1	0.95	496.4		M1	1.9	0.96	510.7
Ahvaz	M0	1.4	0.97	522.4	Mashhad	M0	1.05	0.98	454.3
	M1	1.2	0.99	526		M1	1.01	0.98	458.3
Arak	M0	1.2	0.98	501.7	Orumiyeh	M0	1.05	0.98	472.5
	M1	1.3	0.98	504.7		M1	0.97	0.98	476.5
Anzali	M0	1.4	0.97	612.1	Ramsar	M0	1.2	0.98	679.3
	M1	1.4	0.97	609.1		M1	2.2	0.95	682
Babolsar	M0	1.6	0.97	600.7	Rasht	M0	0.97	0.99	578.4
	M1	1.1	0.98	603.2		M1	1.1	0.98	579.7
Bam	M0	1.5	0.97	413.9	Sabzevar	M0	1.5	0.97	430.7
	M1	0.99	0.99	416.9		M1	1.8	0.97	433.2
Birjand	M0	1.2	0.98	424.4	Sanandaj	M0	1.04	0.98	475.8
	M1	1.5	0.98	427.3		M1	1.2	0.98	477.9
Bandarabbas	M0	2.1	0.96	584.2	Shahrekord	M0	1.5	0.97	496.7
	M1	1.9	0.96	589.3		M1	1.4	0.97	498.8
Bushehr	M0	2.9	0.92	574	Shahrud	M0	1.3	0.98	451.2
	M1	2.8	0.92	575.1		M1	1.5	0.97	454.06
Dezful	M0	1.2	0.98	563.1	Shiraz	M0	1.3	0.98	507.1
	M1	1.3	0.98	566.2		M1	1.1	0.98	511.6
Esfahan	M0	1.6	0.97	449.5	Semnan	M0	1.4	0.97	425
	M1	1.2	0.98	449.3		M1	1.6	0.98	428.3
Ghazvin	M0	0.99	0.98	433.4	Tabriz	M0	2.2	0.95	431.8
	M1	0.98	0.98	436.8		M1	2.1	0.95	433.8
Gorgan	M0	1.7	0.97	510.6	Tehran	M0	1.6	0.97	445.2
	M1	1.6	0.97	513.2		M1	1.4	0.97	446.6
Hamedan	M0	2.7	0.93	468.1	Torbat	M0	0.98	0.99	446.9
	M1	3.2	0.91	469.8		M1	1.4	0.97	446.8
Iranshahr	M0	1.4	0.97	481.9	Yazd	M0	1.5	0.97	415.4
	M1	1.4	0.97	484.6		M1	1.4	0.98	417.9
Kerman	M0	3.4	0.89	498.4	Zahedan	M0	1.3	0.98	414.5
	M1	2.3	0.93	491		M1	0.8	0.99	414.6
Khoy	M0	1.4	0.98	451.2	Zabol	M0	1.2	0.98	440.8
	M1	1.6	0.97	454.3		M1	1.3	0.98	437.8
Kermanshah	M0	1.5	0.97	471.9	Zanjan	M0	1.5	0.97	437.8
	M1	0.94	0.99	474.7		M1	1.7	0.96	439.3

We have two models, M0 and M1. For enhanced clarity and distinction, the M0 model is highlighted with a blue background, and the M1 model is highlighted with a white.

, the non-stationary GEV model, which includes time as a covariate, was the best fit for the Esfahan, Krman, and Anzali stations based on RMSE, NSE, and AIC criteria. In two stations, Khorram Abbad and Zabol, the non-stationary time-covariate model provided a better fit for rainfall extremes than the stationary model based on the AIC. Therefore, according to the AIC, the stationary model performed better, whereas based on RMSE and NSE, the non-stationary model performed better at most stations.

3.3. Return Level Estimate

3.3.1. Stationary Return Levels of Extreme Rainfall

Based on the stationary model, the greatest increase in magnitude of return levels of rainfall extremes occurred in the north, south, and southwest areas, with an average increase of 3 percent from a return period of 2 to 10 years. These areas correspond to the maximum return level of extreme rainfall for different return periods. According to Asakareh and Shahbaee [30], extreme rainfall is most intense in the Zagros Mountains, along the Caspian Sea coast in Iran, and in parts of the southern coast of the Oman Sea and the Persian Gulf. The lowest rainfall amounts occurred in the 50- and 100-year return periods (Figure 5). Stations where the return level of extreme rainfall is higher in the stationary model show an increasing trend in the intensity and frequency of extreme rainfall extremes. Conversely, stations with a lower return level in the stationary model than in the non-stationary model show a decreasing trend.

3.3.2. Non-Stationary Return Level Change

To account for the temporal and spatial changes in extreme rainfall risk due to climate change, the study area’s non-stationary frequency analysis results for extreme rainfall with return periods of 2, 10, 25, 50, and 100 years were interpolated from all stations.

The magnitude of extreme rainfall events occurring within the 10-, 25-, 50-, and 100-year return periods increased under non-stationary conditions across 52.8% of the study area. Meanwhile, the stationary-based method with constant GEV parameters underestimates the risk of extreme rainfall across the study area. It demonstrates the significant uncertainty that arises from ignoring the impact of non-stationary conditions due to climate change on extreme rainfall.

Non-stationary analysis indicated that less than half of the analyzed stations experienced a decreasing trend in extreme rainfall intensity. These stations are distributed across all regions except the west. The trends in the severity and frequency of rainfall extremes show spatial inconsistency. However, the west and southwest areas show consistency in increasing trends in the intensity and frequency of extreme rainfall, while the east and southeast show consistency in decreasing trends (Figure 6). The stations in the western part of the country show an increased trend in rainfall intensity. Changes in atmospheric pressure in the upper atmosphere over the eastern Mediterranean region affect rainfall patterns in the west. The southwestern regions show more variability and heterogeneity in extreme rainfall than other areas. This variability in extreme rainfall intensity is influenced by Sudanese low-pressure systems and the sheltering effect of the mountains, resulting in different weather patterns in this area. In the northern regions, extreme rainfall is increasing from west to east. Anzali and Rasht in the west of the Caspian Sea show a decreasing trend. These findings are consistent with Masoudian and Darand [31], who reported a decreased extreme rainfall intensity in northern Iran. The intensity and frequency of extreme rainfall in all stations in the southeast show a decreasing trend, except for the Kerman station. The behavior of extreme rainfall at Kerman station differs from other stations, likely due to local conditions and the topography. Kerman is situated in a unique geographical location that includes highlands and mountainous terrains. This topography can significantly influence rainfall patterns, leading to more intense rainfall events in certain areas due to orographic lift. Non-stationary analysis indicated that in 47.2% of the study areas, the non-stationary GEV model predicted lower return levels than the stationary model. However, these areas showed an increasing trend in extreme rainfall under non-stationary conditions. Conversely, in the remaining 52.8%, the non-stationary model predicted higher return levels, with a decreasing trend in extreme rainfall events over time. The trends variation in extreme rainfall based on a time-varying non-stationary model indicates that the stationary assumption underestimates the risk of extreme events.

3.3.3. Non-Stationary Return Periods and Return Levels

In this study, we analyzed changes in the intensity and frequency of extreme rainfall events under non-stationary conditions by incorporating time as a covariate to examine the effective return level of rainfall extremes. Figure 7 shows the non-stationary return levels in the 2-, 10-, 25-, 50-, and 100-year return periods of rainfall extremes for the non-stationary model in the stations that have been classified as non-stationary based on the AIC value. As observed, return levels of extreme rainfall vary with time. The results show that return levels for rainfall extremes have increased with time in three stations, including Khoramabbad, Esfahan, and Kerman in the west, center, and southeast, respectively. This suggests that, over time, the frequency and intensity of rainfall at these stations are expected to rise. In other words, rainfall extremes of a certain magnitude that previously occurred in more than a 2-year return period have now reduced to occurring in 2 years or less. The non-stationary framework highlights the significant uncertainty associated with using stationary assumptions to assess the risk of extreme rainfall events. In two stations, Zabol in the east and Anzali in the north, with the highest rainfall in Iran, the return levels of rainfall extremes decreased significantly with time. Figure 7 illustrates that non-stationarity in rainfall extremes has occurred in areas with different rainfall regimes. Anzali station, located in the high-rainfall northern area along the southern Caspian Sea, experiences the highest levels of rainfall. Esfahan station, located in the arid central region, receives less rainfall than the average.

As illustrated in Figure 8, the expected waiting time for extreme rainfall return periods has become shorter in most stations. This reduction in expected waiting time suggests that rainfall extremes of a certain magnitude, historically occurring at longer intervals, are now occurring more frequently. This indicates that the intervals between rainfall extremes in different return periods have become shorter, indicating that these events happen more often than before. The lowest time occurrence of return periods was in Esfahan in the center. As shown in Figure 8, the expected waiting time for a 40 mm rainfall event in the Esfahan station decreased from 100 years to 30 years, highlighting a substantial increase in the likelihood of such an event. Our results showed that in 22% of the stations, the return periods of rainfall extremes have remained constant or increased in expected waiting time. In 27% of the observed stations, despite experiencing a reduction in extreme rainfall return levels, the expected waiting time for the occurrence of extreme rainfall has decreased, indicating an increase in the frequency of rainfall extremes. This reduction has been 20 years across all return periods. For example, an extreme rainfall event that would previously occur once every 100 years instead occurs every 80 years.

3.4. Tail Behavior of Rainfall Extremes

The non-stationary predictive probability density in Figure 9 evaluates extreme rainfall at three different times: 31, 60, and 93 years from the first observation. The 93-year point extends beyond the actual observation period (62 years), meaning we project the observed trend into the future and infer from it.

The probability density function indicates that the intensity and frequency of rainfall extremes have increased at many observation stations. Additionally, prediction models suggest that this trend will likely persist in the future in these stations, indicating that extreme rainfall events may become more intense over time. As illustrated in Figure 8, the results demonstrated that the majority of the maximum or minimum values were found over the 93-year period (see the tails of the distribution). The study revealed significant variability in the trends of extreme rainfall events. Some stations showed decreased intensity with more frequency, while others exhibited the opposite pattern.

4. Discussion

The results indicate that non-stationary conditions have significantly altered the return periods of extreme rainfall. We observed that trend patterns in the intensity and frequency of extreme rainfall are not uniform across the country. These findings are consistent with the results of Masoudian and Darand [31], who reported non-uniform behavior in the trends of extreme rainfall indices in different regions of the country. This spatial difference in the trends and non-stationarity of extreme rainfall is affected by local-scale spatial factors such as latitude, longitude, elevation, slope direction, and gradient [32,33,34]. The occurrence of non-stationary rainfall extremes was not uniform across different regions of the country. This indicates the influence of local factors such as topography and rainfall climatology on extreme rainfall events in each area. Our results show that changes in extreme rainfall do not follow a consistent pattern. The increasing trend in extreme rainfall in the west of Iran has the most consistent pattern. In eastern, southwestern, northeastern, and northwestern regions, a noticeable decreasing trend of rainfall extremes has been observed. The behavior of rainfall-inducing weather systems over Iran has shifted in recent decades due to the influence of global warming and climate change [32,33]. Recent studies indicate an increase in the frequency of extreme rainfall events. Asakereh and Ashrafi [34] report an increasing trend in the west and southwest. Sari Sarraf and Asakereh [35] observe a rise in extreme rainfall in the western region, particularly the west Zagros Mountains. Masoudian and Darand [31] note a heightened intensity and frequency of extreme rainfall events in recent decades in the southwest and west of the country. Nazaripour et al. [36] highlight more intense rainfall in the Zagros Mountains. The most significant decreasing trends in rainfall extremes have been in the east, southeast, northwest, and west of the Caspian Sea. These findings are consistent with Pajouh and Darand [37], who reported a reduction in extreme rainfall in the eastern regions, Masoudian and Darand [32] in part of the northern regions, and the results of Asakereh [38] for the Zanjan and northwest of Iran. However, our results differ from Masoudian and Darand’s [31] findings, which indicate no trend in extreme rainfall in eastern Iran, and Fakur and Osternoll’s [39] findings regarding the increase in extreme rainfall in northwest Iran. The increase in the intensity and magnitude of extreme rainfall events indicates a higher risk of extreme rainfall-related phenomena such as floods and droughts, especially in arid regions that experience extreme rainfall anomalies. In 47.2% of the 36 examined stations, the stationary model estimated higher return levels of extreme rainfall than the non-stationary model. In these stations, the stationary model provides higher estimates of extreme rainfall return levels across various periods, indicating an increase in the intensity of extreme rainfall in these regions.

5. Conclusions

This study examined whether trends in extreme rainfall have changed using a non-stationary time covariate model in a changing climate. We aimed to identify the best model (non-stationary time covariate and stationary model) for frequency analysis of these extremes, using different criteria, and to estimate the uncertainties in expected return levels with a Bayesian approach. We also analyzed the return levels and return periods of extreme rainfall to evaluate changes in the intensity and frequency of these events using both stationary and non-stationary approaches, independent of model performance. Then, we compared the results obtained from both models. The case study included daily rainfall data from 36 meteorological stations across Iran. Although the evaluation of the models showed that the non-stationary and stationary models performed better based on RMSE, NSE, and AIC evaluation criteria, respectively, there was no significant difference between the two models. This study demonstrates that the non-stationary model, with time-varying scale and location parameters, performed better at five of the examined stations based on the AIC metric. Additionally, this accounts for 20 stations according to the RMSE and NSE metrics. There has been a change in the trend of extreme rainfall in Iran, with both decreases and increases observed. Under the impact of human-induced climate change, the frequency and intensity of extreme rainfall events have changed. Our results showed that about 47.2% of stations show a decreasing trend in extreme rainfall time series from 1960 to 2021. These were roughly consistent with the findings in previous studies [30,31]. The uncertainty of extreme rainfall events revealed significant spatial variability in their probabilities, highlighting notable spatial differences in extreme rainfall risk. We observed that the spatial consistency of trends and changes in the intensity and frequency of extreme rainfall varies across different regions [30,40]. The western and southwestern regions show the most spatial consistency in increasing trends in extreme rainfall intensity and frequency, and the eastern and southeastern demonstrate consistency in decreasing trends. The spatial differences in extreme rainfall events can be attributed to local-scale spatial factors such as geographic coordinates (latitude and longitude), elevation, and the orientation and gradient of slopes, leading to spatial heterogeneity in the trends [41,42,43]. The stationary models may underestimate or overestimate the return levels for the stations with significantly increasing or decreasing trends. Our findings reveal that the return levels of rainfall extremes increase in areas with high and low rainfall. This can cause extreme rainfall anomalies and lead to the simultaneous occurrence of floods and droughts in these areas [29]. Changes in extreme rainfall do not follow a consistent pattern; in some stations in the same area, the trend is decreasing, while in others, it is increasing. The spatial differences in extreme rainfall trends can be attributed to the effects of climate change. These impacts are visible in rainfall characteristics and varying spatial and temporal extremes, including the intensity, frequency, and amount of extreme rainfall. We found that the return periods for extreme rainfall events changed over time, indicating variation in extreme rainfall frequency. Further, previous studies discussed the possible mechanisms for the changes in rainfall extremes in different regions in Iran. The reduction in extreme rainfall in northeastern Iran might be due to changes in the Siberian high-pressure system, which plays a key role in creating frontal rainfall [34,43,44]. Changes in atmospheric pressure in the upper atmosphere over the eastern Mediterranean region affect rainfall patterns in the west [34]. The establishment of Rex blocking and moisture convergence in the southern Mediterranean is linked to increased rainfall in western Iran [45]. Sudanese low-pressure systems significantly contribute to rainfall in southwestern regions and the slopes of the Zagros Mountains. Also, a slight increase in the activity of the Sudanese system has led to more rainfall in the higher parts of southwestern Iran [46].