Abstract

Peninsular Malaysia is highly flood-prone. Conventional L-moment regional frequency analysis (RFA) is sensitive to post-2013 extreme monsoon outliers at multiple stations across the peninsula. This study presents the first parallel TL-moment and L-moment RFA for the comprehensive 179-station DID network (1971–2023), simultaneously evaluating GEV, GLO, GPA, and K3D-II distributions across three climatologically distinct regions. Ward's minimum-variance hierarchical clustering was applied to TL-moment site characteristics, optimized by the average silhouette width (ASW) criterion. Discordancy, heterogeneity, goodness-of-fit (Z-test), and quantile estimation were executed in strict parallel under both estimation frameworks. Parametric bootstrap (B = 1,000 replicates) was applied to derive 90% confidence intervals for all regional quantiles. Three acceptably homogeneous regions were delineated: R1 (N = 55; west coast, mean = 115.0 mm), R2 (N = 94; interior, mean = 117.7 mm), and R3 (N = 30; east coast interior, mean = 200.9 mm). Under L-moments, GLO was best for R1 and R2; GPA was the sole passing distribution for R3. Under TL-moments, K3D-II was best for R2, GPA for R3. No standard distribution passed for R1 under TL-moments. TL-moment quantiles were 7–44% lower than L-moment estimates at T ≥ 10 years. Bootstrap 90% CI widths for T = 100 years were 0.076 (R1), 0.037 (R2), and 0.066 (R3) growth-factor units. L-moment and TL-moment quantiles should be used jointly as upper and lower design-rainfall bounds. For life-safety-critical structures, the L-moment estimate is the conservative upper bound.

1. Introduction

Peninsular Malaysia is one of the most flood-prone regions in Southeast Asia ¹. The Titiwangsa Range divides the peninsula into climatologically distinct zones, with the Northeast Monsoon (October–March) and Southwest Monsoon (May–September) producing sharply contrasting rainfall regimes on the east and west coasts, respectively ^{2, 3}. The 2014–2015 floods displaced 200,000 people; the 2006–2007 events caused RM 1.5 billion in infrastructure damage ^{4, 5}. Reliable regional design-rainfall estimates are therefore critical for safe hydraulic infrastructure planning under the Department of Irrigation and Drainage (DID) guidelines and the MSMA 2nd Edition design standard ⁶.

Regional frequency analysis (RFA), formalised by Hosking and Wallis ⁷, addresses the limitation of short at-site records by pooling data from statistically homogeneous networks. The L-moment framework ⁸ is the standard estimation basis. Still, equal weighting of all order statistics renders it sensitive to extreme outliers, a recognised problem for Malaysian annual maximum series, where post-2013 monsoon events substantially inflate regional skewness ^{11, 15}. Trimmed L-moments (TL-moments), introduced by Elamir and Seheult ⁹, assign zero weight to the most extreme order statistics at both tails, improving robustness without discarding data.

Four methodological gaps motivate this study: (i) no Malaysian RFA has simultaneously applied TL-moments to all four candidate distributions, including the Kappa Type-II (K3D-II) ¹⁷; (ii) all prior Malaysian studies use data through 2012 at most, predating the exceptional 2013–2023 monsoon extremes; (iii) no systematic parallel comparison of L-moment and TL-moment goodness-of-fit for GEV, GLO, GPA, and K3D-II exists for any Malaysian dataset; and (iv) no Malaysian RFA study has reported bootstrap confidence intervals for regional quantile estimates, leaving engineering practitioners without formal uncertainty bounds for design-rainfall decisions.

This study pursues four objectives: (i) delineate homogeneous regions from 179 quality-controlled DID stations using Ward's clustering with TL-moment site characteristics; (ii) compare L-moment and TL-moment discordancy, heterogeneity, and goodness-of-fit across regions; (iii) estimate normalized regional growth-factor quantiles at T = 2–200 years under both frameworks and all four distributions; and (iv) derive parametric bootstrap 90% confidence intervals for all regional quantile estimates. Results are directly applicable to DID MSMA design standards and national flood risk management practice.

2. Study Area and Data

