Non-stationarity of extreme rainfalls and its impact assessment in the northeastern part of India

Climate change has impacted almost every region globally, which has resulted in a rise in extreme precipitation occurrences and, thereby, disastrous consequences. Floods, soil erosion, agricultural damage, landslides, water pollution, ecological devastation, and loss of life are all potential consequences of heavy, intense precipitation. Hence, it is crucial to address and prepare for these events for the future, considering the anticipated rise in high rainfall occurrences globally. In recent decades, frequency analysis of rainfalls with consideration of non-stationarity has gained attention and numerous investigations have been conducted in this domain (El Adlouni et al., 2007; Aissaoui-Fqayeh et al., 2009; Villarini et al., 2009; Tramblay et al., 2013; Agilan & Umamahesh, 2017). To account for the non-stationarity, many researchers have commonly assumed that the parameters of rainfall frequency distributions vary based on the explanatory variables (Coles et al. 2001; Villarini et al., 2009; Gilroy & McCuen, 2012). Gilroy & McCuen (2012) proposed a novel approach for conducting flood frequency analysis that incorporates multiple non-stationary factors. However, there are numerous studies that have identified non-stationarity in hydrological extremes with trends (Agilan & Umamahesh, 2017), and frequency analysis with return levels that have the parameters of a given distribution varying with time (Obeysekera & Salas, 2016; Yilmaz et al., 2017; Hajani & Rahman, 2018; Song et al., 2020; Hajani, 2022; Jayaweera et al., 2023). The studies have unanimously verified that the non-stationary models have been able to better characterize the extreme precipitation than the stationary models. However, such studies are rare and limited for the study region of northeast India that covers the meteorological sub-divisions 2b and 2c. The majority of the existing studies in northeast India reported the trends in either rainfalls or temperature with a stationary assumption (Das & Goswami, 2003; Jain et al., 2013; Laskar et al., 2014; Gharphalia et al., 2018; Pradhan et al., 2019; Datta & Bose, 2020).

Northeast India is one of the highest rainfall receiving regions in the country and floods occur regularly every year, causing major hydrological disasters (Goswami et al., 2010). The powerful monsoonal regimes in the Brahmaputra and Barak valley regions of northeast India are responsible for floods and a source of vulnerability. The region is agrarian and is highly sensitive to extreme hydrologic events. Intense rainfalls triggered floods and landslides, negatively impacting the economic development of the region. Extreme rainfall patterns have experienced sudden shifts as a result of the impact of climate change exacerbated by growing human activities. A combination of various factors, including the region's varied topography, high altitude, and prevalent meteorological conditions, impacts the occurrence of extreme precipitation in northeast India. Observations indicate an increasing temperature trend in the region (Dash et al., 2012; Laskar et al., 2014) and a significant decline in summer monsoon rainfalls (Choudhury et al., 2019). Sreekesh & Debnath, 2016 found that the monthly rainfalls and temperature for five stations in northeast India showed a rising pattern. Past research conducted by Guhathakurta & Rajeevan (2006) and Rupa Kumar et al. (2002) did not identify any notable pattern in rainfalls in northeast India. However, a study conducted by Jhajharia et al. (2012) revealed decreasing trends in 24-h maximum rainfalls in Assam from 1951 to 2003. Dash et al. (2009) observed significant increasing trends in the heavy intensity rainfalls in northeast India. Kiran Kumar & Singh (2021) found the future extreme rainfall events over northeast India to be affected by EL Nino. Deka et al. (2013) witnessed a continuous decreasing trend of rainfalls during monsoon and post-monsoon in the Barak basin. Other major studies for trends in rainfalls in northeast India can be found in the studies of (Ravindranath et al., 2011; Laskar et al., 2014; Das et al., 2015; Gharphalia et al., 2018; Pradhan et al., 2019; Datta & Bose, 2020; Singh et al., 2021).

In recent decades, there has been a collateral shift in urban and rural populations in the northeast region, leading to changes in water availability and causing notable environmental and land use changes that increase the impact of local flooding. Hence, the mere assumption of stationarity for any water resource planning and design may be insufficient under climate change due to the region's constantly changing anthropogenic activities. The climate in the region varies within a few kilometres and is more difficult to interpret accurately in the non-stationary regime. Choosing an appropriate distribution is crucial for accurately representing the rainfall patterns in a certain location. Research shows that two-parameter distributions and three-parameter distributions have been utilized globally to represent severe rainfall occurrences. For example, Shabri et al. (2011) found that the Generalized Extreme Value (GEV) and GLO distributions were the most suitable for representing the statistical characteristics of severe rainfalls in Malaysia. In a research conducted by Trefry et al. (2005) and Adamowski et al. (1996), the GEV distribution was found to be the most suitable distribution for analyzing the frequency of annual maximum rainfall time series.

Only a limited number of studies have conducted frequency analysis of extreme rainfall in the region, such as those by Deka et al. (2009), Bora & Borah (2017), Deka et al. (2011), Koutsoyiannis (2004), and Bora et al. (2017). These studies identify the GEV distribution as the most suitable for modelling extreme rainfall events.

