Transport and Surface Complexation in Subsurface Flow-through Systems Open Access

In subsurface environments the cycling of nutrients and the fate and transport of contaminants is controlled by the coupling of physical flow and transport processes with biogeochemical reactions (Steefel et al. 2005; Essaid et al. 2015; Rolle and Le Borgne 2019). Surface complexation reactions occurring at the mineral–water interface in porous media are of key importance in a number of subsurface flow-through systems both under natural conditions and in engineering applications (Davis and Kent 1990; Appelo and Postma 2005; Sposito 2008; Baũuelos et al. 2023). Notable examples include the transport of pollutants in contaminated groundwater flow systems, where the chemical and electrostatic interactions between dissolved contaminants and different subsurface minerals may control the contaminant migration. Other important examples of flow-through systems in which the reaction at solid surfaces significantly impact pore water quality are the interface zones between different environmental compartments. For instance, surface/solution interactions play an important role for the composition of pore water infiltrating in soils and unsaturated zones, in determining the impact on water quality of hyporheic exchange between groundwater and surface water, and at the mixing zone between fresh and salt water in coastal aquifers. Taking into account surface complexation reactions is also important in subsurface engineering applications in which subsurface flow and transport processes are manipulated by human interventions. For instance, different approaches have been proposed for in situ remediation of contaminated groundwater to prevent the spreading and migration of contaminants by controlling the groundwater flow and/or by inducing or enhancing chemical and biological reactions. Other subsurface engineering fields of growing interest are river bank filtration systems and managed aquifer recharge (MAR). In MAR applications, the subsurface serves to store water for future use and typically involve the infiltration or pumping of surface water in subsurface formations. In such context, reactions at mineral–water interfaces are of interest and often of concern since they can control the release and transport of trace elements and pollutants.

Figure 1 provides a schematic illustration of important examples of subsurface processes and applications in which surface complexation reactions at the mineral/water interface are of key interest.

In this chapter, we address the coupling of transport and surface complexation reactions in subsurface flow-through systems. In the first part, we discuss the typical scales at which these processes are investigated and we illustrate some relevant examples at the laboratory and field scales, as well as the reactive transport modeling tools that have been developed to quantitatively describe these systems. In the second part, we take a closer look at the coupling of transport and surface reactions. In particular, we focus on aspects that have recently advanced the understanding of transport processes in saturated porous media. Such aspects include the effect of incomplete mixing, due to the formation of solute concentration gradients within pore channels at high groundwater flow velocities, the impact of dimensionality and the effects of considering transport processes in multiple dimensions (i.e., 2D and fully 3D porous media), the impact of Coulombic interactions on the displacement of charge solutes, and the importance of heterogeneity. The latter is ubiquitous in subsurface formations at all scales and encompasses both physical properties (i.e., hydraulic conductivity, porosity and tortuosity) and chemical and electrostatic properties (i.e., mineral composition, surface charge). We focus on an example of surface complexation reactions on a mineral assemblage, including a quartz and a metal oxide surface, and we illustrate the behavior of such reactive system in different flow-through systems involving the non-trivial coupling between the physical mass transfer processes and the surface reactions.

The study of subsurface processes is performed at different scales. The key scales that are typically considered in the investigation of surface complexation reactions in flow-through systems are the laboratory and the field scales.

Experimental investigations at the laboratory scale have been typically performed in batch and column setups. These systems are complementary and provide the opportunity of exploring different aspects of surface/solution interactions. Batch experiments have been fundamental to illuminate the mechanisms of surface complexation reactions and many important works have been based on these experimental systems (e.g., Schindler et al. 1976; Dzombak and Morel 1990; Hiemstra and van Riemsdijk 1996; Davis et al. 1998; Dixit and Hering 2003). Although this approach probably remains the most widely used, the possibility to investigate surface/solution interactions in flow-through column setups has attracted the attention of a large number of researchers. Column systems offer several advantages to investigate surface complexation in subsurface porous media, such as: (i) higher solid to solution ratios representing more closely natural soils and sediments; (ii) immobile porous media with less abrasion and disturbance of mineral surfaces; (iii) water flows through the porous media, characteristic of natural subsurface environments and avoiding accumulation of reaction products; (iv) the possibility to integrate in a simple setup the effects of hydrodynamics, variable hydrochemistry, mass transfer limitations and surface reactions. Therefore, column experiments have become a fundamental tool to investigate the fate and transport of a wide number of environmental pollutants and major ions. Column experiments have been performed to investigate transport and surface complexations of key environmental contaminants such as uranium (Kohler et al. 1996; Barnett et al. 2000), arsenic (Darland and Inskeep 1997a,b; Radu et al. 2005), chromium (Jardine et al. 1999), molybdenum (Sun and Selim 2019), fluoride (Meeussen et al. 1996), phosphate (Stollenwerk 1996), heavy metals (Papini et al. 1999; Delholme et al. 2004; Hanna et al. 2009), rare earth elements (Iqbal et al. 2023), pH fronts (Scheidegger et al. 1994, McNeece and Hesse 2017; McNeece et al. 2018). Such column experiments have been performed in a wide variety of porous media, including quartz sand, granular media coated and/or mixed with specific mineral phases (e.g., iron, aluminum and manganese oxides, clay minerals, carbonates), as well as natural soils and sediments.

In a study on the mobility of heavy metals, Hanna et al. (2009) investigated sorption of lead and zinc in roadside soils. These soils can be heavily polluted and behave as a source of contaminants potentially affecting other environmental compartments. The authors performed column experiments in which the elution and transport of Pb and Zn were studied in presence of acetic acid and EDTA, at pH 5 and 7, respectively. Figure 2a shows the measured breakthrough curves of major cations and heavy metals at the outlet of the column for the case of transport with acetic acid. These conditions were selected to mimic flushing by infiltrating water containing organic acids used as additives in de-icing salts. The breakthrough of heavy metals was delayed and their peak concentrations were observed after approximately two pore volumes (PVs), whereas the major cations, released by exchange with protons, eluted earlier from the columns (1 PV). The measured Pb and Zn concentrations also show long tailing patterns due to the desorption from the soil surfaces.

