Mathematical modeling of groundwater security under climate change with optimal control

Introduction

Groundwater is defined as the water within underground formations (aquifers) that is critical for drinking and ecosystems, and its movement is simulated through mathematical models (Abd El Ghany et al., 2024) [1] . Climate change refers to global temperature increases and altered precipitation patterns affecting groundwater systems and broader environmental conditions as described by the Intergovernmental Panel on Climate Change (IPCC) (Davamani et al., 2024a) [10] . Globally, groundwater is facing increased pressure due to climate change, which affects its recharge, discharge, quality, and sustainability. In Illinois, United States, a 10% increase in population leads to a 5.5% decrease in groundwater levels (Asadollahi et al., 2024) [4] . Iran is highlighted as a global hot spot for groundwater depletion, with 90% of Iran negatively impacted by drought due to severe precipitation decline (Maghsoudi et al., 2023) [15] . Saudi Arabia faces challenges as groundwater has been affected by a depletion trend of -0.34 ± 0.01 mm/yr from 2002 to 2021 (Mohamed et al., 2022) [17] . Groundwater depletion rates in India are projected to triple by 2041-2080 due to increased withdrawals in response to warming temperatures, with current depletion rates expected to increase from 8.15 cm/year to 36.01 cm/year (Bhattarai et al., 2023) [6] . In Africa, climate change affects groundwater recharge through changes in precipitation, temperature, and sea level rise. In Egypt, maximum drawdown values range from 59 to 112 meters, equating to about 20-40% of the aquifer's saturated thickness (Shalby et al., 2024) [24] . Over 80% of South Africa is semi-arid to arid, with approximately 60% of total water withdrawals used for irrigation (Dalin, 2022) [9] . In northeastern Algeria, the Tebessa-Morsott alluvial aquifer is experiencing significant groundwater vulnerability, with about 97% of the area classified as high or average vulnerability (Mohammed et al., 2023) [18] . In East Africa, Ethiopia faces changes in precipitation ranging from -11.2% to -14.3% under RCP4.5 and from -9.2% to -10.0% under 1 RCP8.5 (Tilahun et al., 2023) [25] . Rwanda has experienced groundwater decline of approximately 16.4%, from 105.58 mm to 88.27 mm (Hakorimana et al., 2025) [13] . In south-west Tanzania, annual groundwater recharge to the shallow aquifer decreases by 26% under SSP 2-4.5 and 18% under SSP 5-8.5 (Mollel et al., 2023) [19] . In Kenya, the Nairobi Aquifer System (NAS) is experiencing groundwater level decline due to excessive abstraction, with an average annual recharge of 73 mm (9.7% of precipitation). Projections indicate a decrease in groundwater levels by 76 m by 2063 if abstraction rates remain constant, and by 14 m by 2030 if abstraction rates increase while recharge remains constant (Nyakundi et al., 2022) [20] . The evolution of groundwater hydrology from a purely descriptive empirical field to a rigorous domain of applied mathematics began with Henry Darcy's seminal formulations in 1856. Darcy's Law established that fluid velocity through porous media is directly proportional to the hydraulic gradient. When combined with the law of conservation of fluid mass, this relationship yielded the foundational partial differential equations (PDEs) governing transient groundwater flow in saturated aquifers (Bear, 2013) [5] . For distributed aquifer systems, these governing models expand into non-linear parabolic and elliptic second-order PDEs (Todd & Mays, 2004) [26] . While distributed approaches provide detailed grid-based spatial mapping, their immense parameter calibration needs and high computational costs make them poorly suited for real-time decision-making, predictive control, or dynamic tracking frameworks. To bypass these limitations, applied mathematicians developed lumped-parameter systems that utilize macroscopic water mass balance loops to collapse spatial dimensions into continuous temporal paths governed by ordinary differential equations (ODEs). To transition these predictive hydrogeological equations into active management tools, control-theoretic state-space principles must be applied. A continuous-time linear time-invariant state-space framework converts higher-order differential models into a coupled system of first-order vector-matrix differential equations structured as x(t) = Ax(t) + Bu(t) + Ed(t), where x(t) represents internal physical states, u(t) is the managed control input, and d(t) is the environmental disturbance vector (Ogata, 2010) [21] . Once a system is proven to be stable, completely controllable, and observable under the Kalman Rank Criterion, dynamic optimization via variational calculus can be introduced. Applying Pontryagin's Minimum Principle allows the formulation of an optimal statefeedback control policy by defining a scalar control Hamiltonian function that appends the system constraints directly to a quadratic cost index (Kirk, 2004) [14] . Current groundwater management is based on passive, static paradigms that do not account for "aquifer memory" due to delayed gravity drainage, leading to systemic over-extraction and irreversible basin collapse. Previous work uses either complex distributed simulations or simplified first-order models, but an important mathematical gap remains. No literature exists that couples a second-order lumped state-space differential layout with dynamic optimization. This study bridges this gap by formulating a continuous-time linear state-space model based on mass balance and applying Pontryagin's Minimum Principle to provide an automated, optimal extraction tracking law under active climate disturbances. The broad objective of this study was to develop a mathematical model to analyze the aggregate temporal dynamics of groundwater systems under climate variability and formulate optimal control strategies that ensure long-term groundwater security. The specific objectives were, to formulate a second-order lumped hydrogeological differential equation that captures climate-driven recharge variability and delayed gravity drainage based on water mass balance principles; to transform the governing second-order differential model into a continuous-time linear state-space architecture and evaluate its intrinsic structural properties using mathematical stability analysis; to determine the optimal groundwater extraction policy trajectory using Pontryagin's Minimum Principle subjected to a quadratic performance cost criteria, state constraints, and control bounds; and to numerically simulate the open-loop and closed-loop aquifer storage level trajectories under time-varying climate disturbance scenarios using MATLAB.

