Dimensionality effects on multicomponent ionic transport and surface complexation in porous media

Coulombic interactions between charged species in pore water and at surface/solution interfaces are of pivotal importance for multicomponent reactive transport in porous media. In this study, we investigate the impact of domain dimensionality on electrostatically coupled dispersion and surface-solution reactions during transport of acidic plumes and major ions in porous media. Column and quasi two-dimensional ﬂow-through experiments were performed, with identical silica porous media and under the same advection-dominated conditions. Equal mass ﬂuxes of diﬀerent electrolyte solutions (i.e., HCl - pH (cid:1) 2.8, NaBr - 100 mM, HCl - pH (cid:1) 2.8 plus NaBr - 100 mM) were continuously injected in the 1-D and 2-D experiments and breakthrough curves of pH and major ions were measured at the outlet of the domains. The presence of pronounced ionic strength gradients in the transverse direction in the 2-D setup caused distinct retardation and transport behaviors of protons and major ions which were not observed in the one-dimensional column experiments. Furthermore, in the cases of salt electrolytes injection, considerably enhanced release of H + (>61%) from the quartz surface was observed in the multidimensional system compared to the one-dimensional setup. Reactive transport modeling was performed to reproduce the experimental outcomes and to analyse the coupling between transport processes, based on the Nernst-Planck formulation of diﬀusive/dispersive ﬂuxes and on surface complexation reactions at the solid-solution interface. Electrostatic interactions between Na + , Br (cid:3) , and H + , and deprotonation of the quartz surface upon the formation of sodium outer-sphere complexes, are the primary controllers of the spatial and temporal features displayed by the pH and major ions measurements. The reactive transport simulations allowed us to interpret the experimental observations, to visualize the distribution and spatio-temporal evolution of dissolved and solid species, to identify a spatially heterogeneous zonation of Coulombic interactions with distinct behavior at the fringe and core of the injected plumes in the multidimensional setup, and to quantify the diﬀerent components of the Nerst-Planck ﬂuxes of the charged solutes. This study demonstrates that the domain dimensionality directly aﬀects electrostatic interactions between charged aqueous species in the pore water and surface complexation reactions at the solid-solution interface. The non-trivial eﬀects of dimensionality on multicomponent ionic transport result in a signiﬁcantly diﬀerent behavior in 1-D and 2-D systems. (cid:1)


