Surface complexation reactions in sandy porous media effects of incomplete mixing and mass-transfer limitations in flow-through systems

Although mixing and surface complexation reactions are key processes for solute transport in porous media, their coupling has not been extensively investigated. In this work, we study the impact of mass-transfer limitations on heterogeneous reactions taking place at the solid-solution interface of a natural sandy porous medium under advection-dominated flow-through conditions. A comprehensive set of 36 column experiments with different grain sizes (0.64, 1.3 and 2.3 mm), seepage velocities (1, 30 and 90 m/day), and hydrochemical conditions were performed. The injection of NaBr solutions of different concentrations (1 – 100 mM) led to the release of protons via deprotonation reactions of the quartz surface. pH and solute concentration breakthrough curves were measured at the outlet of the columns and the propagation of pH fronts in the column setups was tracked inside the porous medium with non-invasive optode sensors. The experimental results show that the deprotonation of the reactive surfaces, resulting from their interactions with the injected ionic species, strongly depends on the hydrodynamic conditions and differs among the tested porous media despite their apparent similar surface properties. Reactive transport modeling was used to quantitatively interpret the experimental results and to analyze the effects of mass-transfer limited physical processes on surface complexation reactions, propagation of pH fronts, transport of major ions and spatio-temporal evolution of surface composition. A dual domain mass transfer formulation (DDMT) combined with a surface complexation model (SCM) allowed capturing the effects of incomplete mixing on the surface reactions and to reproduce the experimental observations collected in the experiments with high flow velocities. The SCM was parameterized with a single set of surface complexation parameters, accounting for the similar surface properties of the porous media, and was capable of describing the surface complexation mechanisms and their impact on the hydrochemistry over the large range of tested ionic strengths.

Numerous experimental and modeling studies have shown that the degree of mixing between reactants is primarily controlled by molecular diffusion and local scale dispersion (Kitanidis, 1994;Chiogna et al., 2011;Cirpka et al., 2012).Incomplete mixing caused by the interplay between advective and diffusive mass transfer, leads to the formation of concentration gradients within single pore channels (Willingham et al., 2008;Rolle et al., 2012).Furthermore, the heterogeneity of the pore space leads to the stretching and deformation of solute plumes and to complex mixing patterns (Willingham et al., 2008;Rolle et al., 2009;Le Borgne et al., 2010;Le Borgne et al., 2011;Liu and Kitanidis, 2012;Hochstetler et al., 2013).
The impact of incomplete mixing propagates across scales (Raje and Kapoor, 2000;Gramling et al., 2002;Bolster et al., 2011;Chiogna et al., 2011;Rolle et al., 2013b;Sole-Mari et al., 2021) and can significantly affect the macroscopic transport behavior of aqueous solutes as well as the extent of reactions taking place at the solid-solution interface (Li et al., 2008;Molins et al., 2012;Molins et al., 2014).An accurate consideration of mixing processes and mass transfer limitations from pore to field scales is therefore of utmost importance to describe the transport of aqueous species in groundwater.Considering the difficulty in obtaining a sufficient resolution of the microscopic properties of natural porous media (e.g., pore geometry) and the high computational expense of pore-scale simulations, it is often necessary to describe porescale mechanisms with macroscopic representations.In particular, macroscopic models, based on continuum formulations, are commonly used to simulate conservative and reactive transport from laboratory to field scale applications (Steefel et al., 2015;Steefel, 2019).Different approaches such as dual domain mass transfer model (DDMT) (Gottschlich, 1963;Coats and Smith, 1964;van Genuchten and Wierenga, 1976), continuous time random walk (Berkowitz et al., 2006;Heidari and Li, 2014;Sole-Mari et al., 2020), and/or hybrid pore scalecontinuum simulations (Meile and Tuncay, 2006;Battiato and Tartakovsky, 2011) have been used to describe transport in presence of mass transfer limitations both in homogeneous and complex heterogeneous porous media.Among these approaches, the DDMT is probably the most popular and the simplest formulation to capture mass-transfer limitations and has been applied in a number of studies on transport in aggregated porous media (van Genuchten and Wierenga, 1977;Valocchi, 1985;Ball and Roberts, 1991;Gerke and van Genuchten, 1993;Fesch et al., 1998;Qafoku et al., 2005), in heterogeneous subsurface systems (Gerke and van Genuchten, 1993;Luo et al., 2005;Liu et al., 2007;Ronayne and Gorelick, 2010;Blackmore et al., 2018;Binder et al., 2021), as well as in homogeneous porous media where high flow velocity results in preferential flow paths and incomplete mixing (Luo et al., 2007;Liu and Kitanidis, 2012;Rolle and Kitanidis, 2014;Muniruzzaman and Rolle, 2017).
Surface complexation reactions at the solid/solution interface play a key role in the transport of aqueous charged species in porous media (Kent et al., 2007;Sharma et al., 2011;Johannesson and Neumann, 2013;Prigiobbe and Bryant, 2014;Rathi et al., 2017;Tournassat and Steefel, 2019) and are commonly included in continuum scale reactive transport models (Mayer, Benner, and Blowes, 2006;Jessen et al., 2012;Ye and Prigiobbe, 2018;Sprocati et al., 2019;Stolze et al., 2019a;Fakhreddine et al., 2020;Wallis et al., 2020).However, the coupling of surface complexation and mixing processes has received limited attention.Only a few contributions have assessed the impact of mass transfer limitations on adsorption processes (James and Rubin, 1979;Valocchi, 1985;van Genuchten and Wierenga, 1976;Darland and Inskeep, 1997).These investigations have focused on the impact of flow velocity and mass transfer processes in flow through systems but have been limited to simple sorption models (isotherms, ion exchange).A rigorous description of the electrostatic properties at the solid-solution interface is, instead, required to accurately capture the sorption behavior and the reactive transport of aqueous charged species, including important inorganic contaminants such as heavy metals and metalloids (Stachowicz and Hiemstra, 2008;Stolze et al., 2019b).Therefore, the coupling of comprehensive surface complexation models and mass-transfer limited transport should be explored and tested.
Silica sand is of particular interest for a variety of groundwater applications considering its ubiquity in the subsurface and its high permeability, which can lead to incomplete mixing in presence of high groundwater flow velocities (Atchley et al., 2014;Jung and Navarre-Sitchler, 2018a;Muniruzzaman and Rolle, 2017;Wen and Li, 2018).Sand also exhibits a strong surface reactivity with protons and aqueous charged species, which impacts pore water hydrochemistry (Kent and Fox, 2004;McNeece et al., 2018;Garcia et al., 2019;Stolze et al., 2020) and, potentially, other solute/mineral interactions.
In this study, we investigate the effects of mass-transfer limitations induced by incomplete mixing on surface complexation under flowthrough conditions by analyzing the propagation of acidic plumes and the transport of major ions in sandy porous media.We performed a series of 36 column experiments by injecting salt solutions with different electrolyte (NaBr) concentrations to trigger the deprotonation of the quartz sand surface.The propagation of the aqueous species was monitored by measuring the concentration of the aqueous species at the outlet of the columns and by performing non-invasive in-situ measurements of the pore water pH.The experiments were repeated using different grain sizes and flow velocities in order to test the influence of the porous medium granulometry and of the flow regime on mixing and mass-transfer limitations.A reactive transport model, based on dual domain transport and on a description of the chemical and electrostatic surface complexation mechanisms in the sandy porous media, was implemented to quantitatively interpret the experimental results.