Research methodology 2.1 Introduction

This study adopted a quantitative research design grounded in continuous-time mathematical modeling, analytical matrix verification, and deterministic optimal control simulations. The design was selected to provide an objective, mathematically rigorous framework capable of tracking and optimizing groundwater storage trajectories under non-stationary environmental anomalies. Rather than relying on computationally heavy distributed partial differential equations (PDEs), this research developed a low-dimensional lumped-parameter architecture that aggregates localized spatial parameters into continuous temporal paths. The structural workflow executed in this research comprised four distinct mathematical milestones: 1. Formulation of a governing second-order hydrogeological ordinary differential equation (ODE) integrating fluid mass conservation loops with vertical drainage delay constraints. 2. Transformation of the structural layout into a continuous-time linear time-invariant (LTI) state-space architecture triplet (A,B,C). 3. Analytical valuation of the plant matrix properties, explicitly checking asymp-4. totic stability, complete state controllability, and complete observability using the Kalman Rank Criterion. 5. Synthesis of an automated feedback extraction schedule policy (u * ) via variational calculus using the Pontryagin Minimum Principle subjected to a multicriteria quadratic performance index.

Nomenclature and System Parameters

To preserve structural consistency across the mathematical derivation, state transformation, and numerical execution stages, the variables and constitutive parameters of the system were explicitly mapped as detailed in Table 2.1. Where S(t) represents the dynamic aquifer storage volume (m 3 ), Inet(t) is the net climate-driven volumetric recharge influx (m 3 /month), and Qout(t) is the anthropogenic volumetric extraction rate (m 3 /month). By normalizing Equation across the total surface basin area Aaq and modeling storage variations as a function of the matrix specific yield (Sy), the instantaneous response of the hydraulic head relative to net vertical boundary fluxes was formulated as:

Sy dy(t)/ dt = I(t)-EV(t)-u(t)

Where I(t) is precipitation infiltration recharge per unit area (m/month), EV(t) is thermal driven evapotranspiration loss (m/month), and u(t) is the managed extraction rate (m/month), which serves as the active system control variable. The net climate driving forces were condensed into a singular, continuous environmental disturbance forcing function:

Incorporation of Aquifer Memory Effects

Traditional first-order models assume an instantaneous drainage response across the water table . To account for realistic temporal lags induced by vertical gravity drainage through the unconfined unsaturated sediment matrix, this study incorporated Boulton's Delayed Yield Theory. According to this hydrogeological framework, the release of water from storage is non-instantaneous and follows a characteristic temporal delay index τ (months). This drainage behavior couples the instantaneous hydraulic head deviation x1(t) to the delayed yield response via a convolution integral layout, which simplifies to the following continuous firstorder differential linkage:

Where Rres represents the vertical hydraulic matrix resistance (month) and Kaq denotes the dimensionless natural restoring equilibrium constant of the aquifer formation.

Formulation of the Governing Second-Order Differential Equation

To capture higher-order temporal storage delays without using partial derivatives, Equation 4.3.4 was differentiated with respect to time t: Where the constitutive lumped physical parameters were explicitly mapped as:  M = τR - res 1 represents the Aquifer Memory Coefficient (month 2 ), quantifying vertical drainage lag.  Bres = R - res 1 represents the Matrix Damping Resistance Factor (month), resisting rapid shifts in drawdown velocity.  K = Kaq represents the Dimensionless Restoring Force Constant, governing the natural steady-state basin recovery potential.

State-Space Transformation

To transform the higher-order differential governing model into an active management framework, a continuous-time linear timeinvariant state-space representation was established. The state vector was defined as x(t) = [x1(t),x2(t)] T ∈R 2 , where the state variable components isolate the primary physical indicators:

𝑑𝑡

Differentiating the water table shift velocity component x2(t) with respect to time yielded the secondary state derivative expression:

Substituting Equations into the second-order model and isolating the state acceleration derivative mapped the direct resource depletion: Because monitoring networks track the physical water table heights directly at the surface, the measurable output tracking configuration was isolated. To avoid identity mapping conflicts, the observation vector was mapped as a scalar field tracking variable Thus, the total continuous system dynamics were compactly expressed via the standard LTI state triplet matrices (A,B,C) alongside the non-homogeneous weather disturbance mapping vector E:

Mathematical Stability Analysis

To guarantee that the underlying hydrogeological basin was inherently stable before introducing external optimization rules, the asymptotic stability of the open-loop system was analytically evaluated. The intrinsic stability properties of any continuous LTI system are dictated by the eigenvalues (λ) of its plant matrix A, requiring all real parts to sit strictly within the negative Left-Half Plane (LHP). The characteristic equation was generated by solving the matrix determinant equation det(λI -A) = 0:

Applying the quadratic formula, the exact system roots (poles λ1,2) were analytically solved as: Bres p4KM -B2 Since ℜ(λ1,2) < 0 holds true across all admissible physical parameter configurations, the unmanaged groundwater system was proven to be globally asymptotically stable.

Controllability and Observability Verification

To verify the structural mathematical authority of the tracking configurations prior to synthesis, the structural layout was subjected to the Kalman Rank Criterion.

Controllability Verification

Complete state controllability determines whether the managed pumping input u(t) possesses the mathematical rank authority to guide the internal aquifer states from any arbitrary initial disturbed level to a safe target within a finite time horizon. The algebraic controllability matrix Co ∈R 2×2 for the matrix pair (A,B) was constructed as:

Substituting the parameter matrix equations yielded:

To evaluate row rank properties, the matrix determinant was computed:

Since the Aquifer Memory Coefficient must be strictly non-zero (M > 0), the determinant is non-vanishing (det(Co) ̸ = 0). Thus, rank (Co) = 2 = n, proving that the system is completely state controllable.

Observability Verification

Complete observability determines whether measuring a single surface output sequence-the water table head deviation deviation x1(t)-provides enough information to completely reconstruct the hidden internal drawdown velocity state x2(t). The algebraic observability matrix Ob ∈R 2×2 for the matrix pair (A,C) was defined as:

Substituting the constitutive matrices yielded:

Evaluating the rank of the resulting identity matrix yielded:

Since det (Ob)̸ = 0, the matrix possesses full column rank (rank(Ob)= 2 = n). The system was verified as completely observable, ensuring tracking surface head deviations provides sufficient structural information to state-reconstruct internal basin storage velocity trajectories.

Optimal Control Formulation via Pontryagin's Minimum Principle

To resolve the groundwater over-extraction crisis under active climate disturbances, a dynamic optimization layer was synthesized using the necessary conditions of the Pontryagin Minimum Principle.

Physical and Economic Rationale of the Cost Functional

A continuous quadratic performance cost index functional J was introduced across a finite operational planning horizon tf:

Where the positive semi-definite state coefficients (q1,q2 ≥ 0) and the strictly positive control penalty factor (r > 0) serve as adjustable design metrics. The physical logic governing this quadratic configuration is defined as follows:

 The expression q1x1 2 (t) quadratically penalizes large deviations of the water table from its target target sustainable baseline, ensuring severe operational structural deficits are penalized more heavily than minor tracking errors.  The expression q2x2 2 (t) penalizes high water table shift velocities, suppressing abrupt drawdown spikes that trigger irreversible land subsidence.  The expression ru 2 (t) minimizes operational extraction pumping expenditures, matching the real-world non-linear energy consumption curves associated with lift requirements.

The Control Hamiltonian and Optimality Tracking Law

To convert this constrained optimization problem into an unconstrained system, a timevarying vector of co-state variables p(t) = (p1(t), p2(t)) T ∈R 2 was introduced. The scalar control Hamiltonian function H was constructed by appending the state-space derivative dynamics directly to the continuous integrand cost function:

Adjoint System and Transversality Conditions

To resolve the temporal paths of the co-state multiplier structures, the adjoint differential equations were derived by taking the negative partial derivatives of the control Hamiltonian with respect to the physical state variables: This coupled auxiliary system was structured in backward vector-matrix notation as: window. To preserve mathematical continuity with the state-space differential equations, raw daily millimeter precipitation fluxes (mm/day) were processed into absolute monthly cumulative meters (m/month). The complete, chronological 60-month climate forcing dataset utilized to drive the simulation models is presented in Table 4

Open-Loop Structural Verification and Complex Plane Mapping

To check the operational response pathways before applying optimal tracking frameworks, a unit impulse transmission scenario was simulated. This verification stage isolates how basic environmental pulses circulate through the coupled ordinary differential structure over a baseline observation timeline.

Synthesized Optimal Pumping Schedule Policy Analysis

The optimization engine generated a dynamic, state-dependent extraction profile by evaluating the real-time interaction between internal storage layers and atmospheric forcing vectors. Figure 4.5 charts the synthesized optimal pumping schedule policy trajectory (u * ) executed over the 60-month management window. The simulated pumping profile demonstrates clear adaptive feedback properties, contrasting sharply with rigid, static traditional water allocation quotas. The controller automatically expands extraction rates up to a peak near 0.18 m/month during extreme thermal anomalies and severe climate-driven recharge deficits to satisfy essential urban consumer demands without compromising localized hydraulic safety thresholds. Conversely, during high-volume precipitation events, the schedule automatically throttles extraction down to conserve mass and maximize natural vertical infiltration. As the timeline reaches the final operational boundary (t → 60 months), the pumping rate executes a smooth convergence to zero, driven by the structural enforcement of the terminal transversality boundary conditions.

Pontryagin Co-State Multiplier Dynamics

The underlying mathematical mechanism driving the feedback control loop was evaluated by tracking the temporal pathways of the auxiliary co-state multipliers. The first co-state variable (p1(t)) registers the marginal shadow cost of hydraulic head displacement drift, while the secondary co-state variable (p2(t)) registers the marginal shadow cost of drawdown velocity variations. Figure 4.6 maps the continuous trajectories of both parameters over the 60-month planning horizon. The numerical simulation shows that both auxiliary states fluctuate in perfect synchronization with the volatile, non-homogeneous NASA POWER climate forcing data sequence. The co-state dynamics serve as the core mathematical driver for the adaptive extraction law, passing real-time resource scarcity indicators straight to the tracking loop via the analytical connection u * (t) = rM 1 p2(t). As the system approaches the terminal planning window, both tracking indices smoothly converge to zero (p1(60) = 0 and p2(60) = 0). This precise convergence confirms complete compliance with the standard transversality conditions required by calculus of variations, verifying the mathematical validity of the closed-loop optimization layer.

