Impact Evaluation Using Nonstationary Parameters for Historical and Projected Extreme Precipitation

Abstract

Recent improvements in time series studies of hydro-climatological variables have led to the belief that the effects of nonstationarity are substantial enough to call the idea of traditional stationary approaches into doubt. The mean and variability of annual and seasonal rainfall in Pakistan are changing due to anthropogenic climate change. With the use of stationary and nonstationary frequency analysis techniques, this study set out to assess the impacts of nonstationarity in Southern Punjab, Pakistan, over the historical period of 1970–2015 and the future periods of 2020–2060 and 2060–2100. Four frequency distributions, namely Generalized Extreme Value (GEV), Gumbel, normal, and lognormal, were used. The findings of the nonstationarity impact across Southern Punjab showed different kinds of impacts, such as an increase or reduction in the return level of extreme precipitation. In comparison to other distributions, GEV provided the finest fit. In Bahawalnagar, Bahawalpur, Multan, Rahim Yar Khan and DG. Khan, the annual nonstationarity impacts for the 100-year return level were increased up to 15.2%, 8.7%, 58.3%, 18.7%, and 20%, respectively. Moreover, extreme precipitation was found to be increasing during the historical and projected periods, which may increase floods, while less water availability appeared at a seasonal scale (summer) during 2061–2100. The increased nonstationarity effects emphasized adapting these nonstationarities induced by climate change into the design of water resource structures.

1. Introduction

A comprehensive understanding of extreme precipitation is an important challenge in hydrology. Probability distributions for frequency analysis (FA) are a common technique used by engineers and hydrologists to understand extreme precipitation regimes throughout the world. In accordance with the IPCC’s Sixth Assessment Report, anthropogenic activities have been altering the climate with increasing intensity, resulting in major changes in regional climate and increasingly frequent extreme weather events [1]. The design of major infrastructures, such as urban and highway drainage, flood barriers, hydraulic structures, sewerage systems, and dam spillways, relies on the magnitude and frequency of heavy precipitation events [2]. As a consequence, it is important to understand the annual and seasonal distribution of extreme precipitation events. Frequency analysis calculates the extreme precipitation for a specific return period; however, the methods now in use presuppose that the peak time series are uniformly distributed and independent [3,4], or in other words, without trends and unexpected deviations [5]. In fact, the idea of stationarity has been and continues to be used to build water-related development and flood prevention projects all over the world. The climate system has been stressed in recent decades owing to natural changes in global temperature, and human activity has a potential impact on regional climate, which is eventually deteriorating the natural cycle of water [6]. The hypothesis of stationarity has become widely questionable in the context of regional and global change. However, many studies suggested an increase in the nonstationary nature of climate patterns. Keeping this in mind, various studies have attempted to investigate the validity of this concept in flood regimes in many different parts of the world, taking into account the influence of natural climatic variability [7,8,9,10,11,12,13] or land use changes [14,15,16]. Therefore, there is convincing evidence that the stationarity assumptions must be revised under the changing climate [17], and the development and planning of water resources and the design of hydraulic structures need to be evaluated by considering the nonstationarity approach.

The regional variance in precipitation is intimately related to reliable and precise estimations of global climate trends [18]. Changes in precipitation volume and pattern might have a direct impact on the water resources and agriculture of the affected regions [19]. As a result, knowing the spatiotemporal variance in precipitation at the regional scale is critical in climate monitoring and hydro-climate studies [20]. The spatiotemporal variations of precipitation for various parts of the world have been recorded by several studies [21,22,23,24]. There is increasing consensus that long-term variations in precipitation may affect agricultural and hydrological systems [25], which provide the basis of our understanding of both global and local climate change [26]. These precise and trustworthy precipitation data serve as the foundation for our understanding of regional and global climate change as well as any potential effects on water supplies [27,28]. Warmer temperatures lead to greater water vapor content, resulting in more severe precipitation in seasons, although annual precipitation is modestly reduced [29]. Extreme precipitation events can cause severe flooding, property damage, and potential loss of life, especially in urban areas with high population density [2]. Unfortunately, a lot of recorded data indicates that occurrences of severe precipitation are becoming increasingly common in many parts of the world. For instance, Donat et al. [30], analyzed the data from over 11,000 worldwide weather stations between the years 1951 and 2010, and it is concluded that severe precipitation occurrences increased in both the dry and wet parts of the world. Hence, there is substantial evidence of severe precipitation variations.