Another interesting example is the study of Fakhreddine et al. (2015), who investigated the impact of water composition on arsenic desorption from sediments of a shallow aquifer in Orange County (USA). Their column experiments mimicked the injection of treated recharge water in MAR operations and showed different extents of arsenic elution from the sediments in the case of unamended recharge water compared to the scenarios in which the recharge water was enhanced in calcium and magnesium content through its amendments with lime, dolomitic lime and gypsum (Fig. 2b). The authors reported that, in presence of divalent cations, arsenic elution decreased due to the electrostatic interactions at the surface of clay minerals. The proposed mechanism of surface/solution interaction was the formation of cation-bridging complexes changing the charge of the clay and resulting in surface complexation and reduced elution of arsenate species.

Column experiments have also been used to study transport in subsurface porous media in presence of mass transfer limitations (e.g., Brusseau et al. 1989; Darland and Inskeep 1997a; Fesch et al. 1998). The occurrence of fast and transient groundwater flow, the heterogeneity in the sediment texture and in the pore space, and the presence of preferential flow paths such as macropores and fractures can lead to mass-transfer limited displacement of solutes from the bulk pore water to reactive surface sites. Significant impact of mass transfer limitations on macroscopic transport was observed during the study of uranium migration in sediments from the Hanford site (e.g., Qafoku et al. 2005; Liu et al. 2008). For instance, the work of Liu et al. (2008) investigated desorption of U(VI) from aquifer sediments packed in columns of different dimensions with fine and coarse sediments, respectively. The imposed water flow was dynamic, with intermittent and stop flow events, and resulted in a kinetic desorption of uranium from the sediments (Fig. 2c). Kinetic effects were observed in the breakthrough curves of the eluted U(VI), both from the fine, sieved sediments in the small column and from the coarse, field-textured sediments packed in the larger column. In the latter setup, mass transfer limitation effects were more pronounced and were observed not only for U(VI) but also for different conservative tracers injected in the setup.

Other interesting aspects that have emerged in column transport experiments are the effects of surface reactions on the breakthrough time and characteristic shapes (e.g., self-sharpening and/or rarefaction/diffusive waves) of eluted fronts of charged species such as metals, protons and major ions. The interpretation of observed spatial profiles and breakthrough curves can illuminate the mechanisms of surface/water interactions controlling solute transport in granular porous media (e.g, Bürgisser et al. 1994; McNeece and Hesse 2016). For instance, the interactions of hydrodynamics with surface chemistry reactions can lead to fast pH and metal pulses observed in goethite-coated glass beads columns (Prigiobbe et al. 2012; Prigiobbe and Bryant 2014). As will be discussed below, aspects that have not been covered extensively in well-controlled laboratory experiments are the study of transport and surface complexation in multidimensional setups, the effects of fast flow and incomplete mixing in permeable granular media, as well as the impact of Coulombic interactions between charged species in the pore water (i.e., multicomponent ionic transport).

Field work at contaminated sites shed light on the importance of microscopic surface/solution interactions on the macroscopic fate and transport of contaminants of primary concern. Extensive work focusing on transport and surface complexation reactions in aquifer systems has been done at several contaminated sites such as Cape Code, USA (Kent at al. 1994; Stollenwerk et al. 1996, 1998; Davis et al. 2000), Naturita, USA (Davis et al. 2004; Curtis et al. 2006), Grindsted, Denmark (Kjoller et al. 2004), and Hanford, USA (Zachara et al. 2016). The main focus has been on reactive transport of trace metals and metalloids (e.g., uranium, zinc, nickel, molybdenum, chromium, selenium etc.), as well as on changes in hydrochemistry and transport of other charged species (e.g., phosphate, protons, fluoride, etc.).

As an example of advances from field investigation in contaminated aquifers, the Cape Code research site highlighted important aspects of transport and surface complexation in sandy aquifers. The site was equipped with a dense and high-resolution network of observation wells and thorough hydrogeological and hydrochemical investigation was performed (Le Blanc et al. 1991; Davis et al. 2000). The conditions were particularly suitable to study reactive transport under variable geochemical conditions since a sewage plume was released in the shallow permeable aquifer thus perturbing the original conditions. Beside the sewage contamination, the site was also used to conduct well-controlled tracer experiments. Different contaminants were injected in the sand and gravel aquifer, typically in mixtures with a tracer and the spatio-temporal distributions of the different solute plumes were monitored in the high-resolution network of observation points. Figure 3 shows the outcomes of a tracer experiment in which nickel, complexed with the organic ligand EDTA, and bromide were released and their displacement was monitored in successive sampling and measurement campaigns (Hess et al. 2002).

The tracer studies at Cape Cod have illuminated the distinct transport behavior of different contaminants, undergoing surface reactions with the aquifer solids, compared to the distribution of the bromide tracer. The field experiments also provided unique datasets for the development and testing of reactive transport modeling incorporating surface complexation formalisms. Detailed geochemical characterization and laboratory experiments performed with aquifer sediments and under conditions mimicking those observed in the sand and gravel aquifer were useful to constrain important properties such as the quantity and composition of the surface coatings of the solid grains (Coston et al. 1995; Fuller et al. 1996; Zhang et al. 2011).

Investigation of geogenic contamination also highlighted the key role of solid/solution interactions in subsurface flow-through systems. A notable example is the one of arsenic, whose release and transport are strongly affected by complexation reactions on mineral surfaces. In particular, arsenic has a high affinity for iron (hydr)oxide surface sites that provide large surface areas and charged surface sites. Many factors can control sorption of arsenic to these minerals, including the oxidation state of As, the pH and pore water composition, the concentration of competing ligands, and the mineralogy of the Fe-oxides (e.g., Dixit and Hering 2003; Stachowicz et al. 2006, 2008). Evidence from field studies has shown that pH variations and competitive sorption with other ions in groundwater are important mechanisms, besides the dissolution of iron (hydr)oxides, of arsenic release from aquifer sediments and responsible for the increase of arsenic concentrations in groundwater (e.g., Biswas et al. 2014; Stolze et al. 2019a). In a forced gradient field experiment in a flood-plain aquifer adjacent to a channel tributary to the Red River (Vietnam), Jessen et al. (2012) observed As(III) desorption from the sediment upon infiltration of the channel water, which perturbed the water chemistry and induced a response from the mineral surface. The measurements of arsenic and multiple solutes, including competing sorption ligands, allowed the authors to compare the performance of different surface complexation models to describe multicomponent adsorption under field conditions.

