Groundwater depletion: A mathematical model incorporating climate and human factors

Kottakkaran Sooppy Nisar; Raghad Mohammed Al-Suliman; Mada Samhoud Al-Qahtani

doi:10.2478/ijmce-2026-0016

Introduction

Groundwater is an essential component of the Earth’s water cycle, providing critical support for ecosystems, agriculture, and human consumption. It serves as a primary source of freshwater, especially in regions where surface water is scarce. However, groundwater resources are increasingly under pressure due to human activities such as over-extraction, pollution, and deforestation [1–5]. These environmental challenges disrupt the natural balance of groundwater replenishment, leading to long-term consequences such as depletion and contamination. Reports have highlighted the urgent need for effective groundwater management strategies, emphasizing the climate change’s effects on water resources. However, existing models often lack detailed representations of the additional stresses imposed by human activities like pollution and deforestation [6–11]. In addition to traditional ODEs, the application of fractional calculus has gained traction in various fields, providing a framework to model procedures that display memory impacts and non-local behavior. Fractional models can capture the complexities of systems that are not well represented by classical integer-order models, allowing for more accurate predictions in scenarios influenced by historical states and environmental factors [12–19]. Recent research has also explored the controllability of such systems using advanced mathematical tools, such as fixed point theorems in Banach spaces, to analyze impulsive fractional integro-differential equations, further extending the utility of fractional models in environmental and engineering applications [20]. Numerous studies have focused on developing mathematical models to describe the dynamics of groundwater systems. Researchers have employed ODE models to capture the interactions between atmospheric water, surface water, and groundwater [21–25]. These models help in understanding how various factors, including precipitation, evaporation, infiltration, and pumping, affect groundwater levels over time. However, existing models often lack detailed representations of the additional stresses imposed by human activities like pollution and deforestation. Previous researchers have made significant strides in modeling groundwater dynamics using ODEs. For instance, studies have incorporated key environmental processes such as evaporation, infiltration, and recharge into models to predict groundwater levels under various scenarios. Additionally, some studies have looked at how urbanization and climate change affect water cycles, shedding light on the potential impacts of these factors on groundwater systems [26–31]. These models serve as a foundation for developing strategies to manage groundwater resources more effectively. Moreover, water pollution stands out as a significant environmental challenge confronting developing nations. Mathematical modeling has proven effective in analyzing the transmission of water pollutants and their impact on ecosystems and public health. Numerical approaches, such as the use of shifted Jacobi polynomials, were also adopted to convert the model into an algebraic form and validate its performance against traditional Runge-Kutta methods [32]. In recent research, a system of ODEs was used to model soluble and insoluble pollutants, with sensitivity analysis performed on the reproduction number to assess intervention strategies [33]. These findings underscore the importance of integrating pollution dynamics into groundwater models to better reflect real-world complexities.

Building upon the existing ODE models, our research introduces additional terms to better capture the complex interactions between groundwater and environmental stressors. Specifically, we modify the standard ODE framework by incorporating terms for pollution, frequent water pumping, and deforestation. These new terms provide a more comprehensive representation of the factors affecting groundwater dynamics in modern contexts. Our work focuses on analyzing equilibrium points and stability while deriving numerical simulations to validate the model’s performance across different environmental scenarios [34, 35]. By refining the ODE model and incorporating real-world factors, our approach offers deeper insights into groundwater behavior, making it a valuable tool for resource management and sustainability planning.

This study is structured as follows: Section 2 outlines the mathematical model’s fundamental assumptions and parameter definitions. Section 3 details the system’s equilibrium points. Section 4 examines the stability of these equilibrium points. In Section 5, we use numerical simulations to test the model and evaluate its behaviour in different environmental circumstances. Section 6 summarises the study’s results and implications for managing groundwater resources. Section 7 completes the paper by presenting novelties of this work.

Model description

