All mineral surfaces in contact with aqueous solution acquire charged, but the charging mechanisms vary among minerals (Dzombak and Morel, 1990; Lyklema, 1991; Lyklema, 1995; Lützenkirchen, 2006; Karamalidis and Dzombak, 2010). In the case of oxides, the surface charge is due to H⁺/OH⁻ uptake or release by the surface exposed oxygen atoms or dissociative water sorption on the surface exposed metal sites (Yates et al., 1974; Hohl and Stumm, 1976; Schindler et al., 1976; Davis et al., 1978). The surface oxygen atoms show amphoteric behavior, being able to accept H⁺ and act as a base or donate H⁺ and thus act as an acid. One or two proton binding reactions typically capture this acid-base behavior of oxide surface (Yates et al., 1974; Davis et al., 1978; Van Riemsdijk et al., 1991). One of the most popular protonation schemes, 2-pK (Yates et al., 1974; Davis et al., 1978), assumes the presence of a charge-neutral surface site, denoted as ≡SOH⁰, where S stands for surface metal/metalloid site, which can either accept or donate protons according to the following reactions: ≡ S O H 0 + H + ⇄ ≡ S O H 2 + , K 1 = ≡ S O H 2 + ≡ S O H 0 H + (1) ≡ S O − + H + ⇄ ≡ S O H 0 , K 2 = ≡ S O H 0 ≡ S O − H + (2) where K₁,K₂ stand for two equilibrium constants (hence the name 2-pK, pK= − log₁₀K). The brackets {} in the above reactions represent the surface concentration.

In reality, the charge assigned to the surface site is rarely an integer. The fractional values reflect the actual bonding environment of the surface exposed atoms and the specificity of a given mineral surface. These partial charges can be assigned using Pauling bond valence rules (Hiemstra et al., 1989a; Hiemstra et al., 1989b; Hiemstra et al., 1996; Hiemstra and Van Riemsdijk, 1996) or the first principle quantum-chemical calculations (Zarzycki, 2007a; Zarzycki, 2007b).

The charged surfaces attract electrolyte ions, which accumulate near the surface and form a spatial charge distribution known as the electrical double layer. In many cases, the co- and counter-ion adsorb non-specifically and remain fully hydrated. In these cases, they form the outer-sphere surface complexes, which are assumed to be located in the layer parallel to the surface, β-plane (see Figure 1). The β-plane is separated from the surface by at least one layer of interfacial water. The ions complexation reactions and corresponding equilibrium constants are given by: ≡ S O − + C + ⇄ ≡ S O C + , K C = ≡ S O C + ≡ S O − C + (3) ≡ S O H 2 + + A − ⇄ ≡ S O H 2 A − , K A = ≡ S O H 2 A − ≡ S O H 2 + A − (4)

The concentration of protons at the surface and ions in the β-plane are not the same as in the bulk solution. The ion concentration in the EDL is affected by the presence of the electric field gradient, which is usually accounted by using the Boltzmann weight, that is (Yates et al., 1974; Davis et al., 1978): X z i = X z i bulk ⁡ exp − z i e ψ x k B T (5) where k _B is the Boltzmann constant, T is the temperature, i stands for ion type, z _i is its formal charge, and ψ(x) is the potential at the surface (x = 0, ψ ₀) in the case of proton sorption or potential in the β-plane (x = β, ψ _β ) in the case of electrolyte ions.

The electrostatic Boltzmann weights were introduced into SCM reactions to account for the variation in chemical potential caused by the variable surface charge (Yates et al., 1974; Hohl and Stumm, 1976; Schindler et al., 1976; Davis et al., 1978). Consequently, the surface reactions are expressed using the law of mass action subjected to the electric field (Brown et al., 1999). For example, eq. (Dzombak and Morel, 1990). can be rewritten as: K 1 = K 1 int ⁡ exp e ψ 0 k B T (6) where K 1 int is the intrinsic constant, which is independent of the protonation state of other surface sites. The separation of the equilibrium constant into an intrinsic (chemical) part and electrostatic interaction part was inspired in the earliest SCM studies (Yates et al., 1974; Hohl and Stumm, 1976; Schindler et al., 1976; Davis et al., 1978) by a similar treatment of the ionization constant of the polyprotic acids (King, 1965).