Surface complexation reactions can play an important role also in engineered field applications in subsurface flow-through systems such as managed aquifer recharge (MAR), which aims to balance water disparities by storing in the subsurface excess water volumes during wet periods and making them available for dry seasons. MAR applications are based on different infiltration (e.g., infiltration basins and riverbank filtration) and injection methods (Fakhreddine et al. 2021). Typically, in MAR operations waters with different composition such as the recharge water and the native groundwater are mixed and interact with the sediments. This can be beneficial since mixing and sorption processes can dilute or remove contaminants and improve the quality of the recharge water, but can also lead to the release of geogenic contaminants (e.g., As, F, Mn, Cr, U, Ni, Sb, etc.) via dissolution and desorption. Transport and surface complexation reactions are key processes controlling the mobility and fate of the contaminants. For instance, studies at MAR injection sites, where the release of arsenic was of primary concern (Wallis et al. 2011; Rathi et al. 2017; Fakhreddine et al. 2020), have highlighted that sorption/desorption onto/from iron hydr(oxide) is one of the key factors controlling the concentration of As in groundwater.

Reactive transport modeling is an essential tool to quantitatively describe water flow, solute transport and reactive processes in porous media (Steefel et al. 2005; Li et al. 2017). A number of codes have been developed to simulate reactive transport and complex biogeochemical reaction networks in the subsurface, including surface complexation reactions (Steefel et al. 2015). Examples of these advanced reactive transport simulators include PHREEQC (Parkhurst and Appelo 2013), PHAST (Parkhurst et al. 2010), CrunchFlow (Steefel and Lasaga 1994), HBGC123D (Gwo et al. 2001), PFLOTRAN (Lichtner et al. 2013), OpenGeoSys (Kolditz et al. 2012), MIN3P (Mayer et al. 2001), PHT3D (Prommer and Post 2010), HYDRUS/HPx (Simunek and van Genuchten 2008), ORCHESTRA (Meussen 2003). In recent years, the release of modules such as IPhreeqc (Charlton and Parkhurst 2011) and PhreeqcRM (Parkhurst and Wiessmeier 2015) has facilitated the coupling of advanced multidimensional and multiphysics simulators with the widely used geochemical code PHREEQC (e.g., Nardi et al. 2014; Muniruzzaman and Rolle 2016; Jara et al. 2017; Sprocati et al. 2019, 2023; Ahmadi et al. 2022a). Surface complexation reactions can be incorporated and solved in the framework of modern reactive transport simulators. Surface complexation models based on different descriptions of chemical and electrostatic interactions at the surface/solution interface, such as the constant capacitance model, the diffuse layer model and the triple layer model (e.g., Davis and Kent 1990; Appelo and Postma 1999; Koretsky 2000; Goldberg et al. 2007), have been implemented within the framework of reactive transport simulators. As described by Davis et al. (1998), applications of surface complexation models for mineral assemblage can be divided into two main categories: the CA approach, considering the interactions of the solute species with distinct minerals composing the assemblage, and the GC approach, considering sorption equilibria with generic surface sites (Dzombak et al. 2025, this volume). Both of these approaches have been successfully implemented in reactive transport simulators of subsurface flow-through systems.

A large number of experimental studies have benefited from the unique capabilities of reactive transport simulators to integrate water flow, solute transport, aqueous speciation and surface complexation reactions and, thus, to allow the quantitative interpretation of breakthrough curves and spatial profiles obtained in column experiments. In the following, we provide three examples of reactive transport modeling applied to laboratory datasets from column setups:

• As an example of successful application of the CA approach, Meeussen et al. (1996) simulated reactive transport of fluoride in a column setup with goethite-coated sand by explicitly considering the goethite and silica surfaces. The mechanistic description of transport and an electrostatic surface complexation model describing the reactions at the two surfaces allowed capturing very well not only the transport of fluoride but also important hydrochemistry dynamics such as the propagation of acidity fronts resulting from the surface reactions of protons in the goethite–silica sand column and from the change in ionic strength between injected and resident solutions.

• Parkhurst et al. (2003) simulated column experiments investigating sorption/desorption of phosphorous in sediments collected at the Cape Code field site. They employed an electrostatic surface complexation model based on the GC approach, describing protonation, deprotonation and phosphate binding reactions occurring at a generic surface site. In addition to the SCM, exchange reactions were defined for the major cations. The model captured the breakthrough curves of phosphorus and pH in column experiments under variable hydrochemistry conditions switching between the average composition of treated sewage effluent and the one of uncontaminated groundwater.

• Mass transfer limitations and their effects on surface complexation reactions have been described with reactive transport models based on different approaches. The experiments of mass transfer limited uranium transport in the columns packed with sediments from the Hanford site described above (Fig. 2c) allowed to develop and successfully test chemical versus physical non-equilibrium models. The former approach describes diffusion-limited surface complexation reactions with a distribution of first-order kinetics of surface complexation reactions (Qafoku et al. 2005, Liu et al. 2008). In the physical non-equilibrium approach, surface complexation reactions are considered instantaneous and the mass transfer limitations are described as kinetic exchange between mobile and immobile pore water. Dual-domain (i.e., one mobile and one immobile continuum) and/or multiple continua have been used to represent mass exchange between the advective region of the porous medium and the surface complexation sites in intragrains, aggregates and coating domains (Fesch et al. 1998; Liu et al. 2008; Greskowiak et al. 2011).

Reactive transport simulations are precious tools also at the field scale. In fact, during solute and contaminant transport in aquifer systems, many physical, geochemical, electrostatic and biological processes occur simultaneously and reactive transport modeling offers a unique framework to integrate these processes. Integrating surface complexation mechanisms in field-scale reactive transport simulators is important to understand and predict the migration of contaminants of concern taking into account spatially and temporally variable hydrochemical conditions and subsurface heterogeneity, which are common conditions in field-scale applications. Many studies have demonstrated the need and added value of multispecies reactive transport models implementing surface complexation reactions and allowing to advance and improve the description of subsurface flow-through systems compared to simplified formulations based on empirical distribution coefficients and sorption isotherms (e.g., Freundlich and Langmuir models). For instance, at the Cape Code field mentioned above, reactive transport modeling was successfully applied to describe the pH change due to the release of the sewage plume and its impact on transport of different contaminants such as zinc (Kent et al. 2000, 2007), molybdenum (Stollenwerk 1996), and phosphate (Stollenwerk 1998, Kent et al. 2007). Reactive transport models implementing surface complexation reactions have also been applied to describe systems impacted by transient surface-groundwater interactions, such as those observed at the Hanford site and affecting uranium transport (Greskowiak et al. 2011; Ma et al. 2014), as well as the river–groundwater exchange important for arsenic transport documented at a field site in Vietnam (Wallis et al. 2020). Furthermore, the model-based description of transport and surface complexation reactions is important in engineering applications such as remediation of contaminated sites and managed aquifer recharge (e.g., Masi et al. 2017; Rathi et al. 2017).