In this part of the paper, we present the model studied as follows 1 ${\begin{array}{l} \frac{d A}{d t} = ρ_{1} B (t) + ρ_{2} C (t) - ϖ A (t) - ı A (t) + κ \frac{B (t)}{C (t) + 1}, \\ \frac{d B}{d t} = ϖ A (t) - v B (t) C (t) - ρ_{1} B (t) + ϑ C (t) B (t) + η C (t) - ζ \frac{C (t)}{A (t) + 1}, \\ \frac{d C}{d t} = v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) - ϰ C (t) + ς \frac{A (t)}{B (t) + 1} . \end{array}$ \left\{ {\matrix{ {{{dA} \over {d}} = {\rho _1}B() + {\rho _2}C() - \varpi A() - \imath A() + \kappa {{B()} \over {C() + 1}},} \hfill \cr {{{dB} \over {d}} = \varpi A() - vB()C() - {\rho _1}B() + \vartheta C()B() + \eta C() - \zeta {{C()} \over {A() + 1}},} \hfill \cr {{{dC} \over {d}} = vC()B() - {\rho _2}C() - \vartheta C()B() - \eta C() - C() + \varsigma {{A()} \over {B() + 1}}.} \hfill \cr } } \right.

The total amount of water N(𝔱) = A(𝔱) + B(𝔱) + C(𝔱). Initial condition A(0) = A₀, B(0) = B₀ and C(0) = C₀. The model parameters of equation (1) are defined by Table 1.

Table 1

Model parameterization.

Parameter	Parameter Description
A(𝔱)	Atmospheric water
B(𝔱)	Surface water
C(𝔱)	Groundwater
ρ₁&ρ₂	Surface and groundwater evaporation rates to atmospheric water, correspondingly
ϖ	Rate of atmospheric precipitation reaching surface water bodies
κ	Influence of surface water on atmospheric water, moderated by groundwater
ι	Rate at which atmospheric water dissipates
V	Surface-to-groundwater infiltration rate
ϑ	The effect of pollution on groundwater and surface water
η	Recurrence rate of groundwater extraction
ζ	Reduction in surface water due to groundwater and atmospheric interaction
ϰ	Pace of deforestation
ς	Contribution of atmospheric water to groundwater recharge, depending on surface water

Theorem 1.

For each non-negative initial condition, ∃ a unique solution of (1).

Proof.

The method utilized by Moustafa [16] is applied. The region $R = {A (t), B (t), C (t) \in ℝ_{+}^{3} : | A |, | B |, | C |}$ R = \left\{ {A({\rm{t}}),B({\rm{t}}),C({\rm{t}}) \in _ + ^3:|A|,|B|,|C|} \right\}. We use the term T = (A, B, C) and $\bar{T} = (\bar{A}, \bar{B}, \bar{C})$ \bar T = (\bar A,\bar B,\bar C), create a mapping $E (T) = (E_{1} (T), E_{2} (T), E_{3} (T))$ E(T) = \left( {{E_1}(T),{E_2}(T),{E_3}(T)} \right) where ${\begin{array}{l} E_{1} (T) & = & ρ_{1} B (t) + ρ_{2} C (t) - ϖ A (t) - ı A (t) + κ \frac{B (t)}{C (t) + 1}, \\ E_{2} (T) & = & ϖ A (t) - v B (t) C (t) - ρ_{1} B (t) + ϑ C (t) B (t) + η C (t) - ζ \frac{C (t)}{A (t) + 1}, \\ E_{3} (T) & = & v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) - ϰ C (t) + ς \frac{A (t)}{B (t) + 1} . \end{array}$ \begin{array}{lcc}E_1(T)&=&\rho_1B(\mathfrak t)+\rho_2C(\mathfrak t)-\varpi A(\mathfrak t)-\i A(\mathfrak t)+\kappa\frac{\displaystyle B(\mathfrak t)}{\displaystyle C(\mathfrak t)+1},\\E_2(T)&=&\varpi A(\mathfrak t)-vB(\mathfrak t)C(\mathfrak t)-\rho_1B(\mathfrak t)+\vartheta C(\mathfrak t)B(\mathfrak t)+\eta C(\mathfrak t)-\zeta\frac{\displaystyle C(\mathfrak t)}{\displaystyle A(\mathfrak t)+1},\\E_3(T)&=&vC(\mathfrak t)B(\mathfrak t)-\rho_2C(\mathfrak t)-\vartheta C(\mathfrak t)B(\mathfrak t)-\eta C(\mathfrak t)-\varkappa C(\mathfrak t)+\varsigma\frac{\displaystyle A(\mathfrak t)}{\displaystyle B(\mathfrak t)+1}.\end{array}

We consider $\begin{array}{l} ‖ E (T) - E (\bar{T}) ‖ & = & | E_{1} (T) - E_{1} (\bar{T}) | + | E_{2} (T) - E_{2} (\bar{T}) | + | E_{3} (T) - E_{3} (\bar{T}) | \\ \leq |_{ρ 1} B (t) + ρ_{2} C (t) - ϖ A (t) - ı A (t) + κ \frac{B (t)}{C (t) + 1} - ρ_{1} \bar{B (t)} - ρ_{2} \bar{C (t)} + ϖ \bar{A (t)} + \\ ı \bar{A (t)} - κ \frac{\bar{B (t)}}{\bar{C (t)} + 1} | + | ϖ A (t) - v B (t) C (t) - ρ_{1} B (t) + ϑ C (t) B (t) + η C (t) - \\ ζ \frac{C (t)}{A (t) + 1} - ϖ \bar{A (t)} + v \bar{B (t) C (t)} + ρ_{1} \bar{B (t)} - ϑ \bar{C (t) B (t)} - η \bar{C (t)} + ζ \frac{\bar{C (t)}}{\bar{A (t)} + 1} | \\ + | v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) - ϰ C (t) + ς \frac{A (t)}{B (t) + 1} - v \bar{C (t) B (t)} \\ - ρ_{2} \bar{C (t)} - ϑ \bar{C (t) B (t)} - η \bar{C (t)} - ϰ \bar{C (t)} + ς \frac{\bar{A (t)}}{\bar{B (t)} + 1} - \frac{\bar{B (t)}}{\bar{C (t)} + 1} | \\ = & | A (t) - \bar{A (t)} | (- ı - \frac{C (t)}{{(A (t) + 1)}^{2}} + \frac{ζ}{B (t) + 1}) + | B (t) - \bar{B (t)} | (\frac{κ}{C (t) + 1} - \\ \frac{A (t)}{{(B (t) + 1)}^{2}}) + | C (t) - \bar{C (t)} | (- \frac{κ B (t)}{{(C (t) + 1)}^{2}} - \frac{ζ}{A (t) + 1} - ϰ - \frac{ζ A (t)}{{(B (t) + 1)}^{2}}) \\ \leq & ℂ_{1} | A (t) - \bar{A (t)} | + ℂ_{2} | B (t) - \bar{B (t)} | + ℂ_{3} | C (t) - \bar{C (t)} | \\ \leq & ℂ_{\max} | Ϝ - \bar{Ϝ} |, \end{array}$ \matrix{ {E(T) - E(\bar T)} \hfill & = \hfill & {\left| {{E_1}(T) - {E_1}(\bar T)} \right| + \left| {{E_2}(T) - {E_2}(\bar T)} \right| + \left| {{E_3}(T) - {E_3}(\bar T)} \right|} \hfill \cr {} \hfill & {} \hfill & { \le {\rho _1}B({\rm{t}}) + {\rho _2}C({\rm{t}}) - \bar \varpi A({\rm{t}}) - \imath A({\rm{t}}) + \kappa {{B({\rm{t}})} \over {C({\rm{t}}) + 1}} - {\rho _1}\overline {B({\rm{t}})} - {\rho _2}\overline {C({\rm{t}})} + \bar \varpi \overline {A({\rm{t}})} + } \hfill \cr {} \hfill & {} \hfill & {\imath \overline {A({\rm{t}})} - \kappa {{\overline {B({\rm{t}})} } \over {\overline {C({\rm{t}})} + 1}}| + |\omega A({\rm{t}}) - vB({\rm{t}})C({\rm{t}}) - {\rho _1}B({\rm{t}}) + \vartheta C({\rm{t}})B({\rm{t}}) + \eta C({\rm{t}}) - } \hfill \cr {} \hfill & {} \hfill & {\zeta {{C({\rm{t}})} \over {A({\rm{t}}) + 1}} - \varpi \overline {A({\rm{t}})} + v\overline {B({\rm{t}})C({\rm{t}})} + {\rho _1}\overline {B({\rm{t}})} - \vartheta \overline {C({\rm{t}})B({\rm{t}})} - \eta \overline {C({\rm{t}})} + \zeta {{\overline {C({\rm{t}})} } \over {\overline {A({\rm{t}})} + 1}}|{\mkern 1mu} } \hfill \cr {} \hfill & {} \hfill & { + |{\mkern 1mu} vC({\rm{t}})B({\rm{t}}) - {\rho _2}C({\rm{t}}) - \vartheta C({\rm{t}})B({\rm{t}}) - \eta C({\rm{t}}) - C({\rm{t}}) + \varsigma {{A({\rm{t}})} \over {B({\rm{t}}) + 1}} - v\overline {C({\rm{t}})B({\rm{t}})} } \hfill \cr {} \hfill & {} \hfill & { - {\rho _2}\overline {C({\rm{t}})} - \vartheta \overline {C({\rm{t}})B({\rm{t}})} - \eta \overline {C({\rm{t}})} - \overline {C({\rm{t}})} + \varsigma {{\overline {A()} } \over {\overline {B({\rm{t}})} + 1}} - {{\overline {B()} } \over {\overline {C({\rm{t}})} + 1}}|{\mkern 1mu} } \hfill \cr {} \hfill & = \hfill & {|A() - \overline {A({\rm{t}})} |\left( { - \imath - {{C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}} + {\zeta \over {B({\rm{t}}) + 1}}} \right) + |B({\rm{t}}) - \overline {B({\rm{t}})} |\left( {{\kappa \over {C({\rm{t}}) + 1}} - } \right.} \hfill \cr {} \hfill & {} \hfill & {\left. {{{A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}}} \right) + |C({\rm{t}}) - \overline {C({\rm{t}})} |\left( { - {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}} - {\zeta \over {A({\rm{t}}) + 1}} - - {{\zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}}} \right)} \hfill \cr {} \hfill & \le \hfill & {{_1}|A({\rm{t}}) - \overline {A({\rm{t}})} | + {_2}|B({\rm{t}}) - \overline {B({\rm{t}})} | + {_3}|C({\rm{t}}) - \overline {C({\rm{t}})} |} \hfill \cr {} \hfill & \le \hfill & {{_{\max }}| - \bar |,} \hfill \cr } where ℂ_max = max{ℂ₁, ℂ₂, ℂ₃}.