Laboratory experiments
Flow-through experiments were performed in cylindrical glass columns (length: 14.4 cm, inner diameter: 1.75 cm) with a natural quartz sand originating from Dorsten, Germany (Aquagrandry, Euroquartz).The columns were wet packed with the porous medium to minimize entrapment of air in the pores (Haberer et al., 2012) and connected to high precision peristaltic pumps (IPC-N24, ColeParmer, United States) using two ports at each end of the columns.Continuous vertical injection from the bottom was applied to establish steady-state laminar flow conditions.The column experiments were initiated by flushing the domain with Milli-Q water until steady pH conditions were reached (i.e., pH ~ 5.9).This allowed equilibrating the quartz sand surface with the initial pore water solution.Subsequently, a NaBr solution was injected as a step change from the bottom of the columns to alter the surface charge of the quartz via electrostatic interactions and, thus, to trigger a change in the quartz sand surface protonation.As illustrated in Fig. 1, the experiments were performed with the same quartz sand but with (i) different average grain sizes (0.64, 1.3 and 2.3 mm), (ii) different concentrations of NaBr electrolyte (0.1, 1, 10 and 100 mM) and, (iii) different seepage velocities (1, 30 and 90 m/day; see Table S1) in a comprehensive set of 36 flow-through columns.The selected high seepage velocities are representative of fast natural flow in permeable aquifers, as well as of engineering applications such as groundwater pumping, forced-gradient tracer tests, artificial recharge, delivery of reactants and/or mixing enhancement for remediation (e.g., Manning and Solomon, 2005;Magal et al., 2010;Kahler and Kabala, 2016).
For all column experiments, pH was measured at the outlet of the columns by connecting one of their outlet ports to 10 mL flow-through vials equipped with pH probes (Hach Intellical PHC281).The pH was monitored with a time interval of 5 min for the experiments conducted with a seepage velocity of 1 m/day and with a time interval of 10 s for the experiments conducted with seepage velocities of 30 and 90 m/day.Subsequently, the effluent pH breakthrough curves were calculated from the pH measurements using the continuous stirred tank reactor (CSTR) to account for dilution effects in the flow-through vial (see Supplementary Material).Furthermore, breakthrough curves of ions were measured for the column experiments operated at 1 m/day by collecting samples from the remaining outlet port with a time interval of 1 h.Sodium concentrations were measured from acidified samples (2% HNO 3 ) analyzed with inductively-coupled-plasma optical emission spectroscopy (ICP-OES, PerkinElmer Avio 200, in axial and/or radial mode).Bromide concentrations were determined from unacidified samples using ion chromatography (IC, Dionex ICS-5000, ThermoScientific).
During the experiments conducted with a seepage velocity of 1 m/ day, non-invasive high spatial resolution (5 mm spacing) monitoring of pore water pH was performed along the axes of the columns using prototype optical pH luminescence-based organic polymer sensors (SF-LG1, PreSens Gmbh, Germany).Although these pH sensors are primarily designed for physiological applications, this technique is here applied to geologic porous media in order to qualitatively track the pH front propagation along the column setup.Similar to the application of noninvasive optode sensors to track the propagation of oxygen fronts (Haberer et al., 2015;Battistel et al., 2019), pH sensor strips were glued onto the inner walls of the columns before packing the columns with sand.pH spatial profiles were acquired at selected time intervals by using an optic fiber on the outer wall of the columns in the experiments performed at 1 m/day with the highest ionic strength (100 mM NaBr), due to the sensitivity of the sensor that is best suited for a minimum ionic strength of 50 mM.The parameters used for converting the optical measurements into pH values were directly calibrated by comparison between the optical sensors response and the measurements with an electrode in a flow-through vial at the column outlet.This was achieved by continuously monitoring the fluorescence lifetime of an optode sensor spot placed in the flow-through vial used for measuring the pH breakthrough curve and referencing this signal to the pH determined with the electrode probe.The optical sensors were also used in the column experiments performed with a seepage velocity of 30 m/day (with 100 mM NaBr), in which pH breakthrough curves were recorded at 3 spatial locations (i.e., 3.6, 7.2, 10.5 cm from the inlet of the columns).
All experiments were conducted at 20 • C and the aqueous solutions were prepared using Milli-Q water.Before each column experiment, the sand was first washed with HCl pH 4 and NaOH pH 10 solutions and, successively, with a pH 2 solution to remove the remaining impurities, and thoroughly rinsed with Milli-Q water before packing the columns (McNeece and Hesse, 2016).The applied flow rates were measured at the end of all column experiments.Porosity was determined for the different quartz sand grain sizes by fitting the breakthrough curves of the conservative anion bromide measured in column experiments performed at a seepage velocity of 1 m/day (Table S2).The quartz sand surface areas were determined from samples collected before and after the column experiments with single and multi-point BET sorption with N 2 at − 195 • C.

