Authors: Patrick Ober, Willem Q. Boon, Marjolein Dijkstra, Ellen H. G. Backus, René van Roij, Mischa Bonn
Categories: Article, Fluid dynamics, Geochemistry, Optical spectroscopy, Surface spectroscopy
Source: Nature Communications
The charging and dissolution of mineral surfaces in contact with flowing liquids are ubiquitous in nature, as most minerals in water spontaneously acquire charge and dissolve. Mineral dissolution has been studied extensively under equilibrium conditions, even though non-equilibrium phenomena are pervasive and substantially affect the mineral-water interface. Here we demonstrate using interface-specific spectroscopy that liquid flow along a calcium fluoride surface creates a reversible spatial charge gradient, with decreasing surface charge downstream of the flow. The surface charge gradient can be quantitatively accounted for by a reaction-diffusion-advection model, which reveals that the charge gradient results from a delicate interplay between diffusion, advection, dissolution, and desorption/adsorption. The underlying mechanism is expected to be valid for a wide variety of systems, including groundwater flows in nature and microfluidic systems.
**Subject ** Geochemistry, Optical spectroscopy, Surface spectroscopy, Fluid dynamics
Reactions at the interface between a charged solid surface and a flowing fluid play a key role on macroscopic scales in geochemical cycles^1–4^, as well as on microscopic scales in micro- and nanofluidic systems^5^. They are also central in technological applications in areas as diverse as froth flotation^6^, electrophoresis^7^, water desalination^8^, soil remediation^9^, and even dentistry^10^. Only a few years ago, Lis et al.^4^ presented the first experimental evidence that the surface potential of mineral surfaces (silica and calcium fluoride) in contact with water changes substantially when liquid flow is applied. This observation implies that fluid flow can directly affect a chemical equilibrium. For the mineral fluorite (CaF2) under acidic conditions, the change of the surface charge upon flow was argued to be due to the dilution of reactive ions partaking in the surface charging reaction^4^. This dilution is caused by the concentration difference between the fresh solution from the reservoir and that in the flow channel. The dilution changes the charging equilibrium and increases the surface charge, which manifests itself by an increase of the v-SFG response. A similar explanation has also recently been given by Xi et al.^11^. Even though the flow-induced disturbance of the charging equilibrium qualitatively explains several experimental features, there are also inconsistencies. For instance, the question remains how a concentration difference between the reservoir and channel is generated. The surface charging reaction itself cannot supply the excess of ions in the channel over the bulk concentration, as the number of reactive ions on the surface is too small to support a steady flux of ions, especially over many flow cycles. Therefore, a fully self-consistent quantitative model is missing. In fact, several different hypotheses have been put forward for the change of the v-SFG signal upon flow, for instance, a surface conduction model by Werkhoven et al.^12^, which would lead to a surface charge gradient and a 1D toy model by Lian et al.^11^ that considers net dissolution as the driving force. Moreover, Schaefer et al.^13^ concluded that a change in surface potential could be explained by a flow-induced change of the concentration of ions that screen the surface charge. A flow-dependent surface charge not only paves the way for novel electrokinetic effects^14^, but it also has important consequences for the interpretation of zeta potential measurements. However, the mechanism by which the flow alters the surface charge is poorly understood as of yet^15^.
One needs an interface-specific method to obtain a fundamental understanding of the events occurring at interfaces. At solid-water interfaces, vibrational sum frequency generation (v-SFG) is a well-established method that is interface-specific and provides molecular information^2,4,13,16–27^. The interface-specificity of v-SFG results from the selection rule that inversion symmetry must be broken for a v-SFG signal to be generated. If the solid–liquid interface is charged, the resulting electric field aligns the dipole moments of the interfacial water molecules. In addition, the electric field also polarizes the water molecules some nanometers into the solution. Since both phenomena break the symmetry, the strength of the observed v-SFG signal increases with the average polarization and orientation of the water molecules. Therefore, the v-SFG response is a measure of the electric field at the surface (see Supplementary Note 1 for details)^2,4^. By Gauss’ law, this field is connected to the surface charge and is directly related to the surface potential.
In this work, we investigate the flow-dependent surface potential by measuring v-SFG spectra at several positions in a channel with a length of 2.48 cm and a radius of 0.24 cm. The local changes in the surface potential upon flow are tracked through the intensity of the OH-stretch band in v-SFG spectra, which reports on the alignment and polarization of interfacial water molecules and, therefore, the surface charge. Using this approach, we show that the surface charge not only depends on the flow rate but also on the position in the channel, along the flow direction. We explain these observations both qualitatively and quantitatively in terms of coupled reaction–diffusion–advection processes that give rise to flow-dependent heterogeneous ion concentrations, causing the surface charge of CaF2 to vary spatially along the flow direction in the channel. Our theoretical model reveals that the macroscopic lateral gradient of the surface charge density (and hence of the v-SFG response) along the mineral can be traced to a gradient in the interfacial fluoride concentration. This fluoride gradient also causes a gradient in the dissolution rate along the mineral surface. Such a gradient of the net dissolution rate in confined media is argued to be important for understanding geochemical processes, for instance, in groundwater flows^28,29^. In addition, our approach of combining quantitative theory with position-resolved v-SFG experiments allows determining whether the surface is charged by desorption or adsorption when the dissolution mechanism is known.
Similar to Lis et al.^4^ and Xi et al.^26^, we conducted flow experiments at the CaF2-water interface under acidic conditions, first at a single position in the channel. Figure 1a illustrates our experiment. Visible (Vis) and infrared (IR) laser pulses overlap in space and time at the CaF2-water interface, such that a nonlinear optical v-SFG process generates a response at a frequency equal to the sum of the two incident frequencies. Due to spectroscopic selection rules, this process is forbidden in centrosymmetric media such as bulk water. In contrast to the bulk, centrosymmetry is broken at the interface, which renders the v-SFG signal interface-specific^16,17^. Moreover, the resonance of the IR pulse with the water stretching mode at around 3300 cm^−1^ enhances the overall process and provides additional molecular specificity^2,4,18–25^.
Fig. 1 Studying the CaF
2-water interface by v-SFG.a Illustration of the measurement method. In vibrational sum frequency generation (v-SFG) spectroscopy, a visible (Vis) and infrared (IR) pulse overlap in space and time at the interface of interest, generating radiation at the sum of the frequencies (SFG). If an IR pulse is in resonance with the OH-stretch vibration of water molecules, the obtained v-SFG spectra provide information on the orientation and polarization of interfacial water molecules. The illustration shows water molecules for the first few hydration layers as well as a schematic distribution of arbitrary ions (anions in yellow, cations in green). At the charged CaF2-water interface, the orientation and polarization of the interfacial water molecules is determined by the surface charge. Thus, the v-SFG signal can be used as an indirect measure for the surface charge. b v-SFG spectra in arbitrary units (arb. units) in the OH-stretch region of the CaF2-water interface at pH 3 (1 mM HCl and 1 mM NaCl) under flow-off conditions (blue) and flow-on conditions (red). c Time trace of integrated SFG spectra (black circles) at the CaF2-aqueous solution interface at pH 3 (1 mM HCl and 1 mM NaCl) with one flow on-off cycle (black curve). Each circle represents one spectrum, integrated between 2650 and 3600 cm^−1^. The solid orange line is a ten-point moving average to guide the eye. The vertical blue and red lines highlight the steady-state regimes for flow-off and flow-on conditions. As the intensities are hardly changing over time, a steady-state can be assumed.
Our measurements are performed in total internal reflection geometry to increase the signal and reduce the acquisition time, allowing time resolution on an approximately second timescale. Figure 1b shows two v-SFG spectra of the CaF2-water interface for a 1 mM NaCl and 1 mM HCl (so pH 3) solution, once under flow-off and once under flow-on conditions. The measured signal is plotted vs. the wavenumber of the IR pulse being in resonance with the OH-stretch vibration of water molecules. The flow-on (red) spectrum in Fig. 1b is 40% higher in intensity compared to the flow-off (blue) spectrum in Fig. 1b. This increase can be correlated to an increase in surface charge (see Supplementary Note 1 for details). Thus, we reproduced that flow increases the surface charge at the CaF2~-water interface under acidic conditions^4,26^. The underlying mechanism for the flow-induced increase in the surface charge will be considered in the section ‘Model of surface chemistry’. When we calculate the intensity of the two spectra of Fig. 1b by integrating over the spectral area, we can track their time-dependence upon switching on and off the flow with a resolution of seconds, shown in Fig. 1c, where each black circle stems from one spectrum. We note that the v-SFG response is constant before the liquid flow is switched on, which therefore represents a steady state. Upon applying flow, another steady state with an increased v-SFG response is reached within about 5 minutes. When the flow is turned off again, the initial steady state is restored but only after about 10 minutes. As explained above, the v-SFG response can be used as a measure for the surface charge. Thus, the flow-induced change in the v-SFG signal can be correlated to a reversible change in surface charge due to flow. The observation of reversible changes in the v-SFG response and thus in surface charge due to liquid flow is consistent with the studies of Lis et al.^4^ and Xi et al.^26^.
In order to gain an understanding of the interfacial events along the mineral surface upon flow, we compare the changes in v-SFG response at different positions along the flow channel. Another approach is to compare the change in v-SFG signal upon flow at one position, which is not the center of the flow cell, and reversing the flow direction, i.e., interchanging the role of the inlet and outlet as illustrated in Fig. 2a. Following the second approach, time traces of the integrated v-SFG spectra at three different positions with flow on/off cycles in clockwise and counter-clockwise flow directions are shown in Fig. 2b. First, it can be seen that the v-SFG response increases upon flow at every position. In addition, we observe that at the center of the flow channel (orange spot and trace in Fig. 2b) essentially identical changes of the v-SFG intensity for both flow directions are observed, which is fully consistent with the study of Lis et al.^4^. Here we do not obtain additional spatial information as reversing the flow direction interchanges the center with itself due to the flow cell’s symmetry. In contrast, at the two opposite sites at equal distance from the center (8 mm), the reversible changes upon flow depend on the flow direction. At a clockwise flow direction, the increase in the v-SFG response at the red position in Fig. 2 is substantially and reproducibly higher than at a counter-clockwise flow direction, while the opposite is the case at the other side of the flow cell (blue spot and blue trace). We also note that the magnitude of the increase in v-SFG response in the center is approximately in between that at the inlet and outlet. Thus, we conclude that the increase in the v-SFG response upon flow decreases in the flow direction, such that our measurements show that along the mineral there is a flow-induced surface charge gradient in the flow direction.
Fig. 2 Position and flow direction-dependent change in the v-SFG response upon flow.a Schematic representation of the flow set-up. The aqueous solution reservoir (~100 mL of a 1 mM NaCl and HCl (pH 3) solution) is pumped via tubes and the flow cell using a peristaltic pump. The flow channel is 24.8 mm long, 4.3 mm deep, and 4.8 mm wide. The flow direction can be changed from clockwise to counterclockwise. The v-SFG spectra were recorded at three spots in the channel. One spot in the center (yellow) and two points separated by 8 mm from the center in the two directions of the flow channel (blue and red, respectively). b Time trace of integrated SFG spectra (black circles) at the interface of the solution and calcium fluoride for several flow on–off cycles (black curve). As indicated by the arrows and the different signs of the flow rate, different flow directions were employed. The colored solid lines are ten-point averages to guide the eye and highlighting at which spot in the flow channel the spectra were recorded. For the time trace measured at the center, the intensity levels of the spectra with (Ion) and without flow (Ioff) are highlighted by arrows.
To further investigate the proposed gradient along the flow channel, we quantitatively compare the flow-induced changes at different positions along the centimeter-sized flow channel for the same flow direction. To do so, we focus on the relative increase of the v-SFG response under flow-on conditions (compared to flow-off) as a function of the distance from the center of the channel. In addition, we study the influence of the flow rate. The relative increase in v-SFG intensity (IOn−IOffIOff) accounts for slight changes in the alignment that necessarily occur when changing the position of the flow cell. The position-dependent relative increase in intensity is determined by averaging the intensity in steady states, as indicated by the marked regimes in Fig. 1c. The obtained data are shown in Fig. 3, from which it is clearly shown that the intensity throughout the channel increases in the flow-on state, the more so at higher flow rates and closer to the inlet, up to 45% at the inlet for the highest flow rate of 6 mL/min (shear rate ~9 s^−1^) that we consider here. Note, however, that the intensity also increases significantly, by about 10%, at the outlet at the lowest flow rate of 1 mL/min (shear rate ~1.5 s^−1^). With these flow rates, we capture a laminar flow with corresponding Reynolds numbers between ~5 and ~25. For each flow rate, the increase of the intensity upon applying flow varies monotonically between inlet and outlet. To the best of our knowledge, other spectroscopic flow experiments conducted so far only varied the flow rate in ranges where no dependency has been observed^4,26^.
Fig. 3 Gradient in the increase of the v-SFG intensity for different flow rates.At different positions, the relative increase in the v-SFG intensity at the CaF
2-water interface is calculated based on averaged steady-state spectra under flow-on and flow-off conditions as marked in Fig. 1c and Supplementary Fig. 1, for the three indicated flow rates. The open circles represent the relative increase upon flow and the filled symbols the mean at their corresponding positions. The error of the individual ratios (open circles) due to small fluctuations within the steady-states of one flow-on/off cycle (see Fig. 1c of the main text and Supplementary Fig. 1) is negligible compared to the spread. The solid lines are guides to the eye. The solution used was a 1 mM NaCl solution in pH 3 HCl with an additional 1 µM NaF. The addition of a very low NaF concentration constrains the number of unknown parameters in our model. Comparing the black data points here with the data points in Supplementary Figure 2 see shows that the added 1 µM NaF does not noticeably influence the interfacial events. The position is given as displacement from the center, with values increasing in the direction of flow, indicated by the arrow.
In summary, the v-SFG data reveal that flow triggers an increase of surface charge at every position in a channel. This increase varies monotonically along the channel length so that it is largest at the inlet and smallest at the outlet. This change in surface charge upon flow is reversible. In addition, we found that the surface charge variation is also a function of the flow rate. To explain our observations quantitatively, we need a model that combines microscopic surface chemistry and macroscopic processes such as advection and diffusion in a flow-system, which goes beyond the qualitative explanation by Lis et al.^4^. In the following section, we will introduce a self-consistent model that extends that of Lis et al.^4^ and describes not only the flow-dependent surface charge but also the dependency on flow rate and position.
Here we will briefly discuss the qualitative features of our model. We incorporate two reactions that take place at the CaF2-water interface (i.e., the wall of the flow channel): (i) The surface charging reaction of Eq. (1) in which only F^−^ desorbs from the surface to leave behind a positively charged CaF^+^ unit on the surface, together with its back-reaction that involves F^−^ adsorption. This reaction is responsible for the surface charging, and has been considered responsible for the flow-induced change of the surface charge in previous work^4^. However, as demonstrated below, a second reaction is necessary to explain the (ii) The dissolution reaction of Eq. (2) in which a neutral CaF2 unit dissolves from the surface in the form of three dissolved ions (one Ca^2+^ and two F^−^), together with its back-reaction, i.e., precipitation. This reaction is the source of F^−^ ions and does not change the surface charge. The reactions are given
Both reactions have a fluoride ion in common, which couples them. Figure 4a illustrates this view of the surface chemistry. While Lis et al.^4^ and Xi et al.^26^ only considered the charging reaction from Eq. (1) to qualitatively explain the flow-induced changes in v-SFG spectra, we argue that a consistent and quantitative description of the experimental results also requires to incorporate the dissolution as described by Eq. (2). The charging reaction of Eq. (1) alone could not serve as a steady-state source for fluoride ions. During the continuous dilution with fresh solution from the reservoir in the flow-on period, the surface charge and thus the v-SFG signal would have to keep increasing, which is not observed experimentally as a steady-state v-SFG signal is reached after about 5 min. Therefore, our model considers the CaF2 surface not only as an adsorption–desorption site for fluoride anions but also as an ion source for charge-neutral dissolution and precipitation of Ca^2+^ and two F^−^ according to Eq. (2). In fact, the dissolution reaction is not only a theoretical necessity but is also supported by the presence of Ca^2+^ cations in flow experiments with CaF2^4^.
Fig. 4 Model of surface chemistry along the flow channel.a Illustration of the two interfacial reactions. CaF
2is not only an adsorption–desorption site for fluoride anions (Eq. (1)), which alters the surface charge, but also experiences charge-neutral dissolution–precipitation processes (Eq. (2)). The shared fluoride anion couples the two reactions. b Illustration of the flow channel. Without flow (upper part), dissolution creates an enhanced F^−^ concentration in steady-state, constant over the dissolving mineral surface, if the concentration reaches the saturation concentration. Dissolution rates and surface charge are also homogeneous along the mineral’s surface. However, upon flow (lower panel), the concentration is lowered everywhere in the channel, the more so closer to the inlet. This leads to an increase in the surface charge and higher net dissolution than in no-flow conditions, both becoming smaller at the outlet.
The appearance of a surface charge gradient upon flow can be explained when considering a flow-induced gradient in fluoride concentration along the mineral surface, as illustrated in Fig. 4b. This gradient forms when the laminar flow transports dissolved fluoride anions along the interface, which leads to a build-up of fluoride ions at the outlet and a positive gradient in fluoride concentration in the flow direction. Moreover, because the solution flowing into the channel from the reservoir has a lower fluoride concentration, the fluoride ions are also diluted by the fresh solution, the more so the closer to the inlet. Thus, flow reversibly lowers the fluoride concentration everywhere in the channel, yet more so closer to the inlet. As the surface charge increases with lower fluoride concentration according to Eq. (1), the surface charge becomes larger at positions closer to the inlet. This qualitatively explains the gradient of the increase in surface charge along the surface upon flow and thus in the v-SFG response.
We model the channel as a cylinder of length 2L and radius R, and introduce the radial and axial coordinate r and x, respectively. The charged and dissolving mineral surface of the channel is located at r=R for x∈[−H,H] (with H<L). The channel inlet and outlet are at x=±L, and ∣x∣≫L denotes the interior of the two reservoirs in which we specify bulk concentrations and bulk pressure. The model geometry is depicted in Fig. 5a. Throughout, we will assume azimuthal symmetry and consider water to be an incompressible fluid with viscosity η and dielectric constant ε at fixed room temperature T. We are interested in finding the areal density σ(x) of charged surface groups on the mineral surface, the fluid velocity u(x,r) in the channel, and the concentrations ρi(x,r) of ion species i=H,Cl,Na,Ca,F of which, however, only calcium and fluoride will turn out to be relevant. The system is solved for steady-state without any explicit time-dependence. The fluid flow is driven by an imposed static pressure drop between the two reservoirs such that the pressure profile P(x,r) is also to be determined. We model the surface chemistry in terms of the dynamic Langmuir equation
that accounts for the charging reaction of Eq. (1). Here kdes and kads denote the desorption and adsorption rate constants of fluoride, respectively, and Γ denotes the surface density of chargeable groups. The surface charge density is thus denoted by eσ(x), with e the proton charge. In the steady-state conditions of interest here, ∂tσ=0, we can solve the dynamic Langmuir Eq. (3) to obtain the Langmuir isotherm
Fig. 5 Channel geometry and numerical solutions to the reaction–diffusion–advection Eqs. (3)–(7).Illustration of the flow channel as used for the numerical calculations in COMSOL with resulting concentration profiles. The center of the flow channel is at the origin (x=0,r=0). The channel length and diameter are 2L=40 mm and 2R=4.8 mm, respectively, and the dissolving surface length is 2H=25 mm. The color represents the fluoride concentration. In a, the result of the numerical calculations in the full three-dimensional model geometry can be seen in flow-off conditions. Panel (b) is a zoom-in on the channel under flow-off conditions (top) and flow-on conditions (bottom) at 1 mL·min^−1^ in the positive x-direction as also indicated by the parabolic Poiseuille flow profile (white arrows). With flow, the concentration in the center of the channel essentially equals that in the reservoir. Nonetheless, a significant diffuse boundary layer remains at the dissolving surface, with a boundary layer thickness of δ(0)≈0.5 mm. The thickness of the boundary layer scales as δ∝(x+H)1/3. The boundary layer, which is much smaller than the channel radius, reduces the effective advection rate by orders of magnitude, as the fluid velocity is low near the channel wall.
Note that when the surface charge is in steady-state, ∂tσ=0, there is no net source of charge directed into the solution, and hence the flux of calcium ions JCa equals twice the flux of fluoride ions, 2JF=JCa. The fluoride concentration at the surface needs to be determined and depends not only on the position x but also (implicitly) on the flow velocity u. The profiles of the ionic concentrations, the pressure, and the fluid velocity follow from the advection–diffusion equations (Eq. (5)) for ionic transport and the incompressible Stokes equation (Eq. (6)) for the fluid velocity u,
Here D is the diffusion constant that is assumed to be equal for all ionic species. It is computationally not feasible at the cm length scales of the experiment to resolve the electric double layers at the solid–liquid interface on the nm-scale. Hence, we neglect a conduction term to Eq. (5) and a coupling to the Poisson equation that would be required to account for electrostatics. Given that the concentration profiles of sodium, hydrogen, and chloride ions have no sources and sinks and are not affected by conduction within this approximation, we can assume them to be spatially constant, with values set by the reservoir. Moreover, since ∂tσ=0 and JCa−2JF=0, we can reduce the number of concentration profiles further by imposing local charge neutrality and hence setting ρF=2ρCa in the simplest version of the theory that we focus on in the main text. However, in Supplementary Discussion 1 we go beyond the assumption of local charge neutrality and allow for electric double layers by solving a similar, but orders of magnitude smaller, system with a conduction term and the Poisson equation; we find that the non-electrostatic version of the theory captures the key physics. To solve for the fluoride concentration, we seek solutions to Eqs. (5) and (6) that satisfy the following boundary conditions. Deep in the two reservoirs, we impose identical ionic bulk concentrations ρi,b for each species in the reservoir, and a pressure that differs by an imposed pressure drop ∆P. On the solid channel wall we apply no-slip boundary conditions for the flow velocity and vanishing normal fluxes for fluoride for H<∣x∣<L. On the dissolving wall, at ∣x∣<H, we impose the dissolution and precipitation reaction (2), resulting in n⋅∇ρF(x,R)=JF(x)/D, where JF(x) is the normal flux of fluoride ions into the solution given by
Here we introduced the surface normal n pointing into the solution, the dissolution rate constant kdis of the CaF2 units, and the precipitation rate constant kprec, in line with the reaction of Eq. (2); the prefactor 2 is a stochiometric constant resulting from the fact that two F^−^ ions go into the solution when one CaF2 unit dissolves. The ratio of the rate constants defines the saturation concentration ρF,sat=(2kdis/kprec)1/3 upon which no further dissolution occurs. Note that the local concentrations govern the net precipitation rate at the surface, which couples to the flow through the advection term in Eq. (5), such that a flow-dependent surface charge is obtained if the dissolution and charging reactions share an ion. For CaF2 the shared ion is F−.
Due to dissolution, the net flux of ions into the flow channel creates, even without flow, a reaction–diffusion steady-state concentration profile that extends over the length of the channel and into the two reservoirs, as illustrated in Fig. 5a. With flow, a lateral and a normal gradient of fluoride ions develop because the fluid flow increases the rate of transport out the channel, the more so on the central axis of the channel where the flow velocity is largest, as is evident from Fig. 5b. The radial concentration profile near the surface is well described by diffuse boundary layer theory (Supplementary Discussion 3), and has a thickness δ(0)≈0.4 mm at the center (x=0) of the channel, and gets thinner/thicker upon approach of the inlet/outlet following approximately the scaling δ(x)∝(x+H)1/3. Because the diffuse boundary layer is localized at the surface where the fluid velocity is much smaller than the channel-averaged velocity, the effective transport rate here is accordingly lower, allowing for a significant concentration to remain close to the surface.
We have numerically calculated the fluoride concentration profiles and fit the experimental v-SFG results using experimental and literature values for the model parameters (Fig. 6). We used kdis=0.027μmol m^−2 ^s^−1 ^^30,31^, ρF,sat=10ρF,b=10μM^32^, kads/kdes=10μM^−1^, and Γ = 10 nm^−2^. The equilibrium constant kads/kdes and the density of chargeable groups Γ are not directly reported in the literature. However, the combination chosen here is not only in agreement with the experimental equilibrium surface potential of about 70 mV (if the electric double layer is taken into account)^32–34^, but is also in agreement with the observed magnitude of typical flow-induced changes of the zeta potential (≈15 mV) as reported in streaming potential measurements^15^. The dissolution rate constant kdis is low and roughly corresponds to one monolayer of the surface dissolving every hour, as expected for a poorly soluble mineral^30,31^. Figure 6a shows the comparison between the experimentally observed changes in v-SFG intensity as a function of flow rate and position and those found in our numerical calculations (in which we tune the pressure drop to fit the experimental flow rate). The overall increase in the v-SFG response with increasing flow rate and its decrease from the inlet to the outlet is very similar to the characteristics of the surface charge (Fig. 6b). The surface charge increases by 50–100% compared to no-flow conditions, and varies by 10–20% laterally. The corresponding change in the surface potential ψ (as calculated by the Gouy–Chapman relation ψ=2kBTesinh−1(2πσλBλD)) is significantly lower in this nonlinear screening regime where ψ exceeds the thermal voltage of 25 mV, differing by only around 20% between the flow-on and flow-off states. Here λB = 0.7 nm is the Bjerrum length of water at room temperature and λD=7 nm the Debye length for a solution of 1 mM of NaCl and pH 3. Finally, to estimate the change in the SFG intensity, we use the approximation (Ion−Ioff)/Ioff=(ψon/ψoff)2−1, as discussed in Supplementary Note 1.
Fig. 6 Results from the numerical model.a Calculated (lines) and experimentally obtained (symbols) relative changes in the v-SFG intensity along the mineral surface for different flow rates. We show only the averaged values from Fig. 3 without the spread of the individual data points for clarity. b Calculated surface charge along the mineral surface for different flow rates, showing enhancements under flow-on conditions by ~50% compared to flow-off conditions. c Calculated fluoride concentration at the surface ρF(x,R) relative to the bulk along the mineral surface for several flow rates. d The calculated net dissolution rate J(x) along the mineral surface for different flow rates. Note the orders of magnitude increase upon turning on the flow. Under flow-off conditions, the dissolution rate is a small fraction of the maximum dissolution rate ∣J∣≈10−5kdis while during flow the net dissolution increases by several orders of magnitude, to approximately half the maximal dissolution rate ∣J∣≈kdis. The position is given as displacement from the center with values increasing in the direction of flow, indicated by the arrow in panel (c).
The validity of the Gouy–Chapman relation is verified in the Supplementary Discussion 1, where we show that the curvature of the channel wall on electrostatics can be neglected, which is expected as the Debye length is much shorter than the channel radius, λD/R≪1. This also justifies approximating the dissolving surface, which is a flat plate in the experiments, by a cylindrical channel wall in our numerical calculations. Similarly, the diffuse boundary layer and the electric double layer do not affect each other as their thicknesses differ by orders of magnitude, λD/δ≈10−4.
The calculated interfacial fluoride concentration (Fig. 6c) decreases with increasing flow rate and increases from inlet to outlet. The reverse relation between fluoride concentration and surface charge stems from the Langmuir isotherm (Eq. (4)). Interestingly, although the fluoride concentration increases in the flow direction, Fig. 6c shows that the flow-induced dilution ensures that the concentration is lower throughout the channel than in the no-flow state. The concentration at the boundary of the dissolving area, at x=−H, is equal to the bulk concentration ρF,b in flow-on conditions, beyond which it rises rapidly and then slowly reaches a maximum at x=H. This rapid increase of the fluoride concentration at x=−H is due to the precipitation being less at the inlet compared to upstream locations near x=H.
In the no-flow state, the fluoride concentration in the center corresponds to the saturation concentration ρF,sat. Without a precipitation reaction, the concentration would even be orders of magnitude larger. In our calculations that describe the experiment, the flow-induced change in fluoride concentration is on the order of ~30% of the saturation concentration (ρF,sat=10 µM), as is evident from Fig. 6c. In an absolute sense, this change of concentration is orders of magnitude lower than the background ionic strength (~2 mM). While the change in F^−^ concentration is sufficient to affect the surface charge, it is negligible for the total ionic strength that governs the screening of the surface charge. Thus, we can conclude from our combined experiments and theory that the flow-induced changes in the v-SFG signal are connected to the surface charge itself rather than changes in charge screening.
Finally, we draw attention to Fig. 6d, which reveals that the net dissolution rate also exhibits a gradient along the mineral surface and increases with flow rate. Similar to the position- and flow-rate-dependent shift in surface charge, the position-dependent net dissolution rate is also a result of the heterogeneous concentration profile. More importantly, however, we can see that the dissolution rate increases by orders of magnitude from flow-off to flow-on conditions. In the flow-off state, the surface concentration is nearly equal to the saturation concentration ρF(x,R)≈ρF,sat, which causes the net dissolution rate in flow-off conditions to be ~10^5^ times smaller at the center than the maximum dissolution rate 2kdis. During flow the dissolution is almost as large as kdis. The net dissolution rate thus increases by ~5 orders of magnitude due to flows of ~mL·min^−1^ (or shear rates 1.5–9 s−1), from effectively 10−7 to 10−2μM·m^−2^·s^−1^.
To investigate the universal scaling behavior of the fluoride concentration with channel geometry, dissolution, and flow rate, we derive in Supplementary Discussion 2 an analytic expression for the radially averaged fluoride concentration profile ρ¯F(x) under the assumption that the concentration profile is independent of the radial coordinate (such that ρF(x,R)=ρ¯F(x)) and upon neglecting precipitation. Key parameters of this expression are the Peclet number Pe=2u¯L/D, where u¯ is the channel-averaged fluid velocity, and the dimensionless concentration difference ΔρF=kdisL2/DRρF,b between the channel center and the reservoir under flow-off conditions. Here Pe captures the importance of advective transport relative to diffusive transport, and ΔρF is a measure for the dissolution rate. In Fig. 7a, we plot our analytic expression of ρ¯F(x)/ρF,b for the geologically relevant case ρF,b=10−6 mol·L^−1^ and ΔρF = 3.6 for several values of Pe, yielding remarkable agreement with our numerically exact solution of Eqs. (5)–(7) for the long-channel limit. Interestingly, the parabola-like low-Pe regime and the linear-like high-Pe regime of the graphs in Fig. 7a are borne out by the limiting cases of our analytic expression, which are given by
Fig. 7 Numerical and analytic results in the long-channel limit.a Numerical (dashed) and analytic (solid) results for the radially averaged fluoride concentration profile in the channel as a function of the distance from the center for Peclet numbers ranging from Pe=0.31 (average flow velocity 3.2×10−8 m·s^−1^, green) to 3100 (average flow velocity 3.2×10−4 m·s^−1^, pink) in steps of a factor 10. b Radially averaged concentration at the center of the channel as a function of the Peclet number, comparing analytic (solid lines) and numerical (symbols) results for geologically realistic values (kdis=10−11 mol·m^−2^·s^−1^, L=10−2 m, 10< average flow velocity <10^−4^ m·s^−1^,D = 10^−9^ m^2^·s^−1^, η = 10^−3^ Pa·s). Note that the considered geometry has length L=1cm and radius R=0.1mm ensuring that the radius is much smaller than the boundary layer thickness, which is not the case for the experimental geometry.
This expression makes explicit that all excess fluoride gets washed out of the channel for Pe≫1, the more so at the inlet (x=−L) than at the outlet (x=L). In Fig. 7b, we plot the fluoride concentration at the center of the channel (x=0) in the full range of Pe, which reveals a diffusion-dominated (Pe≪1) and an advection-dominated (Pe≫1) regime with Pe-independent limiting fluoride concentrations; a significant dependence on Pe is only found in a rather narrow crossover regime at Pe≈1. This narrow Pe regime could explain why previous authors did not observe a change in surface charge when the flow rate was varied^4,7,9,26^, as there is little flow dependency in the high Peclet regime (Pe≫1). Note that if the radial concentration profile is not constant, but rather described by a diffuse boundary layer with thickness δ≪R (as is the case for our experiments), the transition between flow and no-flow states is expected to broaden significantly and to occur at higher flow rates as the Peclet number has to be replaced by the Sherwood number (Sh∝Pe2/3), which we discuss in Supplementary Discussion 3.
How an almost insoluble mineral such as CaF2 can give rise to a significant concentration difference (ΔρF>1) in our experiment can be understood by interpreting ΔρF as a ratio of the timescale τdif=L2/D≈105 s for a fluoride ion to diffuse from center to outlet and the dissolution timescale τdis=ρF,bR/kdis≈10s by which dissolution can significantly change the bulk concentration. While the typical time for the concentration to change by dissolution is slow, the timescale by which diffusion can smooth out this concentration change is even slower in the present case of a long macroscopic channel. For this reason, the slow dissolution is offset by the even longer diffusion time, leading to a significant concentration difference (ΔρF=τdif/τdis>1). Note that the dissolution time is determined by the surface-to-volume ratio of the channel and is hence linear in the radius R, while the diffusion time is quadratic in the channel length L. Thus the dimensionless concentration difference ΔρF, which needs to be of order unity or larger for a flow-dependent surface charge and dissolution rate, is strongly dependent on both the absolute length L as well as the aspect ratio L/R of the channel. In fact, in the limit of large channels, ΔρF inevitably becomes so large that the assumption of a negligible precipitation reaction is invalidated. Moreover, the inverse scaling of ΔρF with the bulk fluoride concentration ρF,b explains the experimental finding of ref. ^4^ that there is no flow-dependent surface charge when an excess of fluoride is added. We experimentally reproduced this behavior for a 4 mM fluoride bulk concentration (Supplementary Fig. 3a) and confirm that under these conditions no surface charge gradient was observed (Supplementary Fig. 3b).
Our theoretical model and v-SFG experiments can be transferred to other interfaces to determine the charging reaction mechanism and dissolution reaction mechanism. In the present case of calcium fluoride, the surface charging reaction (Eq. (1)) is desorptive, which implies that a decrease in the fluoride concentration increases the surface charge. Meanwhile, the dissolution reaction (Eq. (2)) is a source of reactive ions and hence increases the fluoride concentration. As flow brings the concentration closer to the reservoir’s bulk concentration, the surface charge increases, more so at the inlet than at the outlet. The sign of the charges (i.e., whether an anion or cation is released) does not play a role in the surface charge’s flow dependency. Interestingly, the combination of the surface being charged by desorption and dissolution being a source of ions partaking in this charging reaction is only one of several ways to establish a flow-dependent surface charge. For instance, a precipitating surface in contact with an oversaturated solution is a sink for the reactive ions. In this case, the concentration would increase upon flow, which for a desorptive charging mechanism would lead to a decrease, rather than an increase, in surface charge closer to the inlet. Another charging mechanism is the adsorption of a reactive ion, where the surface charge increases with increasing concentration of reactive ions. This causes the surface charge to increase upon flow when the surface is a sink and to decrease upon flow when the surface is a source of reactive ions. Note that whether the surface is a sink or a source of ions is not necessarily determined by whether the mineral is dissolving or precipitating. For instance, protons play a role in the charging of iron oxide surfaces and are consumed when iron oxide dissolves^35^. In this case, the surface is a sink below the saturation concentration and a source above saturation concentration. As protons are often consumed in dissolution and partake in the charging of metal oxides, we expect such counter-intuitive behavior to be relatively common.
We summarize the expected change in surface charge upon flow for different combinations of charging and dissolution reactions in Fig. 8. If either the dissolution or charging reaction is known, the other reaction mechanism can be established by inspecting the sign of the v-SFG change upon flow or the slope of the surface charge gradient. This could be useful for the investigation of other minerals, for instance, the notoriously complex surface chemistry of silica^11,36–40^.
Fig. 8 Overview of charging and dissolution coupling.Indication whether the surface charge increases (+) or decreases (−) upon flow. Column position indicates if the surface is a source or sink of reactive ions, row position indicates if the reactive ion desorbs or adsorbs to generate surface charge. The change of surface charge upon flow is always larger at the inlet than at the outlet.
Next we discuss in what natural setting flow-dependent charge and dissolution could be relevant, considering a common geological system. While we have focused on the position-dependent surface charge, as this is a quantity indirectly measured by v-SFG, the dissolution flux of Eq. (7) also depends on the fluoride concentration. Therefore, our model directly and naturally implies that also the dissolution flux is flow- and position-dependent. This insight brings us to discuss the more general case of mineral dissolution in geological environments, particularly whether and where such a dissolution rate gradient can be expected. As already emphasized, the gradients of interfacial concentrations can only be significant if there is a sufficiently large relative concentration difference (ΔρF≫1) between the reservoir and the channel (or porous material) without flow. At first sight, this requires high dissolution rates. However, also the geometry and the channel dimensions play a role, and we have seen that macroscopic and elongated channels are favorable for a large ΔρF. These conditions are met naturally in porous networks such as rocks and soils. In addition to the condition of a significant concentration difference, also a significant fluid velocity, or Peclet number, is required. Using our analytic results (Eq. (8)), we estimate that in pores of centimeter lengths (L=10−2 m) and sub-millimeter widths (R=10−4 m) and for fluid velocities of u¯= 10^−5^ m·s^−1^ and dissolution rate constants that exceed kdis=10−11 mol·m^−2^ s^−1^) , there will be a significant surface charge gradient. We note that these dissolution rates are typical for silicate minerals^28,29^, and these flow rates are typical for water flow through soils^41,42^. We compare numerical and analytic results for these parameter values at different flow rates in Fig. 7b. The concentration profile changes dramatically between flow-off state and a flow-on state with a flow velocity of u¯=10−5 m·s^−1^, corresponding to Pe≈102. In the case that the concentration in the pore becomes close to the saturation concentration (Δρ≈ρsat/ρF,b), there will also be a flow-rate and position-dependent dissolution rate. Such a flow-dependent dissolution rate has been suggested to explain commonly observed dependencies of dissolution rate with pore size^29^. Therefore, our experimental findings and numerical models can be expected to be of geological relevance. In addition, our approach of combining the model and position-resolved v-SFG experiments itself could be used to investigate charging and dissolution mechanisms, which are important for soil remediation^9^ and industries such as froth flotation^6^.
In conclusion, we demonstrated that the surface charge and dissolution rate of calcium fluoride show macroscopic gradients along the interface when flow is applied. This was shown by the combined use of surface-specific v-SFG spectroscopy and full numerical calculations of reaction–diffusion–advection equations, including the interfacial chemistry at water–mineral interfaces. We extended v-SFG spectroscopy spatially to investigate the surface chemistry of minerals upon flow on a macroscopic scale. The observed gradient in the v-SFG response can be correlated to a gradient in the surface charge upon flow. A dissolution–diffusion–advection process can entirely explain such a gradient. The key physicochemical mechanism is captured by the coupling of a dissolution reaction with flow and surface charge. This dissolution reaction creates a steady-state concentration profile of fluoride ions, which depends on both position and flow rate. Both the surface charge and dissolution rate are influenced by the local fluoride concentration. Thus, dissolution also causes a position-dependent dissolution rate and surface charge during flow. Our reaction–diffusion–advection model can be generalized to arbitrary surface and dissolution reactions, and we expect that the method of position-resolved v-SFG spectroscopy can be transferred to other systems in order to determine charging and dissolution mechanisms. Interfaces exhibiting both dissolution and charging in flowing water occur naturally in a wide variety of systems, and only minor confinement is needed to ensure that both surface charge and dissolution rate depend on position and flow rate. Thus, our finding of a flow-induced gradient in the dissolution reaction is expected to impact geological research of porous materials. Also, scientific disciplines that rely on the use of surface potentials, such as electrokinetics, microfluidics, and nanofluidics, or those involving measurements of surface potentials may be interested in its dependency on the position and flow rate.
The SFG spectra presented in this work are recorded with an experimental setup based on a Ti:sapphire amplifier (Solstice® Ace^TM^, Spectra Physics) that generates 800 nm pulses with 1 kHz repetition rate and ~40 fs duration. For producing broadband infrared pulses (full width at half maximum, FWHM ~400 cm^−1^) with 4 μJ pulse energy, a commercial optical parametric amplifier (TOPAS Prime, Spectra Physics) is used together with a non-collinear difference frequency generation (NDFG) scheme. By guiding the visible upconverting pulse through a Fabry–Perot etalon (SLS Optics Ltd), it is spectrally narrowed to an FWHM of ~20 cm^−1^ with ~15 μJ pulse energy. The spectral resolution of the SFG signal is obtained by using a spectrograph (Acton Spectro Pro^®^ SP-2300, Princeton Instruments). An electron multiplied charge-coupled device (emCCD) camera (ProEM 1600, Roper Scientific) is used for detection. All spectra are recorded in ssp polarization combination (s-polarized SFG, s-polarized visible, and p-polarized IR) and at incident angles θvis≈73∘ and θIR≈76∘ of the visible and infrared pulse, respectively.
Sodium chloride was purchased from Carl Roth (>99.5%, CAS 7647-14-5), sodium fluoride from Merck (99.99%, CAS 7681-49-4), and concentrated hydrochloric acid from VWR (37 w% CAS 7647-01-0). All reagents were used as received without further purification.
The rectangular fluorite prism (dimensions 4 cm × 1 cm × 0.5 cm) purchased from EKSMA optics was baked at 500 °C for at least 2 h to remove organic residues. Before and after baking, it was rinsed with demineralized water, filtered with a Millipore unit (resistivity = 18 MΩ cm). The salt was dissolved in a pH 3 solution of hydrochloric acid prepared from the concentrated acid. If not stated otherwise, a 1 mM NaCl solution at pH 3 was used.
A scheme of the flow set-up and a technical drawings of the flow cell are shown in Supplementary Fig. 4. The design of the flow cell was optimized to ensure laminar flow. The prism is mounted on top of the flow channel. A 0.03-mm Teflon film is used for sealing. The flow rate was between 1 and 6 mL/min, corresponding to Reynold number of the order of ~5−25 and shear rate of ~1.5–9 s^−1^, which ensures laminar flow over the investigated flow cell. A reservoir with ~100 mL of the solution is connected with the flow cell via Tygon® tubes (ISMCSC0832 MHSL 2001; inner diameter 4.8 mm, outer 8.0 mm). The peristaltic pump (Masterflex 77924-70 L/S Ethernet/IP Network-Compatible Pump) is controlled by a self-written LabView program. The same program controls the motors (Newport Spectra-Physics GmbH CONEX-TRA25CC) of the x- and y-stages (Manual Linear Stage, Crossed-Roller Bearings, 25.4 mm Travel, M6).
One stage is used to bring the center of the channel perpendicular to the flow direction into the focus of the incident laser beams. To do so, a prism with 100-nm gold coating in the area of the flow channel is mounted on the flow cell, and the stage is moved to the edges where the gold signal disappears. The center is chosen in the middle of the edges. The center of the flow cell in flow direction was marked with an engraving on the flow cell’s surface. The second stage orthogonal to the first one is used to adjust the laser focus to the center of the flow cell in flow direction by eye with the help of the engraving. Afterward, the flow cell and the tubes are rinsed with MilliQ water. Before spectra are collected, the flow cell with the freshly cleaned prism mounted is flushed 2 min at a high flow rate of ~600 mL/min with the solution, followed by 10 min of rest to ensure equilibration of the interface. The second stage is used to change the position along the flow direction. At each position, spectra were recorded with an acquisition time of 5 s.
The spectra were subtracted by a background obtained by measuring a spectrum with a blocked IR beam for the same time as for the usual spectrum. Cosmic rays were identified by a self-written MatLab script once a tolerance for the difference between the intensities of neighboring pixels was exceeded. The cosmic ray was replaced by the average of the two previous points. Usually, a spectrum (1600 data points) contains <1% of cosmic rays. Integration was done by summation of the background-subtracted intensities in the range of the IR pulse. By inserting a polystyrene foil in the IR beam path prior to a gold sample, a spectrum with dips due to the infrared absorption of foil is obtained. The absorption dips were used for calibrating the wavenumbers. To calculate relative increases of the signal upon flow, spectra over ~5 min prior to starting the flow and spectra over ~5 min prior to turning the flow off are averaged. Several flow on/off cycles as well as single flow on/off cycles were taken into account for the presented data (see Supplementary Fig. 1). From Supplementary Fig. 1, it is also clear that the spectra chosen for averaging represent a steady-state.
We model a system involving fluid flow, diffusion, and two chemical reactions. The reactions are a desorption reaction responsible for the surface charging, and a dissolution reaction, responsible for the flux of reactive ions in the channel. Mathematically, the reaction–diffusion–advection problem is described by Eqs. (3)–(7). We use the approximation of local charge neutrality (ρF(x,r)=2ρCa(x,r) and ρNa(x,r)+ρH(x,r)=ρCl(x,r)=2 mM (a spatial constant)) to reduce the computational complexity. Note that NaF rather than CaF2 was used in the experiments to fix the bulk fluoride concentration, which only introduces a minor discrepancy in the Ca^2+^ concentration. The dissolution of CaF2 is taken into account by the boundary condition for the dissolution equation,
We ignored electrostatics, the Poisson equation, and a conduction term in Eq. (4). All this would require a nm-scale resolution that is computationally intractable on the cm-scale of the channel. Small scale calculations, however, in which electrostatics is fully incorporated through Poisson–Nernst–Planck–Stokes equations, can be found in the Supplementary Discussion 1; they reveal essentially the same key physics as the charge-free case. Nevertheless, we can convert the areal density σ(x) of CaF^+^ groups on the surface to a surface charge eσ(x) and to a surface potential ψ(x) if we make use of the Gouy–Chapman relation (which is also validated in the Supplementary Discussion 1),
Here e is the proton charge, λB=e2(4πϵ0ϵrkBT)−1 is the Bjerrum length, λD=(8πλBρs)−1/2 the Debye length where the total ionic strength ρs=2 mmol/L is determined by the background NaCl and HCl salt concentration in the experiments, e is the electron charge, ϵ0 the electric permittivity of vacuum, ϵr the relative permittivity of water, kB the Boltzmann constant, and T the temperature. The Debye length λD is spatially constant as the concentration of reactive ions is too low to change the total salt concentration.
The differential equations are solved numerically in the commercial finite-element software package COMSOL MULTIPHYSICS. The modeling geometry is described using cylindric coordinates (x,r,θ), as can be seen in Fig. 6 and is assumed to be azimuthally symmetric. The origin is located at the center of the channel. The channel length runs from x=−L to x=L, with L=20 mm, and the reactive surface runs from x=−H to x=H, with H=12.4 mm. The channel is connected to two bulk reservoirs of length Lb=42mm. The radius of the channel is R= 2.4 mm. The diffusion constant D=10^−9^ m^2^·s^−1^ and viscosity η=10^−3^ Pa s were used. A sketch of the model geometry can be seen in Fig. 5a. As boundary conditions, we impose bulk concentrations ρi(−Lb)=ρi(Lb)=ρb,i at the reservoir in- and outlet with ρb,F=1μM, a pressure difference P(Lb)−P(−Lb)=ΔP over the channel, no-slip boundary conditions at the channel wall, u(R)= 0, and no-flux boundary condition at the edge of the channel without a reactive surface n⋅J=0, where n is the surface normal. Dissolution is responsible for a net flux of ions out of the channel wall into the channel in the region −H<x<H as given by Eq. (7).
The steady-state solution for the surface charge is given by the Langmuir isotherm σ(x)=Γ(1+kadsρF(R,x)/kdes)−1 and the solution for the lateral component of the velocity profile is given by Poiseuille flow ux(r)/umax=1−(r/R)2, with umax=R2ΔP/8ηL.
P.O. gratefully acknowledges financial support by the Max Planck Graduate Center with the Johannes Gutenberg University Mainz (MPGC). All authors acknowledge the technical support of Marc-Jan van Zadel particularly the design of the flow cell. All authors thank anonymous Reviewer 2 for pointing us to relevant literature on diffuse boundary layer theory. This work is part of the D-ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).
P.O., E.H.G.B., and M.B. designed the experimental part of the research project and provided a qualitative interpretation of the experimental results. P.O. performed the experiments and analyzed the data. W.Q.B., M.D., and R.v.R. provided a quantitative interpretation of the experimental results. W.Q.B. performed the numerical and analytic calculations and generalized the findings. All authors discussed the results and wrote the manuscript.
Open Access funding enabled and organized by Projekt DEAL.
The data that support the findings of this study are available from the corresponding authors upon reasonable request.
The authors declare no competing interests.
René van Roij, Email: r.vanroij@uu.nl.
Mischa Bonn, Email: bonn@mpip-mainz.mpg.de.
The online version contains supplementary material available at 10.1038/s41467-021-24270-x.
The data that support the findings of this study are available from the corresponding authors upon reasonable request.