As a result, we deduce that the +^vely invariant set induced by the model (1) is the area R. In the region R, the model is well-posed both mathematically and biologically. Thus, the existence criterion of the system (1) is established.

Theorem 2.

The solution of system (1) which start in $R_{+}^{3}$ R_ + ^3 are uniformly bounded and non-negative.

Proof.

Since the total amount of water is N(𝔱) = A(𝔱) + B(𝔱) + C(𝔱), the rate of change of the water content N(𝔱), $\begin{array}{l} \frac{d N}{d t} & = \frac{d A}{d t} + \frac{d B}{d t} + \frac{d C}{d t} \\ = ρ_{1} B (t) + ρ_{2} C (t) - ϖ A (t) - ı A (t) + κ \frac{B (t)}{C (t) + 1} + ϖ A (t) - v B (t) C (t) - ρ_{1} B (t) \\ + ϑ C (t) B (t) + η C (t) - ζ \frac{C (t)}{A (t) + 1} + v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) \\ - ϰ C (t) + ς \frac{A (t)}{B (t) + 1} \\ \leq - ı A (t) - ϰ C (t) + κ \frac{B (t)}{C (t) + 1} - ζ \frac{C (t)}{A (t) + 1} + ς \frac{A (t)}{B (t) + 1} . \end{array}$ \matrix{ {{{dN} \over {d{\rm{t}}}}} \hfill & { = {{dA} \over {d{\rm{t}}}} + {{dB} \over {d{\rm{t}}}} + {{dC} \over {d{\rm{t}}}}} \hfill \cr {} \hfill & { = {\rho _1}B({\rm{t}}) + {\rho _2}C({\rm{t}}) - \varpi A({\rm{t}}) - \imath A({\rm{t}}) + \kappa {{B({\rm{t}})} \over {C({\rm{t}}) + 1}} + \varpi A({\rm{t}}) - vB({\rm{t}})C({\rm{t}}) - {\rho _1}B({\rm{t}})} \hfill \cr {} \hfill & { + \vartheta C({\rm{t}})B({\rm{t}}) + \eta C({\rm{t}}) - \zeta {{C({\rm{t}})} \over {A({\rm{t}}) + 1}} + vC({\rm{t}})B({\rm{t}}) - {\rho _2}C({\rm{t}}) - \vartheta C({\rm{t}})B({\rm{t}}) - \eta C({\rm{t}})} \hfill \cr {} \hfill & { - C({\rm{t}}) + \varsigma {{A({\rm{t}})} \over {B({\rm{t}}) + 1}}} \hfill \cr {} \hfill & { \le - \imath A({\rm{t}}) - C({\rm{t}}) + \kappa {{B({\rm{t}})} \over {C({\rm{t}}) + 1}} - \zeta {{C({\rm{t}})} \over {A({\rm{t}}) + 1}} + \varsigma {{A({\rm{t}})} \over {B({\rm{t}}) + 1}}.} \hfill \cr }

This implies that $\frac{d N}{d t}$ {{dN} \over {d{\rm{t}}}} is bounded above and below by terms involving – (ιA(𝔱) + ϰC(𝔱)). Integrating the above inequality and using initial conditions, we obtain $N (0) e^{(ι + ϰ) t} \leq N (t) \leq N (0) e^{(ι + ϰ) t} .$ N(0){e^{(l + )}} \le N() \le N(0){e^{(l + )}}.

Considering t → ∞, we have $\lim_{t \to \infty} i n f N (t) \leq N (t) \leq \lim_{t \to \infty} s u p N (t) .$ \mathop {\lim }\limits_{ \to \infty } in\;fN({\rm{t}}) \le N({\rm{t}}) \le \mathop {\lim }\limits_{ \to \infty } supN({\rm{t}}).

Hence the feasible region for the system (1) is $R = {A (t), B (t), C (t) \in ℝ_{+}^{3} : 0 < | A |, | B |, | C | \leq N (0) e^{(ı + ϰ) t}} .$ R = \left\{ {A({\rm{t}}),B({\rm{t}}),C({\rm{t}}) \in _ + ^3:0 < |A|,|B|,|C| \le N(0){e^{(\imath + ){\rm{t}}}}} \right\}.

Hence the region R is positive invariant so that no solution path moves beyond the boundary of R. Thus above theorem ensures that the proposed model is feasible both biologically and mathematically.

