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.