Two monsoon systems govern Peninsular Malaysia's climate. The NEM (October–March) delivers heavy rainfall to east-coast states exposed to the South China Sea, while the Southwest Monsoon (May–September) is attenuated on the west coast by the Titiwangsa Range. Based on these relationships, three rainfall regimes are established: west coast, highland interior, and east coast interior, which are associated with three homogeneous regions defined from the clustering analysis (Figure 1).

Annual maximum daily rainfall data were obtained from the Department of Irrigation and Drainage (DID), Ministry of Natural Resources and Environment (MINE), Malaysia. The original network comprised 304 stations; after quality-control screening, 179 stations were retained based on three criteria: (a) record length n ≥ 15 years; (b) ≥ 300 valid daily observations per station-year; and (c) coefficient of variation CV ≤ 0.80. Seven stations were excluded (CV range 0.84–2.02). A sensitivity analysis retaining all 186 stations produced no changes in H₁ (L) or Z-test outcomes for R2 and R3; the H₁ (L) value for R1 increased marginally from −1.631 to −1.412, remaining within the acceptably homogeneous range. The retained dataset spans 1971–2023 (4,182 station-years; record lengths 15–49 years). Descriptive statistics are presented in Table 1.

3. Literature Review

Hosking ⁸ introduced L-moments as linear combinations of order statistics with lower bias than conventional moments for shape-parameter estimation in small, heavy-tailed samples. Hosking and Wallis ⁷ formalised the four-stage RFA framework: discordancy, heterogeneity, goodness-of-fit, and quantile estimation, now universally adopted ^{25, 30, 31, 32, 37}. Equal weighting of all order statistics, however, can inflate parameter estimates when a small number of extreme observations dominate the sample.

TL-moments ⁹ trim the extreme order statistics from each tail, achieving lower mean squared error than L-moments for leptokurtic populations. At t = 1, the RRMSE advantage is 15–30% at T ≤ 10 years ^{22, 34}, largest for sample sizes of 15–50 years, typical of Malaysian DID records. Hosking ¹⁰ analytically demonstrated lower MSE for GEV shape estimation when κ &gt 0, directly applicable to tropical annual maxima. Shabri et al. ¹¹ conducted the first Malaysian TL-moment RFA (40 Selangor stations), confirming GLO as the best-fit distribution but excluding K3D-II and data beyond 2010.

Hosking ¹⁶ introduced the four-parameter Kappa distribution, unifying GEV, GLO, GPA, and Gumbel as limiting cases. K3D-II is obtained by imposing h = −κ, retaining three parameters while covering a wider region of the moment-ratio diagram than any single standard distribution. Noor et al. ¹⁷ demonstrated improved K3D-II fit for Malaysian peninsular rainfall under L-moments; no study has evaluated K3D-II under TL-moment estimation.

Sahrin et al. ¹⁵ applied Ward's clustering to 83 peninsular stations, delineating seven regions, the only prior peninsula-wide objective clustering study, but restricted to pre-2013 data and excluding TL-moments and K3D-II. Ghobadi and Kang ²³ confirmed Ward's superiority over k-means and fuzzy methods across 85 regionalization studies; Hussain and Pasha ²¹ further demonstrated its effectiveness when combined with L-moments for regional flood frequency analysis.

Islam et al. ¹⁹ demonstrated that Bayesian RFA produces wider uncertainty bounds than L-moment point estimates at T = 100–200 years. Rahman et al. ²⁶ showed bootstrap uncertainty bounds are essential for reliability-based design in Australia. A universal limitation of Malaysian RFA studies is the absence of confidence intervals, a gap that the present study closes by deriving parametric bootstrap 90% CIs for all regional quantile estimates.

4. Methodology

The five-stage workflow (Figure 6) was implemented in Python 3.12 (NumPy, SciPy, scikit-learn) and applied in strict parallel under L-moments and TL-moments (t = 1). All code is openly available at: https:// doi.org/ 10.5281/ zenodo.15489203.

L-moments are linear combinations of order statistics ⁸. Dimensionless L-moment ratios are: τ₂ = l₂/l₁ (L-CV); τ₃ = l₃/l₂ (L-skewness); τ₄ = l₄/l₂ (L-kurtosis).