The activities should replace the bulk concentration. Here, the activity coefficients γ _i are calculated using the Davies equation (Davies, 1962): log γ i = − 1.825 × 1 0 6 z i 2 ϵ r T 3 / 2 I 1 + I + 0.3 I (7) where I is the ionic strength of the electrolyte solution and ϵ _r is the relative dielectric constant of the solution.

The surface complexation reactions, eqs. (Dzombak and Morel, 1990; Lyklema, 1991; Lützenkirchen, 2006; Karamalidis and Dzombak, 2010). have been frequently rephrased as multicomponent Langmuir-type isotherms to emphasize that surface charging and ions complexation are treated as localized (chemi)sorption of ions at the specified lattice sites in the presence of the mean-field type electric field (Dzombak and Morel, 1990; Lyklema, 1995). Consequently, the surface concentrations can be replaced by the fractional coverages, θ _i (Yates et al., 1974; Charmas et al., 1995; Rudzinski et al., 1999; Prélot et al., 2002; Zarzycki et al., 2004a; Zarzycki et al., 2004b; Zarzycki et al., 2005a; Zarzycki et al., 2005b; Piasecki, 2006; Zarzycki, 2006; Zarzycki, 2007a; Zarzycki, 2007c; Piasecki et al., 2007; Piasecki et al., 2010): θ i = i N s , a n d ∑ i θ i = 1 (8) where i =SO⁻, SOH, SOH 2 + , SOC⁺, SOH₂A⁻ and N _s is the surface site density in sites/nm².

Note that the assumption that protons are located, or chemisorbed, in the 0-plane and electrolyte ions are localized in the β-plane are consistent with the requirement that the charge distribution in the rigid part of EDLs is governed by the Laplace equation (d ² ψ/dx ² = 0, constant capacitor models).

By using θ-notation, the surface charge is defined as a sum of positive and negative charges at the surface: σ 0 = B s θ SOH 2 + + θ SOH 2 A − − θ SO − − θ SOC + (9) where B _s = eN _s . Similarly, the charge in the β-plane is defined as the average of the charges accumulated in β-plane (Yates et al., 1974; Davis et al., 1978): σ β = B s θ SOC + − θ SOH 2 A − (10) The charge neutrality condition allows us to define the diffuse charge as: σ 0 + σ β = − σ d (11) The surface complexation reactions eqs. (Dzombak and Morel, 1990; Lyklema, 1991; Lützenkirchen, 2006; Karamalidis and Dzombak, 2010), along with the mass balance eq. (Kosmulski, 2002) and charge neutrality eq. (Kosmulski, 2018) form a set of six nonlinear equations with seven unknown variables: ψ ₀, ψ _β and θ _i . In order to solve these equations, one needs to reduce the number of unknowns to become equal to the number of equations. The reduction is accomplished by assuming that i the EDL can be considered as the parallel-plate capacitors connected in series, and i i using the Gouy-Chapman theory that gives the charge-potential relationship for the diffuse part of the EDL.

The first assumption allows us to connect surface, β-plane, and diffuse potentials using the integral capacitances c ₁ and c ₂ and surface and diffuse layer charge densities (Figure 1): c 1 = σ 0 ψ β − ψ 0 a n d c 2 = − σ d ψ β − ψ d (12) Now, the surface and β-layer potentials can be rewritten as functions of charges and ψ _d (Yates et al., 1974): ψ 0 = ψ β + σ 0 c 1 = ψ d − σ d c 2 + σ 0 c 1 (13) The second assumption provides a link between the diffuse layer potential and diffuse charge density, and is known as the Grahame relationship (Lyklema, 1991; Lyklema, 1995): σ d = − f G C ⁡ sinh e ψ d 2 k B T (14) The proportionally coefficient f _GC is given by: f G C = 2 ϵ r ϵ 0 k B T e κ = 8 ϵ ϵ 0 k T I (15) where ϵ ₀ is the vacuum permittivity, and κ is a measure of the EDL thickness and is defined through the inverse of the Debye length (λ _d ) as (Lyklema, 1991; Lyklema, 1995): λ d = κ − 1 = ϵ r ϵ 0 k B T 2 I e 2 (16) Note that eq. (Yates et al., 1974) is back-propagated to the surface and β-plane through capacitances via eq. (Kosmulski, 2001).