Results of equilibrium

In order to identify the model’s equilibrium points (1), we must solve $\frac{d A}{d t} = \frac{d B}{d t} = \frac{d C}{d t} = 0.$ {{dA} \over {d{\rm{t}}}} = {{dB} \over {d{\rm{t}}}} = {{dC} \over {d{\rm{t}}}} = 0.

The system then assumes the subsequent configuration 2 ${\begin{array}{l} ρ_{1} B (t) + ρ_{2} C (t) - ϖ A (t) - ı A (t) + κ \frac{B (t)}{C (t) + 1} = 0, \\ ϖ A (t) - v B (t) C (t) - ρ_{1} B (t) + ϑ C (t) B (t) + η C (t) - ζ \frac{C (t)}{A (t) + 1} = 0, \\ v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) - ϰ C (t) + ς \frac{A (t)}{B (t) + 1} = 0. \end{array}$ \left\{ {\matrix{ {{\rho _1}B({\rm{t}}) + {\rho _2}C({\rm{t}}) - \varpi A({\rm{t}}) - \imath A({\rm{t}}) + \kappa {{B({\rm{t}})} \over {C({\rm{t}}) + 1}} = 0,} \hfill \cr {\varpi A({\rm{t}}) - vB({\rm{t}})C({\rm{t}}) - {\rho _1}B({\rm{t}}) + \vartheta C({\rm{t}})B({\rm{t}}) + \eta C({\rm{t}}) - \zeta {{C({\rm{t}})} \over {A({\rm{t}}) + 1}} = 0,} \hfill \cr {vC({\rm{t}})B({\rm{t}}) - {\rho _2}C({\rm{t}}) - \vartheta C({\rm{t}})B({\rm{t}}) - \eta C({\rm{t}}) - C({\rm{t}}) + \varsigma {{A({\rm{t}})} \over {B({\rm{t}}) + 1}} = 0.} \hfill \cr } } \right.

When it comes to pollution-free and groundwater-free pumping, we consider B = B₀.

Let ${\bar{E}}_{1} (\bar{A (t)}, \bar{B (t)}, \bar{C (t)})$ {{\bar E}_1}(\overline {A({\rm{t}})} ,\overline {B({\rm{t}})} ,\overline {C({\rm{t}})} ) be the pollution free equilibrium point. Using B = B₀ in the system (2), and solving the system of algebraic equation, we obtain $\bar{A (t)} = \frac{ρ_{1} B_{0} (t) + ρ_{2} \bar{C (t)} + \frac{κ B_{0} (t)}{C (t) + 1}}{ϖ + ı},$ \overline {A({\rm{t}})} = {{{\rho _1}{B_0}({\rm{t}}) + {\rho _2}\overline {C({\rm{t}})} + {{\kappa {B_0}({\rm{t}})} \over {C({\rm{t}}) + 1}}} \over {\varpi + \imath }}, and $\begin{array}{l} \bar{C (t)} [v B (t) - ρ_{2} - ϑ B (t) - ϰ - η] = - \frac{ζ \bar{A (t)}}{B (t) + 1}, \\ \bar{C (t)} = \frac{ϰ \bar{A (t)}}{(B (t) + 1) (ρ_{2} + ϑ B_{0} + η + ϰ - v B_{0})} . \end{array}$ \matrix{ {\overline {C({\rm{t}})} \left[ {vB({\rm{t}}) - {\rho _2} - \vartheta B({\rm{t}}) - - \eta } \right] = - {{\zeta \overline {A({\rm{t}})} } \over {B({\rm{t}}) + 1}},} \hfill \cr {\overline {C({\rm{t}})} = {{\overline {A({\rm{t}})} } \over {(B({\rm{t}}) + 1)\left( {{\rho _2} + \vartheta {B_0} + \eta + - v{B_0}} \right)}}.} \hfill \cr }
The Atmospheric-Surface water equilibrium point E₂ = (A(t), B(t), 0) where $\begin{array}{l} ρ_{1} B (t) - ϖ A (t) - ı A (t) + κ B (t) = 0, \\ (ϖ + ı) A (t) - (ρ_{1} + κ) B (t) = 0, \\ A (t) = \frac{(ρ_{1} + κ) B (t)}{ϖ + ı} . \end{array}$ \matrix{ {{\rho _1}B({\rm{t}}) - \varpi A({\rm{t}}) - \imath A({\rm{t}}) + \kappa B({\rm{t}}) = 0,} \hfill \cr {(\varpi + \imath )A({\rm{t}}) - \left( {{\rho _1} + \kappa } \right)B({\rm{t}}) = 0,} \hfill \cr {A({\rm{t}}) = {{\left( {{\rho _1} + \kappa } \right)B({\rm{t}})} \over {\varpi + \imath }}.} \hfill \cr }
The Atmospheric-Ground water equilibrium point E₃ = (A(𝔱), 0, C(𝔱)), $\begin{array}{l} ρ_{2} C (t) - ϖ A (t) - ı A (t) = 0, \\ C (t) = \frac{ϖ + ı}{ρ_{2}} A (t) . \end{array}$ \matrix{ {{\rho _2}C({\rm{t}}) - \varpi A({\rm{t}}) - \imath A({\rm{t}}) = 0,} \hfill \cr {C({\rm{t}}) = {{\varpi + \imath } \over {{\rho _2}}}A({\rm{t}}).} \hfill \cr }
The Surface-Ground water equilibrium point E₄ = (0, B(𝔱), C(𝔱)), $\begin{array}{l} v C (t) B (t) - ρ_{2} C (t) - ϑ C (t) B (t) - η C (t) - ϰ C (t) = 0, \\ B (t) C (t) (v - ϑ) = C (t) (ρ_{2} + η + ϰ), \\ B (t) = \frac{ρ_{2} + η + ϰ}{v - ϑ} . \end{array}$ \matrix{ {vC({\rm{t}})B({\rm{t}}) - {\rho _2}C({\rm{t}}) - \vartheta C({\rm{t}})B({\rm{t}}) - \eta C({\rm{t}}) - C({\rm{t}}) = 0,} \hfill \cr {B({\rm{t}})C({\rm{t}})(v - \vartheta ) = C({\rm{t}})\left( {{\rho _2} + \eta + } \right),} \hfill \cr {B({\rm{t}}) = {{{\rho _2} + \eta + } \over {v - \vartheta }}.} \hfill \cr }

The local stability of steady state points

The stability criteria of the possible stable state points are described as follows.

Theorem 3.

The steady state point E₀ is conditionally locally asymptotically stable if the following condition satisfied ρ₁ < κ [36].

Proof.

At E₀, the Jacobian matrix for (1) is given by $\begin{matrix} J (A, B, C)= (\begin{matrix} - ı - ϖ & ρ_{1} + \frac{κ}{C (t) + 1} & ρ_{2} - \frac{κ B (t)}{{(C (t) + 1)}^{2}} \\ ϖ - \frac{ζ C (t)}{{(A (t) + 1)}^{2}} & - v C (t) - ρ_{1} + ϑ C (t) & - v B (t) + ϑ B (t) + η - \frac{ζ}{A (t) + 1} \\ \frac{ς}{B (t) + 1} & v C (t) - ϑ C (t) - \frac{ς A (t)}{{(B (t) + 1)}^{2}} & v B (t) - ρ_{2} - ϑ B (t) - η - ϰ \end{matrix}) \\ J (E_{0}) = (\begin{matrix} - (ı + ϖ) & ρ_{1} + κ & ρ_{2} - κ B (t) \\ ϖ & - ρ_{1} & - (η + ζ) \\ 0 & 0 & - (ρ_{2} + η + ϰ) \end{matrix}) \end{matrix}$ \matrix{ {J{\rm{(}}A{\rm{, }}B{\rm{, }}C{\rm{) = }}\left( {\matrix{ { - \imath - \varpi } & {{\rho _1} + {\kappa \over {C({\rm{t}}) + 1}}} & {{\rho _2} - {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}}} \cr {\omega - {{\zeta C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}}} & { - vC({\rm{t}}) - {\rho _1} + \vartheta C({\rm{t}})} & { - vB({\rm{t}}) + \vartheta B({\rm{t}}) + \eta - {\zeta \over {A({\rm{t}}) + 1}}} \cr {{\zeta \over {B({\rm{t}}) + 1}}} & {vC({\rm{t}}) - \vartheta C({\rm{t}}) - {{\zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}}} & {vB({\rm{t}}) - {\rho _2} - \vartheta B({\rm{t}}) - \eta - } \cr } } \right)} \cr {J\left( {{E_0}} \right) = \left( {\matrix{ { - (\imath + \varpi )} & {{\rho _1} + \kappa } & {{\rho _2} - \kappa B({\rm{t}})} \cr \varpi & { - {\rho _1}} & { - (\eta + \zeta )} \cr 0 & 0 & { - \left( {{\rho _2} + \eta + } \right)} \cr } } \right)} \cr } (λ₁ – a₃₃)(λ² – (a₁₁ + a₂₂)λ + (a₁₁a₂₂ – a₁₂a₂₁)) = 0.