Modeling approach
Numerical modeling of the column experiments was performed at the continuum scale, considering the physical processes controlling solute transport and the surface complexation reactions at the surface of the sand grains.

Surface complexation reactions
The interactions between the aqueous charged species and the sand surface were simulated with a surface complexation model.The model explicitly considered both the quartz surface and a generic metal oxide surface (i.e., accounting for the low fraction of Fe <0.02 wt% and Al <0.007 wt% coatings) following the approach implemented in our previous study (Stolze et al., 2020).
Specifically, the quartz sand and metal oxides were described with a Basic Stern model (Stern, 1924;Grahame, 1947) and a triple layer model (Yates et al., 1974;Davis, James, and Leckie, 1978), respectively.Both surfaces were represented with a single type of surface site having an amphoteric behavior (i.e., including protonation and deprotonation reactions) and capable to react with major ions (Na + and Br − ) through the formation of outersphere complexes.The latter are defined in the model β-plane and account for ions that electrostatically interact with the surface to compensate the charge deficit or surplus of the model 0-plane which is determined by the surface protonation state.Fig. 1 presents a schematic illustration of the electrostatic description of triple layer models.The reactive surface areas of the generic metal oxide for the three different sand grain sizes considered in the present study were calibrated and the reactive surface area of the quartz was defined as A quartz = A BET − A oxide .