Reactive transport modeling is also essential for upscaling and predicting the large-scale impact of surface complexation reactions in heterogeneous aquifer systems. Only a few studies have used reactive transport models to investigate the coupling of transport and surface reactions in multi-dimensional domains with spatially variable distributions of hydraulic properties and minerals, the latter controlling the extent of surface complexation. For instance, Wang et al. (2018) performed reactive transport simulations of Cr(VI) transport and surface complexation in 2D synthetic scenarios characterized by spatial physical and geochemical heterogeneity and investigated the effects of connectivity and spatial correlation lengths on the surface complexation capacity of the system.

In this section we analyze in detail the coupling of physical processes and surface reactions in flow-through systems. We present an illustrative example of surface complexation reactions at mineral surface/solution interface and we explore different aspects of coupling this reactive system in flow-through porous media with different processes and conditions, including flow, transport, and electrostatic interactions, as well as system dimensionality. We then discuss the impact of physical, chemical and electrostatic heterogeneity on reactive transport in a two-dimensional domain representative of an aquifer cross section.

The considered reactive system focused on the interactions between natural sand and protons and major ions in the pore water solution. Although surface complexation reactions in the subsurface have been mostly studied for iron and aluminum oxides (e.g., Davis et al. 1978; Dzombak and Morel 1986, 1987; Meeussen et al. 1996; Appelo et al. 2002; Dixit and Hering 2003; Catalano et al. 2008; Hiemstra 2018), silicate minerals are the most abundant materials in aquifer systems and their importance and non-trivial surface behavior have been investigated in a number of recent studies (e.g., Sulpizi et al. 2012; Pfeiffer-Laplaud et al. 2015; McNeece and Hesse 2017; Garcia et al. 2019; Neumann et al. 2021). Furthermore, quartz grains are typically coated by a complex mineral assemblage of iron, aluminum and other oxides that, despite representing a minor fraction of the solid, impacts the overall surface chemistry behavior and the interactions with the ions in pore water. The study of the reactive system presented in this section is relevant to understand subsurface systems where the propagation of acidic fronts is of concern (Kjoller et al. 2004) and also for the migration of trace elements since pH and major ions control the surface charge behavior of minerals, which, in turn, can greatly impact the dynamics of trace elements displacement.

The conceptual model for the reactive transport system is illustrated in Figure 4. The figure highlights the interactions between the dissolved ions and the mineral surfaces on the sand grains (i.e., quartz and metal oxide minerals). Sodium ions can bind as outer-sphere complexes to the deprotonated quartz surface sites. The released protons form an acidic front that is transported downgradient by the flowing pore water but their transport is not conservative. The H⁺ ions interact with the surface of the metal oxides, which becomes more positively charged. The natural sand surface was described using a “hybrid” approach, inspired both by the component additive model (Davis et al. 1998; i.e., considering two distinct surfaces: quartz and metal oxide) and by the general composite approach (Davis et al. 1998; i.e., considering a generic metal oxide to represent the natural coating of quartz sand). The generic metal oxide surface was attributed a surface charge behavior representative of iron and aluminum oxides (PZC in the range 6.7–10) and significantly different from quartz (PZC in the range 1.8–3).

Figures 4b,c present schematic illustrations of the surface complexation models used for the quartz and metal oxide surfaces, as well as the computed surface charge behavior with varying pH at different ionic strength of the pore water solution (0.1–100 mM NaBr). The adsorption processes on the quartz surface were simulated by applying the Basic Stern model (BS), whereas the triple layer model (TLM) was adopted to describe sorption on the metal oxides. These models account for the impact of electrostatic potential and have been widely used in previous studies on quartz (McNeece and Hesse 2016, 2017; McNeece et al. 2018; García et al. 2019) and oxide surfaces (Koretsky 2000; Villalobos and Leckie 2001; Peacock and Sherman 2004). They allow capturing the impact of major ions on surface protonation/deprotonation reactions over a wide range of ionic strength by distinguishing inner- and outer-sphere surface complexes (Hayes and Leckie 1987; Hayes et al. 1990). The surface properties and reactions considered in the two selected surface complexation models are summarized in Table 1.

Calibration was performed using the particle swarm optimization (PSO) method implemented in the MATLAB^® environment (e.g., Rawson et al. 2016; Stolze et al. 2019b), with parallelization of the simulations to constrain the parameterization on the entire set of experimental observations.

The governing reactive transport equations for the 1D and 2D flow-through setups discussed in the following are expressed as:

where C_i is the concentration of a solute species, S_i is the concentration of a surface species, θ is the porosity, σ_s is the density of the solid grains, D_ij is the cross-coupled dispersion tensor, γ_i is the activity coefficient of a solute species, R_i is the reaction rate and υ_i represents the stoichiometric coefficient. These general mass balance equations describe multicomponent ionic transport and surface complexation reactions in saturated porous media. Coulombic interactions between charged species in the pore water are included considering a Nernst–Planck based formulation of the diffusive/dispersive fluxes (Rolle et al. 2013a, 2018; Rasouli et al. 2015; Tournassat and Steefel 2015):

where D_i is the self-dispersion tensor of species i.

In a one-dimensional system the cross-coupled dispersion tensor, D_ij, is represented by the longitudinal component, whereas in two-dimensional systems, oriented along the principal flow direction, it has two diagonal entries, $DijL$ and $DijT$⁠, representing the longitudinal and transverse cross-coupled dispersion matrices (Muniruzzaman et al. 2014; Muniruzzaman and Rolle 2015):

where δ_ij is the Kronecker delta, which is equal to 1 when and equal to 0 when $DiL$ and $DiT$ are the longitudinal and transverse components of the hydrodynamic dispersion of species i, which were parameterized according to a linear (e.g., De Carvalho and Delgado 2005) and a nonlinear compound-specific formulation (Chiogna et al. 2010; Rolle et al. 2012), respectively:

where $DiP$ is the pore diffusion coefficient, approximated as the product between θ and the self-diffusion coefficient of the different ionic species $Diaq$⁠, Pe is the Péclet number, $Pei=vd/Diaq$⁠, with d being the characteristic grain size of the porous medium. δ is the ratio between the length of a pore channel and its hydraulic radius, and β is an empirical exponent accounting for the incomplete mixing in the pore channels. δ and β were set to 5.37 and 0.5, respectively, based on a compilation of flow-through experiments (Ye et al. 2015a).

