Bounding copula functions for dependence modeling with applications to deccan plateau groundwater models

Introduction

Picture a monsoon that fails across two adjacent talukas simultaneously: both bore-well clusters draw down below the pump intake, tanker deliveries begin, and a local water crisis unfolds within weeks. Now picture the same rainfall deficit hitting one taluka while the other receives a late burst the crisis is halved. Which scenario is more likely? The answer depends not on the marginal probability of drought at each location but on the copula that ties the two marginal distributions together [1] . Classical Fréchet-Hoeffding bounds place the joint distribution between an upper co-monotone limit and a lower counter-monotone limit, but these bounds are often too wide to be operationally useful [2] . This research tightens those bounds by incorporating auxiliary information the rank correlation and a smoothness constraint derived from the subsurface flow equation into the Fréchet framework. The result is a family of "informed" bounds that shrink the uncertainty band on joint tail probabilities by roughly a third without committing to a specific parametric copula. The Deccan Plateau, stretching across Maharashtra, Karnataka, Madhya Pradesh, and Telangana, sits on Deccan basalt a hard-rock aquifer system with limited primary porosity and highly heterogeneous secondary (fracture) permeability [3] . Groundwater is the primary irrigation source for over 60% of the Plateau's farmland, and well failures during consecutive drought years have become increasingly common [4] . Modelling co-depletion risk requires capturing the spatial dependence between wells, a task for which copulas are well suited because they separate the marginal behaviour of each well from the joint dependence structure. This research fits five copula families to 34 well pairs, derives tightened analytical bounds on the joint survival function, validates the bounds by bootstrap, and applies Lyapunov stability analysis to classify each well as trending toward sustainable equilibrium or irreversible decline.

Theoretical framework

A copula C: [0,1]² → [0,1] satisfies the Fréchet-Hoeffding bounds max(u+v-1, 0) ≤ C(u,v) ≤ min(u,v) for all (u,v) ∈ [0,1]² [1] . The tightened bounds introduced here add a Spearman rank-correlation constraint: given observed ρ_S, the feasible set of copulas is reduced to those satisfying 12∫∫ C(u,v) du dv -3 = ρ_S. A variational argument over this constrained set yields point-wise upper and lower bounds C⁺(u,v;ρ_S) and C⁻(u,v;ρ_S) that are strictly tighter than Fréchet-Hoeffding whenever |ρ_S| > 0. The CW complex framework enters through the multi-well dependence graph: each well is a 0-cell, each pairwise copula a 1-cell, and each triple dependence a 2-cell, forming a simplicial complex whose Euler characteristic summarises the topology the dependence network [5] . Lyapunov stability was assessed by defining V(x) = Σ (h_i(t) -h_i*)², where h_i is the water table at well i and h_i* the equilibrium level, and checking whether dV/dt < 0 along trajectories of the coupled finite-element flow model [6] .

Materials and Methods

Materials

This research was carried out at the Bhopal University of Quantitative Studies, Bhopal, Madhya Pradesh, India, from August 2018 to December 2022. Monthly pre-monsoon (May) and post-monsoon (November) water-table depth records were obtained from the Central Ground Water Board (CGWB) for 34 monitoring wells distributed across five districts of the northern Deccan Plateau: Bhopal, Sehore, Raisen, Vidisha, and Hoshangabad. Records spanned 2008-2022 (15 years, 30 observations per well). Aquifer lithology (basalt, weathered mantle, alluvial patches) and estimated specific yield were sourced from CGWB district-level hydrogeological maps.

Methods

Marginal distributions were fitted to each well's time series using the L-moment method: gamma for alluvial wells (n = 11) and log-normal for basaltic wells (n = 23), selected by Anderson-Darling goodness-of-fit. Pseudo-observations (probability-integral transforms) were computed from the fitted marginals and used to estimate five copula families by maximum pseudo-likelihood [1] . AIC selected the best family for each well pair. Bootstrap resampling (B = 2 000) generated 95% confidence intervals for the copula parameter and the joint exceedance probability P(U > u, V > v) at the 90th-percentile threshold. The finite-element model used a 2D triangular mesh (3 412 elements, mean edge 500 m) with steady-state Darcy flow and heterogeneous transmissivity calibrated from pump-test data [7] . Lyapunov stability was assessed numerically by integrating the coupled flow-recharge ODE forward for 50 years under the 2008-2022 average recharge rate.

Implementation details

Copula fitting was performed in R 4.2 using the VineCopula and copula packages. The tightened Fréchet bounds were computed by a custom optimisation routine (Nelder-Mead in the nloptr package) that maximised and minimised C(u,v) subject to the ρ_S constraint on a 100×100 grid of (u,v) values. The finite-element model was built in MODFLOW-2005 with the FloPy 3.3 Python interface, and Lyapunov trajectories were integrated using SciPy's solve_ivp with the LSODA solver. Total computation time for all 34 wells was approximately 2.1 hours on a 4-core laptop. The CW complex was constructed using the gudhi 3.7 Python library for topological data analysis, and its Euler characteristic was computed from the alternating sum of simplex counts. The Student-t copula was selected at 24 of 34 well pairs, reflecting the heavy-tailed dependence structure expected in drought-prone hard-rock aquifers (Table 1). The estimated degrees-of-freedom parameter ν ranged from 3.2 to 7.8, with a median of 4.7. The tightened Fréchet bounds reduced the uncertainty band on joint 90th-percentile exceedance probabilities by an average of 31.4%, with the Student-t pairs showing the largest tightening (34.7%) because their high ρ_S values carry the most constraining power. Lyapunov analysis classified 21 wells as asymptotically stable (dV/dt < 0 everywhere along the 50-year trajectory) and 13 as unstable (Table 2). The unstable wells were concentrated in the basaltic zone, where low specific yield and limited lateral recharge pathways prevent recovery after consecutive drought years. All 9 stable alluvial wells showed dV/dt < -0.15, indicating rapid convergence to equilibrium.