The study by Monish & Rehana (2020) suggest GEV distribution as the best distribution for representing rainfalls in northeast India. In the paper published by Boulomytis et al. (2018), they preferred gamma function distribution for annual maximum daily rainfall calculation in a probabilistic approach. A common aspect of the conducted studies was that frequency analysis was done based on stationarity assumptions with no attribution to climate change. Therefore, the dependence of the parameters of a distribution to identify the behaviour of extreme rainfalls is not available and remains unexplored. Therefore, extreme value distributions that incorporate the dependence of location, scale, and shape parameters to climate change for the region demands necessary investigation. We shall limit our investigation to the case of dependence only on time. To the best of the author's knowledge, no research has been conducted on the topic of climate change's impact on extreme rainfall estimates for the lower Brahmaputra and Barak basins in northeastern India, specifically considering non-stationarity analysis. This research is, therefore, expected to enhance the understanding of the connection between global warming and the occurrence of extreme rainfalls in a highly climate change vulnerable region of India. The objectives of the study are, therefore, as follows: (i) development of various combination of non-stationary models using the two distributions (GEV and Pearson Type III (PE3)) and determine the optimum non-stationary model using each distribution, (ii) selection of the best non-stationary model for each gauge station among the two distributions based on visual comparison and two numbers of performance criteria, and (iii) quantifying the magnitude of the effect of non-stationarity in each station. The current investigation of rainfall frequency analysis considering the non-stationary of extreme rainfall occurrence shall provide valuable information on the climate change scenario and enhance the study of climate impact research in the northeast region.

The annual daily maximum rainfall for the stations is collected from the Indian Meteorological Department (IMD) for a period of 38 years (1980–2017) and the basic site characteristics are presented in Table 1.

To attain the study objectives, the following steps are followed: (i) assessment of randomness and non-stationarity of the time series; (ii) development of stationary and non-stationary models for GEV and PE3 distributions and determining the best non-stationary model using performance metrics; (iii) selection of the best probability distribution to assess non-stationarity for the stations using visual comparison and two performance metrics; and (iv) assessment of the impact of non-stationarity on the 12 stations. In literature, the most commonly used tests for determining the non-stationarity time series is done by using the Augmented Dickey–Fuller (ADF) test, Philips–Perron (PP) test, and the Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test (Wang & Tomek, 2007; Fedorová 2016). The ADF test and the PP test are among the most commonly used for determining the non-stationarity present in a time series (Yoo, 2007). Hence, in the present study, the non-stationarity is determined using both ADF and PP tests. In the ADF test, the presence of a unit root or time series is non-stationary and is taken as the null hypothesis, while the presence of stationarity in the time series is taken as the alternative hypothesis. Generally, for time series data, the ADF test is more flexible and accurate since the rainfall time series data contain trends or serial correlations. A benefit of the PP unit root test compared to the ADF test is that it can withstand different types of heteroscedasticity and serial correlations in handling error terms. Another benefit of the PP test is that the user need not define a lag period for test regression. The test for both the tests is done using the ‘tseries’ package in R platform.

A common statistical procedure for estimating distribution parameters is the use of a maximum likelihood estimator, since this method can be easily extended to the non-stationary case. The maximum likelihood estimation method is used for estimating the parameters of the distribution. The maximum likelihood estimate (MLE) is an often employed and favoured statistical approach for estimating, since it possesses desirable asymptotic properties. It is well recognized as the method generates accurate and reliable estimates of the parameters (Strupczewski et al., 2001; Hamed & Rao, 2019). By optimizing the likelihood function, the MLE approach provides estimates for the parameters of a probability distribution. Obtaining the parameters of extremal distributions requires minimizing the negative log-likelihood function through an iterative numerical procedure. Genetic algorithm is a powerful evolutionary optimization technique and has the advantage of handling a high degree of complexity and is highly efficient in dealing with larger datasets. It explores the solution space in multiple directions through multiple offsprings in less time (Tabassum & Mathew, 2014). Hence, the genetic algorithm has been adopted to obtain a better fit of the estimate of optimal parameters using the MLE method.

In the stationary PE3 model, the parameters are constant and do not vary with time. The estimation of the parameters in the study is carried out with the help of ‘evd’, ‘GA’, ‘EnvStats’, ‘lmomco’, ‘PearsonDS’, and ‘extRemes’ packages in R.

In the study, in addition to the stationary models (GEV and PE3), we have also examined eight variations of the GEV and PE3 distribution parameters. These variations involve assuming linear and quadratic trends for the location parameter, linear and exponential trends for the scale parameter, and their different combinations are presented in Table 2.

Table 2

Stationary and non-stationary models for GEV and PE3 distributions.