The reactive transport models used to study transport and surface complexation under the different flow-through conditions described in the next sections, were implemented in the geochemical code PHREEQC (Parkhurst and Appelo 2013) coupled with MATLAB^® by using the IPhreeqc module (Charlton and Parkhurst 2011) for the 1D systems. A similar approach was employed to simulate the 2D laboratory and field-scale systems, but in these cases PHREEQC was coupled to a multidimensional and multicomponent transport simulator (Muniruzzaman and Rolle 2019) by using the module PhreeqcRM (Parkhurst and Wiessmeier 2015).

Proton propagation fronts were investigated in a series of 21 flow-through column experiments in different porous media, including two types of natural sand and quartz beads, and considering the injection of both salt and acidic solutions (Stolze et al. 2020). Here, we report on a subset of that study, limiting our attention to one porous medium (i.e, quartz sand from Dorsten, Germany) and on the salt (NaBr) injections performed at different concentrations.

Figure 5 illustrates the observations of dissolved species at the outlet of the column setups. pH was measured with an electrode in a flow-through vial and dilution effects were removed by applying the continuous stirred tank reactor equation (Wright and Kravaris 1989). Major cations and anions were determined by inductively-coupled-plasma optical emission spectroscopy and ion chromatography. Figures 5a-d show the proton breakthrough curves at the outlet of the columns. A sharp drop in pH is observed after 1 PV followed by a slow recovery towards the pH of the injected solution. The amount of protons released from the sand surface differs considerably between the experiments performed at different ionic strengths and increases at higher injected NaBr concentrations. Moreover, a maximum capacity of proton release is observed, as shown by the similar H⁺ breakthrough curves in the case of injection of the 10 and 100 mM NaBr solutions. The impact of the surface reactions on solute displacement is also apparent from the measured major ions breakthrough curves (Fig. 5e–h). While the transport behavior of the anion Br⁻ appears to be little affected by surface/solution interactions, the displacement of sodium is retarded. The retardation of Na⁺ is particularly important at low ionic strengths, whereas its effect is less pronounced when higher salt concentrations are injected. Na⁺ retardation is not visible for the highest value of NaBr injected concentration (100 mM, Fig. 5h) although in that experiment the release of protons is the highest (Fig. 5d). This behavior is due to the fact that, at high concentrations, the amount of Na⁺ in solution largely exceeds the amount adsorbed.

The reactive transport simulations performed with the 1D transport and surface complexation model described above captured the experimental observations. The simulations reproduced the pH propagation fronts observed in the different porous media (Fig. 6 in Stolze et al. 2020) and showed the capability of the “hybrid” approach formulation considering both quartz and metal oxide surfaces to precisely simulate the distinct features of the observed pH fronts (Fig. 5a–d). The multicomponent reactive transport model also allowed describing the behavior of the other charged species in the column setup and the impact of ionic strength of the injected salt solution on their transport and retardation (Fig. 5e–h).

Furthermore, the model allowed exploring the surface behavior. Figure 6 displays simulated spatial profiles of different solute and surface species and their temporal dynamics. The profiles show the evolution of pH and protons (Fig. 6a–b), the change in surface composition (Fig. 6c–f) and the change of surface charge density of quartz (Fig. 6g–i) along the length of the column in the case of 100 mM NaBr injection. As the NaBr propagates through the porous medium, it contacts an increasing amount of quartz surface, thus releasing more protons along the flow path. The change in quartz surface composition indicates that the proton release is primarily due to the formation of sodium outer-sphere complexes and the concomitant decrease of SiOH and, to a lower extent, of surface species, compensating the positive charge brought to the surface by Na⁺ (Fig. 6b-e). The effect of Na⁺ on quartz deprotonation is also confirmed by the decrease of the charge density at the quartz 0-plane mirrored by the surface charge density increase at the surface 1-plane and, to a lesser extent, at the surface 2-plane (Fig. 6g–i).

In this section, we examine the same reactive transport system described above with the aim to investigate the effects of mass-transfer limitations on surface complexation reactions in granular porous media. In particular, we consider the impact of incomplete mixing induced by high pore water velocities. Under these conditions, incomplete mixing results in nonuniform solute distribution and concentration gradients within individual pore channels (e.g., Li et al. 2008; Rolle et al. 2013b; Rolle and Kitanidis 2014; Jimenez-Martinez et al. 2015; Wienkenjohann et al. 2023). The effects of these pore-scale mass-transfer limitations on macroscopic transport and reactions at the mineral-solution interface have been recently investigated in a series of 36 column experiments considering sandy porous media with different grain sizes, different groundwater flow velocities and hydrochemical conditions (Stolze and Rolle 2022). The same column setup and sand described in the previous section were used. The focus was on the investigation of incomplete mixing and its impact on surface reactions and pore water quality through the continuous injection of salt solutions (NaBr at different ionic strengths) triggering changes in protonation of the sand surfaces. Figure 7 illustrates the experimental conditions, the incomplete solute mixing in pore channels, and the surface complexation and physical transport models used to simulate the flow-through experiments.

The experiments showed strong deprotonation of the sand surface and the release of H⁺ in the pore water triggered by the contact of the sand grains with the injected NaBr solution. Distinct elution of protons was observed under the different tested conditions, highlighting the significant impact of the pore water flow velocity and of the porous media grain sizes on the surface reactions under flow-through conditions. Regarding the hydrochemical conditions, the highest concentration and amount of protons released in the column experiments were obtained when injecting more concentrated NaBr solutions. Differences with respect to the grain sizes were observed particularly at higher ionic strengths and the highest proton concentrations were observed for the porous medium with the smallest grain size. Finally, focusing on the flow-through conditions, the experimental results show that the applied flow rates strongly affect the concentration and the total amount of eluted protons. For all porous media and salt injections, higher pore water velocities decreased the release of protons from the reactions occurring at the sand surface.