Reactive transport equations
Reactive transport modeling was performed accounting for masstransfer limitations of aqueous solutes in the high flow velocity experiments (30 and 90 m/day) with a dual-domain approach (Gottschlich, 1963;Coats and Smith, 1964;van Genuchten and Wierenga, 1976).This approach consists in dividing the porous medium into a mobile domain in which advective and dispersive transport occurs and an immobile domain in which the fluid is considered stagnant and well mixed.A schematic illustration of the DDMT approach is provided in the bottom right plot of Fig. 1.The transport of solutes in the saturated porous media was described with the governing mass conservation equation: where C m,i , C im,i [mol/L] S m,i , S im,i [mol/kg], are the aqueous and adsorbed concentrations of a species i with subscripts m and im denoting the mobile and immobile domains, respectively, θ m and θ im [− ] are the porosity of the mobile and immobile domains, q [m/s] is the specific discharge, ρ s [kg/L] is the bulk dry density of the solid phase, θ [− ] is the total porosity, D L,i [m 2 /s] is the longitudinal dispersion coefficient of specie i that includes the contribution of pore diffusion and mechanical dispersion, the latter described as linearly dependent on the grain size and on the seepage velocity(Guedes de Carvalho and Delgado, 2005;Kurotori et al., 2019): where is the average grain size, and v [m/s] the seepage velocity.
The diffusive mass exchange between the mobile and the adjacent immobile cells was expressed as (Parkhurst and Appelo, 1999): where ξ [s − 1 ] is the first-order mass transfer coefficient.For the experiments performed at high flow velocities (i.e., 30 m/day and 90 m/day), the mass transfer coefficient ξ was fitted to reproduce the measured data.The total amount of reactive surface area was divided among the mobile and immobile domains using a calibrated scaling coefficient φ im where φ im = Aim Atot .The dual domain approach was not applied for the experiments performed at 1 m/day where the data did not indicate the need of a multiple-continua formulation and the simpler advectiondispersion equation coupled to the surface complexation model was sufficient to capture the experimental observations.For all simulations, constant flux boundary conditions were applied at the inlet and outlet of the domain to reproduce the flow field applied in the column experiments.

Numerical implementation
Simulations of the column 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;Muniruzzaman and Rolle, 2016).This approach allowed combining the geochemical, multicomponent transport and dualporosity model capabilities of PHREEQC with the computational methods of Matlab®.In particular, the latter was used to parallelize the simulations of multiple experiments, to conduct automated calibration of model parameters and to perform a comprehensive analysis of the data and model outcomes (e.g., Stolze et al., 2019b;Murray et al., 2019).The 1D domain representing the column setup was discretized into 50 grid cells (Δx = 2.88 mm).The aqueous speciation and reaction calculations were performed based on the thermodynamic database WATEQ4F, with the Davies, Debye-Hückel and extended Debye-Hückel equations used to calculate the ion activity coefficients of aqueous species.The Basic Stern surface complexation model was implemented in PHREEQC by setting the second capacitance C 2 to a large number (i. e., 1000 F/m 2 , McNeece and Hesse, 2016).

Calibration strategy
Calibration of the model parameters was performed in two successive steps.First, the surface complexation parameters were calibrated using the laboratory dataset (i.e., H + , Na + and Br − breakthrough curves) collected in the 12 column experiments performed at the slow flow velocity of 1 m/day, including all tested ionic strengths of the injected solutions and all three grain sizes.Subsequently, the DDMT parameters were calibrated based on the eluted H + breakthrough curves measured in the 24 experiments performed at higher flow velocities (i.e., 12 columns at 30 m/day and 12 columns at 90 m/day).In this step, the surface complexation parameters were fixed to the values determined from the slow flow velocity experiments.A visual illustration of the calibration approach is provided in the Supplementary Material (Fig. S3).
Modeling of the flow-through experiments and calibration of the model parameters for a given subset of experiments (i.e., low and fast flow experiments, respectively) were performed simultaneously by running in parallel the reactive transport simulations.Regarding the DDMT parameters, two mass transfer coefficients ξ for the respective seepage velocities used in the experiments (i.e., 30 and 90 m/day) were optimized.The scaling coefficient φ im was calibrated individually among 6 subsets, each composed of 4 experiments performed with equal flow velocities and grain sizes.The immobile porosity θ im was set to 0.02 based on a previous laboratory study in homogeneous porous media with similar granulometry (Muniruzzaman and Rolle, 2017).The global search particle swarm algorithm (PSO) method implemented in Mat-lab® was used to minimize the objective function.