Model .	Definition of the model .	Parameters .
M0	Stationary model where all model parameters are constant in time
M1	Non-stationary model where only location parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠;
M2	Non-stationary model where only scale parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M3	Non-stationary model where both location and scale parameters vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ,
M4	Non-stationary model where only location parameter vary with a quadratic trend with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M5	Non-stationary model where location parameter has a quadratic trend with time and scale parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M6	Non-stationary model where only scale parameter vary with an exponential trend with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M7	Non-stationary model where location parameter vary linearly with time and scale parameter vary exponentially with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M8	Non-stationary model where location parameter vary quadratically with time and scale parameter vary exponentially with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,

Model .	Definition of the model .	Parameters .
M0	Stationary model where all model parameters are constant in time
M1	Non-stationary model where only location parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠;
M2	Non-stationary model where only scale parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M3	Non-stationary model where both location and scale parameters vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ,
M4	Non-stationary model where only location parameter vary with a quadratic trend with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M5	Non-stationary model where location parameter has a quadratic trend with time and scale parameter vary linearly with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M6	Non-stationary model where only scale parameter vary with an exponential trend with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M7	Non-stationary model where location parameter vary linearly with time and scale parameter vary exponentially with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,
M8	Non-stationary model where location parameter vary quadratically with time and scale parameter vary exponentially with time	GEV: ⁠, ⁠, PE3: ⁠, ⁠,

View Large

After non-stationary models' development, it is important to identify which model better represents the original data. To select the best model, we use the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). The AIC (Akaike, 1974) and the BIC are the two most often used statistical criteria for model selection. They aim to strike a balance between conflicting factors by penalizing the inclusion of more parameters in the model. The AIC, which penalizes the minimized negative log-likelihood for the number of parameters estimated, (Katz, 2013) and the non-stationary model with the minimum AIC value can be considered as the best model for representing the non-stationarity of extreme rainfalls for the selected station. The BIC often imposes a higher penalty on complicated models compared to the AIC (Panagoulia et al., 2014). The most favoured model among the candidate models is the one that minimizes the selected criterion; however, consideration should also be given to models with values that are near the minimum.

The time series of annual extreme daily precipitation for the period 1980–2017 were tested for non-stationarity using the ADF test. The critical values of the ADF test at 10% significance levels were obtained. It was checked if the p-value of ADF statistic was greater than the critical value of 0.05, i.e., if the null hypothesis can be accepted. The obtained p-values for the ADF test were greater than 0.05 for all the 12 stations and, thus, rejected the null hypothesis (i.e., data with no trend). Therefore, the annual extreme daily rainfall time series was found non-stationary. The obtained ADF test statistic and the corresponding p-values for all 12 stations are presented in Table 3. The benefit of using the Phillips–Perron (PP) test is that no lagged difference is considered and it involves the estimation of test statistics using only the method of ordinary least squares. The PP test was applied to all the 12 stations and it produced p-values greater than 0.05 at 10% significance level, thereby suggesting the considered stations for the study time period of 1980–2017 to be following non-stationarity. To sum up, considering the aforementioned tests, the annual extreme rainfall series for all the 12 stations is assessed as non-stationary and carried for non-stationary frequency analysis.

Now, to model the non-stationarity of rainfall occurrence for the stations, two important distributions, namely GEV and PE3, are considered. A total of eight models, including seven non-stationary models and one stationary model, were formulated in the study to investigate non-stationarity in each station. The widely applied MLE method was employed to estimate the parameters, and genetic algorithm was utilized to find the optimal value of the parameters. The MLE method and genetic algorithm were utilized together to obtain the optimal model parameters described as MLE-GA in this study. The MLE approach is, significantly, a more straightforward method and convenient to execute in determining the most suitable parameters of a distribution. The models are executed applying the MLE-GA method for each station and the results obtained for the best non-stationary model and the corresponding stationary model of a station are presented in Tables 4 and 5. The estimated parameters and the obtained log-likelihood along with their AIC and BIC values are also presented. Table 4 presents the best non-stationary models among the eight models for each station using the PE3 distribution and the evaluated optimal model parameters. The seven non-stationary models from M1 to M8 are applied to each station and the best non-stationary model is selected based on the least value of the AIC and the BIC. For example, the best non-stationary model is evaluated as M7 for the Golaghat station. The M7 model of the PE3 distribution has a linear trend in the location parameter ξ and an exponential trend in the scale parameter α. The shape parameter β is kept constant. The negative log-likelihood for the stationary model of Golaghat is −78.51 and that for the non-stationary model is −39.11. The model M7 is found best on the basis of AIC and BIC, compared to all other models M1–M8. The stationary model AIC and BIC values are also computed and presented alongside the M7 model. The AIC and BIC values of 76.23 and 82.8, respectively, for the non-stationary model M7 to 141.03 and to 147.95 for the stationary model M0 gives a clear indication of fitment accuracy obtained when non-stationarity is considered. For the Choudhoughat station, the best non-stationary model is selected as M2 among the seven non-stationary models as presented in Table 2. The M2 model has a scale parameter with a linear trend with log-likelihood = −200.07, AIC = 387.77, and BIC = 400.71. For the comparison of log-likelihood of the stationary and non-stationary models for the Choudhoughat station, it can be expressed that the non-stationary model provides a closer acceptance to the observed extreme rainfall series, though not much difference is observed in the values of log-likelihood, but the AIC value of the non-stationary model clearly outperforms the stationary model. Thus, by merely applying a simple linear trend to the scale parameter of the PE3 distribution for the Choudhoughat station, the distribution fit for matching the extreme rainfall series is satisfactorily obtained. As both values of the AIC and the BIC for the Choudhoughat station are found to be lesser than the stationary models’ AIC and BIC, the Choudhoughat station is best represented by a non-stationary model with linear dependence in the scale parameter. Similarly, the best non-stationary models using the PE3 distribution for other stations along with their location parameters, the AIC and BIC values, are available in Table 4.

