In many heterogeneous formations, hydraulic conductivity decreases with depth due to compaction, cementation, or lithologic transitions. This study adopts a simple yet flexible reciprocal-linear depth function to represent such variations, parameterized with a decay coefficient that controls the rate of conductivity decline from the ground surface down to the aquifer base. Under the Dupuit–Forchheimer approximation for steady-state flow in heterogeneous unconfined aquifers, we develop a corresponding discharge potential that leads to a closed-form relationship between hydraulic head and spatial coordinates expressed through the Lambert W function. This transformation converts the original nonlinear head-based equation into a linear Laplace formulation which can be solved by standard analytical or numerical methods. We show that the reciprocal-linear formulation is intrinsically less sensitive to uncertainty in the decay coefficient compared with the conventional exponential model, resulting in improved robustness to parameter uncertainty. The influence of the decay coefficient on flow topology is examined for three representative hydrogeological settings: a circular recharge basin, a pumping well under uniform regional flow, and a wedge-shaped aquifer subject to ambient gradients from boundary streams. For purely radial flow from a recharge basin, the analytical solution agrees closely with numerical counterpart from the Picard-iteration method. The analysis demonstrates that capture-zone geometry, the locus of stagnation point, and well discharge can be highly sensitive to the conductivity-decay coefficient