The characteristic polynomial of J(E₀) is (λ₁ + ρ₂ + η + ϰ)((λ² + λ(ι + ϖ + ρ₁) + (ι + ϖ)ρ₁ – ϖ(ρ₁ + κ)) = 0. It’s observed that all roots of the characteristic polynomial of J(E₀) is equation have a –^ve real part, ρ₁ < κ.

Theorem 4.

The steady point E₁ is conditionally asymptotically stable if the following conditions are satisfied $κ B (t) > ρ_{2}, \frac{ζ}{A (t) + 1} + v B (t) > ϑ B (t) + η$ \kappa B({\rm{t}}) > {\rho _2},{\zeta \over {A({\rm{t}}) + 1}} + vB({\rm{t}}) > \vartheta B({\rm{t}}) + \eta and ϑB(𝔱) + η + ϰ + ρ₂ > vB(𝔱).

Proof.

At E₁, then Jacobian matrix of (1) is, $J (E_{1}) = (\begin{matrix} - (ı + ϖ) & ρ_{1} + κ & ρ_{2} - κ B (t) \\ ϖ & - ρ_{1} & - v B (t) + ϑ B (t) + η - \frac{ζ}{A (t) + 1} \\ \frac{ς}{B (t) + 1} & \frac{- ς A (t)}{{(B (t) + 1)}^{2}} & v B (t) - ρ_{2} - ϑ B (t) - η - ϰ \end{matrix})$ J\left( {{E_1}} \right) = \left( {\matrix{ { - (\imath + \varpi )} & {{\rho _1} + \kappa } & {{\rho _2} - \kappa B({\rm{t}})} \cr \varpi & { - {\rho _1}} & { - vB({\rm{t}}) + \vartheta B({\rm{t}}) + \eta - {\zeta \over {A({\rm{t}}) + 1}}} \cr {{\zeta \over {B({\rm{t}}) + 1}}} & {{{ - \zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}}} & {vB({\rm{t}}) - {\rho _2} - \vartheta B({\rm{t}}) - \eta - } \cr } } \right) λ³ – [a₁₁ + a₂₂ + a₃₃]λ² + [a₂₂a₃₃ – a₂₃a₃₂ + a₁₁a₃₃ – a₃₁a₁₃ + a₁₁ a₂₂ – a₁₂a₂₁]λ –a₁₁[a₂₂a₃₃ – a₂₃a₃₂] + a₁2[a₂₁a₃₃ – a₂₃a₃₁] – a₁₃[a₂₁a₃₂ – a₃₁a₂₂] = 0 when we examine the element of Jacobian matrix, we can see that the signs of elements a₁₁, a₂₂ and a₃₂ are strictly –^ve, while the signs of elements a₁₂, a₂₁ and a₃₁ are strictly +^ve. On the other hand the signs of elements a₁₃, a₂₃ and a₃₃ have –^ve signs provided that $κ B (t) > ρ_{2}, \frac{ζ}{A (t) + 1} + v B (t) > ϑ B (t) + η$ \kappa B({\rm{t}}) > {\rho _2},{\zeta \over {A({\rm{t}}) + 1}} + vB({\rm{t}}) > \vartheta B({\rm{t}}) + \eta and ϑB(𝔱) + η + ϰ + ρ₂ > vB(𝔱) respectively.

Theorem 5.

The steady state point E₂ is conditionally asymptotically stable if the following conditions are satisfied $\frac{ζ C (t)}{{(A (t) + 1)}^{2}} > ϖ, v C (t) + ρ_{1} > ϑ C (t), \frac{ζ}{A (t) + 1} > η$ {{\zeta C()} \over {{{(A({\rm{t}}) + 1)}^2}}} > \varpi ,vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),{\zeta \over {A({\rm{t}}) + 1}} > \eta and ϑC(𝔱) + ςA(𝔱) > vC(𝔱).