INTRODUCTION
Coulombic interactions between charged aqueous species and reactive processes such as dissolutionprecipitation and sorption mechanisms play key roles in the transport of solute plumes in porous media. The study of multicomponent ionic transport is of pivotal importance for many geochemical systems. From the early work of Vinograd and McBain (1941), the role of coupled diffusive fluxes has been increasingly recognized for a variety of natural systems such as seawater ionic mixtures (Felmy and Weare, 1991), fluid-sediment interactions in hydrothermal circulation (Giambalvo et al., 2002), transport of major ions and charged contaminants (Liu et al., 2011;, diffusion in clay rocks including potential repositories for nuclear waste storage (Appelo and Wersin, 2007;Appelo et al., 2008;Appelo et al. 2010;Soler et al., 2019), uranyl diffusion, and displacement of stable and radiogenic isotopes (Druhan et al., 2015).
The coupling of geochemical reactions and mass transport mechanisms is essential to understand the release and transport of solute plumes in subsurface systems. Reactive transport modeling is a well-established approach that quantitatively integrates the joint effects of physical, chemical, electrostatic and biological processes in porous media (Steefel and Lasaga, 1994;Steefel et al., 2015;Li et al., 2016;Tournassat and Steefel, 2019). Bench-scale laboratory experiments are standard methods used to validate reactive transport models and a common approach for subsurface investigation is to perform one-dimensional column experiments. Such approach has the practical advantage of feasibility and relative simplicity of the flow-through experiments but does not allow capturing important processes such as transverse mixing at the interface between solutions with different water quality. Lateral dispersion determines plume dilution (Kitanidis, 1994;Rolle et al., 2013a;Chiogna and Rolle, 2017) and impacts many mixing-controlled biogeochemical reactions at the fringes of groundwater plumes (Rolle et al. 2008;Tartakovsky et al., 2008;Willingham et al., 2008;Bauer et al., 2009;Prommer et al., 2009;Rolle and Le Borgne, 2019;Valocchi et al., 2019). However, little attention has been dedicated to the investigation of the non-intuitive impact of dimensionality and lateral mixing processes on transport of charged species and on their apparently non-mixingcontrolled interactions with charged surfaces.
In this study, we explore the role of system dimensionality for transport of charged species in porous media and its impact on Coulombic interactions and surface complexation reactions. We analyse the transport of acidic plumes and major ions, which is important for a variety of geochemical processes such as soil/groundwater acidification in aquifers with low pH buffering capacity (Hansen and Postma, 1995;Kjøller et al., 2004;Fest et al., 2005), transformation and dissolution-precipitation of reactive minerals (Hansel et al., 2005;Maher et al., 2006;Battistel et al., 2019Battistel et al., , 2021, sorption-desorption of heavy metals onto the aquifer matrix (Zachara et al., 1991;Fuller et al., 1996;Kent et al., 2007), ion exchange (Appelo, 1994), redox reactions (Kocar and Fendorf, 2009;Fakhreddine et al., 2016) and acid mine drainage Muniruzzaman and Pedretti, 2020). We performed a series of flow-through experiments under steady-state flow and transient transport conditions in a 1-D column and a quasi 2-D chamber setups, designed with equal length and run with equal seepage velocity and inlet mass fluxes to allow a direct comparison of the experimental results. Different acidic and electrolyte solutions were continuously injected in the porous media and measurements of pH and major ions' breakthrough curves were performed at the outlet of the setups. Multicomponent reactive transport modeling, based on the Nernst-Planck formulation for the displacement of charged species and on surface complexation models describing the interactions with quartz and metal oxide surfaces, was performed to reproduce and interpret the experimental observations.

Experimental setups
The experimental setups used in this study are shown in Fig. 1 and comprise cylindrical glass columns (length: 30 cm, inner diameter: 1.75 cm) and a quasi twodimensional flow-through chamber. The latter was an acrylic tank with inner dimensions of 30 cm Â 20 cm Â 0. 8 cm (L Â H Â W), equipped with 24 equally spaced ports (5 mm spacing, numbered 1 to 24 from bottom to top throughout the text) at the inlet and outlet. The length of the setup, the injected mass flux, and the seepage velocity were maintained equal between the 1-D and quasi 2-D setups to allow a direct comparison of the experimental results.
All the experiments were conducted in a temperature controlled room (T = 20°C), and laboratory solutions were prepared using Milli-Q water. In both one-dimensional and quasi two-dimensional systems, flow-through experiments were performed by continuously injecting three different electrolyte solutions (i.e., six experiments in total): (i) acidic (HCl, pH $ 2.8), (ii) salt solution (NaBr, 100 mM), and (iii) salt acidic solution (HCl, pH $ 2.8 plus NaBr, 100 mM). Details about the performed flow-through experiments are reported in the Electronic Annex (Table EA1).
The 1-D and quasi 2-D flow-through setups were filled with the porous medium by following a wet-packing procedure, in which the water level was always maintained above the quartz sand to avoid entrapment of air in the pores in the water-saturated domain (Haberer et al., 2012). The porosity resulting from the wet-packing procedure was gravimetrically measured for the one-dimensional setup and further verified by simulating the breakthrough curves of conservative anions.
The inlet solutions were injected by establishing a uniform flow field using high-precision peristaltic pumps (IPC-N24, ColeParmer, United States) connected to the experimental setups by Fluran HCA pump tubing with an inner diameter of 0.64 mm. In all experiments, the seepage velocity was set to 1 m/d, and the flow rates sustained were measured at the end of each experiment.
All the column experiments were performed with continuous vertical injections from bottom to top using two ports at each end of the columns (two needles pierced through the rubber septum). The quasi two-dimensional flow-through experiments were performed by establishing a uniform horizontal flow field at a flow rate of 287.5 mL/min (i.e., 23 ports with 12.5 mL/min each). Five central inlet ports (i.e., ports 10-14) were used to inject the electrolyte solutions in the domain (i.e., corresponding to an injection width of 2.5 cm, marked with a yellow brace in Fig. 1a) with the same mass flux established during the column experiments. For instance, considering the more concentrated electrolyte mixture (i.e., NaBr 100 mM), the solute mass flux was equal to 0.01 mol/d in both setups. The remaining inlet ports were used to inject the background solution (i.e., Milli-Q water for the HCl case and an autoclaved solution of fructose for NaBr and NaBr plus HCl cases to avoid gravity effects).

Porous medium characterization and preparation
The porous medium used for the experiments was a quartz sand originated from Dorsten, Germany (Aquagran -dry, Euroquartz GmbH), with a characteristic grain size diameter equal to 1.1 mm, a solid density equal to 2.65 kg/L and a BET specific surface area of 0.2 m 2 /g sand (i.e., determined by single-and multi-point BET sorption with N 2 at À195°C). Fe and Al coatings were present as small fractions of the sand mass, <0.02 wt% (1.91 mmol/ g) and < 0.007 wt% (0.65 mmol/g), respectively, assuming an oxide molar mass of 100 g/mol (Stolze et al., 2020).
Before each experiment, the quartz porous medium was cleaned according to the method proposed by McNeece and Hesse (2016): first, the sand was washed in an HCl solution (pH $ 4) and thoroughly rinsed with Milli-Q water. The silica medium was then washed in a NaOH solution at pH = 10, cleaned with Milli-Q water until reaching a pH of $ 5.9 and EC < 0.1 mS/cm, and dried at 80°C. Furthermore, the silica sand was autoclaved to minimize microbial activity during the flow-through experiments. Finally, prior to the experiments, the porous media were flushed with an HCl (pH $ 2) solution for approximately 10 pore volumes to dissolve mineral phases present in the quartz sand, which could affect the electrolyte transport (Stolze et al., 2020). Subsequently, the domain was re-equilibrated with Milli-Q water until a pH $ 5.5 was reached.

pH and major ions sampling and analyses
Throughout the experiment, the pH breakthrough was measured by an in-line pH monitoring system. One channel of the peristaltic pump at the outlet was continuously delivering the effluent to a flow-through vial containing a constant volume of solution (i.e., 10 mL) in which a pH electrode performed pH measurements with a time interval of 5 minutes. The effluent H + concentrations were calculated by correcting the measured pH values, removing the dilution effects in the flow-through vial by using the continuous stirred tank reactor (CSTR) equation (Section 2 in the Electronic Annex). Furthermore, samples of the effluent solution were collected at the outlet of the domain at 1hour interval for monitoring the breakthrough of major ions.
During the 1-D flow-through experiments, one outlet port was used for in-line pH breakthrough curve measurement, whereas the other port was used to sample the effluent. Similarly, during the multidimensional experiments, inline monitoring of the pH breakthrough curve at the outlet of the chamber (orange color in Fig. 1) was performed at three individual locations (ports 2, 7 and 12) spaced along the vertical direction and over the upper half of the flowthrough setup by connecting ports 13-23 to one single flow-through vial in order to measure the integrated breakthrough curve. Throughout the text, the data recorded at these four different locations are referred to as bottom, intermediate, central, and integrated, respectively. The breakthrough curves of major ions were monitored from two areas of the chamber outlet: (i) an inner (plume core) location next to the center of the domain (ports 9 and 10), and (ii) an outer (plume fringe) location closer to the bottom of the chamber (ports 4 and 5; red color in Fig. 1).
The in-line pH measurements for both one-dimensional and quasi two-dimensional setups were performed by using Hach Intellical electrodes (PHC 281). The probes were calibrated before each experiment with standard buffer solutions of pH 4, 7 and 10.
Inductively Coupled Plasma Atomic Emission Spectroscopy (ICP-OES, PerkinElmer Avio 200) analyses were conducted on the acidified (2% HNO 3 ) samples to measure the breakthrough curves of Na + . Ion Chromatography (IC, Dionex ICS-5000, ThermoScientific) analyses were carried out on the collected samples to measure the breakthrough curves of the major anions, Br À and Cl À .

MODELING APPROACH
Reactive transport modeling was performed to quantitatively describe the processes governing the transport of ions and pH in the flow-through systems and to evaluate and compare the outcomes of the one-dimensional and quasi two-dimensional flow-through experiments. The simulations considered advective-dispersive transport of the injected electrolyte solutions, including the Coulombic interactions between the transported aqueous species, aqueous speciation, and the surface complexation reactions at the mineral-solution interface. The latter entailed a distinct description for the quartz surface and the metal oxide coatings.

Governing multicomponent ionic transport equations
The transport of charged species in the 1-D and quasi 2-D setups was described with multicomponent ionic transport equations accounting for the diffusive/dispersive properties of each aqueous species and Coulombic interactions between charged species. In the flow-through setups the diffusive/dispersive flux of a charged species i is expressed by the Nernst-Planck equation derived from the gradient of the electrochemical potential (Liu et al., 2011;Rasouli et al., 2015;Rolle et al., 2018): where D i accounts for the components of the hydrodynamic dispersion tensor, C i is the aqueous concentration of a species i, c i is the activity coefficient, R is the ideal gas constant, T is the absolute temperature, z i is the charge number, F is the Faraday's constant, and U is the electrostatic potential. The activity coefficient term can be recast as: and the Nernst-Planck equation can be further expressed as (Appelo and Wersin, 2007): In the absence of an external electric field P N i¼1 z i J i ¼ 0 À Á , the electrical potential gradient can be expressed as (Boudreau et al., 2004): By substituting Eq. (4) in the Nernst-Planck equation, the diffusive flux J i becomes: This equation highlights that the movement of a charged species in pore water depends on concentration gradients, dispersion coefficients, activity coefficients, and charge numbers of that species and all other ions in solution.
Eq. (5) can be formulated in a more compact way as: where D ij is the dispersion tensor, which, for a twodimensional system oriented along the principal flow direction, has two diagonal entries, D L ij and D T ij , representing the longitudinal and transverse cross-coupled dispersion matrices (Muniruzzaman et al., 2014): where d ij is the Kronecker delta function, equal to 1 when i ¼ j and equal to 0 if i-j. D L i and D T i are the longitudinal and transverse components of the hydrodynamic dispersion of species i, which were parameterized according to a linear (Guedes de Carvalho and Delgado, 2005) and nonlinear (Chiogna et al., 2010;Rolle et al., 2012) formulation, respectively: where D P i is the pore diffusion coefficient (i.e., approximated as the product between the porosity h and the self-diffusion coefficient D aq i , values provided in Table EA2), Pe i is the Péclet number defined as, , with d being the characteristic grain size diameter. d is the ratio between the length of a pore channel and its hydraulic radius, and b is an empirical exponent accounting for the incomplete mixing in the pore channels. The values of d and b were set to 5.37 and 0.5, respectively, based on the compilation of multidimensional flow-through experiments (Ye et al., 2015a).
Finally, the governing equation for multicomponent ionic transport in saturated porous media is obtained from a mass balance over a control volume and integrating the fluxes due to advection, multicomponent ionic dispersion, and gradient of activity coefficients: where S i is the adsorbed concentration of a species i, v is the seepage velocity vector, q s is the solid density, h is the porosity, R i is the reactive term accounting for aqueous speciation reactions, and t i is the stochiometric coefficient of the species i.

Surface complexation model
We simulated adsorption processes occurring during the flow-through experiments by applying a surface complexation model that considers both quartz and metal oxide surfaces. In a previous study (Stolze et al., 2020), the model was shown capable of simulating the protonation/deprotonation behavior of the sand surface for a wide range of ionic strength and NaBr electrolyte background concentration as considered in the present work. Indeed, in a multidimensional system it is of critical importance to adopt a model with enough flexibility to account for ionic strength effects on surface complexation due to the generation of steep concentration gradients caused by transverse dispersion.
The Basic Stern (BS) and Triple Layer (TLM) models were applied to simulate adsorption processes on quartz and metal oxides (i.e., assuming a generic surface to account for both Al and Fe coatings), respectively. Both surfaces were represented with an amphoteric behavior by defining one single type of surface site able to protonate and deprotonate. Surface complexation reactions with formation of Na + , Cl À , and Br À outer-sphere complexes were defined to balance the surplus/deficiency of surface charge in the 0-plane resulting from the protonation/deprotonation reactions. Furthermore, carbonate outer-sphere complexes were implemented in the TLM based on the work of Leckie (2000, 2001). All the complexation reactions and SCM parameters for the quartz model are reported in Table EA3.

Numerical implementation for 1-D and 2-D reactive transport modeling
The simulations of the 1-D column flow-through experiments were performed with the geochemical code PHREEQC-3 (Parkhurst and Appelo, 2013) coupled with MATLAB Ò by using the IPhreeqc module (Charlton and Parkhurst, 2011). This method allows the combination of the transport and geochemical capabilities of PHREEQC with the computational methods of MATLAB Ò (Muniruzzaman and Rolle, 2016;Stolze et al., 2019b).
Regarding the simulations of the two-dimensional flowthrough experiments, the multicomponent ionic transport and the chemical reactions were decoupled utilizing a sequential non-iterative operator splitting approach Rolle, 2016, 2019). The multicomponent ionic transport was solved with a finite volume method (FVM) on a streamline-oriented grid constructed with the approach developed by Cirpka et al. (1999). A Picard iterative scheme was utilized to linearize the nonlinear set of transport equations (Eq. (11)), and both the flow and multicomponent transport calculations were implemented in MATLAB Ò (Muniruzzaman et al., 2014). After the advection and dispersion steps, the geochemical reaction calculations were performed with PHREEQC by using the PhreeqcRM module, which allowed great flexibility in accessing PHREEQC's reaction capabilities (Parkhurst and Wissmeier, 2015;Jara et al., 2017;Rolle et al., 2018;Sprocati et al., 2019).
The one-dimensional domain was discretized into 300 cells (i.e., Dx = 1 mm) along the longitudinal direction, whereas the two-dimensional model domain was discretized into 50 (i.e., Dx = 6 mm) and 115 grid-cells (i.e., Dy = 1 mm) along the longitudinal and transverse dimensions, respectively (Table 1). The boundary conditions were determined based on pH and concentration measurements of the inlet solutions and the flow rate values measured prior to the experiments (Table EA1). Similarly, the initial conditions within the model domain were assigned by performing speciation calculations determining pore water and surface species based on the initial pH values measured at the beginning of the experiments at the domains outlet (Table EA1). The aqueous speciation reactions were calculated by using the thermodynamic database wateq4f. As displayed in Fig. 2a, the pH breakthrough curves, measured during the experiment performed by injecting an HCl solution in both 1-D and 2-D setups, showed a monotonic decrease in pH towards a steady-state plateau. Furthermore, the breakthrough curves of H + highlight the importance of sorption reactions at the metal oxides surface (Stolze et al., 2020), as reflected in the breakthrough tailings observed at late times (i.e., after 1.5 PVs, Fig. 2d). The central port of the 2-D chamber displayed the most similar breakthrough curve compared to the column experiment. The integrated and intermediate location measurements yielded approximately the same H + concentration, whereas no relevant proton concentrations reached the bottom location of the 2-D chamber throughout the experiment. Table 1 Domain geometry and flow and transport parameters (h: porosity, v: seepage velocity) in the one-dimensional and two-dimensional setups.   Fig. 2b displays the pH breakthrough curves for the case of the NaBr injection. For both 1-D and 2-D domains, pH dropped sharply at approximately one pore volume, where the minimum was reached, before recovering towards higher pH values. In the 2-D domain, the highest H + concentration was measured at the central location of the chamber (i.e., 2.4 mM), whereas the proton concentrations decreased at the lower locations of the outlet boundary (i.e., 1.8 mM and 1.3 mM at the intermediate and bottom locations, respectively). In the 1-D system, the H + concentration (i.e., 2 mM) was considerably lower than the measured value at the central port of the 2-D system (i.e., 2 mM). Fig. 2c and f show that the pH and proton breakthrough curves for the NaBr plus HCl experiment were similar to the case of NaBr injection. At approximately one pore volume, a fast wave of protons reached the outlet of the domain during both one-dimensional and quasi twodimensional flow-through experiments. After reaching minimum pH values, the breakthrough curves displayed a slow recovery towards a steady-state plateau. For the 1-D domain, the proton concentration approximately reached the value of the injected solution (i.e., corresponding to pH $ 2.8), whereas in the 2-D system, the pH measured at the outlet boundary increased towards the bottom locations of the chamber (pH $ 3.37, 3.44, 3.63, and 3.73 for the central, integrated, intermediate and bottom ports, respectively, Fig. 2c). The difference between 1-D and 2-D experiments at the steady-state plateau, after the observed gradual recovery, is more visible by looking at the breakthrough curves of normalized H + concentration (Fig. 2f).