Experimental results
Fig. 2 displays the proton breakthrough curves measured at the outlet of the columns for all tested NaBr concentrations, grain diameters and seepage velocities.The experimental results show the strong deprotonation of the quartz sand surface and the release of protons in the pore water induced upon the contact of the sand with the injected salt solution.
The injection of NaBr solutions led to a steep increase of the proton concentrations at the outlet of the column after 1 pore volume followed by a more progressive decrease in proton concentrations towards the pH of the inlet solution (~ pH 5.5) in all experiments.Remarkable differences in the elution of protons can be observed between the different tested conditions.In particular, the flow velocity and the grain sizes significantly impact the surface reactions under flow-through conditions.Regarding the hydrochemical conditions, the maximum concentration and amount of protons generated in the column experiments directly depended on the ionic strength of the inlet solutions, with increasing concentrations of eluted protons corresponding to higher concentrations of NaBr.For instance, in the experiments performed with a grain size diameter of 1.3 mm and a seepage velocity of 1 m/day (full line in Fig. 2e-h), the maximum proton concentration was ~20 times higher when injecting a 100 mM NaBr solution (C H+,max = 1.78 mM) compared to a 0.1 mM NaBr solution (C H+,max = 0.09 mM), respectively.This is also observed for the total amount of protons eluted (e.g., n H+,tot = 3.33 × 10 − 2 mmol and n H+,tot = 5.63 × 10 − 3 mmol for these two particular column experiments).
Comparison between proton breakthrough curves measured during L. Stolze and M. Rolle experiments performed under similar transport and ionic strength conditions also showed that the dynamics and extent of the protons released from the sand surfaces were influenced by the sand grain size.This impact was limited under low ionic strength conditions, as shown by Fig. 2a, e and i with similar maximum concentration of eluted protons (e.g., concentrations in the range [0.089-0.093]mM for the three grain sizes under seepage velocity of 1 m/day and 0.01 mM NaBr), but became more pronounced with increasing ionic strengths.In particular, the quartz sand with the smallest grain size led to the highest proton concentrations.For instance, with a seepage velocity of 1 m/day and a NaBr concentration of 100 mM NaBr, C H+,max was equal to 4.01 mM, 1.78 mM and 2.18 mM for 0.64 mm, 1.3 mm and 2.3 mm grain diameters, respectively.These differences are particularly interesting as they do not directly reflect the differences in BET surface area (i.e., A 2.3mm = 0.32 m 2 /g > A 0.64mm = 0.23 m 2 /g > A 1.3mm = 0.15 m 2 /g) and as all column experiments were performed with the same quartz sand.In contrast, the total amount of protons eluted from the column per unit of exposed surface area n H+ , tot was similar at the low seepage velocity of 1 m/day (i.e., 3.12 × 10 − 3 , 3.83 × 10 − 3 , 3.08 × 10 − 3 mmol/m 2 A BET for the experiments performed with 100 mM NaBr with the 0.64 mm, 1.3 mm and 2.3 mm grain diameters, respectively) demonstrating similar surface reactivity (i.e., equal total pool of protons per unit exposed surface area) for the different grain sizes.These contrasting results demonstrate that specific properties of sand grain sizes (e.g., different amount of coatings retarding the proton front and/or different hydrodynamic properties of the considered porous media) directly impacted the dynamics of the reaction between the ionic species and the reactive surface and therefore the shape of the proton breakthrough curves.
Considering the transport conditions, the results show that both the concentration and the total amount of eluted protons were strongly impacted by the applied flow rates.For all grain sizes and under all ionic strength conditions, higher seepage velocities lowered the release of protons (i.e., full, dashed and dash-dotted lines corresponding to velocities of 1, 30 and 90 m/day in Fig. 2).Furthermore, the protons breakthrough curves measured in experiments performed with seepage velocities of 30 and 90 m/day are nearly identical for average grain size diameters of 1.3 mm and 2.3 mm (Fig. 2e-l).In contrast, for grain diameter of 0.64 mm, the proton breakthrough curves are markedly different for the 3 tested seepage velocities (Fig. 2a-d).These results indicate that the grain size exerts a direct and important control on the degree of incomplete mixing of the injected aqueous solutes and their interaction with the reactive quartz surface under higher flow velocity and, consequently, on the release of H + from the quartz surface.
In addition to the pH breakthrough curves, we measured spatial pH profiles at different times using prototypes optical sensor strips.Fig. 3 shows the pH spatial profiles for the different grain sizes during the salt injections performed at high ionic strength (100 mM NaBr) and under a flow velocity of 1 m/day.For all experiments, the propagation of the proton fronts can be clearly identified from the profiles measured at 1 and 2 h (i.e., before flushing 1 PV).During the salt injections, the transport of H + is characterized by a sharp front followed by a gentle slope of the profile.After 3 h, the measurements show the sluggish recovery of pH within the column pore water in agreement with the pH breakthroughs observed at the outlet vials (solid lines in Fig. 2d, h and l).
Non-invasive in-situ pH measurements were also performed in the column experiments carried out with a flow velocity of 30 m/day and inlet solution with an ionic strength of 100 mM NaBr.To track the propagation of the pH fronts under such high flow velocity, optical sensors were placed at three fixed locations along the column (i.e., 3.6, 7.2 and 10.5 cm from the inlet).Fig. 4 displays the pH breakthrough curves measured at these locations within the porous media, using the non-invasive optode technique for the three different tested grain sizes.
The pH breakthrough curve shows a sharp decrease in proton concentrations when the acidic plume reaches the location of the sensor in agreement with the H + breakthrough curve measured at the outlet of the columns with the pH sensor (Fig. 2).However, these measurements provide additional information about the dynamics of the H + plume formation.In fact, Fig. 4 clearly shows a progressive increase in proton concentrations at the front of the plume as the latter progresses within the porous media towards the outlet of the columns.This pattern results from the increasing amount of quartz sand surface that reacts as the NaBr plume propagates along the column leading to an accumulation of protons at the front of the plume where the quartz surface deprotonation reaction is highest.
The impacts of the ionic strength, transport condition and grain size diameter on the release of protons are summarized in Fig. 5.In order to remove the effects of surface area differences among the tested quartz grains, the maximum concentration of protons measured at the outlet of the column (left panels in Fig. 5) and the total amount of protons eluted after 10 PV (right panels in Fig. 5) were normalized by the measured BET surface area exposed in the respective experiments.
The strong impact of flow velocity on the release of protons is clearly visible in Fig. 5.This is both expressed by the substantial decrease in the maximum concentration of protons measured at the columns outlets (Fig. 5a-d) and the total amount of protons released from the porous media for a given number of pore volume flushed (Fig. 5e-h) when increasing the flow velocity.
Fig. 5 highlights significant differences in the interactions between the aqueous species and the reactive surface for the different sand grains used in the experiments.In particular, the maximum concentration of H + per unit surface area is significantly higher for the smaller grain sizes (Fig. 5d) whereas the larger grain size consistently shows lower proton concentration maxima (green bars in Fig. 5a-d).Furthermore, the relative decrease in maximum H + concentration when increasing the flow velocity increases with the grain sizes (e.g., − 58%, − 71% and − 82% for grain sizes 0.64 mm, 1.3 mm and 2.3 mm and an inlet solution with 100 mM NaBr when increasing the flow velocity from 1 to 30 m/day).
Differently from the maximum concentration of protons measured at the column outlets, the total amount of eluted protons (Fig. 5e-h) is very similar for the three grain sizes in the low seepage velocity experiments and, in particular, for the experiments performed with higher concentration of NaBr (i.e., 10 and 100 mM in Fig. 5g and h, respectively).Such similarity reflects the common surface properties (i.e., equal total pool of protons per unit exposed surface area) of the grain sizes, although the intermediate grain size systematically shows slightly higher total amount of eluted protons (red bars in Fig. 5e-h).This slight difference suggests lower buffering capacity of the metal oxide coatings for the intermediate grain size and becomes increasingly more pronounced at higher flow velocity (e.g., case at 100 mM NaBr; Fig. 5h).
The experimental results, obtained under broad ranges of hydrodynamic and hydrochemical conditions, indicate that the dynamics and extent of the reaction between the dissolved ionic species and the solid surface are significantly influenced by a complex interplay between grain size specific properties, flow velocities and incomplete mixing in the pore channels.

