Transport of multicomponent electrolyte solutions in saturated porous media is affected by the electrostatic interactions between charged species. Such Coulombic interactions couple the displacement of the different ions in the pore water and remarkably impact mass transfer not only under diffusion, but also under advection‐dominated flow regimes. To accurately describe charge effects in flow‐through systems, we propose a multidimensional modeling approach based on the Nernst‐Planck formulation of diffusive/dispersive fluxes. The approach is implemented with a COMSOL‐PhreeqcRM coupling allowing us to solve multicomponent ionic conservative and reactive transport problems, in domains with different dimensionality (1‐D, 2‐D, and 3‐D), and in homogeneous and heterogeneous media. The Nernst‐Planck‐based coupling has been benchmarked with analytical solutions, numerical simulations with another code, and high‐resolution experimental data sets. The latter include flow‐through experiments that have been carried out in this study to explore the effects of electrostatic interactions in fully three‐dimensional setups. The results of the simulations show excellent agreement for all the benchmarks problems, which were selected to illustrate the capabilities and the distinct features of the Nernst‐Planck‐based reactive transport code. The outcomes of this study illustrate the importance of Coulombic interactions during conservative and reactive transport of charged species in porous media and allow the quantification and visualization of the specific contributions to the diffusive/dispersive Nernst‐Planck fluxes, including the Fickian component, the term arising from the activity coefficient gradients, and the contribution due to electromigration.

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 subcontinuum (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 subcontinua and to significantly different breakthrough curves.