<p>Accurately predicting regional groundwater flow behavior in layered, heterogeneous aquifers remains a central challenge in hydrogeology. This study presents an analytical model for steady-state groundwater flow in a stratified three-layer aquifer system with anisotropic, depth-dependent hydraulic conductivity. Unlike prior models, this approach incorporates both a water table and surface water bodies as top boundary conditions over a finite vertical domain, addressing key limitations of earlier semi-infinite or topography-constrained frameworks. The model employs an exponentially decaying conductivity profile within each layer and utilizes a separation-of-variables solution technique to derive an analytical expression for hydraulic head. Validations against the existing analytical solution and COMSOL-based numerical simulation demonstrate strong agreement. The framework further enables analysis of recharge/discharge zones, artesian conditions, stagnation points, and groundwater particle travel paths across varied hydrogeological scenarios. The model reveals several novel insights, including the emergence of dual recharge–discharge behavior beneath surface water bodies due to climate variability, fragmentation of the groundwater flow system, and the presence of stagnation zones influenced by boundary conditions and conductivity structures. Key findings highlight the significant influence of surface water body boundaries and vertical heterogeneity on flow paths and artesian conditions. Sensitivity analysis identifies vertical conductivity of upper-layer as key drivers of hydraulic head variation. Additionally, residence time distributions exhibit heavy tails, with the Fréchet distribution providing the best statistical fit. It fills a key knowledge gap by showing how aquifer stratigraphy governs subsurface water connectivity and enables long-distance contamination via preferential flow paths, intensified by climate-induced waterbody shrinkage.</p>