Regional Frequency Analysis of Extreme Rainfall in Haryana Using L-Moments

Z^DIST = [ τ₄^DIST – τ_bar₄ + B₄ ] / S₄

Where τ₄^DIST is the L-kurtosis of the fitted candidate distribution, τ_bar₄ is the regional average sample L-kurtosis, B₄ is the bias, and S₄ is the standard deviation of τ_bar₄ obtained via Monte Carlo simulations (usually 500+ trials). The fit is deemed acceptable if |Z^DIST| ≤ 1.64 (at the 90% confidence level). If multiple models pass, the one with the smallest |Z^DIST| is selected as the best-fit distribution.

Region	Distribution	Z-Statistic Value	Location (ξ)	Scale (α)	Shape (k / γ)	Fit Status
Region I	GLO**	-0.10	0.9324	0.2111	-0.1869	Best Fit
	GEV*	-1.18	0.8099	0.3149	-0.0262	Satisfactory
	GNO*	-1.39	0.9254	0.3727	-0.3856	Satisfactory
Region II	GLO*	1.32	0.9038	0.2813	-0.1985	Satisfactory
	GEV**	-0.72	0.7414	0.4154	-0.0439	Best Fit
	GNO*	-1.19	0.8939	0.4964	-0.4100	Satisfactory
Region III	GEV*	1.20	0.7866	0.3958	0.0398	Satisfactory
	GNO*	1.08	0.9312	0.4523	-0.2974	Satisfactory
	PE3**	0.44	1.0000	0.4808	0.8800	Best Fit

6. Regional Growth Curves and Rainfall Quantile Estimates

Using the **index-flood method**, the regional growth curve q(F) for each homogeneous region is calculated using the best-fit distribution parameters. The site-specific rainfall quantile at return period T (non-exceedance probability F = 1 – 1/T) is then computed as:

Q_i(F) = μ_i × q(F)

Where μ_i is the site-specific mean maximum monthly rainfall (the index flood). Table 6 presents the regional quantiles (in mm) scaled for a site with the average regional mean, and Table 7 details the station-wise estimated rainfall depths for various return periods (T = 2 to 100 years).

Region	T = 2 yrs (F=0.5)	T = 5 yrs (F=0.8)	T = 10 yrs (F=0.9)	T = 20 yrs (F=0.95)	T = 50 yrs (F=0.98)	T = 100 yrs (F=0.99)	T = 200 yrs (F=0.995)
Region I (GLO)	310.38	421.61	501.54	586.46	712.67	822.24	945.79
Region II (GEV)	161.75	250.49	311.58	372.12	453.63	516.71	581.59
Region III (PE3)	207.68	306.16	367.12	422.50	490.61	539.52	586.64

Region	Station	T = 2 yrs	T = 5 yrs	T = 10 yrs	T = 20 yrs	T = 50 yrs	T = 100 yrs
Region I (GLO)	Ambala	277.53	376.99	448.45	524.38	637.24	735.21
	Karnal	259.43	352.40	419.20	490.19	595.68	687.26
	Jagadhari	345.31	469.06	557.98	652.46	792.88	914.77
	Kalka	360.29	489.41	582.19	680.77	827.28	954.47
Region II (GEV)	Sirsa	139.18	214.55	265.42	314.89	379.78	429.26
	Hansi	100.89	155.52	192.39	228.26	275.30	311.16
	Farukhnagar	173.25	267.06	330.38	391.96	472.73	534.32
	Faridabad	231.37	356.66	441.22	523.47	631.34	713.58
	Mahendragarh	152.34	234.84	290.51	344.67	415.70	469.85
	Khol	127.39	196.38	242.94	288.22	347.62	392.91
	Palwal	190.58	293.79	363.44	431.19	520.04	587.79
	Bhiwani	119.03	183.48	226.98	269.30	324.79	367.10
	Tohana	132.24	203.86	252.18	299.19	360.85	407.86
	Sohana	192.26	296.38	366.64	434.99	524.63	592.97
	Dujana	189.01	291.37	360.44	427.63	515.75	582.94
	Salhawas	151.32	233.27	288.57	342.36	412.92	466.71
	Beri	186.22	287.06	355.12	421.31	508.13	574.33
Region III (PE3)	Hisar	156.11	230.14	275.97	317.60	368.80	405.56
	Sonipat	243.18	358.50	429.88	494.73	574.49	631.75
	Rohtak	199.77	294.49	353.14	406.41	471.92	518.96
	Nuh	239.59	353.21	423.54	487.43	566.01	622.43
	Jhajjar	222.96	328.68	394.13	453.59	526.71	579.21
	Bawal	221.12	325.97	390.87	449.84	522.36	574.42
	Panipat	196.19	289.22	346.81	399.13	463.47	509.67
	Kurukshetra	209.64	309.05	370.59	426.49	495.25	544.61
	Narwana	186.02	274.23	328.84	378.44	439.45	483.25
	Kaithal	202.61	298.68	358.16	412.19	478.63	526.34

7. Code Tutorial: Implementing L-Moments in Python and R

To enable researchers to perform these calculations, we provide two ready-to-use snippets demonstrating L-moment estimation and extreme-value fitting.

R Script (Using the `lmom` Package)

# Install and load the lmom library
if (!requireNamespace("lmom", quietly = TRUE)) install.packages("lmom")
library(lmom)

# Example: Maximum monthly rainfall data for a station
rainfall_data <- c(150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250)

# 1. Compute sample L-moments (L1, L2, L3, L4)
sam_lmom <- samlmu(rainfall_data)
cat("Sample L-Moments:\n")
print(sam_lmom)

# 2. Extract L-Cv, L-Cs (t_3), L-Ck (t_4)
# Note: samlmu returns L-location, L-scale, L-skewness (t_3), L-kurtosis (t_4), etc.
l_cv <- sam_lmom[2] / sam_lmom[1]
cat("L-Cv:", l_cv, "\nL-Cs (t_3):", sam_lmom[3], "\nL-Ck (t_4):", sam_lmom[4], "\n")

# 3. Fit a Generalized Extreme Value (GEV) distribution
gev_params <- pelgev(sam_lmom)
cat("\nFitted GEV Parameters:\n")
print(gev_params)

# 4. Estimate quantiles for T = 10, 50, and 100 years
return_periods <- c(10, 50, 100)
probabilities <- 1 - 1 / return_periods
quantiles <- quagev(probabilities, gev_params)

# Display results
results <- data.frame(ReturnPeriod_Yrs = return_periods, Quantile_mm = quantiles)
print(results)

Python Script (Using the `lmoments3` Package)

import numpy as np
# Note: install via: pip install lmoments3
import lmoments3 as lm
from lmoments3 import distr

# Example: Maximum monthly rainfall data
rainfall_data = [150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250]

# 1. Compute sample L-moments and ratios
lmom_ratios = lm.lmom_ratios(rainfall_data, nmom=4)
print("Sample L-moments and Ratios:")
print(f"Mean (L1): {lmom_ratios[0]:.4f}")
print(f"L-scale (L2): {lmom_ratios[1]:.4f}")
print(f"L-skewness (t3): {lmom_ratios[2]:.4f}")
print(f"L-kurtosis (t4): {lmom_ratios[3]:.4f}")

# 2. Fit a Generalized Extreme Value (GEV) distribution
fitted_gev = distr.gev.lmom_fit(rainfall_data)
print(f"\nFitted GEV Parameters: {fitted_gev}")

# 3. Compute return level quantiles for T = 10, 50, and 100 years
return_periods = [10, 50, 100]
for T in return_periods:
    F = 1 - 1 / T
    quantile = distr.gev.ppf(F, **fitted_gev)
    print(f"T = {T:3d} years (F = {F:.2f}) -> Quantile: {quantile:.2f} mm")

8. Agricultural and Engineering Implications

The results of this regional study have critical applications for the development and policy planning of Haryana:

• Hydraulic Structures: For Region I (Wet zone, fitted to GLO), designs must accommodate larger return-period rainfall quantities, where a 100-year event can exceed 950 mm in Kalka.
• Agricultural Drainage: In Region II (Dry zone, GEV) and Region III (Central zone, PE3), drainage infrastructure must cope with 50-year rainfall events ranging from 270 mm to 570 mm depending on the exact location. Over-designing can waste valuable rural infrastructure budget, while under-designing can cause widespread waterlogging of sensitive agricultural crops, ruining seasonal yields.
• Water Harvesting: Estimating return levels helps calculate maximum design inflows for farm ponds, reservoirs, and check dams, helping farmers store surplus rainwater for dry season irrigation.

References