Reactive transport simulations were performed to quantitatively interpret the experimental observations. The reactive transport model was based on the same surface complexation reactive module described above, considering the surfaces of quartz and a metal oxide with a basic Stern model and a triple layer model, respectively. Transport was performed with the classic 1D advection-dispersion equation for the experiment at low flow velocity. To capture the physical non-equilibrium induced by high flow rates, a dual domain mass transfer (DDMT) formulation was adopted (e.g., Valocchi 1985; Liu et al. 2008; Muniruzzaman and Rolle 2017). The surface complexation parameters of the reactive transport models were calibrated using the experimental dataset at slow flow velocity (1 m/day) and assuming local physical equilibrium (i.e., no DDMT included in the transport description). The model allowed reproducing the H⁺ breakthrough curves measured for the different grain diameters and injected solutions with a single set of surface complexation parameters, as well as the transport of the major ions. Na⁺ was strongly retarded at low ionic strengths by the formation of outer-sphere complexes at the quartz surface, whereas Br⁻ showed little interaction with the surfaces of the considered porous media. These simulations revealed that the difference in proton release in the different porous media at the slow flow velocity primarily results from the distinct amounts of oxide coating in the sandy porous media considered. Such oxide coatings and their protonation behavior buffer the decrease in pH caused by the deprotonation of the quartz surface. The surface complexation reactive module developed and validated for the slow flow experiments (i.e., physical equilibrium conditions) was subsequently applied to the high velocity cases. In these experiments the fast flow induced physical non-equilibrium conditions causing incomplete mixing and mass transfer limitations in the pore channels. These conditions resulted in the decrease of emitted proton concentrations both at 30 and 90 m/day. The macroscopic DDMT model coupled to the tested surface complexation reactive module could reproduce the measured proton concentrations for the different grain sizes and ionic strengths tested in the experiments (Fig. 8). Conversely, a classic advection-dispersion formulation of the reactive transport model (dashed lines in Fig. 8) overestimated the H⁺ concentrations measured at the outlet of the columns. This overestimation was particularly pronounced for the cases with larger grain sizes where pore-scale concentration gradients and incomplete mixing are more pronounced.

This example highlights the important interactions of physical flow and mass transfer processes with surface complexation reactions and the capability of a simple macroscopic formulation of reactive transport modeling to capture the key features and concentration patterns induced by microscopic mass transfer limitations and surface/solution interactions.

The reactive transport behavior of the sandy porous medium described above was investigated in another recent study (Cogorno et al. 2022) exploring other aspects of the coupling between physical mass transfer and surface/solution interactions. In that work, the focus was on the impact of system dimensionality and Coulombic interactions between charged species in pore water on surface complexation reactions. To this end, flow-through experiments were designed and performed in a 1D column and a 2D flow-through chamber with the same length and filled with the same sandy porous medium. The experiments were performed with equal flow velocity and inlet mass fluxes to allow a direct comparison of the results. Continuous injection of electrolyte solutions and acidic solutions were performed in both experimental setups. Figure 9a schematically illustrates the experimental setups, their geometry, inlet conditions and outlet measurements of pH and major ion concentrations.

Figure 9b illustrates the conceptual model for the case of continuous injection of the electrolyte solution (100 mM NaBr), highlighting the impact of the salt injection on the proton release from quartz and protonation of the metal oxide surfaces, as well as the Coulombic interactions between the ions in the pore water and on their lateral displacement and gradients at the plume fringe. As discussed above, in the experiments with NaBr injection in the sandy porous medium, sodium sorbs to the quartz surface and displaces protons. The released H⁺ further interact with the metal oxide coatings present on the sand. The electrostatic interactions between the ions in the pore water affect the displacement of Br⁻, which becomes progressively more coupled to H⁺ at the plume fringe, due to the Na⁺ retardation (i.e., formation of outer-sphere complexes at the quartz surfaces). Conversely, the electroneutrality of the pore water is dominated by the ions present in excess concentration (i.e., Na⁺ and Br⁻) in the core of the plume. The effects induced by system dimensionality were quantified by comparing the proton fluxes eluted at the outlet of the 1D column and of the 2D flow-through chamber. Figure 10a shows the flux-weighted boundary-normalized integrated breakthrough curves of H⁺ (i.e., normalized by the sum of the flux-weighted inflow concentrations of H⁺ and Na⁺) at the outlet of both experimental setups. The results are very different and show that the amount of protons released in the 2D chamber is considerably higher than in the 1D column. For instance, the peak proton concentration in the 2D system is almost three times higher than the one observed in the 1D column. Since the injected mass fluxes were identical in both setups, the increased impact of surface complexation reactions in the 2D chamber depends on the greater amount of surface area exposed to the injected Na⁺, which displaces laterally due to the transverse dispersive fluxes in 2D. Modeling was performed by solving the governing multicomponent reactive transport equations in 1D and 2D (Eqs. 1-5) and considering the same component additive surface complexation description of the surface/solution reactions. The results show that the simulations capture the amount and concentration observed in both experimental setups and allowed quantifying the increase in the moles of protons released (+61%) in the 2D chamber.

Figure 10b shows the dynamics of the transverse pH fronts triggered by the continuous injection of the 100 mM NaBr solution in the 2D chamber. The vertical profiles of proton concentration initially increase to reach a maximum after ~7h from the injection of the NaBr solution and then diminish to reach steady-state conditions, at which the effects of surface complexation reactions and proton release from the surface vanish. After reaching the peak concentration, the pH fronts develop a vertically more spread shape with a plateau corresponding to the core of the plume where the surface of the quartz was saturated by the injected Na⁺.

The 2D multicomponent reactive transport model also allows mapping and visualizing the spatial distribution of the key dissolved and adsorbed species and surface properties in the 2D setup. Figure 11 shows the outcomes after 0.6 PV from the NaBr injection. The computed concentrations of the aqueous species are reported in the first row (Fig. 11a–c) and allow visualizing the H⁺ plume induced by the surface reactions and the continuously injected plumes of Na⁺ and Br⁻. The quartz surface generating the acidic plume deprotonates (Fig. 11d) upon sodium outer-sphere surface complexes formation (Fig. 11e), thus retarding Na⁺ displacement (Fig. 11g) and modifying the quartz surface charge density at the 0- and 1-planes (Fig. 11h and 11i). The multicomponent ionic transport features of the system are also apparent from the excess of proton and bromide concentrations (compared to sodium) at the plume fringe due to the Na⁺ delay and the electroneutrality constraint in the pore water solution (Fig. 11j–l).

The impact of Coulombic interactions on the displacement of charged species in the 2D setup is illustrated by the analysis of the transverse dispersive flux components. Figure 12 shows the spatial distribution of the different terms in the Nernst–Planck fluxes (Eqn. 2) for H⁺ after 0.6 PVs. Interesting patterns can be observed: for instance, the migration flux of H⁺ has positive and negative contributions at the inner and outer plume fringe, respectively (Fig. 12c), depending on the electrostatic coupling of the lateral proton displacement with the other ions in solution.