Modeling results at slow flow velocity
The surface complexation parameters were calibrated using the experimental dataset of the column experiments performed with the slow seepage velocity (1 m/day) and assuming local physical equilibrium (i.e., no DDMT included in the model).Fig. 6 displays the comparison between simulated and measured proton breakthrough curves for the 12 experiments performed at 1 m/day.
The model allowed reproducing the proton breakthrough curves measured for all tested grain sizes and injected solutions with a single set of surface complexation parameters (Table 1).In particular the model captures the dynamics of the proton concentrations that peak after 1 PV and the dependency of the maximum proton concentration on the specific ionic strength conditions (i.e., more proton generated under higher ionic strength conditions for all grain sizes).
Additional simulations were performed by individually inhibiting the contribution of some specific mechanisms (i.e., surface complexation reactions with metal oxide coatings) and by testing the impact of the longitudinal hydrodynamic dispersion (Supplementary Material).The outcomes of these simulations revealed that the difference in proton release from the different porous media, during the experiment performed with a seepage velocity of 1 m/day, primarily results from the distinct amount of oxide coatings among the different grain sizes considered (Fig. S3).The presence of these oxides and the protonation of their surface sites buffer the increase in the pore water acidity resulting from the deprotonation of the quartz surface.The lower amount of oxide in the coatings of the sand of intermediate grain size (Table 1) could explain the systematic higher total amount of eluted protons compared to the smaller and larger gain size (Fig. 5).Furthermore, the model shows that, as a result of the increase in the longitudinal hydrodynamic dispersion term with increasing grain size (i.e., Eq. 3), the average grain diameter has a secondary but substantial impact on the maximum concentration of protons.In fact, lower longitudinal hydrodynamic dispersion led to lower spreading of the proton plumes and, thus, sharper proton fronts exhibiting higher H + concentration maxima.
The model could also capture the transport of ionic species as shown by the agreement between the simulation outcomes and the measured breakthrough curves of sodium and bromide in Fig. 7.
The breakthrough curves of ionic species show that the transport  behavior of sodium and bromide differ.Bromide exhibits a conservative pattern in all flow-through experiments whereas sodium is significantly retarded in particular for experiments performed with lower concentrations of injected NaBr (Fig. 7a-b, 7e-f and 7i-j).The retardation of sodium can be explained by the formation of sodium outersphere complexes on the quartz surface.This surface complexation mechanism leads to the release of protons for charge conservation both at the mineral surface (i.e., negative charge of the surface is compensated by the positive charge of sodium) and in the pore water (i.e., negative charge of chloride compensated by the released protons and by the sodium remaining in solution).The apparent lack of retardation of sodium in the breakthrough for the experiment performed with inlet solutions of highest ionic strength (Fig. 7d, h and l) is caused by the presence of sodium in great excess relative to the adsorption capacity of the reactive sand surface.