Table 4

Estimated parameters and performance of stationary (S) and best non-stationary (NS) models using the PE3 distribution.

Station .	Model .	Location .			Scale .		Shape .	Log-Likelihood .	AIC .	BIC .
		.	.	.	.	.	.
Choudhoughat	S	102.67	–	–	42.83	–	1.62	−200.46	398.40	401.11
	NS: M2	102.68	–	–	43.28	0.0012	1.62	−200.07	387.77	402.71
Tezpur	S	66.89	–	–	26.97	–	1.15	−168.37	342.75	347.59
	NS: M2	66.71			27.42	0.0023	1.14	−168.25	341.58	341.02
Golaghat	S	50	–	–	26	–	5.28	−78.51	141.03	147.95
	NS: M7	46.4	0.4267	–	1.85	0.0002	6.26	−39.11	76.23	82.78
Neamatighat	S	−132.63	–	–	2.03	–	110.82	−188.15	382.30	387.21
	NS: M8	−127.30	0.0583	0.075	0.95	0.00085	118.73	−183.81	377.62	385.81
Shillong	S	29.98	–	–	33.07	–	3.90	−148.07	304.14	307.06
	NS: M2	30.44	–	–	34.76	−0.0036	3.90	−147.84	303.68	310.23
Kampur	S	−28.85	–	–	5.56	–	22.4	1,506.78	−3,007.57	−3,002.73
	NS: M8	−32.35	−0.075	−0.077	2.38	−0.0008	23.4	1,375.84	−2,739.68	−2,730.02
North Lakhimpur	S	89.5	–	–	28.4	–	2.1	−176.09	358.19	363.02
	NS: M2	93.5	–	–	28.5	0.0099	2.1	−173.51	355.02	361.47
Cherrapunjee	S	264.99	–	–	111.98	–	2.88	−221.24	558.59	593.32
	NS: M2	267.78	–	–	118.74	0.0182	2.88	−279.17	551.59	580.16
Dibrugarh	S	46.02	–	–	11.80	–	8.93	173.10	−340.21	−335.30
	NS: M2	46.43	–	–	15.84	−0.002	4.9	−68.96	145.93	152.48
Dholai	S	−109.88	–	–	4.34	–	58.47	−194.40	394.81	399.64
	NS: M7	−100.5	1.47	–	1.473	0.0002	58.5	−191.65	391.29	397.74
Imphal	S	12.55	–	–	9.73	–	7.16	78.67	−151.34	−146.42
	NS: M3	9.74	0.05	–	9.8	0.002	7.2	76.78	−143.57	−135.38
Beki Rd Bridge	S	19.76	–	–	10.78	–	13.80	600.26	−1,194.53	−1,189.70
	NS: M8	28.54	0.0020	0.069	2.24	0.00047	10.79	316.88	−621.77	−612.11

Station .	Model .	Location .			Scale .		Shape .	Log-Likelihood .	AIC .	BIC .
		.	.	.	.	.	.
Choudhoughat	S	102.67	–	–	42.83	–	1.62	−200.46	398.40	401.11
	NS: M2	102.68	–	–	43.28	0.0012	1.62	−200.07	387.77	402.71
Tezpur	S	66.89	–	–	26.97	–	1.15	−168.37	342.75	347.59
	NS: M2	66.71			27.42	0.0023	1.14	−168.25	341.58	341.02
Golaghat	S	50	–	–	26	–	5.28	−78.51	141.03	147.95
	NS: M7	46.4	0.4267	–	1.85	0.0002	6.26	−39.11	76.23	82.78
Neamatighat	S	−132.63	–	–	2.03	–	110.82	−188.15	382.30	387.21
	NS: M8	−127.30	0.0583	0.075	0.95	0.00085	118.73	−183.81	377.62	385.81
Shillong	S	29.98	–	–	33.07	–	3.90	−148.07	304.14	307.06
	NS: M2	30.44	–	–	34.76	−0.0036	3.90	−147.84	303.68	310.23
Kampur	S	−28.85	–	–	5.56	–	22.4	1,506.78	−3,007.57	−3,002.73
	NS: M8	−32.35	−0.075	−0.077	2.38	−0.0008	23.4	1,375.84	−2,739.68	−2,730.02
