Soil-based water filtering devices can be described by models of viscous flow in porous media coupled with an advection–diffusion–reaction system modelling the transport of distinct contaminant species within water, and being susceptible to adsorption in the medium that represents soil. Such models are analysed mathematically, and suitable numerical methods for their approximate solution are designed. The governing equations are the Navier–Stokes–Brinkman equations for the flow of the fluid through a porous medium coupled with a convection-diffusion equation for the transport of the contaminants plus a system of ordinary differential equations accounting for the degradation of the adsorption properties of each contaminant. These equations are written in meridional axisymmetric form and the corresponding weak formulation adopts a mixed-primal structure. A second-order, (axisymmetric) divergence-conforming discretisation of this problem is introduced and the solvability, stability, and spatio-temporal convergence of the numerical method are analysed. Some numerical examples illustrate the main features of the problem and the properties of the numerical scheme.