The example of NaBr injection in sandy porous media discussed above reveals interesting effects of multicomponent ionic transport and surface complexation and reveals a different behavior of the same chemical system studied in experimental and modeling setups with different dimensionality.

In this last section, we explore transport and surface complexation in physically, chemically and electrostatically heterogeneous porous media at the field scale. We consider a two-dimensional domain (20 m × 2 m) representing a cross section of an aquifer in which reactive inclusions, with the same properties and reactive surfaces (i.e., quartz and metal oxide, Table 1) considered in the experiments described above, are embedded in an inert matrix. Binary fields are generated with a stochastic approach (Dykaar and Kitanidis 1992) based on a Gaussian covariance model for an auxiliary variable and a thresholding process to obtain 20% of the domain coverage by the reactive inclusions (e.g., Lee et al. 2018; Muniruzzaman and Rolle 2023). We focus on three different scenarios with the following key features:

• Scenario 1: Physically homogeneous domain with the same hydraulic conductivity (K = 6.14 × 10⁻³ m/s) in the inert sandy matrix and in the reactive inclusions;

• Scenario 2: Physically, chemically and electrostatically heterogeneous domain with the same hydraulic conductivity as in the previous case for the reactive inclusions (K = 6.14 × 10⁻³ m/s) and high hydraulic conductivity (K = 6.14 × 10⁻² m/s) assigned to the background inert matrix;

• Scenario 3: Physically, chemically and electrostatically heterogeneous domain with the same hydraulic conductivity as in the previous cases for the reactive inclusions (K = 6.14 × 10⁻³ m/s) and low hydraulic conductivity (K = 6.14 × 10⁻⁴ m/s) assigned to the background inert matrix.

Fixed head boundary conditions are applied at inlet and outlet boundaries to establish values of hydraulic gradient resulting in average flow velocities (~1 m/day) typical of sandy aquifers. The steady-state groundwater flow is solved in a rectangular grid (Δx = 0.1 m and Δz = 0.01 m) and the computed distribution of hydraulic head and stream functions allowed the construction of streamline-oriented grids (Cirpka et al. 1999) used for the solution of the multicomponent ionic transport equations (Eqn. 1). Figure 13 illustrates the three heterogeneous domains, with the spatially-distributed reactive inclusions and the computed streamlines.

In the physically homogeneous domain (Scenario 1), the computed streamlines show a regular pattern and the flow velocity is uniform in the whole domain. In the other two cases the streamline behavior is irregular due to the different hydraulic conductivities of the reactive inclusions causing focusing and defocusing of the groundwater flow (e.g., Rolle et al. 2009; Ye et al. 2015b) and resulting in spatially variable groundwater flow velocities. The considered reactive transport problem was similar to the cases investigated in the laboratory and discussed in the previous sections. A solution of 100 mM NaBr was continuously injected at the inlet boundary of the heterogeneous domains. The surface complexation reactions, occurring in the reactive inclusions, were simulated according to the same reactive module described above and the reactive properties summarized in Table 1. The interplay between the groundwater flow and the spatially distributed reactive zones resulted in different distribution of main chemical species within the heterogeneous domains. Figure 14 shows the spatial distribution of key aqueous and surface species for the different scenarios after flushing 0.4 pore volumes. In the physically homogeneous case (Scenario 1), the fronts of dissolved species are sharp and show that the ions have travelled until ~8 m. The H⁺ distribution reflects the location of the reactive inclusions from where the protons are released (panel a1). The changes in the spatial distribution of surface species are indicative of the effects of surface complexation reactions in the reactive zones, showing a depletion of protonated sites (panel c1) and an increase of sodium surface complexes (panel d1).

In the case with low-permeability reactive inclusions (Scenario 2), the fronts of dissolved species are more irregular as a result of the superposition of the spatially distributed reactive zones and the heterogeneous groundwater flow. The latter cause a farther downgradient propagation of the transported ions compared to the physically homogeneous case due to the faster water flow in the permeable matrix surrounding the low-K zones. The reactive inclusions act as sources of protons, due to the surface reactions triggered by the salt injection, which are slowly released by back-diffusion into the permeable matrix (panel a2). Concentration gradients are also observed for the injected sodium that can slowly enter the low-permeability zones. Mass transfer limitations in this scenario are also noticeable in the spatial distribution of the surface species (panels c2 and d2). Finally, in Scenario 3, the high-permeability inclusions cause a much faster displacement of the solutes in the domain. The protons are completely released from the reactive zones on the left of the domain and a significant amount of H⁺ has already reached the outlet boundary. The main Na⁺ front appears to be located approximately in the middle of the domain, although the ions transported by flow focusing on the main high-K inclusions already reached the domain's outlet. The maps of the surface species also show that surface complexation reactions have already occurred in most of the reactive inclusions present in the cross section (panels c3 and d3). The effect of the different heterogeneity considered in the simulations can be appreciated not only from the spatial distribution of ions and surface species at a given time in the domains but also from the breakthrough of solutes at the outlet.

Figure 15 shows flux-weighted breakthrough curves of pH at the outlet boundary. This integrated measure of proton concentrations also shows important differences between the considered scenarios. In the physically homogeneous scenario, the pH drop due to the deprotonation of the quartz surface is the most pronounced, reaching almost a value of pH = 2.5 at 1 PV, and the breakthrough has a very sharp front followed by a slow recovery. The low-permeability reactive inclusions (Scenario 2) result in an earlier breakthrough and also in higher pH values corresponding to the peak concentrations of protons eluted. Finally, in Scenario 3, the high-K inclusions cause a very early breakthrough and the strong flow focusing events also affect the integrated pH breakthrough curve (Fig. 15). The curve shows a faster recovery of the pore water pH, which tends to recover more rapidly towards the value of the salt solution injected at the inlet boundary.