Proof.

At E₂, then Jacobian matrix of (1) is given by,

J(E₂) = [Jij]_i,j = 1, 2, 3. $\begin{matrix} J (E_{2}) = (\begin{matrix} - (ı + ϖ) & ρ_{1} + \frac{κ}{C (t) + 1} & ρ_{2} \\ ω - \frac{ζ C (t)}{{(A (t) + 1)}^{2}} & - v C (t) - ρ_{1} + ϑ C (t) & η - \frac{ζ}{A (t) + 1} \\ ς & v C (t) - ϑ C (t) - ς A (t) & - (ρ_{2} + η + ϰ) \end{matrix}) \\ A_{1} = - (a_{11} + a_{22} + a_{33}), A_{2} = - a_{22} a_{33} + a_{23} a_{32} - a_{11} a_{33} + a_{31} a_{13} - a_{11} a_{22} + a_{21} a_{12}, \\ A_{3} = - a_{11} a_{22} a_{33} + a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} - a_{31} a_{21} a_{32} + a_{13} a_{31} a_{22}, \end{matrix}$ \matrix{ {J\left( {{E_2}} \right) = \left( {\matrix{ { - (\imath + \varpi )} & {{\rho _1} + {\kappa \over {C({\rm{t}}) + 1}}} & {{\rho _2}} \cr {\omega - {{\zeta C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}}} & { - vC({\rm{t}}) - {\rho _1} + \vartheta C({\rm{t}})} & {\eta - {\zeta \over {A({\rm{t}}) + 1}}} \cr \zeta & {vC({\rm{t}}) - \vartheta C({\rm{t}}) - \zeta A({\rm{t}})} & { - \left( {{\rho _2} + \eta + } \right)} \cr } } \right)} \cr {{{\cal A}_1} = - \left( {{a_{11}} + {a_{22}} + {a_{33}}} \right),{{\cal A}_2} = - {a_{22}}{a_{33}} + {a_{23}}{a_{32}} - {a_{11}}{a_{33}} + {a_{31}}{a_{13}} - {a_{11}}{a_{22}} + {a_{21}}{a_{12}},} \cr {{{\cal A}_3} = - {a_{11}}{a_{22}}{a_{33}} + {a_{11}}{a_{23}}{a_{32}} - {a_{12}}{a_{21}}{a_{33}} + {a_{12}}{a_{23}}{a_{31}} - {a_{31}}{a_{21}}{a_{32}} + {a_{13}}{a_{31}}{a_{22}},} \cr } when we analyze the components of the Jacobian matrix, we notice that a₁₁ and a₃₃ have strictly negative signs, whereas a₂₁, a₂₂, a₂₃ and a₃₂ have negative sign provided that $\frac{ζ C (t)}{{(A (t) + 1)}^{2}} > ϖ, v C (t) + ρ_{1} > ϑ C (t), \frac{ζ}{A (t) + 1} > η$ {{\zeta C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}} > \varpi ,vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),{\zeta \over {A({\rm{t}}) + 1}} > \eta and ϑC(𝔱) + ςA(𝔱) > vC(𝔱) respectively.

Under the specified criteria, A₁, A₃ and A₁A₂ > A₃ are all positive signals.

Now, let $A_{2} (Y) = Y^{3} + A_{1} Y^{2} + A_{2} Y + A_{3} .$ {{\cal A}_2}() = {^3} + {{\cal A}_1}{^2} + {{\cal A}_2} + {{\cal A}_3}.

Thus, this polynomial’s discriminant is $D (A_{2}) = 18 A_{1} A_{2} A_{3} + {(A_{1} A_{2})}^{2} - 4 A_{3} A_{1}^{3} - 4 A_{2}^{3} - 27 A_{3}^{2} .$ D\left( {{{\cal A}_2}} \right) = 18{{\cal A}_1}{{\cal A}_2}{{\cal A}_3} + {\left( {{{\cal A}_1}{{\cal A}_2}} \right)^2} - 4{{\cal A}_3}{\cal A}_1^3 - 4{\cal A}_2^3 - 27{\cal A}_3^2.

Matignon’s criteria and the Routh-Hurwitz criterion state that the equilibrium point E₂ is locally asymptotically stable if $\frac{ζ C (t)}{{(A (t) + 1)}^{2}} > ϖ, v C (t) + ρ_{1} > ϑ C (t), \frac{ζ}{A (t) + 1} > η, ϑ C (t) + ς A (t) > v C (t)$ {{\zeta C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}} > \varpi ,vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),{\zeta \over {A({\rm{t}}) + 1}} > \eta ,\vartheta C({\rm{t}}) + \varsigma A({\rm{t}}) > vC({\rm{t}}) and D₂(A₂) > 0.

Theorem 6.

The steady state point E₃ is conditionally asymptotically stable if the following conditions are satisfied $\frac{κ B (t)}{{(C (t) + 1)}^{2}} > ρ_{2}, v C (t) + ρ_{1} > ϑ C (t), v B (t) + \frac{ζ}{A (t) + 1} > ϑ B (t) + η, ϑ C (t) + \frac{ς A (t)}{{(B (t) + 1)}^{2}} > v C (t)$ {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}} > {\rho _2},vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),vB({\rm{t}}) + {\zeta \over {A({\rm{t}}) + 1}} > \vartheta B({\rm{t}}) + \eta ,\vartheta C({\rm{t}}) + {{\zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}} > vC({\rm{t}}) and ρ₂ + ϑB(𝔱) + η + ϰ > vB(𝔱).

Proof.