Economic losses from weather- and climate-related disasters have increased with large spatial and inter-annual variability. There is a lack of data at the local level on disasters and disaster risk reduction, which can constrain improvements in local vulnerability reduction [31]. Pakistan is listed among the top vulnerable countries to climate change, which is experiencing catastrophic climate extreme events, and these extreme events are predicted to increase in frequency, intensity, and duration. The climate change scenario is alarming because of the existence of climate-change-related risks in Pakistan, including sea level increase, glacial depletion, high average temperatures, floods, and high-frequency droughts, and half of the country’s population is at risk [32]. Moreover, Pakistan has been listed as the eighth most impacted country in the world by extreme occurrences, according to the Global Climate Risk Index (CRI) [33]. In Pakistan, extreme weather events, including heatwaves, intense rainfall, floods, and droughts, have become more frequent and violent during the past few decades [34]. Numerous studies that compared precipitation statistics for various time periods and examined precipitation extremes and mean changes typically used trend analysis (e.g., the Mann–Kendall test and Sen’s slope estimate) [35]. Hence, it is evident that most of the studies remained silent, considering nonstationary approaches in evaluating return levels of extreme precipitation. The Punjab province of Pakistan produces the majority of the nation’s agricultural commodities and is extremely vulnerable to changes in most meteorological parameters due to extreme events [36]. Furthermore, to address the gaps in earlier research, local-scale studies evaluating changes in intense precipitation are necessary. Many investigations have been undertaken over the years to demonstrate how the concept of stationarity (for example, the notion of regular return intervals) may be adversely affected as a result of environmental change [37,38]. For instance, Aziz et al. [39] focuses on the transient inconsistency in yearly as well as seasonal extreme rainfall in Turkey, which was investigated by using stationary and nonstationary frequency approaches. Sertac et al. [40] studied the nonstationary investigation of extreme rainfall and observed that nonstationary (NST) models outperformed the stationary model at the 17 stations. Nashwan et al. [41] investigated how a change in the environment (nonstationarity) affects the average and intensity of annual and seasonal rainfall in Malaysia and the magnitude of maximum rainfall at 2-, 10-, 25-, 50-, and 100-year return periods has increased at most of the stations. In terms of environmental change, the assessment of the recurrence time of exceptional precipitation occurrences according to static assumptions may thus be insignificant. The frequency of extreme temperatures, droughts, floods, and rainfall has risen during the last century [42,43]. However, most studies only used historical and CMIP5 GCMs to assess different aspects of the hydrological cycle, such as atmospheric precipitation [44,45,46,47,48] and surface hydrological processes [49,50]. Therefore, these studies are evident in the validity of historical and future precipitation scenarios and should also be updated with SSP2 and SSP5 scenarios.

Changing climate and land utilization circumstances may modify the probability of hydrodynamically significant occurrences, implying that the characteristics (of location, shape, and scale) of fundamental distributions may vary over time, rendering the view of stationarity unsustainable [39]. The return levels of extreme precipitation can be more accurately predicted for any particular return period from a nonstationary perspective. On the other hand, the predicted return levels are vulnerable to oversights like overestimation or underestimation when nonstationarity is not taken into consideration. This study examines the evaluation of nonstationary impacts for the return levels of extreme precipitation at the local scale.

In this work, we provide an update on nonstationary impacts in return levels of extreme precipitation at the Southern Punjab. Specifically, the aim of this study was (a) to determine the best distribution function under stationary and nonstationary assumptions, and (b) to quantify the nonstationary impacts at historical and projected extreme precipitation. Overall, study findings will aid in predicting the changes in historical and projected extreme precipitation by considering nonstationary parameters in terms of a changing climate.

2. Study Area and Datasets

2.1. Study Area

Pakistan geographically lies between 24°–37° N and 60°–75° E in southwest Asia. The country has an area of 8 × 106 km². The country has a diverse landscape ranging from the Karakoram and Himalayan mountains in the north and northwest, with the agricultural plains of the Indus River basin in the center and the Arabian Sea along the southern coast [51]. The province of Punjab, Pakistan, is the second largest province, and its districts are distributed in the southern, northern, and central parts. This study was carried out in South Punjab, Pakistan and lies between 28.0° to 30.5° North latitude and 69.0° to 73.0° East longitude (Figure 1). The total area of the study region is 205,345 square kilometers. Despite the fact that the region is primarily flat, the southwest includes a few hilly areas. In addition to Cholistan, a desert region of the country, there is a part of the Koh-e-Suleman mountain range that exists towards the west of the study region.

South Punjab experiences a subtropical climate characterized by hot summers and mild winters. The region has two distinct seasons: a hot and dry season from April to June and a monsoon season from July to September. The average maximum and minimum temperature ranges from 38 °C to 45 °C and 25 °C to 30 °C, respectively during the summer period. In the winter, the average maximum and minimum temperatures drop down ranges from 20 °C to 25 °C and 7 °C to 12 °C, respectively. The average annual precipitation varies across different locations within the region. The average annual precipitation in the region ranges between 100 and 400 mm. The majority of rainfall occurs during the monsoon seasons from July to September. During this period, the region receives a significant portion of its annual precipitation. The rest of the year tends to be relatively dry, with sporadic and limited rainfall.

2.2. Datasets

In South Punjab, the impacts of nonstationarities on annual and seasonal min-max precipitation were investigated for the actual period (1970–2015) and for the chosen high-resolution model for the prediction period (2020–2100). Five meteorological stations from the Pakistan Metrological Department (PMD) and five GCM models were used for the analyses that took place in the South Punjab region.