Models geometry
By comparing the proton breakthrough curves of the three experimental scenarios in the 2-D domain, different shapes of plume breakthrough are visible. In particular, the characteristic pattern displayed by the injection of a salt electrolyte (i.e., fast acidic wave at approximately one pore volume, Fig. 2e) dominates the overall breakthrough behavior also when the injection is performed under acidic conditions, resulting in considerably different breakthrough shapes between the HCl (Fig. 2d) and NaBr plus HCl (Fig. 2f) experiments. Fig. 3 displays the comparison between the measured and simulated proton breakthrough curves for the 1-D experiments performed using continuous injections of HCl (Fig. 3a), NaBr (Fig. 3b), and NaBr plus HCl (Fig. 3c).

Model-based interpretation of proton and major ions breakthrough curves
The simulated breakthrough curves are in good agreement with the measured data confirming that the definition of two distinct surfaces (i.e., quartz and metal oxides) allows reproducing the experimental observations (RMSE reported in Table EA4 of Electronic Annex). In the case of the experiment performed by injecting HCl (Fig. 3a), the model can capture the retarded proton breakthrough and the tailing behavior towards the injected proton concentration. In the case of the NaBr injection (Fig. 3b), the wave of protons emitted from the silica sand is wellcaptured by defining the formation of Na + outer-sphere at the quartz interface. Similarly, the experimental results are well captured by the reactive transport model also in the case of NaBr plus HCl (Fig. 3c). Fig. 4 displays the comparison between the measured and simulated proton breakthrough curves of the quasi 2-D flow-through experiments performed with the injection of the same electrolyte solutions: HCl (Fig. 4a-d), NaBr ( Fig. 4e-h), and NaBr plus HCl ( Fig. 4i-n). The outcomes of the simulations are plotted for each location at which the pH was continuously measured. From left to right, the four column panels refer to the integrated, central, intermediate, and bottom ports, respectively. The breakthrough curves of the acidic experiments, HCl (Fig. 4a-d) and NaBr plus HCl ( Fig. 4i-n), are plotted in terms of normalized effluent proton concentrations. The 2-D reactive transport model can reproduce well the H + breakthrough curves in the presence of strong ionic strength gradients in the transverse direction. A good match between the multicomponent ionic transport simulations and the experimental observations was obtained both for the integrated and for the port-resolved outlet observations (RMSE reported in Table EA4 of Electronic Annex). These results show the capability of the transport model to provide an adequate description of the physical and electrostatic masstransfer processes in the pore water, as well as the appropriateness of the surface complexation formulation that allowed accounting for the effect of local ionic strength on the surface protonation/deprotonation behavior. Fig. 5 displays the measured and simulated breakthrough curves of proton and major ions (i.e., Na + , Br À and Cl À ) transported during the NaBr and NaBr plus HCl flow-through experiments. Fig. 5a and d show the out-comes of the 1-D column experiments. The results of the quasi two-dimensional experiments are displayed for the locations at which the effluent samples were collected: ports 9-10 at the plume core ( Fig. 5b and e) and ports 4-5 at the plume fringe ( Fig. 5c and f). The curves are displayed in molar concentrations. Fig. 5a displays the curves of the NaBr injection in the 1-D column. Both sodium and bromide show a nearly conservative transport behavior with concentration values almost  two orders of magnitude higher than H + (i.e., 100 and 2 mM, respectively). Lower concentrations were also measured for the anion Cl À (i.e., $ 1.6 mM) in the case of the NaBr plus HCl 1-D experiment (Fig. 5d). Fig. 5b-c and e-f display the breakthrough curves of H + , Na + , Br À , and Cl À for the cases of NaBr and NaBr plus HCl injection in the 2-D experiments. The difference in concentration between the protons and the highly concentrated NaBr plume tends to diminish towards the fringe of the electrolyte plume (Fig. 5c and f). Furthermore, at the plume fringe in the two-dimensional domain (ports 4-5), Na + is retarded compared to the other ions (retardation factor of $ 1.13 for both experiments) and displays a more pronounced tailing behavior with lower concentrations compared to bromide. In particular, at early times (i.e., $ 1-1.5 PVs), due to the temporary increase of proton concentrations, Br À reaches concentration values higher than Na + , which conversely tends to decrease to counterbalance the local electroneutrality (25% relative difference compared to bromide for both NaBr and NaBr plus HCl experiments at 1.2 PVs). The retardation of Na + and the similar earlier breakthrough of H + and Br À (i.e., $ 0.8 PVs) indicate that the local charge balance at the invading front of the electrolyte plume is mainly ensured by the protons released from the surface and by bromide, which become electrostatically coupled. After one pore volume, H + concentrations decrease and this leads to a progressively more important coupling of Na + and Br À , as reflected in their similar concentration values at late times ($0.9 Â 10 À2 M for both experiments).