TL-moments at trimming level t = 1 assign zero weight to the single most extreme order statistic at each tail. The trimming level t = 1 was selected because it achieves the lowest RRMSE relative to L-moments for leptokurtic populations at record lengths of 15–49 years ⁹; higher trimming (t ≥ 2) imposes efficiency losses that outweigh the additional robustness gains for these sample sizes. TL-moment ratios τ₂⁽¹⁾ (TL-CV), τ₃⁽¹⁾ (TL-skewness), and τ₄⁽¹⁾ (TL-kurtosis) are defined analogously to L-moment ratios.

Four distributions are evaluated: GEV, GLO, GPA, and K3D-II. K3D-II encompasses GEV (κ → 0), GLO (κ = −1), and GPA (κ = 1) as special cases and covers a wider region of the L-moment and TL-moment ratio diagrams than any single three-parameter distribution ¹⁷.

No closed-form TL-moment expressions exist for K3D-II. The theoretical TL-kurtosis τ₄⁽¹⁾ as a function of TL-skewness τ₃⁽¹⁾ was derived numerically over τ₃⁽¹⁾ ∈ [0.05, 0.35]. Within the calibration range: maximum absolute residual = 0.00187, RMSE = 0.00063, Monte Carlo validation R² = 0.9991 (10,000 samples). A sensitivity test showed polynomial coefficient perturbation of ±0.002 shifts K3D-II Z-statistics by at most Δ|Z| = 0.04; no model-selection decision is affected. Extrapolation beyond τ₃⁽¹⁾ = 0.35 is not recommended.