Modeling results at high flow velocity: impact of incomplete mixing and mass-transfer limitations
Fig. 8 presents the comparison of the simulated and measured proton breakthrough curves for column experiments performed at higher flow velocity (i.e., 30 and 90 m/day).The simulations were performed using the set of surface complexation parameters determined in the 1 m/day experiments and by accounting for mass-transfer limitations induced by the fast flow using a DDMT approach.
As shown by the solid lines in Fig. 8, the reactive transport model based on a DDMT approach could satisfyingly reproduce the measured proton breakthrough curves and capture the decrease in H + concentration when increasing the flow velocity for all tested grain sizes and ionic strengths.These results demonstrate the capability of the model to describe the complex transport behavior of the aqueous species observed at fast flow velocities by accounting for mass-transfer limitations and surface complexation.In contrast, the model that did not include the DDMT formulation, thus not accounting for the effect of mass-transfer limitations on the transport of aqueous species and surface reactions, significantly overestimated the concentrations of protons at the column outlets for all experiments (dashed lines in Fig. 8).The differences between proton breakthrough curves for flow velocities of 30 m/day and 90 m/day are particularly clear for the smallest and largest grain sizes (Fig. 8a-d and 8i-l) whereas both the simulations and the experimental dataset did not indicate a decrease in eluted protons between the highest tested flow velocities for the intermediate grain size (Fig. 8e-h).Table 2 lists the values of the DDMT parameters fitted to reproduce the measured proton breakthrough curves.
The calibrated DDMT parameters reflect the specific mass-transfer limited behavior observed in the different porous media when increasing the flow velocity.In particular, the fraction of reactive surface area allocated to the immobile domain (φ im ) increases with increasing grain sizes for the 30 m/day cases and is higher than that for the mobile domain under all tested hydrodynamic conditions.The correlation between φ im and the grain size is not observed for the 90 m/day case as φ im only increases for the smaller grain size when changing the flow velocity from 30 to 90 m/day whereas φ im remains stable for the two other grain sizes.Furthermore, the mass transfer coefficient (ξ) Table 1 Surface complexation model parameters used to simulate the interaction between the aqueous charged solutes and the quartz sand surface.increases with flow velocity which is in agreement with previous studies (e.g., Liu and Kitanidis, 2012).
The calibrated reactive transport models allowed the analysis of the spatio-temporal evolution of the injected solutes and the release of protons from the reactive sand surface.Fig. 9 displays the simulated changes in protons and sodium concentrations in the mobile and immobile domains for the experiments performed with a grain size diameter of 0.64 mm and a 10 mM NaBr inlet solution.
For all tested velocities, proton concentrations progressively increase towards the outlet (Fig. 9a-e).In particular, for the fast flow velocity cases, in the mobile domain the protons released from the sand surface travel as a wave characterized by a sharp front and a long tailing during the first pore volume flushed (Fig. 9a, b and d).This wave of protons forms as H + accumulates at the front of the plume, with the advancement of the injected sodium along the column and its contact with a progressively larger amount of reactive surface (Fig. 9f, g and i).More specifically, the change in surface composition simulated by the surface complexation model indicates that it is the deprotonation of the quartz surface upon adsorption of sodium via electrostatic interaction that leads to the acidification of the pore water (Fig. 9k-t).For instance, Fig. 9k and r show that the uptake of sodium via the formation of outersphere complexes is mirrored by the decrease in protonated sites.
Fig. 9 clearly shows the impact of incomplete mixing on the evolution of the proton plume with lower concentration of released protons in the mobile fraction of the pore water when increasing the seepage velocity (Fig. 9a, b and d).These model results also reveal important differences between the mobile and immobile domains in the evolution of the hydrochemical and surface compositions.In particular, the aqueous H + and Na + concentrations profiles are broader and H + concentration is significantly higher in the immobile domain (e.g., Fig. 9c and e compared to Fig. 9a, b and d for protons).Such model outcome reflects the mass transfer limitations between the mobile and immobile domains and provides insights on the effect of incomplete mixing on the sorption and transport behavior of aqueous charged species.In particular, the low volume/surface ratio in the immobile domain leads to the build-up of high proton concentrations and the slow mass transfer from the immobile to the mobile domain results in more pronounced tailing.Conversely, the higher volume/surface ratio in the mobile domain leads to lower H + concentrations and the continuous replacement of the pore water due to the injection of the inlet solution leads to a rapid elution of the proton wave and to sharper H + profiles.The concentration of H + in the immobile domain decreases with time (Fig. 9c and e) as H + diffuses into the mobile domain and is subsequently eluted from the column.Differently, Na + concentration remains steady throughout the entire domain (i.e., mobile and immobile; Fig. 9f-j) as sodium is continuously injected from the column inlet.Regarding the evolution of the surface composition, the change in surface composition is more regular in the immobile domain (e.g., Fig. 9m) compared to the mobile domain where the highest changes in surface composition occur in the vicinity of the column inlet (e.g., Fig. 9l).Finally, the larger surface concentrations in the immobile domain (Fig. 9m, o, r and t) reflect the low volume/surface ratio in this model region.