The diffuse charge can also be defined as the sum of all charges accumulated between the slip plane and the bulk solution phase (Verwey and Overbeek, 1948; Lyklema, 1991; Lyklema, 1995): σ d = ∫ s ∞ ρ x d x w h e r e ρ x = − ϵ 0 ϵ r d 2 ψ x d x 2 (17) where ψ(x), ρ(x) are the electrostatic potential and charge density at a distance x from the surface, respectively. By integrating eq. (Davis et al., 1978) with assumptions that electric field E = −∇ψ and potential ψ disappear in the bulk electrolyte yields eq. (Yates et al., 1974).

Summing up, the SCM is a combination of the chemisorption models for the Helmholtz part of the EDL, and the physisorption model for the diffuse part. From the mathematical point of view, such division is consistent with applying the Laplace and Poisson differential equation for the charge distribution in the rigid and diffuse parts of the EDL, respectively. Note that the sorption model improves the Poisson-Boltzmann model of the EDL by imposing a limit on surface charge and β-layer charge because protonation occurs at specific coordination sites with known surface site density (Dzombak and Morel, 1990), see eqs. (Kosmulski, 2002; Kosmulski, 2011; Kosmulski, 2020).

The governing equations in the SCMs are not only nonlinear due to the exponential dependence of ions activities on the surface potential, eq. (Lyklema, 1995), but they are also convoluted due to the assumption of the parallel-plate capacitor, eq. (Kosmulski, 2001). The complexity of the equations is not immediately apparent. Here, we rewrote the governing equations to show the convoluted character explicitly.

Let’s start by presenting the complexation reactions, eqs. (Dzombak and Morel, 1990; Lyklema, 1991; Lützenkirchen, 2006; Karamalidis and Dzombak, 2010). in terms of the fractional coverages: θ SOH 2 + θ SOH 0 = ≡ S O H 2 + ≡ S O H 0 = K 1 H + exp − e ψ 0 k B T (18) θ SOH 0 θ SO − = ≡ S O H 0 ≡ S O − = K 2 H + exp − e ψ 0 k B T (19) θ SOC + θ SO − = ≡ S O C + ≡ S O − = K C a C ⁡ exp − e ψ β k B T (20) θ SOH 2 A − θ SOH 2 + = ≡ S O H 2 A − ≡ S O H 2 + = K A a A ⁡ exp e ψ β k B T (21) where a _C , a _A stand for the activity of cations and anions (e.g., a C = γ C [ C + ] bulk ).

Solving equations (Hiemstra et al., 1989a; Hiemstra et al., 1989b; Van Riemsdijk et al., 1991; Hiemstra et al., 1996), supplemented with conservation of the surface site density, eq. (Kosmulski, 2002), gives: θ SOH 2 + = K 1 K 2 H + 2 e − 2 e ψ 0 k B T Θ , θ SOH = K 2 H + e − e ψ 0 k B T Θ , θ SO − = 1 Θ (22) θ SOC + = K C a C e − e ψ β k B T Θ , θ SOH 2 A − = K 1 K 2 K A a A H + 2 e − 2 e ψ 0 k B T + e ψ β k B T Θ (23) where the denominator, Θ, is given by: Θ = 1 + K C a C e − e ψ β k B T + K 2 H + e − e ψ 0 k B T 1 + K 1 H + e − e ψ 0 k B T 1 + K A a A e e ψ β k B T (24)