Five sequential stages: discordancy (Dᵢ statistic), region delineation (Ward's clustering optimized by ASW, CH, DB indices), heterogeneity (H₁ statistic from 500 Monte Carlo simulations), goodness-of-fit (Z-statistic; |Z| ≤ 1.64 = PASS at 90% confidence), and quantile estimation (growth factors QT normalized to regional mean = 1).

Parametric bootstrap 90% CIs derived from B = 1,000 synthetic replicates. The 5th and 95th percentiles of replicated QT values define the lower and upper CI bounds for each return period T = 2–200 years.

5. Results

The ASW criterion nominally maximises at k = 2 (ASW = 0.2557), but k = 2 was rejected for five converging reasons: (1) merging east-coast and highland stations raises H₁ (L) to +0.84, worse than either constituent sub-region; (2) no distribution passes the Z-test for the merged cluster under either framework; (3) CH peaks at k = 3 (84.2 vs 79.1) and DB is minimized at k = 3 (1.21 vs 1.39); (4) the ASW difference is only Δ = 0.0032 (1.3%), within noise range; (5) the two distinct Malaysian monsoon systems independently justify three regions.

All six H₁ values are negative, confirming strong within-region statistical consistency. TL-CV is 30–44% lower than L-CV; TL-skewness is 29–48% lower; TL-kurtosis is 33–75% lower.

TL-moment growth factors are 7–44% below L-moment estimates at T ≥ 10 years. Inter-distributional spread is narrower under TL-moments (R1 at T = 100 yr: spread = 0.109 vs 0.192 under L-moments).

6. Discussion

The failure of all four distributions for R1 under TL-moments was investigated through four systematic diagnostic steps. First, R1's TL-kurtosis (0.1881) at TL-skewness 0.2413 exceeds the theoretical ceiling of all four distributions. Second, splitting R1 into northern and southern subsets along the Titiwangsa divide showed that both sub-groups also plot outside distribution loci, ruling out sub-cluster artefacts. Third, a leave-10-out sensitivity test dropped TL-kurtosis to 0.1693, still outside bounds, showing no single station subset drives the failure. Fourth, post-2013 stations exhibit TL-kurtosis values 0.03–0.07 higher than pre-2013 records, implicating intensified northeast monsoon extremes. The physical explanation is a dual-component rainfall regime: a moderate-intensity southwest monsoon core population coexists with a distinct high-intensity northeast monsoon extreme-event population that exceeds what any evaluated three-parameter family can accommodate.

Three approaches are identified for future resolution: (i) the five-parameter Wakeby distribution; (ii) asymmetric TL-moments trimming only from the upper tail; and (iii) at-site heterogeneity decomposition.

K3D-II's selection as best-fit for R2 under TL-moments (|Z| = 0.396 vs |Z-GLO| = 1.322) demonstrates that distributional preference shifts with estimation method, a finding not previously documented for Malaysian rainfall. K3D-II should be included as a standard candidate in any TL-moment-based RFA for Malaysian rainfall.

CI widths increase monotonically with return period. At T = 5 years, widths are 0.012 (R1), 0.008 (R2), 0.010 (R3); at T = 100 years, they reach 0.076, 0.037, 0.066; at T = 200 years, 0.102, 0.047, 0.088. In practical engineering terms, the lower CI bound suits risk-tolerant applications; the point estimate suits standard infrastructure; the upper CI bound governs life-safety-critical structures.

For R1: upper CI bound at T = 100 yr = 288.8 mm (governing value for dam spillways, primary embankments). For R2: design rainfall range at T = 100 yr = 243.1–247.3 mm (narrowest uncertainty among regions). For R3: upper CI bound at T = 100 yr = 451.1 mm (governing value for all consequence-critical east-coast assets).

The three-region partition is more refined than the seven-region solution of Sahrin et al. ¹⁵. GLO dominance under L-moments for R1 and R2 is consistent with Zin et al. ¹⁸. TL-moment quantile attenuation of 7–44% is consistent with the 15–30% RRMSE advantages reported by Shabri et al. ¹¹ and Ibrahim et al. ²².

Seven stations (CV > 0.80) were excluded. K3D-II TL-moment Z-statistics rely on a numerically approximated polynomial (residuals < 0.002). Mann–Kendall trend tests detected significant upward trends (p &lt 0.05) in 41 of 179 stations (22.9%), concentrated in R1 (30.9%) and R3 (30.0%), consistent with northeast monsoon intensification. Bootstrap CIs assume stationarity and represent a lower bound on total uncertainty in the presence of non-stationarity; an additional margin is recommended for life-safety-critical structures in R1 and R3 pending formal non-stationary analysis ^{20, 24, 35}.

7. Conclusion

This study presents the first parallel TL-moment and L-moment RFA for the 179-station peninsular Malaysian DID network (1971–2023). Ward's clustering (k = 3; ASW = 0.2525; CH = 84.2; DB = 1.21) delineated three acceptably homogeneous regions. TL-moments reduced regional skewness by 29–48% and kurtosis by 33–75%, yielding quantiles 7–44% lower at T ≥ 10 years. GLO was best under L-moments for R1 and R2; GPA for R3. Under TL-moments, K3D-II was best for R2; GPA for R3. No standard distribution passed for R1 under TL-moments. Bootstrap 90% CI widths at T = 100 years were 0.076, 0.037, and 0.066 growth-factor units for R1, R2, and R3, respectively.

Future priorities: (a) derive closed-form K3D-II TL-moment expressions; (b) evaluate the five-parameter Wakeby distribution and asymmetric upper-tail trimming for R1; (c) extend bootstrap CIs to TL-moment estimation frameworks; (d) develop formal non-stationary RFA models for R1 and R3 where Mann–Kendall testing detected upward trends in 30.9% and 30.0% of stations, respectively; (e) evaluate t = 2 trimming; (f) release an open-source Python package implementing the full TL-moment RFA workflow.

ACKNOWLEDGEMENTS

The authors acknowledge the Department of Irrigation and Drainage (DID), Ministry of Natural Resources and Environment Malaysia, for providing annual maximum daily rainfall data. Python scripts are openly archived at https://doi.org/10.5281/zenodo.15489203 (Zenodo, CC-BY 4.0).

Appendix A. K3D-II TL-Moment Polynomial: Coefficients and Validation

Appendix B. Code Availability

All Python 3.12 scripts are openly archived under CC-BY 4.0 at: https://doi.org/10.5281/zenodo.15489203

Key modules: (1) pwm_lmoments.py; (2) ward_clustering.py; (3) rfa_testing.py; (4) quantile_estimation.py; (5) bootstrap_ci.py.

References