At E₃, Jacobian matrix for (1), we have $J (E_{3}) = (\begin{matrix} - (ı + ϖ) & ρ_{1} + \frac{κ}{C (t) + 1} & ρ_{2} - \frac{κ B (t)}{{(C (t) + 1)}^{2}} \\ ϖ - \frac{ζ C (t)}{{(A (t) + 1)}^{2}} & - v C (t) - ρ_{1} + ϑ C (t) & - v B (t) + ϑ B (t) + η - \frac{ζ}{A (t) + 1} \\ \frac{ς}{B (t) + 1} & v C (t) - ϑ C (t) - \frac{ς A (t)}{{(B (t) + 1)}^{2}} & v B (t) - ρ_{2} - ϑ B (t) - η - ϰ \end{matrix})$ J\left( {{E_3}} \right) = \left( {\matrix{ { - (\imath + \varpi )} & {{\rho _1} + {\kappa \over {C({\rm{t}}) + 1}}} & {{\rho _2} - {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}}} \cr {\omega - {{\zeta C({\rm{t}})} \over {{{(A({\rm{t}}) + 1)}^2}}}} & { - vC({\rm{t}}) - {\rho _1} + \vartheta C({\rm{t}})} & { - vB({\rm{t}}) + \vartheta B({\rm{t}}) + \eta - {\zeta \over {A({\rm{t}}) + 1}}} \cr {{\zeta \over {B({\rm{t}}) + 1}}} & {vC({\rm{t}}) - \vartheta C({\rm{t}}) - {{\zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}}} & {vB({\rm{t}}) - {\rho _2} - \vartheta B({\rm{t}}) - \eta - } \cr } } \right) where $\begin{array}{l} Y^{3} + A_{1} Y^{2} + A_{2} Y + A_{3} = 0 \\ A_{1} = - (a_{11} + a_{22} + a_{33}), A_{2} = - a_{11} a_{33} + a_{23} a_{32} - a_{11} a_{33} + a_{31} a_{13} - a_{11} a_{22} + a_{21} a_{12}, \\ A_{3} = - a_{11} a_{22} a_{33} + a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} - a_{31} a_{21} a_{32} + a_{13} a_{31} a_{22} . \end{array}$ \matrix{ {{^3} + {{\cal A}_1}{^2} + {{\cal A}_2} + {{\cal A}_3} = 0} \hfill \cr {{{\cal A}_1} = - \left( {{a_{11}} + {a_{22}} + {a_{33}}} \right),{{\cal A}_2} = - {a_{11}}{a_{33}} + {a_{23}}{a_{32}} - {a_{11}}{a_{33}} + {a_{31}}{a_{13}} - {a_{11}}{a_{22}} + {a_{21}}{a_{12}},} \hfill \cr {{{\cal A}_3} = - {a_{11}}{a_{22}}{a_{33}} + {a_{11}}{a_{23}}{a_{32}} - {a_{12}}{a_{21}}{a_{33}} + {a_{12}}{a_{23}}{a_{31}} - {a_{31}}{a_{21}}{a_{32}} + {a_{13}}{a_{31}}{a_{22}}.} \hfill \cr }

Sign of element a₁₁ strictly negative, while the sign of elements a₁₂, a₂₁ and a₃₁ strictly positive. The sign of elements a₁₃, a₂₂, a₂₃, a₃₂ and a₃₃ are negative provided that $\frac{κ B (t)}{{(C (t) + 1)}^{2}} > ρ_{2}, v C (t) + ρ_{1} > ϑ C (t), v B (t) + \frac{ζ}{A (t) + 1} > ϑ B (t) + η, ϑ C (t) + \frac{ς A (t)}{{(B (t) + 1)}^{2}} > v C (t)$ {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}} > {\rho _2},vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),vB({\rm{t}}) + {\zeta \over {A({\rm{t}}) + 1}} > \vartheta B({\rm{t}}) + \eta ,\vartheta C({\rm{t}}) + {{\zeta A({\rm{t}})} \over {{{(B({\rm{t}}) + 1)}^2}}} > vC({\rm{t}}) and ρ₂ + ϑB(𝔱) + η + ϰ > vB(𝔱), respectively.

Theorem 7.

The steady state point E₄ is conditionally asymptotically stable if the following conditions are satisfied $\frac{κ B (t)}{{(C (t) + 1)}^{2}} > ρ_{2}, ζ C (t) > ϖ, v C (t) + ρ_{1} > ϑ C (t), v B (t) + ζ > ϑ B (t) + η, ϑ C (t) > v C (t)$ {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}} > {\rho _2},\zeta C({\rm{t}}) > \varpi ,vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),vB({\rm{t}}) + \zeta > \vartheta B({\rm{t}}) + \eta ,\vartheta C({\rm{t}}) > vC({\rm{t}}) and ρ₂ + ϑB(𝔱) + η + ϰ > vB(𝔱).

Proof.

At E₄, Jacobian matrix for (1), we have $J (E_{4}) = (\begin{matrix} - (ı + ϖ) & ρ_{1} + \frac{κ}{C (t) + 1} & ρ_{2} - \frac{κ B (t)}{{(C (t) + 1)}^{2}} \\ ϖ - ζ C (t) & - v C (t) - ρ_{1} + ϑ C (t) & - v B (t) + ϑ B (t) + η - ζ \\ \frac{ς}{B (t) + 1} & v C (t) - ϑ C (t) & v B (t) - ρ_{2} - ϑ B (t) - η - ϰ \end{matrix})$ J\left( {{E_4}} \right) = \left( {\matrix{ { - (\imath + \varpi )} & {{\rho _1} + {\kappa \over {C({\rm{t}}) + 1}}} & {{\rho _2} - {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}}} \cr {\varpi - \zeta C({\rm{t}})} & { - vC({\rm{t}}) - {\rho _1} + \vartheta C({\rm{t}})} & { - vB({\rm{t}}) + \vartheta B({\rm{t}}) + \eta - \zeta } \cr {{\zeta \over {B({\rm{t}}) + 1}}} & {vC({\rm{t}}) - \vartheta C({\rm{t}})} & {vB({\rm{t}}) - {\rho _2} - \vartheta B({\rm{t}}) - \eta - } \cr } } \right) $Y^{3} + A_{1} Y^{2} + A_{2} Y + A_{3} = 0$ {^3} + {{\cal A}_1}{^2} + {{\cal A}_2} + {{\cal A}_3} = 0 where $A_{1} = - (a_{11} + a_{22} + a_{33}), A_{2} = - a_{11} a_{33} + a_{23} a_{32} - a_{11} a_{33} + a_{31} a_{13} - a_{11} a_{22} + a_{21} a_{12},$ {{\cal A}_1} = - \left( {{a_{11}} + {a_{22}} + {a_{33}}} \right),{{\cal A}_2} = - {a_{11}}{a_{33}} + {a_{23}}{a_{32}} - {a_{11}}{a_{33}} + {a_{31}}{a_{13}} - {a_{11}}{a_{22}} + {a_{21}}{a_{12}}, $A_{3} = - a_{11} a_{22} a_{33} + a_{11} a_{23} a_{32} - a_{12} a_{21} a_{33} + a_{12} a_{23} a_{31} - a_{31} a_{21} a_{32} + a_{13} a_{31} a_{22} .$ {{\cal A}_3} = - {a_{11}}{a_{22}}{a_{33}} + {a_{11}}{a_{23}}{a_{32}} - {a_{12}}{a_{21}}{a_{33}} + {a_{12}}{a_{23}}{a_{31}} - {a_{31}}{a_{21}}{a_{32}} + {a_{13}}{a_{31}}{a_{22}}.