Interpretation

The empirical findings from the continuous-time linear time-invariant state-space simulation reveal critical structural characteristics of the Nairobi Aquifer System. The open-loop plant matrix eigenvalues (λ1 =-0.011767, λ2 =-0.623167) provide definitive mathematical insight into the intrinsic behavior of the physical formation. Because both calculated poles sit strictly on the real axis within the negative Left-Half Plane, the system operates within a globally asymptotically stable, non-oscillatory regime. The resulting damping ratio (ζ = 3.71) indicates a heavily overdamped physical environment, which confirms that the aquifer naturally absorbs sudden pressure impulses without experiencing water table bouncing or fluid pressure shocks. However, the large gap between the fast pole (λ2) and the slow pole

Outcomes

The deployment of the continuous-time linear state-space architecture and variational optimization layer yielded four distinct analytical and empirical outcomes:

Summary of Findings

The quantitative investigations executed across this study yielded several critical structural insights aligned with the specific research objectives: 1. Model Formulation: By incorporating Boulton's Delayed Yield Theory into standard mass balance continuity equations, a low-dimensional governing secondorder lumped ordinary differential equation was successfully derived. This model effectively captures vertical gravity drainage delay ("aquifer memory") without incurring the high parameter calibration burdens typical of distributed partial differential architectures. 2. Structural Properties: The transformation of the hydrogeological governing model into a continuous-time linear timeinvariant (LTI) state-space triplet matrix formulation allowed for rigorous algebraic validation. The plant matrix was proven to be globally asymptotically stable, completely state controllable 3. (rank(Co) = 2), and completely state observable (rank(Ob) = 2) via the Kalman Rank Criterion. 4. Eigenvalue Mapping: Parameterization using site-specific Nairobi Aquifer data mapped the system poles strictly inside the stable Left-Half Plane (λ1 =-0.011767, λ2 = -0.623167). The open-loop system exhibits an overdamped physical behavior (ζ = 3.71), proving that basic impulse inputs and environmental pulses decay smoothly without generating high-frequency oscillations. 5. Dynamic Optimization: Applying Pontryagin's Minimum Principle to a multicriteria quadratic cost functional successfully resolved the groundwater overextraction crisis. The explicit tracking control law (u * (t) = rM 1 p2(t)) establishes that sustainable real-time resource allocation is strictly governed by tracking the marginal shadow cost dynamics of the secondary aquifer pressure co-state vari-able.

Conclusions

Based on the empirical calibration and analytical matrix verifications, the following academic conclusions are established:  Insufficiency of Memoryless Paradigms: Traditional first-order groundwater models assume instantaneous aquifer drainage responses. This study conclusively demonstrates that neglecting the vertical drainage delay parameter (M) leads to systemic parameter bias, resulting in over-extraction patterns that accelerate water table collapse under climate shocks.  The Role of Pumping Penalties: The synthesis of the automated state-feedback tracking law proves that long-term basin stabilization is highly sensitive to the control expenditure weight (r). Tuning the penalty factor forms a critical mathematical lever to regularize mechanical extraction rates against real-world energylift consumption curves.  Active Closed-Loop Superiority: While open-loop unmanaged forecasting simulations show a steep, persistent structural drift toward aquifer failure thresholds, the closed-loop optimization framework dynamically balances current agricultural abstractions with future storage security.

Recommendations

To transition these mathematical proofs into actionable regional resource frameworks, the following recommendations are proposed:

Policy and Management Interventions

 Transition to Dynamic Allocation: Water regulatory authorities within the Nairobi Aquifer boundary should phase out static, historical fixed pumping quotas. They should adopt dynamic, state-dependent extraction paths that automatically scale based on real-time observations of the water table displacement deviation (x1).  Institutional Integration of Co-State Multipliers: Operational extraction schedules for regional wellfields should be calibrated using the derived secondary aquifer pressure co-state trajectory (p2(t)). This ensures that pumping commands are regularized to protect the system against irreversible land subsidence during prolonged dry spells.

Recommendations for

References

1. Abd El Ghany MF, et al. Groundwater modeling for sustainable management. Hydrological Sciences Journal. 2024.
2. Afrifa GA, et al. Mathematical models for groundwater level prediction under climate change. Journal of Hydrology. 2022.
3. Afrifa GA, et al. Groundwater security and climate change adaptation strategies. Water Resources Management. 2022.
4. Asadollahi M, et al. Impacts of climate change and urbanization on groundwater levels in Illinois. Journal of the American Water Resources Association. 2024.
5. Bear J. Dynamics of Fluids in Porous Media. New York: Dover Publications; 2013.
6. Bhattarai N, et al. Projected tripling of groundwater depletion rates in India under warming temperatures. Nature Climate Change. 2023.
7. Boufekane A, et al. Groundwater quality decline in Khartoum State due to over-pumping. Environmental Earth Sciences. 2022.
8. Chepkemoi J, et al. Groundwater sustainability in Kericho County under population growth and climate change. African Journal of Environmental Science and Technology. 2023.
9. Dalin C. Irrigation water withdrawals and groundwater pressure in South Africa. Environmental Research Letters. 2022.
10. Davamani V, et al. Climate change impacts on groundwater systems: IPCC perspective. Current Climate Change Reports. 2024.
11. Davamani V, et al. Policy frameworks for climate change and groundwater management. Environmental Science & Policy. 2024.
12. Foster S, Tuinhof A. Groundwater Management in Arid and Semi-Arid Regions. Washington (DC): World Bank Publications; 2005.
13. Hakorimana E, et al. Groundwater decline due to urban expansion and agricultural intensification in Rwanda. Journal of African Earth Sciences. 2025.
14. Kirk DE. Optimal Control Theory: An Introduction. Mineola (NY): Dover Publications; 2004.
15. Maghsoudi M, et al. Groundwater depletion in the Yasouj aquifer under drought conditions. Iranian Journal of Hydrogeology. 2023.
16. Menichini M, et al. Groundwater modeling for future water supply security under climate change. Hydrogeology Journal. 2022.
17. Mohamed A, et al. Groundwater depletion trends in Saudi Arabia from GRACE satellite data. Remote Sensing of Environment. 2022.
18. Mohammed S, et al. Groundwater vulnerability assessment in the Tebessa-Morsott alluvial aquifer, Algeria. Journal of African Earth Sciences. 2023.
19. Mollel LL, et al. Climate change impacts on groundwater recharge in south-west Tanzania under SSP scenarios. Journal of Hydrology: Regional Studies. 2023.
20. Nyakundi RM, et al. Groundwater level projections for the Nairobi aquifer system under changing abstraction rates. African Journal of Science, Technology, Innovation and Development. 2022.
21. Ogata K. Modern Control Engineering. 5th ed. Upper Saddle River (NJ): Prentice Hall; 2010.
22. Perez M, et al. Groundwater depletion, reduced freshwater access, and increased inequity. Nature Water. 2024.
23. Perez M, et al. Policy interventions for sustainable groundwater management under climate change. Water Policy. 2024.
24. Shalby A, et al. Groundwater depletion due to uncontrolled agricultural pumping in Egypt. Journal of Hydrology. 2024.
25. Tilahun SA, et al. Climate change impacts on groundwater in Ethiopia under RCP scenarios. Journal of Water and Climate Change. 2023.
26. Todd DK, Mays LW. Groundwater Hydrology. 3rd ed. Hoboken (NJ): John Wiley & Sons; 2004.