• Babu, V. B. and Hooda B. K. (2018). Fuzzy Majority Approach for Modeling Spatial and Temporal Distributions of Daily Rainfall in Western Zone of Haryana. International Journal of Agricultural and Statistical Sciences, 14(1), 57-67.
• Greenwood, J. A., Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form. Water Resources Research, 15(5), 1049-1054.
• Hosking, J. R. M. (1990). L-moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. Journal of the Royal Statistical Society (Series B), 52(1), 105-124.
• Hosking, J. R. M. and Wallis, J. R. (1993). Some statistics useful in regional frequency analysis. Water Resources Research, 29(2), 271-281.
• Hosking, J. R. M. and Wallis, J. R. (1997). Regional frequency analysis: An approach based on L-Moments. Cambridge University Press, United Kingdom.
• Hooda, B. K. (2006). Probability Analysis of Monthly Rainfall for Agricultural Planning At Hisar. Indian Journal of Soil Conservation, 34(1), 12-14.
• Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability-weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resources Research, 15, 1055-1064.
• Malekinezhad, H. and Garizi, A. Z. (2014). Regional frequency analysis of daily rainfall extremes using L-moments approach. Atmosfera, 27(4), 411-427.
• Majumder A., Patil S. G., Noman M. D., and Biswas S. (2015). Application of L-moments for regional frequency analysis of maximum monthly rainfall in West Bengal, India. Mausam, 66(2), 273-280.
• Nain, M. and Hooda B. K. (2019). Probability and Trend Analysis of Monthly Rainfall in Haryana. International Journal of Agricultural and Statistical Sciences, 15(1), 221-229.
• Sahrin S., Ismail N., and Alias N. E. (2018). Regional frequency analysis on peninsular Malaysia using L-moments. Far East Journal of Mathematical Sciences (FJMS), 103(8), 1379-1398.

b_r = n^-1 ∑_i=1ⁿ [ (^i-1_r) / (^n-1_r) ] X_i:n

L-Moments Definition

The first four L-moments (λ_r) are linear combinations of the PWMs:

• L-Location (Mean): λ₁ = β₀
• L-Scale: λ₂ = 2β₁ - β₀
• L-Skewness measure: λ₃ = 6β₂ - 6β₁ + β₀
• L-Kurtosis measure: λ₄ = 20β₃ - 30β₂ + 12β₁ - β₀

To characterize distributions independently of their scale, we define dimensionless L-moment ratios:

• L-coefficient of variation (L-Cv, τ): τ = λ₂ / λ₁
• L-coefficient of skewness (L-Cs, τ₃): τ₃ = λ₃ / λ₂
• L-coefficient of kurtosis (L-Ck, τ₄): τ₄ = λ₄ / λ₂

2. Database and Initial Data Screening

The study utilizes maximum monthly rainfall data for the 48-year period (1970–2017) obtained from the National Data Centre, Indian Meteorological Department (IMD), Pune, covering 27 rain gauge stations in Haryana.

Before executing RFA, the assumptions of stationarity, randomness, and independence must be verified for all stations:

• Stationarity: Tested using the Mann-Kendall trend test. Results showed that only 3 out of 27 sites (Karnal, Kaithal, and Bhiwani) had a statistically significant trend, meaning the regional maximum rainfall series can be treated as stationary.
• Randomness: Tested using the Run test. Except for Rohtak, the rainfall series across all other sites were random.
• Independence: Evaluated using the Autocorrelation Function (ACF). Only Rohtak and Kurukshetra showed significant autocorrelation at lag-1. Overall, it is highly reasonable to treat the data as time-independent and suitable for regional frequency analysis.

Station Name	Mann-Kendall Trend (Tau)	MK P-value	Interpretation	No. of Runs	Run P-value
Sirsa	-0.024	0.810	No Trend	24	0.771
Narwana	-0.108	0.282	No Trend	20	0.145
Hisar	-0.126	0.210	No Trend	24	0.770
Karnal	0.285	0.005	Trend	18	0.054
Ambala	-0.162	0.106	No Trend	22	0.381
Jhajjar	-0.155	0.126	No Trend	20	0.233
Hansi	-0.143	0.060	No Trend	25	1.000
Sonipat	-0.160	0.112	No Trend	22	0.243
Rohtak	-0.109	0.074	No Trend	16	0.008 (Not Random)
Panipat	-0.121	0.230	No Trend	19	0.080
Farukhnagar	-0.114	0.259	No Trend	20	0.145
Faridabad	0.101	0.315	No Trend	22	0.381
Kurukshetra	-0.005	0.810	No Trend	18	0.074
Mahendragarh	-0.141	0.160	No Trend	26	0.770
Kaithal	-0.296	0.003	Trend	21	0.243
Khol	0.029	0.776	No Trend	24	0.780
Palwal	-0.164	0.110	No Trend	15	0.710
Bhiwani	-0.246	0.017	Trend	30	0.074
Tohana	-0.078	0.439	No Trend	25	1.000
Sohana	-0.153	0.129	No Trend	26	0.770
Bawal	-0.025	0.810	No Trend	20	0.145
Jagadhari	0.114	0.255	No Trend	21	0.243
Dujana	-0.041	0.693	No Trend	18	0.074
Salhawas	-0.176	0.080	No Trend	22	0.381
Nuh	-0.031	0.763	No Trend	20	0.145
Kalka	-0.196	0.053	No Trend	19	0.136
Beri	-0.012	0.915	No Trend	22	0.381

3. L-Moments and Station-Wise Characteristics

For each of the 27 sites, sample L-moments and L-moment ratios were computed. The values represent the mean maximum monthly rainfall (in mm), L-coefficient of variation (L-Cv), L-skewness (L-Cs), and L-kurtosis (L-Ck).

Station	Mean (mm)	L-Cv (τ)	L-Cs (τ3)	L-Ck (τ4)
Sirsa	154.33	0.270	0.153	0.177
Narwana	200.021	0.290	0.153	0.084
Hisar	167.863	0.249	0.115	0.176
Karnal	278.356	0.224	0.097	0.071
Ambala	297.777	0.185	0.208	0.250
Jhajjar	239.739	0.271	0.130	0.047
Hansi	111.727	0.312	0.268	0.137
Sonipat	261.487	0.254	0.166	0.132
Rohtak	214.803	0.276	0.100	0.076
Panipat	210.958	0.262	0.143	0.099
Farukhnagar	191.856	0.344	0.197	0.217
Faridabad	256.224	0.244	0.205	0.216
Kurukshetra	225.420	0.297	0.222	0.116
Mahendragarh	168.708	0.303	0.236	0.238
Kaithal	217.858	0.274	0.141	0.067
Khol	141.079	0.356	0.169	0.161
Palwal	211.058	0.257	0.318	0.309
Bhiwani	145.872	0.270	0.172	0.096
Tohana	151.323	0.283	0.115	0.140
Sohana	225.158	0.261	0.149	0.100
Bawal	237.758	0.213	0.112	0.217
Jagadhari	370.504	0.251	0.289	0.241
Dujana	209.315	0.308	0.161	0.132
Salhawas	167.579	0.347	0.295	0.263
Nuh	257.627	0.264	0.167	0.215
Kalka	386.581	0.234	0.153	0.212
Beri	206.223	0.322	0.183	0.109

4. Formation and Validation of Homogeneous Regions