DISCUSSION
The outcomes of the one-dimensional and quasi twodimensional flow-through experiments show distinct behaviors of proton and major ions transport. A schematic conceptualization of the mechanisms (i.e., electrostatic interactions in the pore water and at the solid-solution interface) controlling reactive transport for the injection of the acid and salt plumes is presented in Fig. 6. In case of HCl injection, the main mechanisms are the surface complexation of protons onto the metal oxide coatings and their electrostatic coupling with chloride ions in the pore water (Fig. 6a). In the experiments with NaBr injection, Na + sorbs to the quartz surface and displaces the protons, which further interact with the oxide coatings. In the pore water, at the fringe of the plume, the displacement of bromide becomes progressively more coupled to the one of the protons as a consequence of Na + delay. On the contrary, within the core of the plume, the electroneutrality is dominated by the ions in excess concentration (i.e., Na + and Br À ). These effects are more apparent for the multidimensional setup as is further analyzed and discussed in the following.

Impact of dimensionality on proton breakthrough
The impact of the system dimensionality on the extent of surface and aqueous charged species interactions was investigated by quantifying and comparing the total amounts of protons eluted during the flow-through experiments in both one-dimensional and quasi two-dimensional domains. This was done by calculating the ''boundary-normalized integrated breakthrough curves" at the outlet of the domains (Muniruzzaman and Rolle, 2017). The integrated BTCs are generated by vertically integrating the mass flux at the outlet and normalizing by the integrated mass flux at the inlet. For a two-dimensional system, it reads as: where Cðx; y; tÞ is the concentration of the transported solute at a location (x,y) and specific time step t, qð0; yÞ is the specific discharge in the longitudinal direction, C in ð0; y; tÞ is the inflow concentration and y* is the height of the saturated porous medium (i.e., y* = 11.5 cm). Here, the solute concentration, Cðx; y; tÞ is calculated at the outlet of the domain (i.e., x = 30 cm), whereas C in is constant for the entire experimental time. Fig. 7 displays the comparison between the onedimensional and two-dimensional proton integrated BTCs for the six experiments with the three different injections: HCl (Fig. 7a), NaBr (Fig. 7b), and NaBr plus HCl (Fig. 7c). The integrated proton BTCs in the case of HCl injection are normalized to the injected proton mass flux, whereas for the two experiments performed with NaBr as background electrolyte, the integrated BTCs are normalized to the sum of the injected sodium and proton mass fluxes.
The proton integrated BTCs in the 1-D and 2-D domains are almost identical for the experiments performed with HCl (Fig. 7a). The protonation of the metal oxides surface, which is the primary reactive mechanism controlling the transport of the HCl plume, only results in negligible discrepancies between the two domains, suggesting that the lateral spreading of H + was less pronounced and restricted to the inner location of the chamber, exposing comparable amounts of porous medium for sorption reactions between the 1-D and 2-D systems. Conversely, the integrated BTCs of the experiments performed with NaBr ( Fig. 7b) show that the proton flux-weighted concentrations are more significant in the 2-D domain, indicating that the extents at which the pH front propagation is affected by the surface-solute interactions vary with the system dimensions. In particular, after 0.6 PVs from the injection, the amount of aqueous H + generated through the adsorption of Na + onto the quartz surface (surface complexation mechanism illustrated in Fig. 6b) is more than double compared to the one released in the 1-D system (Table 2) and the peak concentration of the proton breakthrough in the 2-D setup is almost three times than in the 1-D column (Fig. 7b). Sim-ilarly, the change in sodium outer-sphere complexes at the quartz surface referred to the initial conditions of the simulation displays a 35% increase in the 2-D domain ( Table 2). The full temporal evolution of the total moles of H + and sodium outer-sphere complexes generated within the 1-D and 2-D domains is shown in Fig. EA1. The curves further confirm the enhancement of the surface complexation reactions in the multidimensional domain. As the inlet mass fluxes were identical between the one-dimensional and two-dimensional setups, the enhancement of the surface complexation reactions is due to the greater amount of quartz surface area exposed to Na + in the presence of transverse dispersive fluxes. Reactive transport simulations indicate that the transverse dispersive fluxes alone contribute to a 61% increase of the moles of H + released from the quartz surface (Table EA5). Similarly, a 34% increase of the sodium outer-sphere complexes is a direct consequence of the transverse displacement of Na + outside the central location of the chamber (Table EA5). As a consequence, the volume of acidic solution (i.e., pH < 3) flushed through the porous medium during the entire duration of the NaBr injection was more than double in the 2-D setup compared to the analogous 1-D case (i.e., 84.4 mL and 37.5 mL, respectively). These findings also show that the quartz surface was fully saturated with Na + reaching its maximum capacity of proton release within the core of the injected  Table 2 Overview of the maximum proton concentrations, minimum pH values, and total amounts (expressed in mol units) of aqueous species and sodium outer-sphere complexes present within the domains at t = 0.6 PVs after injection for the NaBr and NaBr plus HCl flow-through simulations.