North Lakhimpur	S	89.5	–	–	28.4	–	2.1	−176.09	358.19	363.02
	NS: M2	93.5	–	–	28.5	0.0099	2.1	−173.51	355.02	361.47
Cherrapunjee	S	264.99	–	–	111.98	–	2.88	−221.24	558.59	593.32
	NS: M2	267.78	–	–	118.74	0.0182	2.88	−279.17	551.59	580.16
Dibrugarh	S	46.02	–	–	11.80	–	8.93	173.10	−340.21	−335.30
	NS: M2	46.43	–	–	15.84	−0.002	4.9	−68.96	145.93	152.48
Dholai	S	−109.88	–	–	4.34	–	58.47	−194.40	394.81	399.64
	NS: M7	−100.5	1.47	–	1.473	0.0002	58.5	−191.65	391.29	397.74
Imphal	S	12.55	–	–	9.73	–	7.16	78.67	−151.34	−146.42
	NS: M3	9.74	0.05	–	9.8	0.002	7.2	76.78	−143.57	−135.38
Beki Rd Bridge	S	19.76	–	–	10.78	–	13.80	600.26	−1,194.53	−1,189.70
	NS: M8	28.54	0.0020	0.069	2.24	0.00047	10.79	316.88	−621.77	−612.11

View Large

Table 5

Estimated parameters and performance of stationary (S) and best non-stationary (NS) models using the GEV distribution.

Station .	Model .	Location .			Scale .		Shape .	Log-likelihood .	AIC .	BIC .
		.	.	.	.	.	.
Choudhoughat	S	153.68	–	–	41.81	–	−0.103	−199.64	405.28	410.20
	NS: M7	153.60	3.23	–	2.996	0.0006	0.3417	−193.26	396.53	404.72
Tezpur	S	86.85	–	–	19.88	–	−0.0103	−175.36	356.72	361.60
	NS: M6	84.08	–	–	2.94	−0.0003	0.2857	−172.93	353.85	360.39
Golaghat	S	85.55	–	–	26.53	–	−0.10	−181.31	368.61	373.5
	NS: M7	80.08	3.111	–	3.092	−0.0002	0.141	−175.27	360.53	368.72
Neamatighat	S	69.97	–	–	25.05	–	0.222	−180.07	366.15	371.06
	NS: M6	79.65	–	–	3.077	−0.0003	0.23	−172.23	352.48	359.03
Shillong	S	134.16	–	–	70.03	–	−0.0417	−212.4	430.85	435.76
	NS: M6	128.10	–	–	3.994	−0.0001	0.022	−209.79	427.59	434.14
Kampur	S	80.39	–	–	22.62	–	0.201	−179.72	365.41	370.32
	NS: M3	80.17	0.0075	–	21.94	0.0026	−0.089	−177.51	365.02	373.21
North Lakhimpur	S	126.63	–	–	27.4	–	−0.101	−184.81	375.71	380.62
	NS: M7	125.67	3.134	–	3.275	−0.0005	0.369	−181.06	372.11	380.31
Cherrapunjee	S	518.37	–	–	149.95	–	−0.100	−251.6	509.20	514.11
	NS: M4	515.36	0.278	−0.0097	140.22	–	0.0844	−247.47	504.94	513.12
Dibrugarh	S	69.80	–	–	29.63	–	0.1075	−201.56	409.12	414.02
	NS: M8	87.89	3.25	0.0001	2.93	0.00049	0.142	−182.66	377.32	387.15
Dholai	S	130.12	–	–	35.50	–	0.20	−186.82	379.64	384.48
	NS: M3	126.81	3.343	–	3.298	0.0003	0.063	−184.06	378.12	386.18
Imphal	S	71.59	–	–	25.47	–	0.114	−178.5	362.10	367.02
	NS: M3	63.40	0.0139	–	18.63	0.005	0.0192	−175.50	361.01	369.20
Beki Rd Bridge	S	137.66	–	–	37.76	–	0.20	−198.39	402.77	407.68
	NS: M6	138.62	–	–	3.725	−0.0003	0.103	−196.83	401.66	408.22

Station .	Model .	Location .			Scale .		Shape .	Log-likelihood .	AIC .	BIC .
		.	.	.	.	.	.
Choudhoughat	S	153.68	–	–	41.81	–	−0.103	−199.64	405.28	410.20
	NS: M7	153.60	3.23	–	2.996	0.0006	0.3417	−193.26	396.53	404.72
Tezpur	S	86.85	–	–	19.88	–	−0.0103	−175.36	356.72	361.60
	NS: M6	84.08	–	–	2.94	−0.0003	0.2857	−172.93	353.85	360.39
Golaghat	S	85.55	–	–	26.53	–	−0.10	−181.31	368.61	373.5
	NS: M7	80.08	3.111	–	3.092	−0.0002	0.141	−175.27	360.53	368.72
Neamatighat	S	69.97	–	–	25.05	–	0.222	−180.07	366.15	371.06
	NS: M6	79.65	–	–	3.077	−0.0003	0.23	−172.23	352.48	359.03