Subsurface environments are dynamic systems in which water flow, physical transport processes and biogeochemical reactions interact and determine contaminant fate and transport. Reactions occurring at the interface between mineral surfaces and pore water in porous media control the attenuation and release mechanisms of contaminants and major ions and, ultimately, the groundwater quality. Therefore, adopting a flow-through perspective on the study of surface complexation reactions is of pivotal importance for a quantitative understanding of natural and engineered subsurface environments. The inherent heterogeneity of subsurface systems and the spatial and temporal variability of hydrochemical conditions limit the applicability of (over)simplified models of sorption processes (e.g., Bethke and Brady 2002; Goldberg et al. 2007) and require rigorous approaches for the description of transport processes and surface reactions. In this work, we have provided examples highlighting the importance of surface complexation in flow-through porous media, as well as the typical scales and investigation approaches to study reactive transport in subsurface systems. For a selected case of surface complexation reactions on quartz and metal oxide, we have shown the results of coupling such reactive system, considering both chemical and electrostatic interactions, with flow and transport processes. The latter include fast flow velocities leading to incomplete mixing, dimensionality effects and the impact of physical, chemical and electrostatic heterogeneity.

Many aspects of surface complexation reactions in flow-through systems still need to be explored and/or require further investigation. In the following we outline some research directions that, in our view, are of great interest and may be promising avenues for future investigation.

There has been an evolution in the study of surface complexation mechanisms from selected contaminants and selected mineral surfaces towards considering more complex hydrochemical conditions including multicomponent adsorption and ligand competition reactions. An aspect that has not received broad attention yet is the study of such complex hydrochemical conditions in system of mineral aggregates in which the pore water composition can differently affect the charge and surface behavior of different minerals, which can also interact and change, in turn, the pore water chemistry. Another aspect that has been explored (e.g., Darland and Inskeep 1997a, Liu et al. 2008) but is not commonly investigated is the impact of varying hydrodynamic conditions, which in laboratory setups could be mimicked with stop flow events, abrupt changes, and cycles of different conditions. Testing such systems will allow representing more realistically dynamic natural environments and engineered systems and capturing key aspects of the coupling between mass transfer processes and surface reactions in flow-through systems.

Furthermore, despite decades of research, datasets on surface complexation are still limited to a few minerals and contaminants. The latter are most often heavy metals and metalloids; however, surface complexation approaches could be extended to other charged species whose fate and transport in the subsurface is of growing interest and concern such as rare earth elements (Neumann et al. 2023) and persistent and mobile organic compounds like PFAS, pharmaceuticals and pesticides (e.g., Luo et al. 2022; Cogorno and Rolle 2024a,b).

Transport in natural porous media is inherently multidimensional. Therefore, experiments in multidimensional setups will allow exploring the effects of flow and transport in 2D and fully 3D porous media on surface complexation reactions. Multidimensional setups will enable capturing important feedbacks and coupling effects between mass transfer processes and hydro-and surface chemistry that cannot be observed in classic 1D columns, such as the impact of lateral acidity and salinity gradients and the effects of electrostatically coupled displacement of ionic species on surface reactions. Multidimensional setups will also allow investigating surface complexation reactions in presence of complex patterns of water flow and solute transport induced by porous media heterogeneity and anisotropy (e.g., Chiogna et al. 2014, 2015; Ye et al. 2016), as well as chemical (e.g., Battistel et al. 2019, 2021; Salehikhoo et al. 2013; Stolze et al. 2022) and electrostatic (Muniruzzaman and Rolle 2021) heterogeneity.

The investigation of surface complexation both at laboratory and field scales has been mostly focused on saturated porous media. Future studies can explore multiphase systems, such as the vadose zone (e.g., Hopmans and van Genuchten 2005; Jimenez-Martinez et al. 2020; Stolze et al. 2024), unsaturated mine wastes (e.g., Acero et al. 2009; Muniruzzaman et al. 2020 and 2021; Muniruzzaman and Pedretti 2021; Seigneur et al. 2021) and evaporative and interface environments (e.g., Or et al. 2013; Haberer et al. 2015; Ahmadi et al. 2022b, 2024). In these systems, water chemistry and contaminant fate and transport can be controlled by two fluid phases and exchange of volatile components at the gas/water interface can affect reactions at the water/mineral interface. Surface reactions, in turn, could impact the distribution and exchange of species between the aqueous and gas phases in unsaturated porous media. Another important aspect that has received limited attention is the effect of temperature fluctuations on surface complexation reactions. Indeed, temperature can significantly impact the sorption and retention of contaminants and major ions by modifying the prevailing thermodynamic conditions, as well as the surface complexation parameters at the mineral–water interface (e.g., Lützenkirchen and Finck 2019; Machesky et al. 2023). These temperature dependent effects are of growing interest for several subsurface applications, such as hydrothermal systems and nuclear waste repository (e.g., Liu et al. 2015; Montavon et al. 2023).

Field studies on contaminant transport and surface complexation have also evolved in the last years. The investigation of groundwater plumes, which was instrumental to develop the understanding of surface/solution interactions in aquifer systems (e.g., Davis et al. 2000; Stollenwerk 1998; Curtis et al. 2006), has been accompanied by studies focusing on the interfaces between different environmental compartments (e.g., Zachara et al. 2016; Wallis et al. 2020). Besides transport and surface complexation reactions under natural flow conditions, engineered setups are also of interest. For instance, managed aquifer recharge systems, increasingly used to mitigate water scarcity, represent a unique opportunity for future understanding of the role of surface/solution interactions in aquifer systems. In fact, the controlled hydrodynamic and hydrochemical forcing of subsurface conditions, the extensive network of monitoring locations, and the multiple cycles of injection/extraction can illuminate spatial and temporal evolution of surface reactions under flow-through conditions.

Finally, the development of reactive transport simulators has greatly contributed to the quantitative understanding of surface complexation in flow-through systems. Envisioned directions of further modeling development include: (i) the coupling of pore scale and hybrid scale flow and transport models with surface complexation reactions (Maes and Geiger 2018), (ii) the extension and unification of databases and adopted formulations of surface reactions (Lützenkirchen et al. 2015), (iii) the improvement of the efficiency of multi-physics-chemistry codes and their capability to perform high-resolution simulations in multiple dimensions, (iv) the development of upscaling approaches and stochastic simulations to transfer detailed small scale process knowledge of transport and surface/solution reactions to larger field scales and to rigorously quantify uncertainty of model predictions, (v) the applications of data assimilation and inverse modeling techniques, often used in subsurface and contaminant hydrology, to image subsurface physical and chemical heterogeneity (Fakhreddine et al. 2016) and to investigate multi-solutes and multi-minerals transport and surface complexation in physically, chemically and electrostatically heterogeneous flow-through systems.

These (and many other) directions of future investigation can contribute to shed new light on the interplay of flow and transport processes with surface complexation reactions in flow-through porous media, thus highlighting and quantitatively capturing the critical role of these coupled reactive transport processes for the understanding, management and remediation of subsurface systems.