Results

Comprehensive interpretation

The CW complex of the 34-well dependence network had Euler characteristic χ = 7, with 34 vertices (0-cells), 34 copula edges (1-cells), and 7 triple-dependence triangles (2-cells).

The positive χ indicates that the network is topologically simple dominated by tree-like chains of dependence rather than closed loops which aligns with the dendritic fracture geometry of the Deccan basalt. Bootstrap confidence intervals for the Student-t copula parameter ν were narrow (median width 1.4 degrees of freedom), confirming that 15 years of biannual data provide adequate power for copula estimation in this setting.

Discussion

The tightened Fréchet bounds represent the paper's main methodological contribution. By conditioning on the observed Spearman ρ_S, the bounds convert a model-free but wide interval into a semi-parametric and operationally useful one. The 31.4% average tightening means that a water manager using the informed bounds would face roughly one-third less ambiguity when estimating the probability of simultaneous well failure. This improvement comes at zero computational cost once ρ_S is known, making it attractive for resourceconstrained agencies.

The prevalence of the Student-t copula reflects the known behaviour of hard-rock aquifers: extreme drawdowns at one well tend to co-occur with extreme drawdowns at nearby wells because fracture-connected clusters deplete together during prolonged dry spells [3] . The Gaussian copula, which lacks tail dependence, was selected only at alluvial well pairs where drought co-occurrence is milder.

The Lyapunov stability results have policy implications. The 13 unstable wells are candidates for managed aquifer recharge (MAR) interventions: artificial recharge during the monsoon could shift dV/dt into the negative (stable) regime and prevent irreversible depletion. A limitation is that the 50-year projection assumes stationary recharge, which climate projections for the Deccan suggest may not hold under a warming scenario.

Conclusion

This research derived tightened analytical bounds on copula functions by conditioning on Spearman rank correlation and applied them to model groundwater co-depletion across 34 wells on the Deccan Plateau. The Student-t copula was selected as the best-fit family at 24 of 34 pairs, and the informed bounds reduced joint exceedance probability uncertainty by 31.4% relative to classical Fréchet-Hoeffding limits.

Lyapunov stability analysis classified 13 wells as unstable, all in the basaltic hard-rock zone, flagging them as priorities for managed aquifer recharge. The CW complex topology of the dependence network confirmed a tree-like fractureconnectivity structure consistent with Deccan basalt hydrogeology. Future work should incorporate non-stationary copula models that allow the dependence parameter to evolve with climate indices and extend the tightened bounds to the multivariate case using vine copulas.

Acknowledgement

Acknowledgements Funding sourcesThis research was funded through internal departmental allocations. No external grant numbers are applicable. Institutional supportThe authors thank the Bhopal University of Quantitative Studies for computational resources and the Central Ground Water Board for supplying well-monitoring records.

References

1. Nelsen RB. An introduction to copulas. 2nd ed. New York: Springer; 2006. p. 1-269.
2. Joe H. Dependence modeling with copulas. Boca Raton (FL): CRC Press; 2014. p. 1-462.
3. Kulkarni H, Deolankar SB, Lalwani A, Joseph B, Pawar S. Hydrogeological framework of the Deccan basalt groundwater systems. Hydrogeology Journal. 2000;8(4):368-378.
4. Shah T. Taming the anarchy: groundwater governance in South Asia. Washington (DC): Resources for the Future Press; 2009. p. 1-310.
5. Hatcher A. Algebraic topology. Cambridge: Cambridge University Press; 2002. p. 1-544.
6. Khalil HK. Nonlinear systems. 3rd ed. Upper Saddle River (NJ): Prentice Hall; 2002. p. 1-750.
7. Harbaugh AW. MODFLOW-2005, the US Geological Survey modular ground-water model. US Geological Survey Techniques and Methods. 2005;6-A16:1-253.
8. Genest C, Favre AC. Everything you always wanted to know about copula modeling but were afraid to ask. Journal of Hydrologic Engineering. 2007;12(4):347-368.
9. Salvadori G, De Michele C, Kottegoda NT, Rosso R. Extremes in nature: an approach using copulas. Water Science and Technology Library. Vol. 56. Dordrecht: Springer; 2007. p. 1-292.
10. Marechal JC, Dewandel B, Ahmed S, Galeazzi L, Zaidi FK. Combined estimation of specific yield and natural recharge in a semi-arid groundwater basin. Journal of Hydrology. 2006;329(1-2):281-293.
11. Sklar M. Fonctions de répartition à n dimensions et leurs marges. Publications de l'Institut de Statistique de l'Université de Paris. 1959;8:229-231.
12. Rodell M, Velicogna I, Famiglietti JS. Satellite-based estimates of groundwater depletion in India. Nature. 2009;460(7258):999-1002.
13. Efron B, Tibshirani RJ. An introduction to the bootstrap. New York: Chapman and Hall; 1993. p. 1-436.
14. Deshpande RD, Bhatt SB, Gupta SK. Groundwater in the Deccan basalts of India. Memoir of the Geological Society of India. 1998;48:1-24.
15. Demarta S, McNeil AJ. The t copula and related copulas. International Statistical Review. 2005;73(1):111-129.
16. Foster S, Chilton J, Moench M, Cardy F, Schiffler M. Groundwater in rural development: facing the challenges of supply and resource sustainability. World Bank Technical Paper. 2000;463:1-98.