Shillong	S	134.16	–	–	70.03	–	−0.0417	−212.4	430.85	435.76
	NS: M6	128.10	–	–	3.994	−0.0001	0.022	−209.79	427.59	434.14
Kampur	S	80.39	–	–	22.62	–	0.201	−179.72	365.41	370.32
	NS: M3	80.17	0.0075	–	21.94	0.0026	−0.089	−177.51	365.02	373.21
North Lakhimpur	S	126.63	–	–	27.4	–	−0.101	−184.81	375.71	380.62
	NS: M7	125.67	3.134	–	3.275	−0.0005	0.369	−181.06	372.11	380.31
Cherrapunjee	S	518.37	–	–	149.95	–	−0.100	−251.6	509.20	514.11
	NS: M4	515.36	0.278	−0.0097	140.22	–	0.0844	−247.47	504.94	513.12
Dibrugarh	S	69.80	–	–	29.63	–	0.1075	−201.56	409.12	414.02
	NS: M8	87.89	3.25	0.0001	2.93	0.00049	0.142	−182.66	377.32	387.15
Dholai	S	130.12	–	–	35.50	–	0.20	−186.82	379.64	384.48
	NS: M3	126.81	3.343	–	3.298	0.0003	0.063	−184.06	378.12	386.18
Imphal	S	71.59	–	–	25.47	–	0.114	−178.5	362.10	367.02
	NS: M3	63.40	0.0139	–	18.63	0.005	0.0192	−175.50	361.01	369.20
Beki Rd Bridge	S	137.66	–	–	37.76	–	0.20	−198.39	402.77	407.68
	NS: M6	138.62	–	–	3.725	−0.0003	0.103	−196.83	401.66	408.22

View Large

Non-stationary distributions have been found to overcome the necessity to update assessments by the stationary models and provide closer and more accurate fits. The two information criteria, AIC and BIC, are used in the study to select the most appropriate model among the non-stationary models and further comparison to the stationary models. The models with the lowest values of AIC and BIC are regarded as the best candidate models.

As both the values of AIC and BIC for the Choudhoughat station are found to be lesser than the stationary models’ AIC and BIC, the Choudhoughat station is best represented by the non-stationary model having linear dependence in the scale parameter. The values of log-likelihood are also seen such that the best non-stationary models presented in Table 5 have comparatively lower values of the two information criteria. Similarly, using the PE3 distribution, the 12 stations are evaluated for seven non-stationary models and the best non-stationary model is selected and presented in Table 4. Compared to the PE3 distribution, the non-stationary GEV models for the 12 stations had one extra parameter to define the rainfall distribution. For example, Choudhoughat, Tezpur, Shillong, North Lakhimpur, Dibrugarh, and Cherrapunjee were defined by a change in only the linear relation in the scale parameter. The non-stationary model by PE3 distributions provides a better fit to stationary models with lower values of log-likelihood, AIC and BIC, whereas for GEV distributions, all the 12 stations required a more complex relation of the distribution parameters with time. Tezpur and Shillong had an exponential trend in the scale parameter along with a linear trend in the location parameter, whereas Dibrugarh had an exponential trend in the scale parameter and a quadratic trend in the location parameter with time. The complexity in the non-stationary models was increased. However, the non-stationary models outperformed the respective stationary models for the stations with lower values of log-likelihood, AIC and BIC. An important observation of the shape parameter values obtained for both the GEV and PE3 distributions indicate that the distributions fitted are heavy tailed, thereby suggesting the higher probability of obtaining larger values of extreme rainfalls in the region due to the climate changes in future.

After the estimation of the best parameters for both the non-stationary and stationary models, the rainfall quantiles for various return periods are estimated. This step is carried out to compare and determine the appropriate non-stationary model suitable for each station. A regular method that is adopted is to compare the AIC and the BIC obtained by the non-stationary models of the two distributions and select the distribution with the lowest AIC and BIC values. Generally, the resulting AIC and BIC values obtained in any analysis depend on the number of parameters involved in the model and the type of distribution. As the numbers of parameters involved in the selected best non-stationary model for each station are different and come from different distributions, finalizing the best rainfall distribution seems to be less robust. Moreover, the fit by any non-stationary model of a distribution may over or under-estimate some of the observed values that may be highly significant and thus require a stronger agreement. An over-estimation would produce unnecessary larger rainfall values for higher return periods whereas an under-estimation would give us cause to consider highly unsafe lower design values. Thus, a combined judicious selection by relying on both visual comparison and the calculated AIC and BIC values is more robust and reliable. In the study, five selected years, i.e. 1980, 1990, 2000, 2010, and 2017 are considered for determining the quantiles for each station. The fit by the best non-stationary model alongside the stationary model is plotted in the same plot of quantiles for various return periods computed in those years. The plot of quantiles showcases the effect of climate change on extreme rainfalls as we proceed from the past (1980) to the present (2017). The quantile plots for GEV and PE3 distributions are put together aside to compare the fit of observed data to account for non-stationarity during the study years, and are available in Figures (3)–(6).