Daily rainfall measurements (1971–2015) were taken from the Pakistan Met Department, and future projections (2020–2100) were taken from the GCM models, which are in phase 6 (CMIP6). (https://esgfnode.llnl.gov/projects/cmip6/ (accessed on 15 January 2023). Numerous researchers analyzed precipitation and maximum and minimum temperature variables using multiple GCM (more than one) because of computational limitations, resolution of data, availability of data against selected emission scenarios, and to compare different models’ performance [52,53,54,55]. For instance, the AWI-CM-1-1-MR model provides both maximum and minimum temperature data at a 100 km resolution, while CESM2, CESM-WACCM, CESM-WACCM-FV2, and CMCC-CM2-SR5 only provide precipitation data at the same resolution. EC-Earth3-AerChem includes all the necessary input variables for only one SSP scenario, and EC-Earth3-CC, HadGEM-GC31-MH, and GFDL-CM4 provide data for two SSP scenarios on a daily basis at a 100 km resolution. Additionally, some of the GCMs, such as FGOALS-3-H, GFDL-C4C192, and SAMO-UNION, do not provide future projections of numerous meteorological inputs. Based on the data required for further hydrological application, we chose the models listed in

Table 1. Selected CMIP6 Model for Precipitation Data of SSP2 and SSP 5 Scenarios.

Institution ID	Model Name	Resolution (Long, Lat)
NOAA-GFDL	GFDL	1.3° × 1.0°
MIROC	MIROC6	1.4° × 1.4°
MPIESM1-2HR	MPI-ESM1-2HR	0.9° × 0.9°
MPI-M	MPIESM1-2LR	1.9° × 1.9°
MRIESM2-0	MRIESM2-0	1.12° × 1.12°

.

The climate projections of the five selected models were downloaded under shared socioeconomic pathways (SSP) scenarios 2 and 5. In the current era of industrialization, an immediate drastic decrease in the emission of greenhouse gases is difficult, and it is unlikely to meet this scenario [56]. The scenario SSP2 with medium stabilization forcing scenario and the scenario SSP5 with very high radiative forcing scenario are considered in this study.

3. Methodology

Both stationary and nonstationary probability distributions are applied to estimate the effects of nonstationarities. Precipitation analyses incorporate the GEV, Gumbel, normal, and lognormal distributions. Global Circulation Models (GCMs) have been utilized for future forecast periods, while historical analysis has been performed using observed data. Figure 2 provides an explanatory flowchart for the approach used in this study.

3.1. Statistical Analysis

The initial analysis of studied meteorological stations was calculated, including mean, standard deviation (STD), coefficient of variance (Cv), coefficient of skewness (Cs), and coefficient of kurtosis (Ck) for annual precipitation data. The summary of these statistics is given in

Table 2. Mean annual precipitation statistics of the study area.

Stations	Latitude	Longitude	Elevation	Mean	STD	Cv	Cs	Ck
Bahawalnagar	29°20′	73°51′	161.05	254.34	126.55	0.50	0.42	−0.59
Bahawalpur	29°20′	71°47′	110	177.58	119.34	0.67	1.67	4.08
Multan	30°12′	71°26′	121.95	197.48	104.59	0.53	0.46	0.90
Rahim Yar Khan	28°26′	70°19′	82.93	112.35	78.70	0.70	0.86	0.39
DG Khan	30°03′	70°38′	148.1	170.12	109.83	0.65	0.85	0.56

.

The maximum annual mean value of precipitation was found 254 mm at Bahawalnagar in the northeast of the study area. The values of the coefficient of kurtosis and coefficient of skewness varied between −0.59 to 4.08 and 0.42 to 1.67, respectively, of annual precipitation. The values of the coefficient of skewness and kurtosis should be in the range of 0 and 3 for normal distribution. Cs indicated that data on precipitation are positively distributed but with comparable differences among stations. The coefficient of variance was calculated to find out the spatial variability of precipitation at each station. The values of the annual precipitation coefficient of variances varied between 50% and 70%.

3.2. Stationarity in Data Series

The Mann–Whitney U (MWU) test is a non-parametric test developed by Mann and Whitney and is commonly used in frequency analysis of extreme events for the stationarity assumption [57]. Mann–Whitney U Test was performed for the annual and seasonal maximum extreme precipitation. The results stated that data exhibited stationary as significance level alpha = 5% is lesser than the estimated p value as shown in

Table 3. p value (two-tailed) for analyzed series by Mann–Whitney U Test.

Station	Annual	Winter	Spring	Summer
Bahawalnagar	0.09	0.15	0.65	0.15
Bahawalpur	0.88	0.03	0.58	1.00
DG Khan	0.79	0.71	0.84	0.58
Multan	0.38	0.17	0.49	0.38
Rahim Yar Khan	0.30	0.2	0.76	0.31

. Furthermore, the details about the Mann–Whitney U Test can be searched out in the studies of [58,59].

3.3. Stationary and Nonstationary Frequency Analysis

The effects of nonstationarities are computed using both stationary and nonstationary probability distributions. Extreme value distributions are used in hydrology and water resources to examine the probability behavior of severe occurrences like droughts and flooding. The four distributions of generalized extreme value, Gumbel, normal, and lognormal, have been used in this study to carry out the return level for a specific return period. Numerous researchers used these probability distributions for the analysis of return levels of hydro-meteorological variables [39,40,41,60,61,62,63,64]. The analysis of precipitation was carried out for the historical and projected period. The parameters of the distributions are rendered time-dependent by including time as a covariate, allowing for their application in nonstationarity settings [39,65,66]. Coles et al. [67], expressed (the return level is defined as a value that is anticipated to equal or be exceeded with a probability of 1/T once per unit of time (T)) of how return level and return period might alter in nonstationary climate circumstances. Therefore, both stationary and nonstationary GEV, Gumbel, normal, and lognormal probability distributions are applied. The GEV probability distribution function of the nonstationarity form is presented in Equation (1) [65,66].

$$F \left(z, {\theta }_t\right)=exp \left\{-{\left[1+\epsilon \left(\frac{z-{\mu }_t }{{\sigma }_t}\right)\right]}^{\frac{-1}{\epsilon }}\right\}$$

(1)

where ${\theta }_t$ is a time-dependent set of GEV parameters, which comprises a constant shape parameter ($\epsilon$) and a time-dependent location ${(\mu }_t)$ and scale ${(\sigma }_t)$ parameters. The shape parameter ($\epsilon$) is typically not modeled as a time-variant parameter since it is more challenging to predict properly [66,67]. Therefore, the location parameter turns µ_t = µ₁+ µ_2t, and the scale parameter turns σ_t = σ₁+ σ_2t in the case of nonstationary distributions. Moreover, the location and scale parameters under nonstationary conditions are represented by intercept (slope), respectively. Additionally, “t” is the temporal explanatory covariate that causes (µ_t) and (σ_t) to be nonstationary for a studied period ranges from one to the entire number of data periods.

The GEV distribution transforms itself into the Gumbel distribution after the shape parameter gets closer to zero. Furthermore, the other two parameter distributions, normal and lognormal, were also applied under both stationary and nonstationary conditions. By reducing the negative log-likelihood function, distributional parameters were evaluated using the maximum likelihood estimation method. The log-likelihood functions for a stationary and nonstationary variant of the GEV distribution are presented in Equations (2) and (3), respectively.

$$l (\mu , \sigma , \epsilon , x) =-n log \sigma -\left(\frac{1}{\epsilon }+1\right) {\sum }_{i=1}^n{log \left[1+\epsilon \left(\frac{x_i-\mu }{\sigma }\right)\right]-{\sum }_{i=1}^{n}log \left[1+\epsilon \left(\frac{x_i-\mu }{\sigma }\right)\right] }^{\frac{-1}{\epsilon }}$$

(2)

$$l ({\mu }_t, {\sigma }_t, \epsilon ) =-n log \sigma -\left(\frac{1}{\epsilon }+1\right) \sum _{i=1}^{n}log \left[1+\epsilon \left(\frac{x_i-{\mu }_t}{{\sigma }_t}\right)\right] {-\sum _{i=1}^{n}log\left[1+\epsilon \left(\frac{x_i-{\mu }_t}{{\sigma }_t}\right)\right] }^{\frac{-1}{\epsilon }}$$

(3)

where i represents the total number of years from 1 to n and sigma $(\sigma )$ represents the scale parameter of the probability distribution. Equations (4) and (5) can be used to represent the log-likelihood functions for stationary and nonstationary normal distributions, respectively. For both stationary and nonstationary lognormal distributions, the log-likelihood function may be expressed as Equations (6) and (7).

$$l (\mu , \sigma , x) =-\frac{n}{2} ln (2 \pi )-\frac{n}{2} ln (2 \pi {\sigma }^2)-\frac{1}{2 {\sigma }^2} {\sum }_{i=1}^m({x_j-\mu )}^2$$

(4)

$$l ({\mu }_t, {\sigma }_t, x) =-\frac{n}{2} ln (2 \pi )-\frac{n}{2} ln (2 \pi {\sigma }_t^2)-\frac{1}{2 {\sigma }_t^2} {\sum }_{i=1}^n({x_j-{\mu }_t)}^{`2}$$

(5)

$$l (\mu , \sigma , x) = \frac{n}{2} ln (2 \pi {\sigma }^2)-{\sum }_{i=1}^{m}ln (x_i)-\frac{{\sum }_{i-1}^{m}lnln (x_i{)}^2 }{2{\sigma }^2}+ \frac{{\sum }_{i-1}^{m}lnln \left(x_i\right) \mu }{{\sigma }^2}-\frac{n{\mu }^2}{2 {\sigma }^2}$$

(6)

$$l ({\mu }_t, {\sigma }_t, x) = \frac{n}{2} ln (2 \pi {\sigma }_t^2)-{\sum }_{i=1}^{m}ln (x_i)-\frac{{\sum }_{i-1}^{m}lnln (x_i{)}^2 }{2{\sigma }_t^2}+ \frac{{\sum }_{i-1}^{m}lnln \left(x_i\right) {\mu }_t }{2{\sigma }_t^2}-\frac{n{\mu }_t^2}{2 {\sigma }_t^2}$$

(7)

Parameter estimation of these distribution functions (with and without the nonstationarity assumption) can be performed using different packages available in the R-Programming. To name a few, these packages include ‘GEVcdn’ by [68], “ismev” by [69], “extremes” such as [70,71]. We opted for “ismev” package, which enables the estimation of distribution parameters using the “Nelder–Mead” algorithm for optimization under stationarity and nonstationarity assumptions. The temporal evolution of parameters is performed by using a time vector as a covariate. Hence, in the case of nonstationary GEV, the estimated location parameter (μt = μ₀ + μ_{1 *t}) has two components. The μ₀ (first component) is constant (which is analogous to intercept), and μ₁ (second component) is time-multiple (analogous to the slope of the time-variant parameter). This makes the whole parameter μ_t time-variant or nonstationary. Furthermore, the details about the distribution fitting with and without nonstationarity assumptions can be found in [69]. A similar procedure was followed for normal and lognormal distribution fitting under the nonstationarity assumption.

The probability distribution with the lowest negative log-likelihood (NLLH) is projected to offer greater compatibility since the parameters of probability distributions were assessed by minimizing the (NLLH) functions. The return levels for the 100-year return period were evaluated after finalizing the parameters of the distributions for stationary and nonstationary assumptions. The equation used to explain the effects of nonstationarities is expressed as percentage changes between 100-year stationary and nonstationary return values.

$$\mathrm{Non}-\mathrm{stationarity} \mathrm{impacts}= \frac{[non stationary return value-stationary return value]}{stationary return level}\times 100$$

(8)

Nonstationarity has an increasing (or decreasing) influence, which means that greater (or lower) return levels develop when nonstationarity were taken under consideration. Greater (or lower) return levels show that nonstationarity has increased or decreased the likelihood of an extreme event recurring within a certain return period.

4. Results

4.1. Selection of Best Probability Distribution Function

The parameters of each distribution were determined by lowering the negative log-likelihood (NLLH) and akaike information criteria (AIC) values. It has been proposed that the distribution with the lowest NLLH and AIC value is the best. Figure 3 provides bar charts of the NLLH and AIC values of each distribution for stationary and nonstationary conditions for yearly maximum precipitation (MP). Each bar chart shows five NLLH values (left) and AIC values (right), one for each of the five stations investigated in this study. The NLLH and AIC statistics of each distribution for nonstationary conditions are frequently lower than those for stationary conditions. This implies that when time was included as a covariate, each distribution fit comparatively well. This demonstrates that the assumption of nonstationarity improves the fit of distributions relatively. In all stationary and nonstationary conditions, NLLH and AIC values for normal distributions were shown to be bigger than NLLH and AIC values for each of the three probability functions. As indicated by bar charts in almost every station, the NLLH and AIC values of the GEV function were lower than those of the other probability functions in most of the precipitation indices under both stationary and nonstationary conditions. For most precipitation indices, the Gumbel distribution has lower NLLH and AIC values than the other two distributions (normal and lognormal) in both stationary and nonstationary conditions. Furthermore, the statistics reveal that the NLLH and AIC values of lognormal distributions are not much higher than those of Gumbel and GEV distributions. Mainly, the GEV distribution will be utilized for future periods since it has the lowest NLLH and AIC values and has traditionally been used more commonly in literature to undertake frequency analysis of extremes because of its tendency to exhibit the best match for tails.

4.2. Impacts of Nonstationarity for the Historic Period (1970–2015)

4.2.1. Yearly Maximum Precipitation (MP)

Figure 4 displays the percentage fluctuations for the 100-year stationary and nonstationary return levels obtained using four probability distributions for yearly MP. With a few notable differences, a visual examination of the maps in Figure 3 demonstrates that all four distributions have proposed analogous consequences for annual MP. There are a number of comparable effect objectives across the stations that are supported to varying degrees by all four distributions. The findings demonstrated that during yearly MP, all four distributions had positive effects on all stations. The level of effects generated by the GEV and lognormal distributions is shown to be greater when compared to the other two distributions. Overall, all of the stations have positive effects from all four probability distributions.

4.2.2. Seasonal MP

The impacts of nonstationarity on seasonal maximum precipitation are presented in Figure 5. Overall, a pattern of nonstationarity impact from all four distributions is seen in the study of seasonal MP. However, a few locations have very minor variations. Seasonal patterns show the importance of seasonality on the return level of the 100-year precipitation over the course of historical time under changing climatic circumstances. While larger return level values are projected throughout the research region in the autumn, the impact type differs between stations in the winter, spring, and summer.

All four distributions have shown more or less comparable effects from nonstationarity during the span of the winter, according to a study of the data for winter MP. Analyzing the four probability distributions, the findings indicate that all of the stations have positive effects. Positive effects’ nonstationarity variability (up to 68%) is seen, especially with the Gumbel distribution over the Multan and DG-Khan stations. Additionally, consequences of nonstationarity (up to 54%) with a normal distribution were seen at the stations in Bahawalpur, Multan, and DG Khan.

During spring, similar nonstationarity impacts have been observed for four probability distributions. There are certain consistent trends in impact types across all stations. Except for GEV, the results show that the three probability distributions have a negative impact on Rahim Yar Khan. All other stations, on the other hand, exhibit a positive impact. The results suggest that adopting GEV has a positive impact (10% to 50%) across all stations. Positive impacts are found with all other probability distributions (6% to 98%) across all stations except Rahim Yar Khan, where negative impacts (9% to 15%) were identified with the three previously mentioned probability distributions.

The results of the nonstationary impact study show that positive impacts were observed during summer MP. The nonstationary impacts of summer MP are very comparable to those of annual MP. The explanation for this is that the majority of the year’s extreme precipitation occurred during the summer season. During the summer MP, all four distributions had a positive impact on all stations, according to the results.

A nonstationary impacts study from all distributions indicated evidence of positive impacts throughout the autumn season. GEV and lognormal distributions, in particular, emphasize this property. However, there were certain locations with a lesser nonstationary impact from the Gumbel and normal distributions. The results demonstrate that utilizing GEV has a positive impact (10 to 65%) across all stations. Positive impacts (up to 63%) were identified across all stations except Bahawalnagar, where negative impacts (9%) were determined using the Gumbel distribution. Moreover, at Bahawalpur, negative impacts (up to 18%) were determined using the Gumbel and normal distributions.

4.3. Impacts of Nonstationarity for Annual Projected Period (2020–2100)

Figure 6 shows the classified maps of the projected models’ (ssp2, ssp5) nonstationary impacts for 100-year return levels using GEV distribution for yearly MP. The estimated time span is divided into two parts: (2020–2060) and (2061–2100). At ssp2 (2020–2060), all projected models show positive nonstationary impacts across all stations. Nonstationary influences are greater in Bahawalnagar and Rahim Yar Khan than in the other three stations. The projected models indicate both types of nonstationary impacts (positive and negative) at ssp2 (2061–2100). Positive effects can be observed in DG-Khan, Multan, and Bahawalpur, where the nonstationary return level of extreme precipitation has increased. Furthermore, except for the GFDL-ESM4 projected model, all of the projected models in Bahawalnagar exhibit negative nonstationary impacts. In Rahim Yar Khan, three of the five projected models exhibit negative nonstationary impacts. Positive nonstationary impacts have been found at all stations during the ssp5 (2020–2060), with the exception of DG-Khan, where three of the five projected models show negative nonstationary impacts. The negative nonstationary impacts dominate the positive nonstationary impacts across all stations in the ssp5 (2061–2100). Except for the GFDL-ESM4 projected model, all other projected models have been found to have negative nonstationary impacts in Rahim Yar Khan, Bahawalpur, and Bahawalnagar. All of the projected models represent a positive impact on the DG-Khan and Multan stations. The nonstationary impacts of the projected model are mostly varying due to variations in the extreme values of the GCMs models.

4.4. Projected Seasonally MP

Table 4. Projected nonstationary impacts (%) for 100-year return levels using GEV probability distribution during 2020–2060.

Station	GCM Model	GFDL		MIROC6		MPI-ESM1-2HR		MPI-ESM1-2LR		MRISEM1-0
	Season	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5
Bahawalnagar	Winter	25.8	−5.3	23.4	−3.5	22.7	1.9	21.5	1.8	24.1	−4.4
	Spring	33.0	2.4	30.3	−2.9	30.2	−2.6	29.4	−3.8	31.3	2.0
	Summer	33.3	9.5	31.1	2.5	32.7	3.0	32.3	2.6	32.4	2.8
	Autumn	13.9	3.4	14.0	11.3	16.1	5.2	13.7	12.9	14.5	6.9
Bahawalpur	Winter	−4.0	12.6	−9.9	13.5	2.6	13.2	2.1	12.6	−3.3	13.8
	Spring	23.5	7.8	20.8	5.2	20.7	3.0	19.9	5.0	21.8	5.4
	Summer	3.7	3.2	4.1	2.7	2.7	2.1	3.7	3.6	2.9	3.4
	Autumn	5.6	9.9	7.5	−5.2	10.4	−3.1	8.3	−7.1	4.4	−3.7
DG Khan	Winter	5.6	16.0	5.2	15.9	5.0	−5.4	4.7	−4.5	5.2	−6.4
	Spring	−11.3	3.3	−15.3	−11.0	3.5	2.9	2.5	−10.6	3.9	−11.3
	Summer	3.9	1.5	3.8	5.6	3.3	7.3	3.7	6.5	2.6	8.2
	Autumn	8.7	17.1	8.3	13.2	8.0	11.1	12.5	12.8	8.3	13.8
Multan	Winter	−2.5	12.6	−1.9	11.8	2.4	11.5	−6.5	10.8	3.2	12.2
	Spring	19.2	5.3	17.6	15.3	17.5	14.9	12.4	14.5	18.2	15.4
	Summer	2.2	−6.4	2.5	2.4	3.9	2.5	5.4	−5.5	2.0	3.4
	Autumn	12.8	8.4	14.6	12.9	16.8	12.4	11.5	11.7	19.4	4.3
Rahim Yar khan	Winter	11.1	2.5	9.8	4.3	9.4	−4.2	8.8	−5.9	10.1	−4.5
	Spring	−10.5	3.7	−4.2	2.5	−6.3	−9.9	3.9	−9.2	4.1	−11.2
	Summer	12.9	15.5	13.2	5.8	13.5	8.4	12.4	9.5	11.2	10.3
	Autumn	6.9	6.4	6.5	8.8	8.5	5.5	10.2	9.0	6.6	9.8

depicts the projected nonstationary impacts (%) for 100-year return levels generated by the GEV distribution at a seasonal scale. During winter, all the projected model shows higher positive nonstationary impacts in Bahawalpur, Rahim Yar Khan, as compared to the other three stations during the ssp2 (2020–2060). The negative nonstationary impacts (up to 6%) were observed in the DG-Khan and Bahawalpur for all the projected models except MPI-ESM1-2hr, MPI-ESM1-2lr, and MRIESM1-0, which show little positive nonstationary impacts. The Positive nonstationary impacts (up to 35%) have been identified for all projected models in Bahawalpur, Bahawalnagar, and DG-Khan during the spring at ssp2 (2020–2060). Except for MPI-ESM1-2HR, MPI-ESM1-2LR, and Mriesm1-0, two of the three forecast models indicated lower extreme precipitation in Multan and Rahim Yar Khan. Positive nonstationary impacts (up to 15%) have been observed in Bahawalpur and DG-Khan during the ssp5 (2020–2060). All of the forecast models revealed positive (up to 4%) and negative (up to 10%) nonstationary impacts at the remaining three stations. Positive nonstationary impacts (up to 35%) have been found for all forecast models across all stations in the ssp2 (2020–2060). During the ssp5 (2020–2060), all stations had positive nonstationary impacts (up to 20%) for all projected models, with the exception of the DG-Khan station, which had two of the five projected models have negative nonstationary impacts (up to 5%). Positive nonstationary impacts were observed throughout all the stations (up to 25%) for all the projected models in the SSP2 (2020–2060). Higher nonstationary return levels (positive nonstationary impacts) of extreme precipitation were found across all stations for the five projected models in the SSP5 (2020–2060), with the exception of Bahawalpur, where one of the five projected models represent positive nonstationary impacts (up to 8%), and the remaining models form negative nonstationary impacts (up to 5%).

During the time period ssp2 (2061–2100), positive nonstationary impacts (up to 35%) were observed in the DG-Khan, Multan, and Rahim Yar Khan generated from the GEV probability distribution using five projected models. A mix type of impacts (positive and negative) was found in the Bahawalnagar and Bahawalpur stations from all the projected models. Two of three projected models, GFDL-ESM4 and MPIESM1-2LR, show positive nonstationary impacts (up to 5%) in the Bahawalnagar and Bahawalpur stations, while the remaining three projected models exhibit negative nonstationary impacts (up to 3%) in the above-mentioned particular stations. Positive nonstationary impacts (up to 15%) are observed in Bahawalnagar and DG-Khan in the ssp2 (2061–2100). The mix of impacts (positive and negative) was observed in the other three stations for all the projected models as shown in

Table 5. Projected nonstationary impacts (%) for 100-year return levels using GEV probability distribution during 2061–2100.

Station	GCM Model	GFDL		MIROC6		MPI-ESM1-2HR		MPI ESM1-2LR		MRISEM1-0
	Season	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5	SSP2	SSP5
Bahawalnagar	Winter	3.8	7.0	−2.9	−9.7	−3.5	−10.6	7.8	8.3	−4.5	−10.3
	Spring	12.6	2.3	12.2	7.9	12.1	6.1	8.9	6.9	10.1	8.4
	Summer	2.1	8.0	−11.2	−11.7	−9.5	−13.3	−9.4	−11.0	−10.5	−12.6
	Autumn	2.8	3.9	−7.5	3.1	−5.2	3.3	−7.5	3.2	2.3	3.5
Bahawalpur	Winter	4.6	7.1	−4.8	4.1	−4.6	6.8	2.7	4.8	−4.6	6.9
	Spring	2.2	3.1	−1.8	4.0	−2.8	5.2	3.8	4.8	4.5	5.8
	Summer	4.1	4.7	3.3	−12.9	4.2	−14.3	4.3	−11.2	3.9	−9.3
	Autumn	3.7	2.6	−5.5	−6.5	−12.5	−2.5	−7.5	2.2	−10.8	−1.5
DG Khan	Winter	25.4	11.6	27.2	8.5	26.5	9.8	27.1	8.6	26.7	9.1
	Spring	3.5	5.0	−13.4	−14.9	−12.2	−15.1	5.3	7.5	−12.2	−14.6
	Summer	30.1	17.5	28.4	22.4	27.5	16.8	25.5	12.6	22.5	24.1
	Autumn	12.6	2.0	6.4	−1.6	7.9	−1.3	9.2	1.7	13.9	−1.9
Multan	Winter	18.9	7.7	19.5	5.9	18.3	6.2	18.9	6.1	18.3	6.5
	Spring	5.5	18.8	12.2	12.7	9.4	13.1	15.2	12.8	7.4	12.7
	Summer	4.2	9.4	3.2	14.5	4.5	14.9	4.9	13.9	5.9	17.4
	Autumn	11.2	5.1	7.5	−3.3	9.9	−7.2	13.5	3.2	8.3	−8.5
Rahim Yar khan	Winter	1.5	1.4	12.6	2.6	15.5	3.2	11.7	3.6	9.5	4.6
	Spring	5.0	5.6	4.9	8.8	−4.7	8.8	−4.9	8.7	−7.7	9.0
	Summer	3.5	6.0	−10.5	−11.1	−10.4	−13.7	8.4	−9.4	−2.8	−10.3
	Autumn	2.7	6.3	5.1	5.1	8.8	3.3	11.6	4.7	4.8	5.0

. All of the stations in the ssp5 (2061–2100) demonstrate positive nonstationary impacts (up to 25%), with the exception of Multan, where two of the three models show positive nonstationary impacts (up to 5%), while the remaining three models exhibit negative nonstationary impacts (up to 20%). Positive nonstationary impacts (up to 30%) have been observed at all stations in the ssp2 (2061–2100) except Bahawalnagar and Rahim Yar Khan, where both positive and negative impacts were noted for all forecast models. Negative nonstationary impacts (up to 25%) of severe precipitation were seen in the far future ssp5 (2061–2100) across all stations except DG-Khan and Multan, where higher nonstationary return levels (positive impacts) of extreme precipitation occurred for all projected models. In the SSP2 (2061–2100), all stations formed positive nonstationary impacts except Bahawalnagar and Bahawalpur, where both types of impacts (positive and negative) are found in all the projected models. At ssp5 (2061–2100), the Bahawalpur and Rahim Yar Khan stations have been identified as having higher positive nonstationary impacts (up to 10%), whereas the other three stations have both positive (up to 7%) and negative nonstationary impacts (up to 10%) for all forecast models.

5. Discussion

GEV has the best-fitting properties in both stationary and nonstationary settings out of the four probability distributions. Additionally, GEV offered a better fit for nonstationary distributions compared to stationary distributions. In the research examining nonstationary impacts over the historical and projected eras by Aziz et al. [39], the same conclusions were drawn. This fact was brought about by the use of the same methodology for finding the NLLH using the maximum likelihood estimate. The analysis showed that all distributions exhibit positive nonstationarity impacts at the annual scale; however, it is not necessary for all distributions to have positive nonstationarity effects. For instance, the same four distributions were used for the assessment of nonstationary impacts by Aziz et al. [39] and demonstrated both positive and negative nonstationary impacts at an annual scale. This is because the return period of 30 years used by [39] might be the cause of this inconsistency. Another distinction is that the climate of the study area varies from one place to another or sub-basin to another throughout the world [72]. Hence, the positive impacts of nonstationarity are more likely linked with high quantiles extremes of precipitation [73,74].

Figure 4 and Figure 5 explain the nonstationary impacts raised up to 60% during the historical period for annual and seasonal maximum precipitation. Considering these results in a broader aspect to any region around the globe, it emphasizes the climatic change adoption in the water development sector since more frequent and intensified extreme precipitation events are expected under continuous climate change in the future [30,75,76,77]. Such nonstationary changes in extreme precipitation cause more significant damage or even the failure of existing water conservation facilities. Moreover, enhanced impacts of nonstationary also affect the service life of hydraulic construction projects [78] and other urban infrastructure [79,80] or increase the risk of bridge and highway damage or even failure [81].

Furthermore, nonstationary conditions lead to an increase or decrease in the return level of extreme precipitation that may have a positive or negative nonstationary impact [39]. Since the number of hot and rainy days, land use patterns, and soil moisture levels fluctuate on a seasonal basis, the return levels of extreme precipitation increase in a complex way [30]. According to the current study, the return level of extreme precipitation increased over the winter season, indicating that there were more instances of extreme precipitation.

The results also showed that the return level of extreme precipitation under nonstationary conditions is inconsistent for both the historical and future periods. The complexity of the climate system under nonstationary conditions is explained by these tendencies in these results [73]. Furthermore, the historical data in time-dependent variations of the return level for a certain return period would not remain true for a longer length of time in the future. Nonstationary impacts in terms of nonstationary return level for extreme precipitation may increase in the twenty-first century. Moreover, Aziz et al., [61] showed more consistent nonstationarity effects for temperature extremes at seasonal and annual scales from historical to future times. The investigation of historic-based analysis may only be helpful for short-run strategy and decision-making, while analysis based upon both historical and projected data may be more suitable for long-term strategy planning.

6. Conclusions

In this study, the return level of extreme precipitation was estimated under nonstationary and stationary assumptions to evaluate nonstationary impacts for historical and projected extreme precipitation. The conclusions of the study are as follows:

-

Although GEV (initially having three components) is a widely used probability distribution. The findings of this study reveal that alternative distributions are also capable of comparing nonstationary impacts. The less complicated distributions (having two parameters) might prove advantageous at a particular station.

-

The increase in the return level (magnitude) of extreme precipitation in winter and spring showed causes of flood events, and the reduction in return level of extreme precipitation in summer and autumn may cause less water availability.

-

The projected increase in nonstationarity impacts (up to 50%) distinguished the climate change in the region and emphasized the nonstationarity in the design of hydraulic structures (Reservoirs, Barrages, and others).

No doubt, nonstationary exists in the data series of hydro-meteorological parameters. This study was conducted in a flat area that is arid where a nonstationary impact of up to 60% is observed during the historical period for annual and seasonal maximum precipitation. The findings emphasize the water development sector to adopt the climatic change in terms of nonstationarity impacts. However, the results based on historical data analysis can also be evaluated with other parameter estimation methods of probability distributions, such as linear moments (L-moments). Moreover, the results based on projected data by GCM show uncertainties because, in general, it suffers from great deficiencies and errors, being indicative in nature. The CMIP6 precipitation projections are quite uncertain, thus requiring further assessment and correction [82,83].