Conclusions
The impact of mass-transfer limitations is of fundamental importance for the transport of aqueous solutes as it determines the degree of reaction not only between dissolved species but also for solutes undergoing chemical and electrostatic interactions at the solid-solution interface.In this study, we have considered the deprotonation of the quartz surface as an example of surface complexation reaction in saturated porous media.We performed flow-through experiments and reactive transport simulations to investigate the effects of mass-transfer limitations on surface deprotonation induced by the injection of electrolyte solutions.The experiments, were carried out with different grain sizes, seepage velocities and electrolyte concentrations to explore the physical and hydrochemical factors controlling the interactions between aqueous solutes and reactive mineral surfaces.A dual domain mass  transfer formulation combined with a surface complexation model was adopted to reproduce the experiments performed at the high seepage velocities and to capture the effects of mass transfer limitations on the transport of protons and injected ions.
The experimental and modeling results show that incomplete mixing can affect both the extent and the dynamics of surface complexation reactions during the transport of aqueous charged species under advection-dominated conditions.Consequently, such physical mass transfer limitations exert a substantial control on the reactive transport behavior in the flow-through systems.In particular, an increase in seepage velocity substantially hampers surface complexation reactions (i.e., lowering the maxima of eluted proton concentration) due to mass transfer limited transport of the aqueous solutes from the pore water to the reactive surface.Porous media of larger grain size and presenting larger pore throats lead to a higher degree of incomplete mixing at fast flow velocities and, consequently, a more pronounced decrease in surface reactivity.Furthermore, the experimental and modeling results reveal that the difference in the amount of coatings in natural porous media as well as in the hydrodynamic dispersion among grain sizes can lead to a surface reactivity that does not linearly correlate with the measured surface area of the grains.
The comprehensive experimental dataset reported in this study can be used to test different formulations of solute transport affected by heterogeneous reactions, which only apparently are not mixingcontrolled.In fact, the outcomes of these experiments show that surface complexation reactions under fast flow regimes become mass transfer limited due to incomplete mixing in the pore channels.The collected dataset could be of interest to test more complex models, including multiple immobile continua and/or first principle pore-scale simulations solving Stokes flow and advection-diffusion-reaction equations in solid-liquid domains.

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.

Fig. 1 .
Fig. 1.Diagram summarizing the laboratory flow-through experiments performed with three different grain sizes, four different electrolyte concentrations, and three different flow velocities.The fast flow leads to incomplete mixing that impacts the chemical and electrostatic interactions between the aqueous solutes and the reactive sand surface.Such incomplete mixing phenomena are described in the numerical model with a dual domain mass transfer approach, schematically illustrated in the bottom right sketch.

Fig. 2 .
Fig. 2. Proton breakthrough curves measured at the outlet of the columns during experiments performed under different flow velocities, different ionic strengths, and with quartz sand of different grain sizes: d = 0.64 mm (a-d), d = 1.3 mm (e-h), and d = 2.3 mm (i-l).

Fig. 3 .
Fig. 3. pH spatial profiles measured during the column experiments performed with a flow velocity of 1 m/day and the injection of a NaBr 100 mM electrolyte for the three considered grain sizes.

Fig. 4 .
Fig. 4. pH breakthrough curves measured with the non-invasive optode technique at 3 different spatial locations (3.6, 7.2 and 10.5 cm from the column inlet) during the experiments performed with a flow velocity of 30 m/day and the injection of a NaBr 100 mM electrolyte for the three considered grain sizes.

Fig. 5 .
Fig. 5. Maximum measured H + concentration at the outlet of the columns (a-d) and total amount of protons eluted after 10 PV (e-h) as a function of the seepage velocity for the three different grain sizes.The values are normalized by the measured BET surface area.

Fig. 6 .
Fig. 6.Comparison between measured (symbols) and simulated (lines) breakthrough curves of protons at the outlet of the columns for the experiments performed with a seepage velocity of 1 m/day.

Fig. 7 .
Fig. 7. Comparison between the measured and simulated breakthrough curves of ions at the outlet of the columns for the experiments performed with a flow velocity of 1 m/day.

Fig. 8 .
Fig. 8.Comparison between measured and simulated breakthrough curves of protons at the outlet of the columns for the experiments performed with flow velocities of 30 and 90 m/day.The solid lines represent the calibrated DDMT model, whereas the dashed lines represent the outcomes of the simulations without mass transfer limitations.

Fig. 9 .
Fig. 9. Simulated spatial profiles of the change in aqueous species (H + and Na + ) and reactive quartz surface composition for the mobile and immobile model domain during the column experiments performed with grain size diameter of 0.64 mm and an injection of 10 mM NaBr solution.

Table 2
Summary of the calibrated DDMT parameters.