For a 100-year return level, the increase or impact due to non-stationarity is found to increase from station to station with the value increasing from 0.05 to 38.6%. For a 20-year return level, the increase in non-stationarity impact follows similarly for the six stations with slightly lesser values. For example, at Neamatighat station, the increase in non-stationarity impact is from 29.03% in 20-year return level to 38.6% for 100-year return level of extreme rainfalls. The study conducted for the period of 1980–2017 suggests that all the 12 rainfall stations show a varying amount of non-stationarity and, thus, the evolution of change in rainfall extremes. Almost the entire northeastern region of India receives increased rainfall during the summer monsoon and the increased impacts of non-stationarity during the study period depict the inclusion of climate change effects prior to any water resources planning and disaster mitigation for the region. The decreasing trend of rainfall extremes in some stations and the increasing trend in some other stations clearly suggest the influence of climate change and other factors triggering the occurrence of extreme rainfalls for the study period.

The study carried out is the first in its domain in northeast India to identify the effect of climate change and to express the non-stationarity of extreme rainfall regimes. The effect of climate change on rainfalls in the region has long been affected, severe occurrence of rainfalls is observed in non-monsoon periods, and trends are highly influenced due to climate change. The growing anthropogenic activities in the region have largely been regarded for the non-stationarity behaviour of extreme rainfall occurrence. The study of non-stationary modelling of extreme rainfalls using time as an explanatory variable clearly in the study suggests non-stationarity to be significantly different from the stationary model. Various degrees of dependence of non-stationarity with the distribution parameters are established and the results obtained demonstrate non-stationarity in extreme rainfalls to significantly vary across space and time in the region.

Many studies carried out previously for the region have investigated the trends and frequency of extreme rainfalls in the region but the application of non-stationary analysis was not studied. The terrain in northeast India is very complex and diversified with different climate types. Extreme precipitation is the leading factor for the occurrence of floods and landslides annually in the region leading to catastrophic damage to life and property. Prokop & Walanus (2015) studied the stations Cherrapunjee and Shillong in the region and suggested an increase in the frequency of extreme events based on 100 years of data. Other studies on increasing trends in the frequency of extreme rainfall events, seasonal rainfalls, and mean annual rainfalls in northeast India were reported by Sahany et al. (2018); Prathipati et al. (2019); Varikoden & Revadekar (2020); and Das et al. (2015). Similar studies reported by Dash et al. (2009) discovered an increase in the number of heavy rainfall events in northeast India; while Sreekesh & Debnath (2016) discovered that there was an increasing trend in the amount of monthly rainfalls in northeast India, many studies reported the decreasing nature of heavy rainfall events in the region. The study of Begum & Mahanta (2022) found decreasing trends in extreme rainfalls considering nine important stations in Assam, whereas the studies by Mohapatra et al. (2021) and Choudhury et al. (2019) reported weakening trends in northeast India in summer monsoonal rainfalls that bring maximum rain to the region. A study by Gharphalia et al. (2018) also found decreasing trends in annual and seasonal rainfalls for the period 1986–2015 in northeastern India. Viswambharan (2019) and Kothawale & Rajeevan (2017) found decreasing rainfall trends over northeast India in the study of historical periods. Our study presented also finds a mixed result of increasing and decreasing trends in extreme rainfall occurrence for the spatially varied stations in the region. The previous studies do not incorporate the effect of climate change on the rainfall occurrence in the northeastern region of India and its variation with years. The investigation of non-stationarity for the extreme rainfalls is thus highly crucial and is an important factor to relate the climate change effect indications in the region.

Unlike stationary model distributions, the non-stationary model for rainfall occurrence for the study period is not defined with constant parameters of distribution. The present study differs from the works carried out by previous authors by employing a new novel non-stationarity investigation as well as comparison to the stationary distributions. Further selection of a distribution for the stations was done by comparing two extreme value distributions for obtaining a more robust result. Of all the stations studied, the stations Neamatighat, Kampur, Dibrugarh, Golaghat, and Kampur required a higher number of distribution parameters forming more complex models with time and show the most gross changing behaviour in density plot from 1980 to 2017. The stationary models for Neamatighat and Kampur are found to shift significantly from the observed records fitted with the GEV distribution while application with the PE3 distribution for the best non-stationary models provides a closer fit with lesser values of AIC and BIC. Stations that required a smaller number of distribution parameters formed less complex models and the shifting nature of the non-stationary models is found to be almost similar to their stationary counterparts. The stations Tezpur, Imphal, North Lakhimpur, and Cherrapunjee required less complex non-stationary models and their performance was superior to the lower values of AIC and BIC. The non-stationarity behaviour of rainfalls in these stations is not very different from that of a stationary one; however, it provides a more accurate fitment of rainfall distribution. Here, the best selected distributions GEV and PE3 required only the dependence of one parameter with time. The three stations Dibrugarh, Golaghat, and Choudhoughat had linear to quadratic dependence of the location parameter and the scale parameter depending exponentially with time. No single distribution was found suitable for all the stations and both the distributions formed complex models in location and scale parameters with time. The present results are in agreement with the studies of earlier researchers Monish & Rehana (2020), Bora et al. (2017); Deka et al. (2009), though the studies were carried out using only stationary distributions in the region. For example, Deka et al. (2011) and Bora et al. (2017) found the Generalized Pareto distribution to perform better than the GEV distribution in describing the annual extreme rainfalls in northeast India. The study carried out by Monish & Rehana (2020) suggests that the distributions for extreme rainfalls are not simply fitted by the GEV distribution. Also, the study by Deka et al. (2009) suggests the GEV distribution to perform better than the PE3 distribution when parameters are estimated using the L-Moments method for extreme rainfalls in northeast India. Overall, the study puts forward the existence of non-stationarity in annual extreme rainfall regimes in the region. Thus, the comparative study of GEV and PE3 distributions for non-stationary modelling of rainfalls suggests that both distributions have time varying distribution parameters in the region and the presence of non-stationarity cannot be ignored.