Finally, we can rewrite the eq. (Kosmulski, 2001). into two master equations that connect surface and β-layer potentials with the surface charge and diffuse layer potential: ψ 0 ψ 0 , ψ d = ψ d + f G C c 2 sinh e ψ d 2 k B T + σ 0 ψ 0 , ψ d c 1 (25) and ψ β ψ d = ψ d + f G C c 2 sinh e ψ d 2 k B T (26) The convolution is given through the surface charge dependence on the surface and diffuse layer potentials σ 0 ψ 0 , ψ d = B s Q Θ (27) where Q = K 1 K 2 H + 2 e − 2 e ψ 0 k B T + K 1 K 2 K A a A H + 2 e − 2 e ψ 0 k B T + e ψ β k B T − K 2 H + e − e ψ 0 k B T − K C a C e − e ψ β k B T (28) Even in a simple 2-pK TLM model, the nonlinear and convoluted character of the governing equations makes the numerical solution challenging.

The performance of machine learning (ML) models relies on the quantity and representation quality of the data on which these models are trained. In the case of the mineral/electrolyte interface, the available experimental observations are limited, sparse, and often obtained in a non-consistent fashion (Kosmulski, 2001; Lützenkirchen, 2006; Kosmulski, 2009a). On the other hand, the existing analytical surface complexation models (SCMs) can be conveniently used to generate sufficient data to successfully train ML models. To allow ML-models to be applicable for a wide range of the oxide/electrolyte interfaces, the dataset needs to cover wide range of the possible values of the SCM model parameters and environmental conditions–in other words explore the SCMs parameter and conditions space. Figure 2 illustrates the workflow of constructing an AI surrogate for predicting the SCM model parameters values from the input (e.g., experimental datapoints).

Here, we used the 2-pK Triple Layer model of the oxide/electrolyte interface because it an archetypal and the most generic model among all SCMs. We used this analytical model to generate a large dataset of titration curves by randomly sampling the parameter space. We focused on the surface charge/electrokinetic potential as a function of pH, because they are the most commonly reported experimental observations (Figure 3A). We sampled the parameter space composed of the model parameters (ion affinities, capacitances) and environmental conditions (pH, electrolyte concentration, surface site density).

Figure 3 illustrates our approach of data extraction from synthetic dataset to create the AI-learnable dataset. Next, this dataset is normalized and split into training, validation, and test subsets. The dataset generated in this study and preprocesses AI-learnable inputs are available in our GitHub repositories (see Section Data Availability Statement).

The random forest (RF) algorithm consists of several decision tree regressors whose predictions are combined together in an averaging scheme to improve the overall predictive power of the weak learners–decision trees (Géron, 2019).

At the moment, the gradient boosting algorithms are one of the most successful classes of regression algorithms. Similar to the Random Forest, the gradient boosting relies on an ensemble of weak learners (decision trees), but here they are not independently trained on the subset of data with averaged predictions but build sequentially on top of each other to reduce progressively the errors of each weak learner (Géron, 2019).

The RF and XGBoost methods described above were trained using a set of vectors, each consisting of the seven titration points (pH,σ ₀ (pH), ζ) and the global system descriptors (Point of Zero Charge, surface site density, and electrolyte concentration)–see Figure 3B. We used the XGBoost and RF algorithms as implemented in the scikit-learn package (Pedregosa et al., 2011).

Figure 4 shows the model performance on the training set, using five model parameters as the output (target). Figure 5 shows the model performance on the test dataset. The overall accuracy of the both models is incredible with R² score around 0.993.

The high value of the R² score indicates that RF has successfully learned the underlying relationship between the input information and the corresponding SCM model parameters during the training. The correlation between RF prediction and the ground truth is about 0.95 on the test set (see Figure 5). Concluding, RF-based surrogate can predict accurately the SCM model parameters.

The comparison of the performance of XGBoost on training and test sets are shown in Figures 6, 7. Surprisingly, the XGBoost method was less accurate than RF, even though it is usually outperforming RF algorithms (Géron, 2019). Specifically, RF achieved higher prediction accuracy for all model parameters on training set. XGBoost was the least accurate in predicting logK ₁ parameter value (see Table 1). Although XGBoost-based SCM surrogate is inferior to RF-based surrogate in terms of the predictive power, it learn the mapping from input to output about 100 times faster than RF (Table 1).