Coulombic surface-ion interactions induce non-linear and chemistry-specific charging kinetics

While important for many industrial applications, chemical reactions responsible for charging of solids in water are often poorly understood. We theoretically investigate the charging kinetics of solid-liquid interfaces, and find that the time-dependent equilibration of surface charge contains key information not only on the reaction mechanism, but also on the valency of the reacting ions. We construct a non-linear differential equation describing surface charging by combining chemical Langmuir kinetics and electrostatic Poisson-Boltzmann theory. Our results reveal a clear distinction between late-time (near-equilibrium) and short-time (far-from-equilibrium) relaxation rates, the ratio of which contains information on the charge valency and ad- or desorption mechanism of the charging process. Similarly, we find that single-ion reactions can be distinguished from two-ion reactions as the latter show an inflection point during equilibration. Interestingly, such inflection points are characteristic of autocatalytic reactions, and we conclude that the Coulombic ion-surface interaction is an autocatalytic feedback mechanism.

Charged solid-liquid interfaces play a central role in a wide variety of industries such as food and coating production [1][2][3], mining [4][5][6], medicine [7][8][9], soil remediation [10][11][12] and even carbon capture [13]. With the advent of nanoscale fluidics one expects that charged surfaces become ever more important [14,15]. In water and other polar solvents chemical reactions are a common mechanism by which surfaces obtain their charge. For ionic solids the de-or adsorption of a dissolved ionic compound is often preferred over the sorption of its own counterion [16][17][18]; for covalent solids such as polymers and metal oxides the acidic nature of surface groups ensures that the surface (de)protonates in polar solvents and hence becomes charged [17][18][19][20][21]. However, for many processes of industrial and environmental importance relatively little is known about the surface chemistry [17][18][19] as the electrolytes in realistic applications contain a large variety of ions that can all undergo multiple reactions [17,18,22]. Due to experimental limitations the majority of studies investigating surface charging are performed at (quasi)-equilibrium conditions [17,18], with the notable exception of pressurejump experiments [23,24]. Only recently, however, it has been shown that the kinetics of chemical surface reactions can strongly couple to electrokinetic fluid flows, thereby affecting the physical surface properties on macroscopic scales [18,[25][26][27][28][29]. Furthermore, with the recent advent of fast and surface specific non-linear spectroscopy the dynamic measurement of surface charge has become feasible [29][30][31][32][33][34]. In this context it has been explicitly stated that there is an urgent need for theoretical models to describe such experiments [35]. Traditionally, sorption kinetics is typically described by (pseudo)-first-order reactions [12,36,37] that exhibit singleexponential relaxation towards equilibrium; the influence of a time-dependent surface charge is usually neglected entirely [35,38,39]. We are aware of one theoretical work [40] and associated review [41] that considers a surface charge that affects the rate constants of ion-association, which, however, does not consider the (chemistry-specific) non-linear dynamics induced by the electrostatic feedback as we do here.
In this Letter we present a theory for the charging dynamics of solid surfaces. We include the Coulombic ion-surface interactions and reveal an intricate dependence on the reaction mechanism and the valency of the reactive ions already present in a mean-field description. The Coulomb interactions not only affect the time constant of the late-time exponential decay of the surface charge towards equilibrium after an ion concentration (or pH) shock, but they also induce strongly nonlinear dynamics at early times far from equilibrium. Combined with the present-day capability to experimentally measure the time-dependent surface charge density, our theory forms a first step to unveil the surface chemistry of technologically important but ill-understood materials [18,35], such as silica [22,42] and graphene [43], and of processes such as the clean-up of radioactive and heavy metals [10,16,44,45].
Surfaces, for instance silica, in water commonly charge either by desorption of ionic species from neutral surface groups or by adsorption of ionic species onto neutral surfaces. While the exact charging mechanism of the silica-water interface is complex, there is support for charging by desorption of protons at high pH and adsorption of protons at low pH [22,42,46,47], where SiOH (s) is a neutral silanol group that is covalently bound to the (solid) glass and where SiO − (s) and SiOH 2 + (s) denote a silanol group with a proton desorbed or adsorbed in Eqs. (1a) and (1b), respectively. Here ρ denotes the proton density at the solid surface, and the dissociation and association rate k d and k a will be discussed below for the charging kinetics of a single desorptive and a single adsorptive reaction, not only for monovalent reactive ions as in Eqs. (1a) and (1b) but for general valency z. While adsorption isotherms of real materials can rarely be described by just a single charging reaction [16,47], we show in Supplemental Material I (SM I [48]) that charging by multiple reactions can actually be well-approximated by the single-reaction kinetics presented in this Letter for a wide range of experimental conditions. We consider a macroscopic surface with a density Γ of identical surface groups. A group can only be in either a neutral or a charged state. The charging is assumed to take place either by desorption (labeled by −) of a cation of charge ze, or by adsorption (labeled by +) of a cation of charge ze, with z ≥ 0 and e the proton charge. The surface densities of charged and neutral groups are denoted by σ ± > 0 and Γ − σ ± > 0 respectively, and the surface charge density is given by ±zeσ ± . Note that the charging dynamics is invariant under the sign of the reacting ions, and without loss of generality we can restrict attention to reactive cations of (strictly positive) valency z. Assuming the chargeable surface sites to be independent, we can describe the reaction kinetics in terms of the time-dependent surface density σ ± (t) > 0 which satisfies Langmuir kinetics described by [49,50] for desorptive charging reaction (1a), and for adsorptive charging reaction (1b). Here k d and k a are the rate constants of the dissociation and association of the reactive ion and ρ(σ ± ) is the volumetric concentration of reactive ions at the surface, which is defined at the position where the rate-limiting step for the reaction occurs [50,51]. We consider this surface to be impermeable to non-reacting ions and therefore do not account for any Stern layer other than the charged surface groups [27]. The equilibrium surface charge follows from ∂ t σ ± = 0 and is given by σ ±,eq = Γ(1 + (k a ρ eq /k d ) ∓1 ) −1 , which reduces to an explicit "Langmuir isotherm" in the case that the equilibrium concentration of the reactive ions ρ eq ≡ ρ(σ ±,eq ) is a constant independent of σ ±,eq . In general, however, this Langmuir isotherm is a self-consistency equation for σ ±,eq that requires an additional "closure" relation ρ(σ ± ) for an explicit equilibrium solution σ ±,eq . Without (Coulombic) interactions between surface and ions, the local concentration ρ(σ ± ) of reactive species in the vicinity of the surface would be equal to the bulk concentration ρ b of the reactive ions far from the surface (which is independent of σ ± and hence also independent from the reaction mechanism), such that Eqs. (2a)-(2b) would be linear differential equations whose solution can be written as s with the dimensionless charge s ± = σ ± /σ ±,eq such that s ±,eq = 1; here s ± (0) − 1 is the integration constant and denotes the relative deviation from equilibrium at the initial time t = 0. Note that the condition that 0 ≤ σ ± (t) ≤ Γ implies that 0 < s ± (0) < Γ/σ ±,eq , where the lower bound corresponds to an initially neutral surface whereas the upper bound can be as large as O(10 − 100), since typical equilibrium conditions have a charge occupancy of only a few percent of the total number of chargeable groups [52]. Thus from measurements of σ ± (t) at various concentrations of reactive (dissolved) species both k d and k a could in this non-interacting case be determined. However, as the charged surface attracts or repels reactive ions, Eqs. (2a) and (2b) are complicated by a nontrivial rela-  for desorptive reactions when σ −,eq Γ, in (b) s + (t) for adsorptive reactions when σ +,eq Γ, and in (c) s + (t) for adsorptive reactions when σ +,eq Γ. Insets denote semi-logarithmic representations of |s ± (t) − 1|. The case σ −,eq Γ (not shown) is trivial with singleexponential decay for all z.
tion ρ(σ ± ), which causes a charge-dependent decay rate and introduces deviations from purely single-exponential relaxation of σ ± (t). In fact, an explicit function ρ(σ ± ) is needed to investigate and solve the dynamics, which we will develop here. We consider the planar and homogeneous chargeable solid surface discussed above in contact with a bulk solvent with permittivity ε and temperature T with a three-component 1 : For convenience we assume trace amounts of reactive ions and therefore set ρ b ρ s , where ρ s is the bulk salt concentration. We also assume the electrolyte volume to be macroscopically large such that ρ b and ρ s do not change due to surface charging. Furthermore, we assume the charging timescale τ ± , which remains to be derived, to be the slowest timescale of the system. Given that the typical timescale for electric double layer (EDL) equilibration is around 10 −9 − 10 −6 s and that the (geometry and flow dependent) transport timescale for ions in stirred reactors can be as short as 10 −4 s [53], we find a large window τ ± 10 −4 s for reactions to be well-described by our (reaction-limited) theory [54]: for example phosphate desorption shows characteristic reaction timescales of hours [55] and adsorption of transition metals can occur on millisecond timescales [23,24]. The slow-reaction assumption allows us to describe the EDL within an equilibrium theory, for which we take the Gouy-Chapman solution of Poisson-Boltzmann (PB) theory for simplicity [20,56]. Although PB theory is based on a mean-field assumption for a system of point ions, it is known that for all but the highest salt concentrations this theory is quite accurate for 1 : 1 and even 1 : 2 aqueous electrolytes [57], and we expect a similar accuracy for 1 : 1 : z electrolytes in the limit ρ b ρ s of our interest. Within these assumptions the concentration of reactive ions at the surface is determined by a Boltzmann distribution ρ(σ ± ) = ρ b exp[−zφ (σ ± )], where k B T φ (σ ± )/e is the electric potential at the surface with a surface charge ±zeσ ± , with k B the Boltzmann constant. For desorptive charging the surface and ions have opposite charge and hence zφ (σ − ) < 0, while for adsorptive charging ions and surface have the same sign yielding zφ (σ + ) > 0. With this observation the Gouy-Chapman solution for a 1 : 1 electrolyte, which is relevant here as ρ s ρ b , gives φ (σ ± ) = ±2 sinh −1 (zσ ± /σ * ) [20,56], where σ * = (2πλ B λ D ) −1 with λ B = e 2 (4πεk B T ) −1 the Bjerrum length of the solvent and λ D = (8πλ B ρ s ) − 1 2 the Debye screening length. Substituting the Gouy-Chapman potential in the Boltzmann factor yields where the exponent is positive for desorptive charging and negative for adsorptive charging. Because Eq. (3) is reactionmechanism dependent, explicit information on the charging mechanism can be deduced from the reaction kinetics as described by combining Eq.
In order to investigate the influence of the Coulombic ionsurface interactions on the charging dynamics, we numerically solve σ − (t) from the kinetic Langmuir-Gouy-Chapman equations Eqs. (2a)-(3). The symbols in Fig.1(a) present the resulting relative deviations from equilibrium, s − (t) − 1, for a desorptive reaction in the experimentally common case of low equilibrium saturation σ −,eq Γ, both for s − (t = 0) = 2 and 0.01 corresponding to a surface with double the charge compared to equilibrium and an initially almost uncharged surface, respectively, for equilibrium surface potentials of 50 mV (|φ eq | = 2, circles) and 100 mV (|φ eq | = 4, crosses) and for valencies z = 0, 1, 2, 3 indicated by the different colours. Fig.1(a) shows that a desorptive surface that is overcharged (s − > 1) decays to equilibrium faster than one that is undercharged (s − < 1), the more so for larger valencies z. Interestingly, the sorption of uncharged species (z = 0, black symbols) reveals perfect symmetry between the two cases as expected for first order kinetics, which is also manifest in the semi-logarithmic representation of |s − (t) − 1| in the inset of Fig.1(a) that shows a data collapse and a single-exponential decay for z = 0. For z ≥ 1 the inset reveals a non-exponential time dependence with an initially slower decay for undercharged surfaces and an initially faster decay for overcharged surfaces, the difference becoming more pronounced for higher valencies. Fig.  1(b) shows the deviation s + (t) − 1 from numerical solutions of Eqs. (2b)-(3) for an adsorptive charging reaction and the same low equilibrium surface density σ +,eq Γ and the same surface potentials and valencies as in (a). Interestingly, for this reaction the relaxation of an initially undercharged surface to equilibrium is faster, rather than slower as we found for desorptive undercharged surfaces in (a). Hence the two mechanisms can be distinguished by inspecting a single timetrace of the surface charge. We do not plot the dynamics of a desorptive surface that is saturated in equilibrium σ −,eq Γ (in which case k d k a ρ eq ) as the equilibration (dissociation) rate for such a surface ∂ t s −k d (s − 1) is linear and equilibration occurs through trivial single-exponential decay. The lack of non-linearity for such a surface stems from the fact that the dissociation process is unaffected by the electrostatic surface-ion interaction. However, as can be seen in Fig.1(c) the dynamics of an adsorptively charged surface with a saturated charge density σ +,eq Γ is markedly non-linear. As was the case in Fig.1(b) we see that an undercharged surface equilibrates faster than a single-exponential. Clearly, these rather distinctive features of the time-dependent surface charge contain explicit information on not only the reaction mechanism but also the valency of reacting ions. Interestingly, such deviations from single-exponential decay have historically been observed in pressure-jump experiments [23,24]. In SM II [48] we show that these experiments are well described by our theory, alleviating the need of introducing multiple reactions to describe such experiments.
In Figs.1(a)-(b) the dimensionless time on the horizontal axes contains a factor (2z + 1), which as we will show now, is convenient as it leads to a data collapse in the asymptotic nonlinear-screening regime |φ eq | 1 where s ± (t) only depends on the valency, the reaction mechanism, and the initial charge state. To see why the near-equilibrium decay rate includes a factor (2z + 1) in Figs.1(a) and (b) but not (c) we simplify the Langmuir-Gouy-Chapman Eqs. (2a)-(2b) and (3) in the important and common case of large equilibrium surface potentials where zσ /σ * > 1, say beyond 50 mV where |φ eq | ≥ 2. In this limit Eqs. (2a)-(2b) can be rewritten as a single polynomial (Chini [58]) differential equation for which a closed form solution can be obtained by separation of variables only for an adsorptively charged surface with k d k a ρ eq in which case s 2z+1 Near-equilibrium, s ± 1, Eq. (4) simplifies to the linear differential equation ∂ t s ± −(s ± − 1)/τ ± with the near-equilibrium decay rate for desorptive and adsorptive charging given, respectively, by As announced, this timescale τ ± shows that electrostatic attraction can alter the linear, near-equilibrium, decay rate by a factor (2z+1) for z ≥ 1 compared to the neutral case (z = 0) in the experimentally common regime (k a ρ eq /k d ) ∓1 1 where σ ±,eq Γ. As for the majority of surfaces the equilibrium charge is much lower than saturation, σ eq Γ [52], we expect the correction by a factor (2z + 1) to be common and we note in passing that the only other work focusing on the influence of Coulombic ion-surface interactions on kinetics [40,41] does not mention this factor. As already observed in Fig.1, our simple Eq. (4) shows that far-from-equilibrium the dynamics becomes non-linear and importantly the ± sign of the reaction mechanism breaks the near-equilibrium symmetry of dynamics with regard to a charge excess or a charge deficit. As can be seen from Eq. (3), desorptively charged surfaces which are overcharged, s − > 1, will initially attract an excess of reactive ions to the surface, thereby having an increased reaction rate compared to uncharged equilibration. Hence the equilibration for large overcharging is faster than expected from uncharged Langmuir kinetics. Conversely, an initially undercharged surface, s − < 1, will have a shortage of reactive ions and thus equilibration will be slower. For  adsorptively charged surfaces equilibration will be non-linear regardless of σ eq /Γ and here undercharging leads to a shortage of reactive ions compared to equilibrium and hence faster equilibration, as can be seen in Fig.1(b). The rate changing during equilibration is reminiscent of autocatalytic reactions where the equilibration rate changes because a catalyst speeding up the reaction is produced simultaneously with a reaction product [59][60][61], and in SM III [48] we demonstrate the similarity between Eq. (4) and autocatalytic kinetics. A characteristic feature of such autocatalytic reactions is an increasing decay rate up to a maximum and a corresponding inflection point in the time-dependent decay. Interestingly, for so-called ion displacement reactions in which ions are involved in both the forward-and the back-reaction, inflection points are easily realisable as there are now two ions attracted or repelled from the charged surface such that the Coulombic feedback is strengthened: inflection points are hence a smoking gun that multiple ions are involved in a charging reaction. An example of a two-ion reaction where all reacting ions are repelled from the charged surface is the calcium charging of silica [42] of Eq. (6a), while an example of a reaction where all the reacting ions are attracted to the charged surface is the fluoride charging of the biomineral carbonato-apatite [62] of Eq. (6b), where X = Ca 10 (PO 4 ) 6 , k f is the forward (charging) reaction rate, k b is the backward (decharging) rate and ρ i is the concentration of ion species i at the charged surface. In SM III [48] we derive under the same Gouy-Chapman and large surface potential conditions of the main text that the charge equilibration for Eq. (6a) and Eq. (6b) are respectively described by with the resulting near-equilibrium decay constant τ = k f ρ i,eq + 3k b ρ i,eq for both reactions. Comparing Eqs. (7a)-(7b) to Eq. (4) we see that now the time evolution is given by the difference of two polynomials of (non-zero) unequal degree, ensuring that there is always a maximum in the decay rate and hence an inflection point. We find that reactions of the form (6b) have two physically realizable inflection points located at s = 2/3 (if σ eq Γ) and s = (2/5) 1/3 0.75 (if σ eq Γ): the reaction (6a) has only one accessible inflection point s = 4 1/3 1.6 (if σ eq Γ), while its second inflection point s = 4/3 is inaccessible for a saturated surface with σ eq Γ.
We plot the dynamics resulting from the Eqs. (7a)-(7b) for a variety of starting conditions s(0) in Fig.2 in the experimentally common limit σ eq Γ. In Fig.2(a) excellent agreement between the asymptotic Eq. (7a) and full numerical results can be seen, while in Fig.2(b) for large undercharging s(0) < 0.2 a discrepancy between Eq. (7b) and the full numeric solution is found. However in general Eqs. (7a)-(7b) predict the location of the inflection point accurately for a range of common surface potentials (k B T /e)|φ 0 | ∈ [50,100] mV. For ion-displacement reactions of the form (6a) and (6b) involving ions with higher valencies but monovalent surface charge the inflection point will lie even closer to equilibrium. Thus surfaces that are initially undercharged by only 25% or overcharged by 60% will generally equilibrate along sigmoidal curves, which is a distinguishing feature that cannot be observed for the single-ion reactions Eqs. (1a)-(1b). Finally we note that ion-displacement reactions offer a simple explanation for the recently observed sigmoidal equilibration of the surface charge at an aqueous silica interface [61] using only a single charging-reaction of the form Eq. (6a) rather than the proposed autocatalytic dissolution cycle involving multiple steps.
To summarize, in this Letter we present a model for the non-linear dynamics of reaction-limited surface charging, combining Langmuir dynamics with Poisson-Boltzmann theory. The model captures how the screened electrostatic surface-ion interaction affects the reaction rate near and far from equilibrium in terms of a non-linear differential equation, where the electrostatic interaction is described by only using the charge valency of the reactive ion. The Coulombic ion-surface interaction leads to a charge-dependent decay rate, which can be used to gain information on the valency of reacting ions, initial charge, and reaction rate. Deand adsorptive reactions can be distinguished by inspecting whether far-from-equilibrium decay is slower or faster than near-equilibrium decay, while an inflection point is characteristic for two-ion reactions. Interestingly, we note that inflection points are a characteristic feature of autocatalytic reactions and that the electrostatic ion-surface interaction can be seen as a catalytic feedback loop. Hence electrostatics offers a straightforward explanation for the recently measured autocatalytic charging of silica [61].
The single-reaction (one-pKa) charging reaction presented in this Letter is a simplified representation for many real liquid-solid interfaces, which often require multiple surface charging reactions for the reproduction of measured equilibrium Langmuir isotherms [1][2][3][4]. In this Supplemental Material we show that for a wide range of conditions a two-reaction model is well approximated by the single-reaction model presented in this Letter. For simplicity we will consider a tworeaction (two-pKa) system with two distinct surface-sites labeled (SiOH) 1 and (SiOH) 2 , charged by similar desorptive surface reactions but with different ad-and desorption rates and hence the sites have different equilibrium areal densities σ eq,1 = σ eq,2 if the site densities Γ i or equilibrium constants K i = k a,i ρ/k d,i are unequal. The Langmuir equation describing their charging kinetics is with the resulting areal charge density σ e = σ 1 + σ 2 coupling the two reactions. As both reactions have the same reactive ion (H + ) a concentration change ∆ρ at the surface would cause the equilibrium to shift for both reactions, and the equilibration kinetics resulting after this concentration shift would in principle need to be described by the two coupled nonlinear differential equations (S2a) and (S2b) as they are coupled by the Gouy-Chapman "closure", given by Eq. (3) in the Letter. However, here we show that if the two reactions have dissimilar equilibrium constants K i or site densities Γ i one of the two differential equations can almost always be neglected, as when the equilibrium constants differ so does the shift in equilibrium charge density ∆σ i = σ eq,i (ρ + ∆ρ) − σ eq,i (ρ) for a given concentration shock ∆ρ. When the difference in equilibrium charge densities of the two sites is very unequal (either ∆σ 1 ∆σ 2 or ∆σ 1 ∆σ 2 ) the change in surface charge can be described using a single reaction model (as ∂ t σ e ∂ t σ 1 or ∂ t σ e ∂ t σ 2 ). Solving for the coupled-Langmuir kinetics Eq. (2) in steady-state, ∂ t σ i = 0, we find for σ i,eq Γ i that the relative change in the equilibrium surface charge ∆σ 1 /∆σ 2 ∝ (Γ 1 K 1 )/(Γ 2 K 2 ) while for nearly saturated  Figure S1. Replotting of Fig.1(a) showing charge relaxation for ion valencies z = 0, 1, 2, 3 (respectively black, blue, red, green) calculated with curves not originating from Eq. (4) in the main Letter but instead plotted using Eq. (S5) with dimensionless static backgroundcharge x ∈ [0, 1]. The outermost curves with x = 0 revert exactly to those described by Eq. (4) in the main Letter, while curves with x = 1 lie significantly closer to the black curve showing single-exponential relaxation, with intermediate x lying in the shaded region. Note that the near-equilibrium decay rate has changed from (2z + 1)k a ρ eq to (z + 1 + x)(1 + x) 2z−1 k a ρ eq which for z = 3 and x = 1 yields a factor 32 difference in the near-equilibrium decay rate. In general however for small x the charging kinetics is well described by the single reaction of Eq. (4) in the main Letter.
surfaces with σ eq ≈ Γ the relative shift scales as ∆σ 1 /∆σ 2 ∝ (Γ 1 K 2 )/(Γ 2 K 1 ). As long as ∆σ 1 /∆σ 2 deviates significantly from unity, the shift in one of the two equilibrium-densities can be neglected for the shift in the total charge density. To find an explicit maximum bound for which this approximation is valid we calculate the maximum concentration shock ∆ρ m ∈ ρ[−1, ∞] for which this approximation holds by solving for ∆σ 1 /∆σ 2 = 1. This has a trivial solution ∆ρ = 0 when Γ 1 = Γ 2 and K 1 = K 2 and a non-trivial solution which simplifies to ∆ρ m /ρ = (K 1 K 2 ρ 2 ) −1 − 1 in the case that Γ 1 Γ 2 . In general, the range of validity of Eq. (S3) is large when K 1 and K 2 are very unequal, except close to a concentration where ρ 2 K 1 K 2 = 1 where ∆ρ m tends to zero. Interestingly, in this case the charge density is exactly half-occupied σ e = (Γ 1 + Γ 2 )/2, which is rare for most experimental conditions. For all other concentrations ρ a very large range ∆ρ m /ρ of concentration shocks remains over which a two-reaction system essentially equilibrates through a single charging reaction. However, while the dynamics will be governed by a single reaction, there will be a static background charge due to which σ 1 = σ e such that the single reaction is still not exactly equal to the single-reaction kinetics in the Letter. That the deviation from the single-reaction kinetics in the Letter due to this static background charge is minor will be shown in the next paragraph.
Interestingly, the near-equilibrium decay rate τ − is altered from (2z + 1)k a ρ eq to (2z + 1 + x)(1 + x) 2z−1 k a ρ eq and hence the deviation from uncharged Langmuir kinetics becomes even larger if the background surface charge is included. Clearly both the timescale and dynamics revert back to the single-reaction kinetics of Eq. (4) in the Letter when x 1. To check how much the dynamics is affected at intermediate x we replot Fig.1(a) in the Letter but now with Eq. (S5) with z = 0, 1, 2, 3 (black, blue, red, green) and s = 1 and s = −0.99 instead of Eq. (4). The resulting shaded regions in Fig.S1 represent curves with different x ∈ [0, 1]. We observe that for increasing x the dynamics moves closer to single-exponential decay (black line), with this shift being more pronounced for increasing z. However, even for x = 1 there is no over-dramatic difference from the single-reaction kinetics explored in the Letter and in many cases a tworeaction system is well approximated by the single-reaction model in the Letter. However one should be cautious when extracting the ion valency from dynamics around x 1 as here the dynamics closely resembles that of a single ion with z − 1. We expect this may occur in processes such as the adsorption of heavy metal ions from ground water which occurs on pre-charged substrates [3]. Under these conditions surface charging can be readily described using Eq. (S5).
To summarize, here we have shown that the one-pKa charging reaction is a valid approximation for more complex sys-tems, involving multiple charging reactions, when (i) the ratio of the equilibrium constant K 1 /K 2 is not close to unity and (ii) the concentration shock ∆ρ is constrained within a range ∆ρ m which we show to be generally large.

II. SURFACE CHARGING KINETICS FROM PRESSURE-JUMP EXPERIMENTS
While experimental investigations of surface charging kinetics are rare, several kinetic studies employing a pressurejump technique exist [5,6]. In such an experiment a mixture of colloidal particles and aqueous electrolytes is slowly pressurized to more than 1 MPa, thereby shifting the surface reactions at the colloidal surface to a high-pressure equilibrium. When this pressure is suddenly released, the solution pressure converges to atmospheric pressure in 0.1 ms. After this jump, the colloidal surface charge must revert from its high-pressure equilibrium to its atmospheric equilibrium. The change in surface charge is measured indirectly, by using the solution conductance as a proxy for the surface charge. While it is unclear what the exact relation between charge and conductance is, a linear relation is often assumed [6].
Of particular interest is a set of experiments where the adsorption of divalent transition metals such as Cu 2+ , Pb 2+ , Mn 2+ , and Co 2+ onto γ-alumina (Al 2 O 3 ) particles is studied. Here the authors reject simplest adsorption mechanism, for Cu 2+ given by because this single-step reaction does not show singleexponential decay [5]. Their theoretical model used for analysis fixes the surface (zeta) potential, and hence they find that multiple reactions are needed to explain the observed charge equilibration [5][6][7]. While their analysis is valid for fixed surface potential, for a potential varying with surface charge we actually expect deviations from single-exponential decay, as shown by Eq. (4) in our Letter. The Langmuir-Gouy-Chapman dynamics of reaction (S5) with z = 2 for surface potentials larger than 50 mV is given by To test whether the non-linear dynamics observed in the pressure-jump experiment are explained by Eq. (S6), we extract the experimental data from the relevant pressure-jump experiment for copper adsorption (Fig.1 in Ref. [5]). In this experiment it is found that the conductivity after the pressure jump decreases and reaches a constant value within 200 ms. In Fig.S2 we compare the experimental data (symbols) with equilibration expected from Eq. (S6) (black line), where we assume the experimentally common case where σ eq Γ and extract the reaction time τ + = (5k d ) −1 27 ms from the single-exponential, late time, relaxation (green line). As the final equilibrium charge density is not measured in a pressure-jump experiment we use the degree of undercharging s(0) − 1 −0.6 as a fit parameter, where we note that an undercharged surface naturally explains why conductivity decreases: during equilibration mobile copper ions are taken out of solution. It can be seen that the difference between the single exponential decay (green) and experimental data (symbols) is large, but that Eq. (S6) (black) can naturally explain a large part of the deviation from single-exponential decay, without needing to introduce a second reaction with a different timescale.
While further analysis is required to reinstate reaction mechanism (S5), our analysis shows the importance of taking Coulombic surface-ion interactions into account when considering kinetics. Not only do we show that non-linear decay is readily captured by our model, also the desorption time k −1 d is five times slower than the late-time decay τ + , in this divalent case. Furthermore, we have demonstrated that chemical information can easily be extracted from a single equilibration curve, even when the initial degree of undercharging is unknown.

III. SIMILARITY BETWEEN SURFACE CHARGING AND AUTOCATALYTIC KINETICS
Here we we will show the similarity between classical autocatalytic kinetics and surface charging as described by Eq. (4) in the Letter. First we derive an approximate solution for this Chini differential Eq. (4) in the main text, by expanding it up to second order around s = 1 obtaining a Bernoulli differential equation (Supplementary Ref. [8]). Solving this equation by standard methods we obtain s ± (t) − 1 = (s ± (0) − 1)e −t/τ ± 1 ∓ z 1 − e −t/τ ± s ± (0) − 1 ± s ± (0) − 1 (2z + 1)(k d /k a ρ eq ) ±1 + 1 , where τ ± is given by Eq. (5) in the Letter. Our Eq. (S7) reverts to single exponential decay when (s ± (0) − 1)z 1, and is valid only as long as ∓z(s ± (0) − 1) −1. For most practical purposes Eq. (S7) is not of much use, however it is interesting to note that a similar solution exists for autocatalytic equations thereby substantiating the claim that surface charging is autocatalytic. To make this comparison explicit we consider the simplest possible autocatalytic reaction (Eq. (1a') from Ref. [9]), where the aqueous reactant Y (aq) together with reactant A (aq) produces q copies of itself. The reaction is autocatalytic when the autocatalytic order q ≥ 2. When the concentration ρ A is constant the production rate of Y (aq) is given by the Chini differential equation with y = ρ Y /ρ Y,eq , which already shows similarities to Eq.
with τ y = (q − 1)k f ρ A . While Eq. (S10) is already very similar to Eq. (S7), the similarity becomes even more apparent when comparing the autocatalytic dynamics to the dynamics of an adsorptively charged surface with σ +,eq Γ (where k d k a ρ eq ) in which case Eq. (S7) simplifies to s + (t) − 1 = (s + (0) − 1)e −t/τ + 1 − z(s + (0) − 1)(1 − e −t/τ + ) . (S11) Comparing Eq. (S10) with Eq. (S11) we find the only difference is the definition of the timescale τ i and that the ionvalency z replaces the autocatalytic order q/2. This correspondence between the ion valency z and autocatalytic order q supports the interpretation that the Coulombic ion-surface interactions acts autocatalytically.

IV. CHARGING DYNAMICS OF ION DISPLACEMENT REACTIONS
Here we will generalize the derivation of Eq. (4) in the main text from single-ion reactions to two-ion reaction also known