Furthermore, the current study also ascertains the impacts of non-stationarity on extreme rainfall occurrence in the stations. The impacts of non-stationarities are quantified in terms of percentage differences for return periods of 20 and 100 years and are evaluated in the study. Analyzing all the 12 stations gives an impression that station rainfalls in the study region are under the influence of positive and negative non-stationary impacts. The non-stationary impacts are high for stations lying in Assam as compared to Meghalaya and Manipur. All the stations in the study region displayed significant differences in quantile estimates between stationary and non-stationary models. Thus, the adaptation to climate change under non-stationary conditions is found critically important for the northeastern region of India as the extreme rainfall magnitudes show significant fluctuations in quantities in time and space.

The current study depicted the influence of climate change on the annual extreme rainfall series with time, whereas the influence of other important factors and atmospheric variables affecting extreme precipitation is not studied. The inclusion of other variables affecting extreme precipitation in the region can be investigated as a future study.

This study investigated the first look for non-stationary frequency analysis of extreme rainfalls and examined its non-stationary impact considering both GEV and PE3 distributions. The non-stationary frequency analysis of extreme daily rainfalls without resorting to stationary assumptions was investigated in the southern part of the Brahmaputra basin and the Barak basin in northeast India for the determination of non-stationary rainfall regimes. The study findings reveal that annual daily extreme rainfalls exhibit varying degrees of non-stationarities with time during the whole study period of 1980–2017. The investigation was done using both stationary and non-stationary models applying GEV and PE3 distributions and the parameters were obtained by combining genetic algorithm in the MLE estimation method. Our findings can be summarized as follows:

• 1. The application of non-stationary models with time as an explanatory variable significantly provided more correct relationships in describing the occurrence of extreme rainfalls regime under climate change. The stationary models performed poorly compared to the non-stationary models in capturing the extreme rainfall behaviour of the 12 climate stations in the study area. Therefore, it will be misleading to assume the stationarity behaviour of extreme daily rainfalls in the region for design and planning purposes.

• 2. The comparison of the two distributions investigated in this study suggests dependence on a scale parameter with a simple linear trend. The scale parameters of the 12 stations vary linearly to exponentially, suggesting that the values of rain are more dispersed and exhibit more variations, except for Cherrapunjee, where a location parameter with a quadratic trend gave a better fit. The location parameter is also found to vary with time from linear to quadratic among all 12 stations. The genetic algorithm was effective in providing optimal parameters for the MLE estimation method to GEV and PE3 distributions.

• 3. The values of AIC and BIC of the stationary models are also significantly larger compared to their non-stationary counterparts. Neamatighat, Kampur, Dibrugarh, Golaghat, and Beki Road Bridge stations required more complex models to represent non-stationarity. For the 100-year return level, the increase in impact due to non-stationarity is found to range from 0.05 to 38.6%.

• 4. Both GEV and PE3 distributions were tested for fitting non-stationary and stationary models and compared to each other in selecting the best distribution to all the 12 stations under the effect of climate change, thus ensuring the robustness of the results. The process of determining the best non-stationary model, based solely on the AIC and BIC values, was not accurate. It was found resourceful to combine this step with the visual comparison of the obtained rainfall quantiles.

• 5. No clear-cut superiority of performance of GEV over PE3 distribution or, otherwise, in the fitting of non-stationary models or stationary models is achieved in the study of the stations. The stationary models provided poorer fits to the observed records for all stations and could not fit the maximum values of extreme rainfall records for most of the stations.

The authors thank the Regional Meteorological Centre, Guwahati for its cooperation and providing the climatic data. This study was possible because of the availability, to the scientific community, of various free software packages (‘tseries’, ‘evd’, ‘GA’, ‘EnvStats’, ‘lmomco’, ‘PearsonDS’, and ‘extRemes’) in R statistical software.

No funding was obtained for this study.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.