A Computational Framework for Groundwater Flow and Solute Migration in Layered Aquifers

Abstract

Understanding the coupled dynamics of groundwater flow and salinity transport is crucial for the sustainable management of aquifer systems, particularly in irrigated and semi-arid regions where evaporation, recharge variability, and groundwater abstraction significantly alter hydrogeological regimes. In multilayer porous media, groundwater-level fluctuations and salt migration are strongly interdependent, as hydraulic gradients govern solute transport and salinity variations may influence flow behavior through density-dependent mechanisms. This study develops a nonlinear mathematical model for describing groundwater-level evolution and salinity transport in a two-layer porous medium composed of an upper phreatic aquifer and an underlying confined aquifer. The model incorporates filtration processes, interlayer hydraulic exchange, density-dependent effects, and external forcing factors, including surface recharge, evaporation, and pumping. For numerical implementation, the governing equations are discretized using a finite-difference scheme with central spatial approximations and an implicit Crank–Nicolson-type temporal formulation. A hybrid second-order time approximation is introduced for the phreatic-layer equation to enhance numerical smoothness, stability, and convergence behavior. The resulting tridiagonal algebraic systems are solved by the Thomas algorithm within an iterative quasi-linearization framework, ensuring computational efficiency and robustness for nonlinear coupled simulations. The simulation results demonstrate distinct dynamic responses in the two aquifer layers. The phreatic aquifer responds rapidly and with higher amplitude to external forcing, whereas the confined aquifer exhibits slower, smoother, and more inertial hydraulic and geochemical adjustments. Sensitivity analysis shows that the filtration coefficient, transmissivity, porosity, density-dependent parameters, surface fluxes, and pumping intensity are the key factors controlling groundwater-level variations and salinity redistribution. The proposed modeling framework provides a reliable computational tool for analyzing coupled groundwater–salinity processes and offers a scientifically grounded basis for groundwater monitoring, salinization risk assessment, and sustainable aquifer management in multilayer hydrogeological systems.