Sign of element a₁₁ strictly negative, while the sign of elements a₁₂ and a₃₁ strictly positive. The sign of elements a₁₃, a₂₁, a₂₂, a₂₃, a₃₂ and a₃₃ are negative provided that $\frac{κ B (t)}{{(C (t) + 1)}^{2}} > ρ_{2}, ζ C (t) > ϖ, v C (t) + ρ_{1} > ϑ C (t), v B (t) + ζ > ϑ B (t) + η, ϑ C (t) > v C (t)$ {{\kappa B({\rm{t}})} \over {{{(C({\rm{t}}) + 1)}^2}}} > {\rho _2},\zeta C({\rm{t}}) > \sigma ,vC({\rm{t}}) + {\rho _1} > \vartheta C({\rm{t}}),vB({\rm{t}}) + \zeta > \vartheta B({\rm{t}}) + \eta ,\vartheta C({\rm{t}}) > vC({\rm{t}}) and ρ₂ + ϑB(𝔱) + η + ϰ > vB(𝔱), respectively.

Numerical simulation

In this section, we implement a numerical simulation using RK4 to solve a system of ODEs representing the interactions between atmospheric water (A), surface water (B), and groundwater (C). The parameter values are taken from [1] and the values of κ, ζ and ς scale proportionally without causing unbounded effects, resulting in a stable and realistic model. The model incorporates key environmental factors such as precipitation, evaporation, infiltration, deforestation, pollution, and water pumping. Additionally, new terms are added to capture the dynamic interactions between these water bodies.

Through the numerical simulation, we validate the model by observing how these parameters influence the stability, equilibrium points, and overall behavior of the water system under various scenarios.

The proposed model is solved using the classical RK4 method due to its accuracy and stability for nonlinear systems. The model equations are implemented in MATLAB with a fixed time step h=0.01. To ensure robustness, selected results are cross-validated using a shifted Jacobi spectral method, confirming consistency with existing approaches such as those by Ebrahimzadeh et al. [32]. The parameters used in this simulation are summarized in Table 2, with values derived from the data presented in Table 1.

Table 2

Model parameterization with values.

Parameter	Parameter Description	Parameter Value
ρ₁&ρ₂	Surface and groundwater evaporation rates to atmospheric water, correspondingly	0.09 & 0.028
ϖ	Precipitation rate from atmosphere to surface water	0.5
κ	Influence of surface water on atmospheric water, modereted by groundwater	0.05
ι	Rate at which atmospheric water dissipates	0.01
v	Surface-to-groundwater infiltration rate	0.3005
ϑ	The effect of pollution on groundwater and surface water	0.300
η	Recurrence rate of groundwater extraction	0.02
ζ	Reduction in surface water due to groundwater and atmospheric interaction	0.07
ϰ	Pace of deforestation	0.06
ς	Contribution of atmospheric water to groundwater recharge, depending on surface water	0.04

Figure 1 shows how precipitation, dissipation, and evaporation from surface and groundwater transform atmospheric water levels throughout time. Figure 2 shows how precipitation, groundwater penetration, evaporation, and pollution affect surface water levels. Associations with groundwater have had a significant role in surface water changes. Figure 3 illustrates groundwater dynamics, including effects from infiltration, evaporation, deforestation, and water extraction. The mathematical model shows that, while air and surface water levels are declining owing to evaporation and pollution, groundwater levels remain steady. This stability can serve as a buffer against sudden changes in surface and atmospheric water. However, persistent losses in surface and atmospheric water might jeopardise long-term groundwater recharge and sustainability, especially with rising pollution and deforestation.

Results and discussion

The numerical simulations conducted for the proposed ODE model reveal key dynamics in the interaction between atmospheric water (A), surface water (S), and groundwater (C) compartments. The solutions show that under baseline conditions, the groundwater level tends to stabilize after initial fluctuations, indicating the system’s inherent stability. However, the inclusion of anthropogenic factors such as deforestation and frequent pumping shifts the equilibrium, resulting in long-term groundwater depletion. Among all parameters, the groundwater recharge rate and deforestation parameter had the most significant influence on groundwater dynamics, as demonstrated by the sharper decline in C over time when these values increased. These results emphasize the critical importance of preserving forest cover and regulating water extraction.

Conclusion

Groundwater is essential to preserving natural equilibrium and enabling human activity. However, it is increasingly threatened by factors like pollution, overconsumption, and deforestation. To address these challenges, we have developed an advanced ODE model that integrates new terms to represent the complex interactions between atmospheric water, surface water, and groundwater. Through the analysis of equilibrium points, stability, and numerical simulations, we validated the model’s accuracy across diverse scenarios. The numerical simulations, which are confirmed by the figures, disclose numerous important results. They illustrated that surface and atmospheric water levels decrease over time caused by pollution and evaporation, although groundwater remains rather steady in the near term. However, continuous environmental stressors gradually reduce groundwater recharge, highlighting underground water systems’ delayed but vital vulnerability. Compared to prior models that address various water sources in isolation or without environmental coupling, our model provides a more integrated and adaptable framework. This integration broadens its application in measuring water sustainability and implementing more effective groundwater management techniques.

The novelty and contribution of this paper are the following:

We extended an existing ODE-based water cycle model by introducing new nonlinear terms that reflect moderated interactions between atmospheric, surface, and groundwater.
The extended model captures complex interdependencies within the hydrological system more realistically than previous formulations.
We conducted equilibrium and stability analysis, and validated the model using numerical simulations under varying environmental conditions.

The model provides a foundational tool for assessing the long-term impacts of environmental changes and can be extended in future work for scenario analysis under varying climate and land-use conditions.

Groundwater depletion: A mathematical model incorporating climate and human factors

Full Article

Paradigm

My account