To define homogeneous regions, the mean monthly rainfall values were subjected to hierarchical cluster analysis (Ward's Method). The Elbow Method (analyzing the within-cluster sum of squares) indicated that the optimal number of regions is three.

• Region I (Wet/Semi-humid zone - 4 stations): Ambala, Karnal, Jagadhari, and Kalka.
• Region II (Dry/Semi-arid zone - 13 stations): Sirsa, Hansi, Farukhnagar, Faridabad, Mahendragarh, Khol, Palwal, Bhiwani, Tohana, Sohana, Dujana, Salhawas, and Beri.
• Region III (Central/Transition zone - 10 stations): Hisar, Sonipat, Rohtak, Nuh, Jhajjar, Bawal, Panipat, Kurukshetra, Narwana, and Kaithal.

Discordancy Test (D_i)

The discordancy measure D_i (Hosking and Wallis, 1993) is a scaled Mahalanobis distance in a 3D space of L-moments (L-Cv, L-Cs, and L-Ck). A site is considered discordant if its D_i exceeds the critical value (which is 3.0 for regions with ≥15 sites, and smaller for smaller regions, as shown in the table below).

No. of Sites (N)	Critical Di	No. of Sites (N)	Critical Di
5	1.33	10	2.49
6	1.65	11	2.63
7	1.92	12	2.76
8	2.14	13	2.87
9	2.33	14	2.97
≥15	3.00

Applying the discordancy test to our 3 homogeneous regions yielded the following site-specific discordancy values and regional average L-moments:

Region	Station Name	Discordancy Di	Regional L-Moments
Region I (N = 4)	Ambala	1.00	L-Cv (τ) = 0.2237 L-Cs (τ3) = 0.1869 L-Ck (τ4) = 0.1935
	Karnal	1.00
	Jagadhari	1.00
	Kalka	1.00
Region II (N = 13)	Sirsa	0.64	L-Cv (τ) = 0.3004 L-Cs (τ3) = 0.1985 L-Ck (τ4) = 0.1724
	Hansi	1.97
	Farukhnagar	0.97
	Faridabad	0.96
	Mahendragarh	0.32
	Khol	1.12
	Palwal	2.15
	Bhiwani	0.89
	Tohana	0.94
	Sohana	0.82
	Dujana	0.22
	Salhawas	1.36
	Beri	0.64
Region III (N = 10)	Hisar	0.72	L-Cv (τ) = 0.2648 L-Cs (τ3) = 0.1446 L-Ck (τ4) = 0.1232
	Sonipat	0.77
	Rohtak	1.40
	Nuh	1.53
	Jhajjar	0.79
	Bawal	1.85
	Panipat	0.24
	Kurukshetra	1.79
	Narwana	0.54
	Kaithal	0.36

Since all computed D_i values are strictly less than their respective regional critical bounds, no stations were flagged as discordant. This confirms that the regional clustering is robust and mathematically valid.

5. Regional Distribution Selection: Z-Statistic Goodness-of-Fit

Five candidate probability distributions were evaluated for each region using L-moment ratio diagrams and the Z-statistic goodness-of-fit measure (Z^DIST). The candidate distributions were: Generalized Logistic (GLO), Generalized Extreme Value (GEV), Generalized Pareto (GPA), Generalized Normal (GNO), and Pearson Type-3 (PE3).

The goodness-of-fit measure is defined as:

Z^DIST = [ τ₄^DIST - τ_bar₄ + B₄ ] / S₄

Where τ₄^DIST is the L-kurtosis of the fitted candidate distribution, τ_bar₄ is the regional average sample L-kurtosis, B₄ is the bias, and S₄ is the standard deviation of τ_bar₄ obtained via Monte Carlo simulations (usually 500+ trials). The fit is deemed acceptable if |Z^DIST| ≤ 1.64 (at the 90% confidence level). If multiple models pass, the one with the smallest |Z^DIST| is selected as the best-fit distribution.

Region	Distribution	Z-Statistic Value	Location (ξ)	Scale (α)	Shape (k / γ)	Fit Status
Region I	GLO**	-0.10	0.9324	0.2111	-0.1869	Best Fit
	GEV*	-1.18	0.8099	0.3149	-0.0262	Satisfactory
	GNO*	-1.39	0.9254	0.3727	-0.3856	Satisfactory
Region II	GLO*	1.32	0.9038	0.2813	-0.1985	Satisfactory
	GEV**	-0.72	0.7414	0.4154	-0.0439	Best Fit
	GNO*	-1.19	0.8939	0.4964	-0.4100	Satisfactory
Region III	GEV*	1.20	0.7866	0.3958	0.0398	Satisfactory
	GNO*	1.08	0.9312	0.4523	-0.2974	Satisfactory
	PE3**	0.44	1.0000	0.4808	0.8800	Best Fit

6. Regional Growth Curves and Rainfall Quantile Estimates

Using the **index-flood method**, the regional growth curve q(F) for each homogeneous region is calculated using the best-fit distribution parameters. The site-specific rainfall quantile at return period T (non-exceedance probability F = 1 - 1/T) is then computed as:

Q_i(F) = μ_i × q(F)

Where μ_i is the site-specific mean maximum monthly rainfall (the index flood). Table 6 presents the regional quantiles (in mm) scaled for a site with the average regional mean, and Table 7 details the station-wise estimated rainfall depths for various return periods (T = 2 to 100 years).

Region	T = 2 yrs (F=0.5)	T = 5 yrs (F=0.8)	T = 10 yrs (F=0.9)	T = 20 yrs (F=0.95)	T = 50 yrs (F=0.98)	T = 100 yrs (F=0.99)	T = 200 yrs (F=0.995)
Region I (GLO)	310.38	421.61	501.54	586.46	712.67	822.24	945.79
Region II (GEV)	161.75	250.49	311.58	372.12	453.63	516.71	581.59
Region III (PE3)	207.68	306.16	367.12	422.50	490.61	539.52	586.64

Region	Station	T = 2 yrs	T = 5 yrs	T = 10 yrs	T = 20 yrs	T = 50 yrs	T = 100 yrs
Region I (GLO)	Ambala	277.53	376.99	448.45	524.38	637.24	735.21
	Karnal	259.43	352.40	419.20	490.19	595.68	687.26
	Jagadhari	345.31	469.06	557.98	652.46	792.88	914.77
	Kalka	360.29	489.41	582.19	680.77	827.28	954.47
Region II (GEV)	Sirsa	139.18	214.55	265.42	314.89	379.78	429.26
	Hansi	100.89	155.52	192.39	228.26	275.30	311.16
	Farukhnagar	173.25	267.06	330.38	391.96	472.73	534.32
	Faridabad	231.37	356.66	441.22	523.47	631.34	713.58
	Mahendragarh	152.34	234.84	290.51	344.67	415.70	469.85
	Khol	127.39	196.38	242.94	288.22	347.62	392.91
	Palwal	190.58	293.79	363.44	431.19	520.04	587.79
	Bhiwani	119.03	183.48	226.98	269.30	324.79	367.10
	Tohana	132.24	203.86	252.18	299.19	360.85	407.86
	Sohana	192.26	296.38	366.64	434.99	524.63	592.97
	Dujana	189.01	291.37	360.44	427.63	515.75	582.94
	Salhawas	151.32	233.27	288.57	342.36	412.92	466.71
	Beri	186.22	287.06	355.12	421.31	508.13	574.33
Region III (PE3)	Hisar	156.11	230.14	275.97	317.60	368.80	405.56
	Sonipat	243.18	358.50	429.88	494.73	574.49	631.75
	Rohtak	199.77	294.49	353.14	406.41	471.92	518.96
	Nuh	239.59	353.21	423.54	487.43	566.01	622.43
	Jhajjar	222.96	328.68	394.13	453.59	526.71	579.21
	Bawal	221.12	325.97	390.87	449.84	522.36	574.42
	Panipat	196.19	289.22	346.81	399.13	463.47	509.67
	Kurukshetra	209.64	309.05	370.59	426.49	495.25	544.61
	Narwana	186.02	274.23	328.84	378.44	439.45	483.25
	Kaithal	202.61	298.68	358.16	412.19	478.63	526.34

7. Code Tutorial: Implementing L-Moments in Python and R

To enable researchers to perform these calculations, we provide two ready-to-use snippets demonstrating L-moment estimation and extreme-value fitting.

R Script (Using the `lmom` Package)

# Install and load the lmom library
if (!requireNamespace("lmom", quietly = TRUE)) install.packages("lmom")
library(lmom)

# Example: Maximum monthly rainfall data for a station
rainfall_data <- c(150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250)

# 1. Compute sample L-moments (L1, L2, L3, L4)
sam_lmom <- samlmu(rainfall_data)
cat("Sample L-Moments:\n")
print(sam_lmom)

# 2. Extract L-Cv, L-Cs (t_3), L-Ck (t_4)
# Note: samlmu returns L-location, L-scale, L-skewness (t_3), L-kurtosis (t_4), etc.
l_cv <- sam_lmom[2] / sam_lmom[1]
cat("L-Cv:", l_cv, "\nL-Cs (t_3):", sam_lmom[3], "\nL-Ck (t_4):", sam_lmom[4], "\n")

# 3. Fit a Generalized Extreme Value (GEV) distribution
gev_params <- pelgev(sam_lmom)
cat("\nFitted GEV Parameters:\n")
print(gev_params)

# 4. Estimate quantiles for T = 10, 50, and 100 years
return_periods <- c(10, 50, 100)
probabilities <- 1 - 1 / return_periods
quantiles <- quagev(probabilities, gev_params)

# Display results
results <- data.frame(ReturnPeriod_Yrs = return_periods, Quantile_mm = quantiles)
print(results)

Python Script (Using the `lmoments3` Package)

import numpy as np
# Note: install via: pip install lmoments3
import lmoments3 as lm
from lmoments3 import distr

# Example: Maximum monthly rainfall data
rainfall_data = [150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250]

# 1. Compute sample L-moments and ratios
lmom_ratios = lm.lmom_ratios(rainfall_data, nmom=4)
print("Sample L-moments and Ratios:")
print(f"Mean (L1): {lmom_ratios[0]:.4f}")
print(f"L-scale (L2): {lmom_ratios[1]:.4f}")
print(f"L-skewness (t3): {lmom_ratios[2]:.4f}")
print(f"L-kurtosis (t4): {lmom_ratios[3]:.4f}")

# 2. Fit a Generalized Extreme Value (GEV) distribution
fitted_gev = distr.gev.lmom_fit(rainfall_data)
print(f"\nFitted GEV Parameters: {fitted_gev}")

# 3. Compute return level quantiles for T = 10, 50, and 100 years
return_periods = [10, 50, 100]
for T in return_periods:
    F = 1 - 1 / T
    quantile = distr.gev.ppf(F, **fitted_gev)
    print(f"T = {T:3d} years (F = {F:.2f}) -> Quantile: {quantile:.2f} mm")

8. Agricultural and Engineering Implications

The results of this regional study have critical applications for the development and policy planning of Haryana:

• Hydraulic Structures: For Region I (Wet zone, fitted to GLO), designs must accommodate larger return-period rainfall quantities, where a 100-year event can exceed 950 mm in Kalka.
• Agricultural Drainage: In Region II (Dry zone, GEV) and Region III (Central zone, PE3), drainage infrastructure must cope with 50-year rainfall events ranging from 270 mm to 570 mm depending on the exact location. Over-designing can waste valuable rural infrastructure budget, while under-designing can cause widespread waterlogging of sensitive agricultural crops, ruining seasonal yields.
• Water Harvesting: Estimating return levels helps calculate maximum design inflows for farm ponds, reservoirs, and check dams, helping farmers store surplus rainwater for dry season irrigation.

References

• Babu, V. B. and Hooda B. K. (2018). Fuzzy Majority Approach for Modeling Spatial and Temporal Distributions of Daily Rainfall in Western Zone of Haryana. International Journal of Agricultural and Statistical Sciences, 14(1), 57-67.
• Greenwood, J. A., Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form. Water Resources Research, 15(5), 1049-1054.
• Hosking, J. R. M. (1990). L-moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. Journal of the Royal Statistical Society (Series B), 52(1), 105-124.
• Hosking, J. R. M. and Wallis, J. R. (1993). Some statistics useful in regional frequency analysis. Water Resources Research, 29(2), 271-281.
• Hosking, J. R. M. and Wallis, J. R. (1997). Regional frequency analysis: An approach based on L-Moments. Cambridge University Press, United Kingdom.
• Hooda, B. K. (2006). Probability Analysis of Monthly Rainfall for Agricultural Planning At Hisar. Indian Journal of Soil Conservation, 34(1), 12-14.
• Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability-weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resources Research, 15, 1055-1064.
• Malekinezhad, H. and Garizi, A. Z. (2014). Regional frequency analysis of daily rainfall extremes using L-moments approach. Atmosfera, 27(4), 411-427.
• Majumder A., Patil S. G., Noman M. D., and Biswas S. (2015). Application of L-moments for regional frequency analysis of maximum monthly rainfall in West Bengal, India. Mausam, 66(2), 273-280.
• Nain, M. and Hooda B. K. (2019). Probability and Trend Analysis of Monthly Rainfall in Haryana. International Journal of Agricultural and Statistical Sciences, 15(1), 221-229.
• Sahrin S., Ismail N., and Alias N. E. (2018). Regional frequency analysis on peninsular Malaysia using L-moments. Far East Journal of Mathematical Sciences (FJMS), 103(8), 1379-1398.

β_r = E[X {F(X)}^r] = ∫₀¹ x(F) F^r dF

Where x(F) is the inverse cumulative distribution (quantile) function, and r is a non-negative integer. Unbiased sample estimators b_r of β_r are computed from an ordered sample X_1:n ≤ X_2:n ≤ ... ≤ X_n:n using:

b_r = n^-1 ∑_i=1ⁿ [ (^i-1_r) / (^n-1_r) ] X_i:n

L-Moments Definition

The first four L-moments (λ_r) are linear combinations of the PWMs:

• L-Location (Mean): λ₁ = β₀
• L-Scale: λ₂ = 2β₁ - β₀
• L-Skewness measure: λ₃ = 6β₂ - 6β₁ + β₀
• L-Kurtosis measure: λ₄ = 20β₃ - 30β₂ + 12β₁ - β₀

To characterize distributions independently of their scale, we define dimensionless L-moment ratios:

• L-coefficient of variation (L-Cv, τ): τ = λ₂ / λ₁
• L-coefficient of skewness (L-Cs, τ₃): τ₃ = λ₃ / λ₂
• L-coefficient of kurtosis (L-Ck, τ₄): τ₄ = λ₄ / λ₂

2. Database and Initial Data Screening

The study utilizes maximum monthly rainfall data for the 48-year period (1970–2017) obtained from the National Data Centre, Indian Meteorological Department (IMD), Pune, covering 27 rain gauge stations in Haryana.

Before executing RFA, the assumptions of stationarity, randomness, and independence must be verified for all stations:

• Stationarity: Tested using the Mann-Kendall trend test. Results showed that only 3 out of 27 sites (Karnal, Kaithal, and Bhiwani) had a statistically significant trend, meaning the regional maximum rainfall series can be treated as stationary.
• Randomness: Tested using the Run test. Except for Rohtak, the rainfall series across all other sites were random.
• Independence: Evaluated using the Autocorrelation Function (ACF). Only Rohtak and Kurukshetra showed significant autocorrelation at lag-1. Overall, it is highly reasonable to treat the data as time-independent and suitable for regional frequency analysis.

Station Name	Mann-Kendall Trend (Tau)	MK P-value	Interpretation	No. of Runs	Run P-value
Sirsa	-0.024	0.810	No Trend	24	0.771
Narwana	-0.108	0.282	No Trend	20	0.145
Hisar	-0.126	0.210	No Trend	24	0.770
Karnal	0.285	0.005	Trend	18	0.054
Ambala	-0.162	0.106	No Trend	22	0.381
Jhajjar	-0.155	0.126	No Trend	20	0.233
Hansi	-0.143	0.060	No Trend	25	1.000
Sonipat	-0.160	0.112	No Trend	22	0.243
Rohtak	-0.109	0.074	No Trend	16	0.008 (Not Random)
Panipat	-0.121	0.230	No Trend	19	0.080
Farukhnagar	-0.114	0.259	No Trend	20	0.145
Faridabad	0.101	0.315	No Trend	22	0.381
Kurukshetra	-0.005	0.810	No Trend	18	0.074
Mahendragarh	-0.141	0.160	No Trend	26	0.770
Kaithal	-0.296	0.003	Trend	21	0.243
Khol	0.029	0.776	No Trend	24	0.780
Palwal	-0.164	0.110	No Trend	15	0.710
Bhiwani	-0.246	0.017	Trend	30	0.074
Tohana	-0.078	0.439	No Trend	25	1.000
Sohana	-0.153	0.129	No Trend	26	0.770
Bawal	-0.025	0.810	No Trend	20	0.145
Jagadhari	0.114	0.255	No Trend	21	0.243
Dujana	-0.041	0.693	No Trend	18	0.074
Salhawas	-0.176	0.080	No Trend	22	0.381
Nuh	-0.031	0.763	No Trend	20	0.145
Kalka	-0.196	0.053	No Trend	19	0.136
Beri	-0.012	0.915	No Trend	22	0.381

3. L-Moments and Station-Wise Characteristics

For each of the 27 sites, sample L-moments and L-moment ratios were computed. The values represent the mean maximum monthly rainfall (in mm), L-coefficient of variation (L-Cv), L-skewness (L-Cs), and L-kurtosis (L-Ck).

Station	Mean (mm)	L-Cv (τ)	L-Cs (τ3)	L-Ck (τ4)
Sirsa	154.33	0.270	0.153	0.177
Narwana	200.021	0.290	0.153	0.084
Hisar	167.863	0.249	0.115	0.176
Karnal	278.356	0.224	0.097	0.071
Ambala	297.777	0.185	0.208	0.250
Jhajjar	239.739	0.271	0.130	0.047
Hansi	111.727	0.312	0.268	0.137
Sonipat	261.487	0.254	0.166	0.132
Rohtak	214.803	0.276	0.100	0.076
Panipat	210.958	0.262	0.143	0.099
Farukhnagar	191.856	0.344	0.197	0.217
Faridabad	256.224	0.244	0.205	0.216
Kurukshetra	225.420	0.297	0.222	0.116
Mahendragarh	168.708	0.303	0.236	0.238
Kaithal	217.858	0.274	0.141	0.067
Khol	141.079	0.356	0.169	0.161
Palwal	211.058	0.257	0.318	0.309
Bhiwani	145.872	0.270	0.172	0.096
Tohana	151.323	0.283	0.115	0.140
Sohana	225.158	0.261	0.149	0.100
Bawal	237.758	0.213	0.112	0.217
Jagadhari	370.504	0.251	0.289	0.241
Dujana	209.315	0.308	0.161	0.132
Salhawas	167.579	0.347	0.295	0.263
Nuh	257.627	0.264	0.167	0.215
Kalka	386.581	0.234	0.153	0.212
Beri	206.223	0.322	0.183	0.109

4. Formation and Validation of Homogeneous Regions

To define homogeneous regions, the mean monthly rainfall values were subjected to hierarchical cluster analysis (Ward's Method). The Elbow Method (analyzing the within-cluster sum of squares) indicated that the optimal number of regions is three.

• Region I (Wet/Semi-humid zone - 4 stations): Ambala, Karnal, Jagadhari, and Kalka.
• Region II (Dry/Semi-arid zone - 13 stations): Sirsa, Hansi, Farukhnagar, Faridabad, Mahendragarh, Khol, Palwal, Bhiwani, Tohana, Sohana, Dujana, Salhawas, and Beri.
• Region III (Central/Transition zone - 10 stations): Hisar, Sonipat, Rohtak, Nuh, Jhajjar, Bawal, Panipat, Kurukshetra, Narwana, and Kaithal.

Discordancy Test (D_i)

The discordancy measure D_i (Hosking and Wallis, 1993) is a scaled Mahalanobis distance in a 3D space of L-moments (L-Cv, L-Cs, and L-Ck). A site is considered discordant if its D_i exceeds the critical value (which is 3.0 for regions with ≥15 sites, and smaller for smaller regions, as shown in the table below).

No. of Sites (N)	Critical Di	No. of Sites (N)	Critical Di
5	1.33	10	2.49
6	1.65	11	2.63
7	1.92	12	2.76
8	2.14	13	2.87
9	2.33	14	2.97
≥15	3.00

Applying the discordancy test to our 3 homogeneous regions yielded the following site-specific discordancy values and regional average L-moments:

Region	Station Name	Discordancy Di	Regional L-Moments
Region I (N = 4)	Ambala	1.00	L-Cv (τ) = 0.2237 L-Cs (τ3) = 0.1869 L-Ck (τ4) = 0.1935
	Karnal	1.00
	Jagadhari	1.00
	Kalka	1.00
Region II (N = 13)	Sirsa	0.64	L-Cv (τ) = 0.3004 L-Cs (τ3) = 0.1985 L-Ck (τ4) = 0.1724
	Hansi	1.97
	Farukhnagar	0.97
	Faridabad	0.96
	Mahendragarh	0.32
	Khol	1.12
	Palwal	2.15
	Bhiwani	0.89
	Tohana	0.94
	Sohana	0.82
	Dujana	0.22
	Salhawas	1.36
	Beri	0.64
Region III (N = 10)	Hisar	0.72	L-Cv (τ) = 0.2648 L-Cs (τ3) = 0.1446 L-Ck (τ4) = 0.1232
	Sonipat	0.77
	Rohtak	1.40
	Nuh	1.53
	Jhajjar	0.79
	Bawal	1.85
	Panipat	0.24
	Kurukshetra	1.79
	Narwana	0.54
	Kaithal	0.36

Since all computed D_i values are strictly less than their respective regional critical bounds, no stations were flagged as discordant. This confirms that the regional clustering is robust and mathematically valid.

5. Regional Distribution Selection: Z-Statistic Goodness-of-Fit

Five candidate probability distributions were evaluated for each region using L-moment ratio diagrams and the Z-statistic goodness-of-fit measure (Z^DIST). The candidate distributions were: Generalized Logistic (GLO), Generalized Extreme Value (GEV), Generalized Pareto (GPA), Generalized Normal (GNO), and Pearson Type-3 (PE3).

The goodness-of-fit measure is defined as:

Z^DIST = [ τ₄^DIST - τ_bar₄ + B₄ ] / S₄

Where τ₄^DIST is the L-kurtosis of the fitted candidate distribution, τ_bar₄ is the regional average sample L-kurtosis, B₄ is the bias, and S₄ is the standard deviation of τ_bar₄ obtained via Monte Carlo simulations (usually 500+ trials). The fit is deemed acceptable if |Z^DIST| ≤ 1.64 (at the 90% confidence level). If multiple models pass, the one with the smallest |Z^DIST| is selected as the best-fit distribution.

Region	Distribution	Z-Statistic Value	Location (ξ)	Scale (α)	Shape (k / γ)	Fit Status
Region I	GLO**	-0.10	0.9324	0.2111	-0.1869	Best Fit
	GEV*	-1.18	0.8099	0.3149	-0.0262	Satisfactory
	GNO*	-1.39	0.9254	0.3727	-0.3856	Satisfactory
Region II	GLO*	1.32	0.9038	0.2813	-0.1985	Satisfactory
	GEV**	-0.72	0.7414	0.4154	-0.0439	Best Fit
	GNO*	-1.19	0.8939	0.4964	-0.4100	Satisfactory
Region III	GEV*	1.20	0.7866	0.3958	0.0398	Satisfactory
	GNO*	1.08	0.9312	0.4523	-0.2974	Satisfactory
	PE3**	0.44	1.0000	0.4808	0.8800	Best Fit

6. Regional Growth Curves and Rainfall Quantile Estimates

Using the **index-flood method**, the regional growth curve q(F) for each homogeneous region is calculated using the best-fit distribution parameters. The site-specific rainfall quantile at return period T (non-exceedance probability F = 1 - 1/T) is then computed as:

Q_i(F) = μ_i × q(F)

Where μ_i is the site-specific mean maximum monthly rainfall (the index flood). Table 6 presents the regional quantiles (in mm) scaled for a site with the average regional mean, and Table 7 details the station-wise estimated rainfall depths for various return periods (T = 2 to 100 years).

Region	T = 2 yrs (F=0.5)	T = 5 yrs (F=0.8)	T = 10 yrs (F=0.9)	T = 20 yrs (F=0.95)	T = 50 yrs (F=0.98)	T = 100 yrs (F=0.99)	T = 200 yrs (F=0.995)
Region I (GLO)	310.38	421.61	501.54	586.46	712.67	822.24	945.79
Region II (GEV)	161.75	250.49	311.58	372.12	453.63	516.71	581.59
Region III (PE3)	207.68	306.16	367.12	422.50	490.61	539.52	586.64

Region	Station	T = 2 yrs	T = 5 yrs	T = 10 yrs	T = 20 yrs	T = 50 yrs	T = 100 yrs
Region I (GLO)	Ambala	277.53	376.99	448.45	524.38	637.24	735.21
	Karnal	259.43	352.40	419.20	490.19	595.68	687.26
	Jagadhari	345.31	469.06	557.98	652.46	792.88	914.77
	Kalka	360.29	489.41	582.19	680.77	827.28	954.47
Region II (GEV)	Sirsa	139.18	214.55	265.42	314.89	379.78	429.26
	Hansi	100.89	155.52	192.39	228.26	275.30	311.16
	Farukhnagar	173.25	267.06	330.38	391.96	472.73	534.32
	Faridabad	231.37	356.66	441.22	523.47	631.34	713.58
	Mahendragarh	152.34	234.84	290.51	344.67	415.70	469.85
	Khol	127.39	196.38	242.94	288.22	347.62	392.91
	Palwal	190.58	293.79	363.44	431.19	520.04	587.79
	Bhiwani	119.03	183.48	226.98	269.30	324.79	367.10
	Tohana	132.24	203.86	252.18	299.19	360.85	407.86
	Sohana	192.26	296.38	366.64	434.99	524.63	592.97
	Dujana	189.01	291.37	360.44	427.63	515.75	582.94
	Salhawas	151.32	233.27	288.57	342.36	412.92	466.71
	Beri	186.22	287.06	355.12	421.31	508.13	574.33
Region III (PE3)	Hisar	156.11	230.14	275.97	317.60	368.80	405.56
	Sonipat	243.18	358.50	429.88	494.73	574.49	631.75
	Rohtak	199.77	294.49	353.14	406.41	471.92	518.96
	Nuh	239.59	353.21	423.54	487.43	566.01	622.43
	Jhajjar	222.96	328.68	394.13	453.59	526.71	579.21
	Bawal	221.12	325.97	390.87	449.84	522.36	574.42
	Panipat	196.19	289.22	346.81	399.13	463.47	509.67
	Kurukshetra	209.64	309.05	370.59	426.49	495.25	544.61
	Narwana	186.02	274.23	328.84	378.44	439.45	483.25
	Kaithal	202.61	298.68	358.16	412.19	478.63	526.34

7. Code Tutorial: Implementing L-Moments in Python and R

To enable researchers to perform these calculations, we provide two ready-to-use snippets demonstrating L-moment estimation and extreme-value fitting.

R Script (Using the `lmom` Package)

# Install and load the lmom library
if (!requireNamespace("lmom", quietly = TRUE)) install.packages("lmom")
library(lmom)

# Example: Maximum monthly rainfall data for a station
rainfall_data <- c(150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250)

# 1. Compute sample L-moments (L1, L2, L3, L4)
sam_lmom <- samlmu(rainfall_data)
cat("Sample L-Moments:\n")
print(sam_lmom)

# 2. Extract L-Cv, L-Cs (t_3), L-Ck (t_4)
# Note: samlmu returns L-location, L-scale, L-skewness (t_3), L-kurtosis (t_4), etc.
l_cv <- sam_lmom[2] / sam_lmom[1]
cat("L-Cv:", l_cv, "\nL-Cs (t_3):", sam_lmom[3], "\nL-Ck (t_4):", sam_lmom[4], "\n")

# 3. Fit a Generalized Extreme Value (GEV) distribution
gev_params <- pelgev(sam_lmom)
cat("\nFitted GEV Parameters:\n")
print(gev_params)

# 4. Estimate quantiles for T = 10, 50, and 100 years
return_periods <- c(10, 50, 100)
probabilities <- 1 - 1 / return_periods
quantiles <- quagev(probabilities, gev_params)

# Display results
results <- data.frame(ReturnPeriod_Yrs = return_periods, Quantile_mm = quantiles)
print(results)

Python Script (Using the `lmoments3` Package)

import numpy as np
# Note: install via: pip install lmoments3
import lmoments3 as lm
from lmoments3 import distr

# Example: Maximum monthly rainfall data
rainfall_data = [150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250]

# 1. Compute sample L-moments and ratios
lmom_ratios = lm.lmom_ratios(rainfall_data, nmom=4)
print("Sample L-moments and Ratios:")
print(f"Mean (L1): {lmom_ratios[0]:.4f}")
print(f"L-scale (L2): {lmom_ratios[1]:.4f}")
print(f"L-skewness (t3): {lmom_ratios[2]:.4f}")
print(f"L-kurtosis (t4): {lmom_ratios[3]:.4f}")

# 2. Fit a Generalized Extreme Value (GEV) distribution
fitted_gev = distr.gev.lmom_fit(rainfall_data)
print(f"\nFitted GEV Parameters: {fitted_gev}")

# 3. Compute return level quantiles for T = 10, 50, and 100 years
return_periods = [10, 50, 100]
for T in return_periods:
    F = 1 - 1 / T
    quantile = distr.gev.ppf(F, **fitted_gev)
    print(f"T = {T:3d} years (F = {F:.2f}) -> Quantile: {quantile:.2f} mm")

8. Agricultural and Engineering Implications

The results of this regional study have critical applications for the development and policy planning of Haryana:

• Hydraulic Structures: For Region I (Wet zone, fitted to GLO), designs must accommodate larger return-period rainfall quantities, where a 100-year event can exceed 950 mm in Kalka.
• Agricultural Drainage: In Region II (Dry zone, GEV) and Region III (Central zone, PE3), drainage infrastructure must cope with 50-year rainfall events ranging from 270 mm to 570 mm depending on the exact location. Over-designing can waste valuable rural infrastructure budget, while under-designing can cause widespread waterlogging of sensitive agricultural crops, ruining seasonal yields.
• Water Harvesting: Estimating return levels helps calculate maximum design inflows for farm ponds, reservoirs, and check dams, helping farmers store surplus rainwater for dry season irrigation.

References

• Babu, V. B. and Hooda B. K. (2018). Fuzzy Majority Approach for Modeling Spatial and Temporal Distributions of Daily Rainfall in Western Zone of Haryana. International Journal of Agricultural and Statistical Sciences, 14(1), 57-67.
• Greenwood, J. A., Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form. Water Resources Research, 15(5), 1049-1054.
• Hosking, J. R. M. (1990). L-moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. Journal of the Royal Statistical Society (Series B), 52(1), 105-124.
• Hosking, J. R. M. and Wallis, J. R. (1993). Some statistics useful in regional frequency analysis. Water Resources Research, 29(2), 271-281.
• Hosking, J. R. M. and Wallis, J. R. (1997). Regional frequency analysis: An approach based on L-Moments. Cambridge University Press, United Kingdom.
• Hooda, B. K. (2006). Probability Analysis of Monthly Rainfall for Agricultural Planning At Hisar. Indian Journal of Soil Conservation, 34(1), 12-14.
• Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability-weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resources Research, 15, 1055-1064.
• Malekinezhad, H. and Garizi, A. Z. (2014). Regional frequency analysis of daily rainfall extremes using L-moments approach. Atmosfera, 27(4), 411-427.
• Majumder A., Patil S. G., Noman M. D., and Biswas S. (2015). Application of L-moments for regional frequency analysis of maximum monthly rainfall in West Bengal, India. Mausam, 66(2), 273-280.
• Nain, M. and Hooda B. K. (2019). Probability and Trend Analysis of Monthly Rainfall in Haryana. International Journal of Agricultural and Statistical Sciences, 15(1), 221-229.
• Sahrin S., Ismail N., and Alias N. E. (2018). Regional frequency analysis on peninsular Malaysia using L-moments. Far East Journal of Mathematical Sciences (FJMS), 103(8), 1379-1398.

Extreme meteorological events, particularly intense rainfall, pose significant risks to human life, agriculture, and infrastructure. In states like Haryana, which is heavily reliant on agriculture, predicting the return periods of extreme rainfall is essential for designing resilient drainage systems, dams, highways, and bridges, as well as planning effective water resource management strategies. Without robust statistical modeling, infrastructure design can either be dangerously under-engineered (leading to failures) or uneconomically over-engineered.

In hydrological research, at-site frequency analysis often suffers from high sampling variability and instability, particularly when estimating return periods that exceed the available record lengths. To overcome these limitations, Regional Frequency Analysis (RFA) pools data from multiple hydrologically and climatologically homogeneous stations. This study reviews and applies an advanced L-moment-based RFA framework to 48 years (1970–2017) of maximum monthly rainfall data across 27 rain gauge stations in Haryana, India.

1. Core Methodology: The Power of L-Moments

Introduced by Hosking (1990), L-moments are linear combinations of ordered sample values. Unlike conventional moments (which square or cube data values, leading to extreme sensitivity to outliers and sample size bias), L-moments use linear weights. This makes them:

• Highly robust to outliers and extreme values.
• Unbiased and stable even for small sample sizes.
• Exceptionally reliable for selecting and fitting extreme value distributions.

Probability Weighted Moments (PWMs)

L-moments are derived from Probability Weighted Moments (PWMs), defined by Greenwood et al. (1979). For a random variable X with cumulative distribution function F(x), the r-th PWM (denoted as β_r) is given by:

β_r = E[X {F(X)}^r] = ∫₀¹ x(F) F^r dF

Where x(F) is the inverse cumulative distribution (quantile) function, and r is a non-negative integer. Unbiased sample estimators b_r of β_r are computed from an ordered sample X_1:n ≤ X_2:n ≤ ... ≤ X_n:n using:

b_r = n^-1 ∑_i=1ⁿ [ (^i-1_r) / (^n-1_r) ] X_i:n

L-Moments Definition

The first four L-moments (λ_r) are linear combinations of the PWMs:

• L-Location (Mean): λ₁ = β₀
• L-Scale: λ₂ = 2β₁ - β₀
• L-Skewness measure: λ₃ = 6β₂ - 6β₁ + β₀
• L-Kurtosis measure: λ₄ = 20β₃ - 30β₂ + 12β₁ - β₀

To characterize distributions independently of their scale, we define dimensionless L-moment ratios:

• L-coefficient of variation (L-Cv, τ): τ = λ₂ / λ₁
• L-coefficient of skewness (L-Cs, τ₃): τ₃ = λ₃ / λ₂
• L-coefficient of kurtosis (L-Ck, τ₄): τ₄ = λ₄ / λ₂

2. Database and Initial Data Screening

The study utilizes maximum monthly rainfall data for the 48-year period (1970–2017) obtained from the National Data Centre, Indian Meteorological Department (IMD), Pune, covering 27 rain gauge stations in Haryana.

Before executing RFA, the assumptions of stationarity, randomness, and independence must be verified for all stations:

• Stationarity: Tested using the Mann-Kendall trend test. Results showed that only 3 out of 27 sites (Karnal, Kaithal, and Bhiwani) had a statistically significant trend, meaning the regional maximum rainfall series can be treated as stationary.
• Randomness: Tested using the Run test. Except for Rohtak, the rainfall series across all other sites were random.
• Independence: Evaluated using the Autocorrelation Function (ACF). Only Rohtak and Kurukshetra showed significant autocorrelation at lag-1. Overall, it is highly reasonable to treat the data as time-independent and suitable for regional frequency analysis.

Station Name	Mann-Kendall Trend (Tau)	MK P-value	Interpretation	No. of Runs	Run P-value
Sirsa	-0.024	0.810	No Trend	24	0.771
Narwana	-0.108	0.282	No Trend	20	0.145
Hisar	-0.126	0.210	No Trend	24	0.770
Karnal	0.285	0.005	Trend	18	0.054
Ambala	-0.162	0.106	No Trend	22	0.381
Jhajjar	-0.155	0.126	No Trend	20	0.233
Hansi	-0.143	0.060	No Trend	25	1.000
Sonipat	-0.160	0.112	No Trend	22	0.243
Rohtak	-0.109	0.074	No Trend	16	0.008 (Not Random)
Panipat	-0.121	0.230	No Trend	19	0.080
Farukhnagar	-0.114	0.259	No Trend	20	0.145
Faridabad	0.101	0.315	No Trend	22	0.381
Kurukshetra	-0.005	0.810	No Trend	18	0.074
Mahendragarh	-0.141	0.160	No Trend	26	0.770
Kaithal	-0.296	0.003	Trend	21	0.243
Khol	0.029	0.776	No Trend	24	0.780
Palwal	-0.164	0.110	No Trend	15	0.710
Bhiwani	-0.246	0.017	Trend	30	0.074
Tohana	-0.078	0.439	No Trend	25	1.000
Sohana	-0.153	0.129	No Trend	26	0.770
Bawal	-0.025	0.810	No Trend	20	0.145
Jagadhari	0.114	0.255	No Trend	21	0.243
Dujana	-0.041	0.693	No Trend	18	0.074
Salhawas	-0.176	0.080	No Trend	22	0.381
Nuh	-0.031	0.763	No Trend	20	0.145
Kalka	-0.196	0.053	No Trend	19	0.136
Beri	-0.012	0.915	No Trend	22	0.381

3. L-Moments and Station-Wise Characteristics

For each of the 27 sites, sample L-moments and L-moment ratios were computed. The values represent the mean maximum monthly rainfall (in mm), L-coefficient of variation (L-Cv), L-skewness (L-Cs), and L-kurtosis (L-Ck).

Station	Mean (mm)	L-Cv (τ)	L-Cs (τ3)	L-Ck (τ4)
Sirsa	154.33	0.270	0.153	0.177
Narwana	200.021	0.290	0.153	0.084
Hisar	167.863	0.249	0.115	0.176
Karnal	278.356	0.224	0.097	0.071
Ambala	297.777	0.185	0.208	0.250
Jhajjar	239.739	0.271	0.130	0.047
Hansi	111.727	0.312	0.268	0.137
Sonipat	261.487	0.254	0.166	0.132
Rohtak	214.803	0.276	0.100	0.076
Panipat	210.958	0.262	0.143	0.099
Farukhnagar	191.856	0.344	0.197	0.217
Faridabad	256.224	0.244	0.205	0.216
Kurukshetra	225.420	0.297	0.222	0.116
Mahendragarh	168.708	0.303	0.236	0.238
Kaithal	217.858	0.274	0.141	0.067
Khol	141.079	0.356	0.169	0.161
Palwal	211.058	0.257	0.318	0.309
Bhiwani	145.872	0.270	0.172	0.096
Tohana	151.323	0.283	0.115	0.140
Sohana	225.158	0.261	0.149	0.100
Bawal	237.758	0.213	0.112	0.217
Jagadhari	370.504	0.251	0.289	0.241
Dujana	209.315	0.308	0.161	0.132
Salhawas	167.579	0.347	0.295	0.263
Nuh	257.627	0.264	0.167	0.215
Kalka	386.581	0.234	0.153	0.212
Beri	206.223	0.322	0.183	0.109

4. Formation and Validation of Homogeneous Regions

To define homogeneous regions, the mean monthly rainfall values were subjected to hierarchical cluster analysis (Ward's Method). The Elbow Method (analyzing the within-cluster sum of squares) indicated that the optimal number of regions is three.

• Region I (Wet/Semi-humid zone - 4 stations): Ambala, Karnal, Jagadhari, and Kalka.
• Region II (Dry/Semi-arid zone - 13 stations): Sirsa, Hansi, Farukhnagar, Faridabad, Mahendragarh, Khol, Palwal, Bhiwani, Tohana, Sohana, Dujana, Salhawas, and Beri.
• Region III (Central/Transition zone - 10 stations): Hisar, Sonipat, Rohtak, Nuh, Jhajjar, Bawal, Panipat, Kurukshetra, Narwana, and Kaithal.

Discordancy Test (D_i)

The discordancy measure D_i (Hosking and Wallis, 1993) is a scaled Mahalanobis distance in a 3D space of L-moments (L-Cv, L-Cs, and L-Ck). A site is considered discordant if its D_i exceeds the critical value (which is 3.0 for regions with ≥15 sites, and smaller for smaller regions, as shown in the table below).

No. of Sites (N)	Critical Di	No. of Sites (N)	Critical Di
5	1.33	10	2.49
6	1.65	11	2.63
7	1.92	12	2.76
8	2.14	13	2.87
9	2.33	14	2.97
≥15	3.00

Applying the discordancy test to our 3 homogeneous regions yielded the following site-specific discordancy values and regional average L-moments:

Region	Station Name	Discordancy Di	Regional L-Moments
Region I (N = 4)	Ambala	1.00	L-Cv (τ) = 0.2237 L-Cs (τ3) = 0.1869 L-Ck (τ4) = 0.1935
	Karnal	1.00
	Jagadhari	1.00
	Kalka	1.00
Region II (N = 13)	Sirsa	0.64	L-Cv (τ) = 0.3004 L-Cs (τ3) = 0.1985 L-Ck (τ4) = 0.1724
	Hansi	1.97
	Farukhnagar	0.97
	Faridabad	0.96
	Mahendragarh	0.32
	Khol	1.12
	Palwal	2.15
	Bhiwani	0.89
	Tohana	0.94
	Sohana	0.82
	Dujana	0.22
	Salhawas	1.36
	Beri	0.64
Region III (N = 10)	Hisar	0.72	L-Cv (τ) = 0.2648 L-Cs (τ3) = 0.1446 L-Ck (τ4) = 0.1232
	Sonipat	0.77
	Rohtak	1.40
	Nuh	1.53
	Jhajjar	0.79
	Bawal	1.85
	Panipat	0.24
	Kurukshetra	1.79
	Narwana	0.54
	Kaithal	0.36

Since all computed D_i values are strictly less than their respective regional critical bounds, no stations were flagged as discordant. This confirms that the regional clustering is robust and mathematically valid.

5. Regional Distribution Selection: Z-Statistic Goodness-of-Fit

Five candidate probability distributions were evaluated for each region using L-moment ratio diagrams and the Z-statistic goodness-of-fit measure (Z^DIST). The candidate distributions were: Generalized Logistic (GLO), Generalized Extreme Value (GEV), Generalized Pareto (GPA), Generalized Normal (GNO), and Pearson Type-3 (PE3).

The goodness-of-fit measure is defined as:

Z^DIST = [ τ₄^DIST - τ_bar₄ + B₄ ] / S₄

Where τ₄^DIST is the L-kurtosis of the fitted candidate distribution, τ_bar₄ is the regional average sample L-kurtosis, B₄ is the bias, and S₄ is the standard deviation of τ_bar₄ obtained via Monte Carlo simulations (usually 500+ trials). The fit is deemed acceptable if |Z^DIST| ≤ 1.64 (at the 90% confidence level). If multiple models pass, the one with the smallest |Z^DIST| is selected as the best-fit distribution.

Region	Distribution	Z-Statistic Value	Location (ξ)	Scale (α)	Shape (k / γ)	Fit Status
Region I	GLO**	-0.10	0.9324	0.2111	-0.1869	Best Fit
	GEV*	-1.18	0.8099	0.3149	-0.0262	Satisfactory
	GNO*	-1.39	0.9254	0.3727	-0.3856	Satisfactory
Region II	GLO*	1.32	0.9038	0.2813	-0.1985	Satisfactory
	GEV**	-0.72	0.7414	0.4154	-0.0439	Best Fit
	GNO*	-1.19	0.8939	0.4964	-0.4100	Satisfactory
Region III	GEV*	1.20	0.7866	0.3958	0.0398	Satisfactory
	GNO*	1.08	0.9312	0.4523	-0.2974	Satisfactory
	PE3**	0.44	1.0000	0.4808	0.8800	Best Fit

6. Regional Growth Curves and Rainfall Quantile Estimates

Using the **index-flood method**, the regional growth curve q(F) for each homogeneous region is calculated using the best-fit distribution parameters. The site-specific rainfall quantile at return period T (non-exceedance probability F = 1 - 1/T) is then computed as:

Q_i(F) = μ_i × q(F)

Where μ_i is the site-specific mean maximum monthly rainfall (the index flood). Table 6 presents the regional quantiles (in mm) scaled for a site with the average regional mean, and Table 7 details the station-wise estimated rainfall depths for various return periods (T = 2 to 100 years).

Region	T = 2 yrs (F=0.5)	T = 5 yrs (F=0.8)	T = 10 yrs (F=0.9)	T = 20 yrs (F=0.95)	T = 50 yrs (F=0.98)	T = 100 yrs (F=0.99)	T = 200 yrs (F=0.995)
Region I (GLO)	310.38	421.61	501.54	586.46	712.67	822.24	945.79
Region II (GEV)	161.75	250.49	311.58	372.12	453.63	516.71	581.59
Region III (PE3)	207.68	306.16	367.12	422.50	490.61	539.52	586.64

Region	Station	T = 2 yrs	T = 5 yrs	T = 10 yrs	T = 20 yrs	T = 50 yrs	T = 100 yrs
Region I (GLO)	Ambala	277.53	376.99	448.45	524.38	637.24	735.21
	Karnal	259.43	352.40	419.20	490.19	595.68	687.26
	Jagadhari	345.31	469.06	557.98	652.46	792.88	914.77
	Kalka	360.29	489.41	582.19	680.77	827.28	954.47
Region II (GEV)	Sirsa	139.18	214.55	265.42	314.89	379.78	429.26
	Hansi	100.89	155.52	192.39	228.26	275.30	311.16
	Farukhnagar	173.25	267.06	330.38	391.96	472.73	534.32
	Faridabad	231.37	356.66	441.22	523.47	631.34	713.58
	Mahendragarh	152.34	234.84	290.51	344.67	415.70	469.85
	Khol	127.39	196.38	242.94	288.22	347.62	392.91
	Palwal	190.58	293.79	363.44	431.19	520.04	587.79
	Bhiwani	119.03	183.48	226.98	269.30	324.79	367.10
	Tohana	132.24	203.86	252.18	299.19	360.85	407.86
	Sohana	192.26	296.38	366.64	434.99	524.63	592.97
	Dujana	189.01	291.37	360.44	427.63	515.75	582.94
	Salhawas	151.32	233.27	288.57	342.36	412.92	466.71
	Beri	186.22	287.06	355.12	421.31	508.13	574.33
Region III (PE3)	Hisar	156.11	230.14	275.97	317.60	368.80	405.56
	Sonipat	243.18	358.50	429.88	494.73	574.49	631.75
	Rohtak	199.77	294.49	353.14	406.41	471.92	518.96
	Nuh	239.59	353.21	423.54	487.43	566.01	622.43
	Jhajjar	222.96	328.68	394.13	453.59	526.71	579.21
	Bawal	221.12	325.97	390.87	449.84	522.36	574.42
	Panipat	196.19	289.22	346.81	399.13	463.47	509.67
	Kurukshetra	209.64	309.05	370.59	426.49	495.25	544.61
	Narwana	186.02	274.23	328.84	378.44	439.45	483.25
	Kaithal	202.61	298.68	358.16	412.19	478.63	526.34

7. Code Tutorial: Implementing L-Moments in Python and R

To enable researchers to perform these calculations, we provide two ready-to-use snippets demonstrating L-moment estimation and extreme-value fitting.

R Script (Using the `lmom` Package)

# Install and load the lmom library
if (!requireNamespace("lmom", quietly = TRUE)) install.packages("lmom")
library(lmom)

# Example: Maximum monthly rainfall data for a station
rainfall_data <- c(150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250)

# 1. Compute sample L-moments (L1, L2, L3, L4)
sam_lmom <- samlmu(rainfall_data)
cat("Sample L-Moments:\n")
print(sam_lmom)

# 2. Extract L-Cv, L-Cs (t_3), L-Ck (t_4)
# Note: samlmu returns L-location, L-scale, L-skewness (t_3), L-kurtosis (t_4), etc.
l_cv <- sam_lmom[2] / sam_lmom[1]
cat("L-Cv:", l_cv, "\nL-Cs (t_3):", sam_lmom[3], "\nL-Ck (t_4):", sam_lmom[4], "\n")

# 3. Fit a Generalized Extreme Value (GEV) distribution
gev_params <- pelgev(sam_lmom)
cat("\nFitted GEV Parameters:\n")
print(gev_params)

# 4. Estimate quantiles for T = 10, 50, and 100 years
return_periods <- c(10, 50, 100)
probabilities <- 1 - 1 / return_periods
quantiles <- quagev(probabilities, gev_params)

# Display results
results <- data.frame(ReturnPeriod_Yrs = return_periods, Quantile_mm = quantiles)
print(results)

Python Script (Using the `lmoments3` Package)

import numpy as np
# Note: install via: pip install lmoments3
import lmoments3 as lm
from lmoments3 import distr

# Example: Maximum monthly rainfall data
rainfall_data = [150, 220, 180, 290, 110, 310, 420, 95, 130, 210, 175, 250]

# 1. Compute sample L-moments and ratios
lmom_ratios = lm.lmom_ratios(rainfall_data, nmom=4)
print("Sample L-moments and Ratios:")
print(f"Mean (L1): {lmom_ratios[0]:.4f}")
print(f"L-scale (L2): {lmom_ratios[1]:.4f}")
print(f"L-skewness (t3): {lmom_ratios[2]:.4f}")
print(f"L-kurtosis (t4): {lmom_ratios[3]:.4f}")

# 2. Fit a Generalized Extreme Value (GEV) distribution
fitted_gev = distr.gev.lmom_fit(rainfall_data)
print(f"\nFitted GEV Parameters: {fitted_gev}")

# 3. Compute return level quantiles for T = 10, 50, and 100 years
return_periods = [10, 50, 100]
for T in return_periods:
    F = 1 - 1 / T
    quantile = distr.gev.ppf(F, **fitted_gev)
    print(f"T = {T:3d} years (F = {F:.2f}) -> Quantile: {quantile:.2f} mm")

8. Agricultural and Engineering Implications

The results of this regional study have critical applications for the development and policy planning of Haryana:

• Hydraulic Structures: For Region I (Wet zone, fitted to GLO), designs must accommodate larger return-period rainfall quantities, where a 100-year event can exceed 950 mm in Kalka.
• Agricultural Drainage: In Region II (Dry zone, GEV) and Region III (Central zone, PE3), drainage infrastructure must cope with 50-year rainfall events ranging from 270 mm to 570 mm depending on the exact location. Over-designing can waste valuable rural infrastructure budget, while under-designing can cause widespread waterlogging of sensitive agricultural crops, ruining seasonal yields.
• Water Harvesting: Estimating return levels helps calculate maximum design inflows for farm ponds, reservoirs, and check dams, helping farmers store surplus rainwater for dry season irrigation.

References

• Babu, V. B. and Hooda B. K. (2018). Fuzzy Majority Approach for Modeling Spatial and Temporal Distributions of Daily Rainfall in Western Zone of Haryana. International Journal of Agricultural and Statistical Sciences, 14(1), 57-67.
• Greenwood, J. A., Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability weighted moments: Definition and relation to parameters of several distributions expressible in inverse form. Water Resources Research, 15(5), 1049-1054.
• Hosking, J. R. M. (1990). L-moments: Analysis and Estimation of Distributions Using Linear Combinations of Order Statistics. Journal of the Royal Statistical Society (Series B), 52(1), 105-124.
• Hosking, J. R. M. and Wallis, J. R. (1993). Some statistics useful in regional frequency analysis. Water Resources Research, 29(2), 271-281.
• Hosking, J. R. M. and Wallis, J. R. (1997). Regional frequency analysis: An approach based on L-Moments. Cambridge University Press, United Kingdom.
• Hooda, B. K. (2006). Probability Analysis of Monthly Rainfall for Agricultural Planning At Hisar. Indian Journal of Soil Conservation, 34(1), 12-14.
• Landwehr, J. M., Matalas, N. C., and Wallis, J. R. (1979). Probability-weighted moments compared with some traditional techniques in estimating Gumbel parameters and quantiles. Water Resources Research, 15, 1055-1064.
• Malekinezhad, H. and Garizi, A. Z. (2014). Regional frequency analysis of daily rainfall extremes using L-moments approach. Atmosfera, 27(4), 411-427.
• Majumder A., Patil S. G., Noman M. D., and Biswas S. (2015). Application of L-moments for regional frequency analysis of maximum monthly rainfall in West Bengal, India. Mausam, 66(2), 273-280.
• Nain, M. and Hooda B. K. (2019). Probability and Trend Analysis of Monthly Rainfall in Haryana. International Journal of Agricultural and Statistical Sciences, 15(1), 221-229.
• Sahrin S., Ismail N., and Alias N. E. (2018). Regional frequency analysis on peninsular Malaysia using L-moments. Far East Journal of Mathematical Sciences (FJMS), 103(8), 1379-1398.