Low‐permeability aquitards, such as clay layers and inclusions, are of utmost importance for contaminant transport in groundwater systems. Although most dissolved species, contaminants, and clay surfaces are charged, the role of electrostatic interactions in subsurface flow‐through systems has not been extensively investigated. This study presents a two‐dimensional multicomponent reactive transport investigation of diffusive/dispersive and electrostatic processes in homogeneous and heterogeneous clay systems. The proposed approach is based on multiple continua and is capable to accurately describe charge interactions during ionic transport in the free water, diffuse layer, and interlayer water of charged porous media. The diffuse layer composition is simulated by considering a mean electrostatic potential following Donnan approach, whereas the interlayer composition is calculated by adopting the Gaines‐Thomas convention. Diffusive/dispersive fluxes within each sub‐continuum (free water, diffuse layer, and interlayer) are calculated solving the Nernst‐Planck equation while maintaining a net zero‐charge flux. Furthermore, the multidimensional flow and transport model is coupled with the geochemical code PHREEQC, by utilizing the PhreeqcRM module, thus enabling great flexibility to access all PHREEQC's reaction capabilities. The code is benchmarked in 1‐D systems against other software and previously published experimental data. Successively, reactive transport simulations are performed in 2‐D clayey‐sandy flow‐through domains with spatially variable physical and electrostatic properties at both laboratory and field scales. The results reveal that different properties of surface charge, diffuse layer, and Coulombic interactions impact the transport of charged species and lead to distinct spatial distribution of the ions in the different sub‐continua and to significantly different breakthrough curves.