Water pollution is a significant threat for human health, particularly in developed countries. This study advances the mathematical understanding of WP transmission dynamics by developing a fractional–fractal derivative framework with non-singular kernels and the Mittage–Leffler function, which successfully preserves the non-local behavior of pollutants. The fractional–fractal derivatives in sense of the Atangana–Baleanu–Caputo formulation inherently captures the non-local and memory-dependent behavior of pollutant diffusion, addressing limitations of classical differential operators. A novel parameter, $\gamma$, is introduced to represent the recovery rate of water systems through treatment processes, explicitly modeling the bridge between natural purification mechanisms and engineered remediation efforts. Furthermore, this study establishes stability analysis, and the existence and uniqueness of the solution are established through fixed-point theory to ensure the mathematical stability of the system. Moreover, a numerical scheme based on the Newton polynomial is formulated, by obtaining significant simulations of pollution dynamics under various conditions. Graphical results show the effect of important parameters on pollutant evolution, providing useful information about the behavior of the system.

This research develops a novel coupled system of nonlinear hybrid stochastic fractional differential equations that integrates neutral effects, stochastic perturbations, and hybrid switching mechanisms. The system is formulated using the Atangana–Baleanu–Caputo fractional operator with a non-singular Mittag–Leffler kernel, which enables accurate representation of memory effects without singularities. Unlike existing approaches, which are limited to either neutral or hybrid stochastic structures, the proposed framework unifies both features within a fractional setting, capturing the joint influence of randomness, history, and abrupt transitions in real-world processes. We establish the existence and uniqueness of mild solutions via the Picard approximation method under generalized Carathéodory-type conditions, allowing for non-Lipschitz nonlinearities. In addition, mean-square Mittag–Leffler stability is analyzed to characterize the boundedness and decay properties of solutions under stochastic fluctuations. Several illustrative examples are provided to validate the theoretical findings and demonstrate their applicability.