NaBr
NaBr plus HCl 1.30 Â 10 À3 1.30 Â 10 À3 1.30 Â 10 À3 1.30 Â 10 À3 DSiONa 3.14 Â 10 À5 4.23 Â 10 À5 2.72 Â 10 À5 2.77 Â 10 À5 plume in the 2-D chamber, and larger amounts of H + were generated by exposing more quartz surface area to the transported electrolyte. Similar to the case of NaBr, the proton flux-weighted concentrations for the experiments performed with NaBr plus HCl (Fig. 7c) are more prominent in the 2-D system. Interestingly, the percentage increase of moles of H + generated as a consequence of the transverse displacement of Na + is lower compared to the NaBr injection (i.e., 44% against 61%). Furthermore, the amount of sodium outersphere complexes is almost identical between 1-D and 2-D setups. These discrepancies with the NaBr case originate from the competition for sorption processes onto the quartz surface between H + and Na + . Indeed, during the NaBr plus HCl experiment, both proton and sodium were transported within the core of the electrolyte plume and consequently they competed for complexation reactions at the quartz surface. This resulted in lower amounts of sodium outer-sphere complexes and in an enhancement of quartz protonation.

Transverse mixing and spatial distribution of solute and surface species in the 2-D setup
In order to analyse the transverse displacement of the pH front in the 2-D flow-through experiments, Fig. 8 shows the vertical profiles at the outlet boundary and at specific times from early breakthrough ($6-7 hours) to steady state.
The proton breakthrough occurs earlier for the experiments performed with NaBr as background electrolyte where significant concentrations are already detected after 6 hours from the injection. After the breakthrough, the H + concentration increases monotonically, assuming a Gaussian distribution for the HCl experiment (Fig. 8a). On the contrary, in the presence of NaBr (Fig. 8b-c), the plumes spread more along the vertical direction, developing a non-Gaussian shape with a plateau of maximum concentration values within the central location of the chamber where the quartz surface was fully saturated with Na + . The concentration profiles display an interesting temporal behavior with an initial increase, until approximately 7 hours from injection, followed by a reversal of the trend and a decreasing behavior towards steady state. In the case of the NaBr experiment (Fig. 8b), the H + concentration at steady state is negligible compared to the NaBr plus HCl injection (Fig. 8c). In the latter experiment the steadystate proton profile has a Gaussian distribution similar to the case of only HCl injection. However, the H + steadystate profile in the presence of the highly concentrated NaBr is more spread and displays lower peak concentration. This effect is driven by the Coulombic interactions in the pore water since the presence of NaBr with higher excess (concentration values approximately 100 higher than HCl) ensures the electroneutrality of the pore water mainly by sodium and bromide ions. This allows H + to migrate with a mobility close to its self-diffusive/dispersive properties. Conversely, in the case of pure HCl injection (Fig. 8a) the H + displacement is dampened by the electrostatic coupling with Cl À , in order to maintain the electroneutrality of the pore water solution, resulting in a less spread transverse profile.
In Fig. 9, we illustrate the 2-D maps of the aqueous species concentration and the change in quartz surface composition with respect to the initial conditions for the experiment performed with NaBr injection. The 2-D maps of the NaBr plus HCl experiment and a hypothetical conservative case where the surface reactions are deactivated during the experiment performed with NaBr plus HCl are shown in Fig. EA2 and EA3, respectively.
The H + plume travels at the front of the injected electrolyte (Fig. 9a) and is more spread along the transverse direction compared to Na + and Br À plumes (Fig. 9b-c). The wave of protons observed in Fig. 8b (6 h vertical profile) is highlighted by the warmest colors (darkest red area) in Fig. 9a (H + $ 2.37 mM corresponding to pH $ 2.62). The acidic front is generated by the quartz surface, which undergoes deprotonation (Fig. 9d) to allow the formation of Na + outer-sphere complexes (Fig. 9e). The surplus of positive charge is balanced by bromide both in the solid (Fig. 9f) and aqueous (Fig. 9g) phases. The difference in aqueous concentration between Na + and Br À (Fig. 9g) confirms the retardation of sodium and the movement of bromide at the front of the main electrolyte plume. The increase of sodium outer-sphere complexes concentration is mirrored by the decrease of surface charge at the 0plane of the quartz surface (Fig. 9h). Furthermore, Fig. 9l and m show that the excess of H + and Br À relative to Na + becomes significant at the fringe of the plume due to the retardation of Na + and the need for H + and Br À to ensure the electroneutrality. Finally, Fig. 9n shows the zones where both H + and Br À are in excess relative to Na + . Different electromigration fluxes of H + and Br À are expected in these zones as the ions are fully electrostatically coupled compared to the core of the plume, where Na + and Br À dominate electroneutrality (as illustrated in Fig. 6b, within the orange fringe area, representing the zone of excess of H + and Br À relative to Na + , proton and bromide are electrostatically coupled, whereas within the core of the NaBr plume, Na + and Br À dominate the electroneutrality).
As expected from the experimental breakthrough curves a similar behavior is displayed by the case of NaBr plus HCl injection (Fig. EA2). In this case, the adsorption of Na + is partially inhibited by the concomitant presence of H + in the core of the electrolyte plume. Proton and sodium outer-sphere complexes concentration and surface charge at the 0-plane are weakened, as also reported in Table 2.
The 2-D multicomponent ionic transport model allows us to understand and visualize the Coulombic interactions between charged aqueous species by mapping the spatial distributions of the different components of the Nernst-Planck fluxes. According to Eqs. (1), (8), and (10), the total transverse dispersive flux can be written as: , representing the contribution due to the electrostatic coupling between charged aqueous species, and an activity gradient flux J Act i À Á . In Fig. 10 we illustrate the components of the proton total dispersive flux at t = 0.6 PVs for the three performed experiments in the 2-D flow-through chamber: HCl (Fig. 10a-c), NaBr (Fig. 10d-f), and NaBr plus HCl (Fig. 10g-i). Furthermore, a fourth case (Fig. 10l-n) is shown where the surface reactions are deactivated during the injection of NaBr plus HCl.
The 2-D maps for the case of the HCl injection ( Fig. 10a-d) highlight a negative proton migration flux, leading to a reduction of the total transverse spreading of the acidic plume. During the transport in pure water, H + is attracted by Cl À to ensure the electroneutrality of the solution. Since the self-aqueous diffusion coefficient of chloride is considerably smaller than the one of H + (i.e., 1.81 Â 10 À9 and 8.64 Â 10 À9 m 2 /s, respectively), the transverse displacement Fig. 9. Simulated two-dimensional maps at t = 0.6 PVs for the case of the NaBr injection: (a-c) aqueous H + , Na + , Br À concentrations, (d-f) change in silica surface composition, (g) difference in aqueous concentration between Na + and Br À , (h) silica surface charge at the 0-plane and (i) 1-plane, (l) aqueous ratios between H + and Na + , (m) Br À and Na + , and (n) their product. of the proton plume is limited (Fig. 8a). Conversely, in the case of NaBr (Fig. 10e-h) and NaBr plus HCl injections ( Fig. 10i-n), the negative migration flux becomes positive within the core of the plume as Na + is present at a high concentration. Indeed, in the presence of a mixture of two electrolytes in different concentrations, the background solution (i.e., in this case, the core of the plume) is charged balanced by the ions present in excess concentration (i.e., Na + and Br À ), whereas those in smaller amounts (i.e., H + and Cl À ) can travel according to their self-diffusive/ dispersive properties or can even accelerate compared to their liberated condition . This explains the wider transverse displacement of H + observed in Fig. 8b-c compared to the case of the HCl injection (Fig. 8a). At the fringe of the invading front of the electrolyte plume, where H + and Br À ensure the electroneutrality of the system as they are present in excess concentration compared to Na + (Fig. 9n), the electromigration flux ( Fig. 10g and m) contributes negatively to the total dispersive flux due to the lower self-aqueous diffusion coefficient of Br À compared to H + (i.e., 1.86 Â 10 À9 m 2 /s).
It is interesting to note that the multidimensionality of the domain contributes to the formation of distinct zones of Coulombic interactions, corresponding to different H + sources: for example, considering the NaBr plus HCl experiment, the core of the plume hosts protons injected from the inlet boundary, whereas the protons released from the silica porous medium are located at the fringe of the plume. When the surface complexation reactions are deactivated ( Fig. 10o-r), the zone with the negative migration flux is absent, and the proton front accelerates due to the positive contribution of J Mig i . The activity gradient fluxes have a negative contribution in all cases because the activity coefficient gradients typically occur in the opposite direction of the concentration gradient (Eq. (1)) ( Fig. 10d, h, n, and r). For H + , these flux components are two orders of magnitude smaller compared to the purely diffusive/dispersive fluxes in all experiments. However, it should be noted that, for the species with higher concentrations (Na + and Br À ), the magnitude of the computed activity gradient fluxes is comparable to the one of the dispersive and migration components (results not shown).
The above discussion focused on proton transport and the propagation of pH fronts. However, we believe that the role of dimensionality should also be addressed for different geochemical systems, including reactive transport of charged organic and inorganic contaminants (Prigiobbe and Bryant, 2014;Lyu et al., 2018;Ye and Prigiobbe, 2018;Brusseau et al., 2019), transport in presence of reactive minerals and mineral assemblages (Li et al., 2014;Stolze et al., 2019a,b;Jung and Navarre-Sitchler, 2018), and electrokinetic transport processes (Alizadeh and Wang, 2019;Sprocati and Rolle, 2020;Sprocati et al., Fig. 10. Maps of different transverse flux components of H + at t = 0.6 PVs: (a-d) HCl, (e-h) NaBr, (i-n) NaBr plus HCl, and (o-r) NaBr plus HCl with the surface complexation reactions deactivated. The column panels display the total, purely dispersive, electromigration, and activity gradient fluxes from left to right, respectively. Note that the direction from the core to the fringe of the plume is considered positive for the calculated fluxes. 2021). Furthermore, we expect that the effects induced by transverse mixing would be even more pronounced within fully three-dimensional systems (Ye et al., 2015a). Physical heterogeneities and anisotropies, as well as complex flow topologies in 3-D porous media, can lead to different mixing mechanisms (Chiogna et al., 2014;Ye et al., 2015b,c; with significant impact on the fate and transport of charged solutes and on their interactions with solid surfaces.

CONCLUSIONS
In this study, we evaluated dimensionality effects on multicomponent ionic transport and surface complexation by comparing one-dimensional and quasi two-dimensional flow-through experiments performed in the same silica sand porous medium, under the same advection-dominated flow regimes, hydrochemistry, and injection conditions. Three different electrolyte solutions with sodium bromide and hydrochloric acid were continuously injected through the domain, and measurements of pH and major ions breakthrough curves were performed at the outlet of the setups. The laboratory outcomes were quantitatively interpreted by performing one-and two-dimensional simulations using a multicomponent reactive transport model which included the electrochemical cross-coupling of charged species in solution and a surface complexation formulation of the silica sand surface.
Our experimental and modeling investigation highlights distinct features in the 1-D and 2-D setups. In particular, we demonstrate that the electrostatic coupling between transverse dispersive fluxes of charged aqueous species and their interaction with reactive mineral surfaces profoundly impacts the evolution of acidic plumes in multidimensional flow-through systems. Specifically, the transverse displacement of Na + exposes more quartz surface area to complexation reactions resulting in higher amounts of H + released from the silica porous medium upon adsorption of sodium (more than double amounts of H + transported in the 2-D domain compared to the equivalent 1-D case). In fact, the surface saturation in Na + occurring within one-dimensional systems leads to a considerably lower volume of acidic solution (e.g., pH < 3) that can be transported in analogous multidimensional domains (i.e., relative difference up to 125%). Furthermore, the acidic plumes generated through the release of protons from the silica porous medium display wide transverse concentration profiles with non-Gaussian shapes and a non-monotonic temporal development, with a rapid increase of acidity at early breakthrough times followed by a slow decrease of proton concentration towards steady-state conditions. The outcomes of the numerical simulations also highlight that the generation of strong ionic strength gradients in the multidimensional domain contributes to a spatially heterogeneous zonation of Coulombic interactions. The charge balance of the domain is fulfilled differently in the fringe and core of the electrolyte plume depending on how adsorption reactions affect relative concentrations of aqueous charged species.
This investigation quantitatively showed the importance of dimensionality for subsurface reactive transport. Onedimensional column experiments may mask the role of microscopic Coulombic interactions in both aqueous and solid phases, which instead becomes important in multidimensional domains. The findings of this study are relevant for a wide variety of geochemical systems. For instance, mixing of groundwater with different chemical composition in subsurface critical zones and at the fresh/salt water interface typically entail pH fluctuations that play a key role in controlling the release and mobility of trace metals and metalloids such as nickel, cadmium, zinc and arsenic (e.g., Zachara et al., 1991;Kent et al., 2007;Stachowicz et al., 2008). Moreover, redox reactions and the propagation of redox fronts in different subsurface environments (e.g., Engesgaard and Kipp 1992;Haberer et al., 2015;Battistel et al., 2019), as well as the dissolution rates of many reactive minerals (e.g., Chen and Brantley, 1997;Sjö berg and Rickard, 1984) are dynamically coupled to the groundwater pH and buffering capacity. The outcomes of this study highlight the need to develop approaches capable of integrating the complexity of physically based, multidimensional transport processes with complex geochemical reactions for an improved description of subsurface geochemical systems.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.