Charge regulation and surface complexation modeling in nanoscale 2D geometries: benchmarking and test cases of a novel code (CRESCENDO)

Abstract

Mineral surfaces in contact with aqueous solutions develop an electric double layer (EDL) through surface (de-)protonation reactions and adsorption of ions, diffusion, and electrostatic forces, resulting in a Stern- and a diffuse layer of ions. Most current models used for surface speciation calculations do not consider changes in surface chemistry caused by charge regulation effects, i.e. effects of interacting EDLs of surfaces in close proximity. Charge regulation modeling requires equilibrium calculation of every involved surface simultaneously, while also solving the Poisson–Boltzmann equation (PBE) to quantify electrostatic interaction. Since analytical solutions of the PBE for complex geometries do not exist it becomes necessary to solve such problems numerically. A Python code is presented that combines a general chemical speciation code, Three Plane Surface Complexation Model, and a Finite Element solution of the PBE on two-dimensional domains. The Finite Element PBE solver is benchmarked against analytical solutions and the speciation code is benchmarked against a PHREEQC model as well as an existing 1D charge regulation code. A test case involving charge regulation in a corner of two perpendicular surfaces is modeled. Charge regulation modeling on a nanoscale enables simulations of the electrostatic environment and surface chemistry in nano-confined systems and interactions of nanoparticles. This may also improve simulations of environmental and biological systems, cementitious materials and modeling of the electrostatic environment and sorption on nanoporous clay materials. Such information can be vital for the in depth understanding of natural and engineered barrier systems of nuclear waste repositories or other environmental scenarios.

---

Introduction

Surface chemistry plays a critical role in a wide range of fields, particularly due to its influence on colloidal stability and sorption processes. Gaining deeper understanding and creating accurate predictive modeling capabilities for the effects of surface chemistry is vital for environmental and technological applications from remediation of polluted wastewater,¹ (trace) element mobility in natural and engineered systems,² and colloid stability³ to radioactive waste disposal in deep geologic repositories.⁴

Electrostatic surface complexation models (SCMs) are commonly used to model chemistry at the solid–liquid interface, determining the chemical speciation of the surface and the adjacent solution. Common SCM variants include the Constant Capacitance Model (CCM) and the Diffuse Layer Model (DLM), along with more detailed models, such as the Basic Stern Model (BSM) or the Triple Layer Model (TLM), Three Plane Model (TPM), and Four Layer Model (FLM), which subdivide the interface into distinct regions to account to some extent for finite ion size effects.⁵ The approaches that are most closely related to the Stern-Grahame model feature at least one layer of specifically bound ions at the surface (Stern Layer) and a diffuse layer of mobile ions screening the remaining charge, which can be modeled using the Poisson–Boltzmann equation (PBE) or simplified versions such as the Debye–Hückel or Grahame equations. For certain cases (geometries and electrolyte compositions) analytical solutions of the PBE exist.^6,7 A variety of SCM approaches can be implemented in common (open source or freeware) geochemical modeling software like PHREEQC,⁸ Visual MINTEQ,⁹ FITEQL¹⁰ and others. Additionally, various commercial software solutions with surface complexation modeling capabilities and some custom-made codes exist, but these are usually not freely available for use in research applications, unlike the aforementioned options.

One aspect of particular interest in nanoporous media is the behavior of solid/liquid interfaces in situations where electric double layers of two or more surfaces overlap.¹¹ Such overlap leads to complex electrostatic interactions that are challenging to model, like spillover of potential on mineral faces and edges (like on mica and clay platelets),¹² inhomogeneous surfaces,¹³ and confined geometries. In these cases, the surface chemistry deviates from the equilibrium obtained by standard calculations of a single surface in contact with solution as the electrostatic interactions need to be taken into account. Two classical boundary conditions for electrostatic simulations are Constant Charge (CC) or Constant Potential (CP) that represent an upper and a lower bound of the osmotic pressure and particle interaction forces between interacting surfaces. Another approach was introduced by Ninham and Parsegian in 1971¹⁴ and later expanded by Chan et al.¹⁵ as well as Behrens and Borkovec,¹⁶ called Charge Regulation (CR). In this approach the interacting surfaces adjust their chemistry, surface charges, and potentials in response to the overlapping EDLs, utilizing intrinsic equilibrium constants for the chemical speciation reactions. The resulting interfacial charges and potentials are in between the upper and lower bounds set by the CC and CP approaches. Approximation methods are suggested in literature, e.g. by combining the CC and CP solutions with a regulation parameter.¹⁶ While such approximations deliver good results for interaction forces, surface charges, and – potentials, they do not involve the calculation of the surface speciation, which is a particular focus of our model developments.

Current popular modeling software for chemical speciation does not include an option for charge regulated surface interaction and most of the rather few publications dealing with CR calculations are focused on interpreting data in the context of DLVO theory,^11,17,18 with the full codes usually notopenly available. The exception is the Python code SINFONIA, which was created by Gil-Díaz et al.¹⁹ for calculating the chemical speciation, surface charges and potentials as well as interparticle forces for the case of two approaching surfaces. It includes a four-layer SCM that can be simplified to the other electrostatic models mentioned above and a 1-D numerical solution of the PBE that is solved using an ordinary differential equation (ODE) solver as included in SciPy.

In this study, we developed a comprehensive SCM which includes CR, titled CRESCENDO (Charge Regulation, Electrostatic Surface Complexation and Equilibria on (2D) Nanoscale DOmains). It is capable of simulating overlapping EDLs in two dimensions on the Poisson–Boltzmann level, applicable to arbitrary geometries with flat surfaces. Using the Finite Element method, a numerical solution of the nonlinear PBE is obtained and incorporated into a full SCM, simulating the chemical speciation of all involved surface segments and adjusting it depending on the mutual influences through diffuse layer potentials. For this purpose, a custom, open-source Python code was developed, utilizing the FEniCSx library with a Newton solver for numerically solving the PBE and iterating it with a custom up to three-plane SCM.

The code is thoroughly benchmarked. First the numerical solution of the PBE is compared against analytical solutions for various electrolyte types, ionic strengths and diffuse layer potentials. The results of the SCM for a simple 1D – single surface system are compared to the implementation of an equivalent system in PHREEQC, both the speciation of a simple NaCl solution as well as monodentate Am-sorption onto a quartz surface using the parameters of García et al.²⁰ These are for the purpose of code validation as CR is not relevant for this calculation. The code is then compared to two benchmarks implemented in SINFONIA,¹⁹ validating 1D CR as well as asymmetric interactions in the case of surfaces with dissimilar chemistry.

Finally, the code is used to model a simple test case of the effect of overlapping EDLs on Am(III) sorption in a rectangular corner of two equivalent quartz surfaces, as it might exist on the rough surface of a mineral or as an approximation for a wall segment in an angular nanopore. This case demonstrates the capability of the algorithm to calculate the chemical speciation and electrostatic interactions in complex geometries with general electrolytes and nanoscale potential overlap.

Theoretical and technical framework

Surface charge and potential

Contact between a mineral surface and an aqueous solution leads to the development of a surface charge due to the presence of broken bonds and subsequent chemical reactions at the surface. The charge created by (de-)protonation of oxidic functional groups at the surface and specifically sorbed ions is balanced by an opposite charge of counter-ions in a diffuse layer at the solution side of the interface. This structure is called the electric double layer (EDL). The inner part of the interface and the specifically sorbed ions are referred to as the Stern layer, the outer part of the EDL is called the diffuse layer. The plane separating the Stern layer from the diffuse layer is called the d-plane. This plane is often assumed to coincide with the slip plane, and the d-plane potential is then approximated with the ζ-potential, but this assumption is frequently inaccurate.²¹

Multiple numerical models for the EDL are available that treat the diffuse and Stern layers differently. The simplest model is the Constant Capacitance Model (CCM) that assumes a linear potential drop from the surface into the solution, effectively excluding a diffuse layer and assuming that the entire surface charge is shielded by specifically adsorbed ions. This model can be rather accurate in high ionic strength solutions as such conditions suppress the extent of the diffuse layer. The other extreme, relevant for low ionic strength, is the Diffuse Layer Model (DLM) which assumes no Stern layer, resulting in the entire surface charge being balanced by a diffuse layer of counterions and a potential drop that follows the PBE. The simplest model that assumes both a Stern layer and a diffuse layer is the Basic Stern Model (BSM). Various other models exist that split the Stern layer into multiple smaller sections like the three-layer model (TLM), three-plane model (TPM) or the four-layer model (FLM) (Fig. 1).

Fig. 1 Schematic graphs of the electrical potential of a surface in contact with an aqueous solution for different surface complexation models. The basic Constant Capacitance Model (CCM) and Diffuse Layer Model (DLM) are shown in the top row and the three-plane model (TPM) in the bottom. The relevant equations for the capacities and charge-potential relationships are shown where applicable.

The electric charge of ions is located in hypothetical planes within the Stern layer, resulting in a linear variation of the potential (ψ) between these planes. Only in the diffuse layer a continuous change in charge density and correspondingly an exponential potential drop according to the PBE is considered. The charge-potential relationship within the Stern layer equals that of a flat plate capacitor

σ₀ = C₁(ψ₀ − ψ_d) (1)

with the electric potential at the surface ψ₀ and at the d-plane ψ_d (V), the surface charge density σ₀ (C m⁻²), Stern layer thickness d (m), capacity C₁ (F m⁻²), vacuum permittivity ε₀ (F m⁻¹), and the relative permittivity of the Stern layer ε₁ (dimensionless).

In a model that subdivides the Stern layer, each section acts as a flat plate capacitor with its own relative permittivity, capacity and charged plane(s) (e.g., TPM in Fig. 1).

The remaining unshielded charge at the head end of the Stern layer (d-plane) is balanced within the diffuse layer. The charge-potential relationship in the diffuse layer is typically calculated using the PBE or its simplifications for certain geometries, electrolytes and potential ranges like the Gouy–Chapman or Grahame equations. The PBE can be obtained from Poisson's equation of electrostatics:

with the relative permittivity ε_r (dimensionless) of the medium and the volume charge density ρ (C m⁻³), which can be calculated as the sum over all ion concentrations and their respective charges throughwhere c_i and z_i are the concentration (m⁻³) and the charge (dimensionless) of the i-th ionic species in solution and e is the elementary charge (C). Assuming that the distribution of ions follows Boltzmann statistics, due to the electric potential present in the diffuse layer caused by surface and layer charges, the expression for the local electric charge density changes towith c⁰_i being the concentration in the bulk solution and the Boltzmann constant k_B (J K⁻¹) as well as the absolute temperature T (K). Inserting the expression for the volume charge density into Poisson's equation results in a nonlinear partial differential equation for the potential profile in the diffuse layer, the so-called Poisson–Boltzmann equation:with ε_r being the relative permittivity of the diffuse layer. When working with molar concentrations a conversion factor of 1000 × N_A converts the concentrations from mol L⁻¹ to m⁻³, where N_A is Avogadro's constant (mol⁻¹). The general form of the PBE is valid for arbitrary electrolytes within the bounds of the mean-field approximation of the Poisson–Boltzmann theory.

Surface complexation models

The concept of surface complexation can be explained using a simple system of a surface exposing the surface functional group [double bond splayed left]SOH to a simple NaCl solution, described with a TPM. The surface site can deprotonate, leading to the reaction:

[double bond splayed left]SOH ⇌ [double bond splayed left]SO⁻ + H⁺,  K^app₁ (7)

where K^app₁ is the apparent equilibrium constant for the deprotonation of the hydroxyl group bound to the surface site [double bond splayed left]S, calculated using the mass law equation:with square brackets representing concentrations and curly brackets representing activities.

The apparent equilibrium constants K^app_i are calculated this way for every relevant surface reaction.²² Cations are able to adsorb to the surface following the reaction

[double bond splayed left]SOH + Na⁺ ⇌ [double bond splayed left]SONa + H⁺, K^app₂. (9)

Alternatively, the surface site might also be protonated

[double bond splayed left]SOH + H⁺ ⇌ [double bond splayed left]SOH₂⁺, K^app₃ (10)

and this positively charged site is then able to adsorb anions

[double bond splayed left]SOH₂⁺ + Cl⁻ ⇌ [double bond splayed left]SOH₂Cl, K^app₄. (11)

The equilibrium constants are called apparent because they depend on the surface potential. To obtain a true constant for the modelling of surface reactions they have to be split into an intrinsic term and a coulombic correction factor that includes the energy necessary to move an ion from solution through an electric potential towards the surface.²² The total free energy is therefore split into an intrinsic and a coulombic term

ΔG⁰_tot = ΔG⁰_int + ΔG⁰_coul. (12)

The coulombic term ΔG⁰_coul can be expressed as

ΔG⁰_coul = ΔzFψ (13)

where Δz is the amount of charge transferred to the surface by the reaction, F is Faraday's constant (C mol⁻¹) and ψ is the potential value of the plane the ion is moved to, in the case of the surface ψ = ψ₀. When combining the equations this results inwhere K^int_i is the so-called intrinsic equilibrium constant of a given reaction.²²

The reactions along with their equilibrium constants and coulombic correction factors can be used to calculate the chemical equilibrium of the surface and solution.

Following the method of Westall,²³ a linearly independent set of components is defined in a way that each species can be written as the product of a reaction of components and no component can be written as the reaction product of other components. Then the reactions forming all species present in the chemical system are formulated.

For the calculation of the chemical equilibrium the coulombic correction factors are treated as components of the system. The concentrations of species and components are gathered in vectors and vectors of activity coefficients are also calculated. The activity coefficients in this study were calculated using the Davies equation, an empirical formula only dependent on the ionic strength (IS) and ion charge

In the present example the species (C) and components (X) vectors would be
In the next step all reactions with their respective equilibrium constants and stoichiometries are defined, the equilibrium constants are gathered in another vector while the stoichiometries are written in matrix form, the so-called A-matrix. For convenience it is usual to write the reactions in log₁₀ scale. The resulting full description of the chemical system is shown in Table 1.

Table 1 Vectors and matrix used for calculation of mass and material balance for a 2-pK triple layer surface model with (de-)protonation of the surface site as well as specific sorption of anions and cations. Setup of the vectors and matrix follow the method of Westall²³

Species			Components					Eq. const.
log10 C + log10 γs	=		A · (log10 X + log10 γc)					−log10 K
	Na^+	Cl^−	H^+	[double bond splayed left]SOH	coul ψ0	coul ψβ	coul ψd
log10[Na^+]	1	0	0	0	0	0	0	0
log10[Cl^−]	0	1	0	0	0	0	0	0
log10[H^+]	0	0	1	0	0	0	0	0
log10[OH^−]	0	0	−1	0	0	0	0	−log10 KOH
log10[[double bond splayed left]SOH]	0	0	0	1	0	0	0	0
log10[[double bond splayed left]SO^−]	0	0	−1	1	−1	0	0	−log10 K^int1
log10[[double bond splayed left]SOH2^+]	0	0	1	1	1	0	0	−log10 K^int3
log10[[double bond splayed left]SONa]	1	0	−1	1	−1	1	0	−log10 K^int2
log10[[double bond splayed left]SOH2Cl]	0	1	1	1	1	−1	0	−log10 K^int4

---

The components coul ψ_i represent the decadic logarithm of the coulombic correction factors for the relevant planes. These vectors and the matrix define the chemical system. The entire system can be solved by

(log₁₀ C + log₁₀ γ_s) = log₁₀ K + A · (log₁₀ X + log₁₀ γ_c). (17)

To solve the chemical equilibrium problem a vector of the total analytical concentrations of the components, T, is also needed. The error in the material balance is calculated by obtaining the difference between the predefined total amount of components and the ones obtained by an initial guess using the equation

Y = A^T · C − T (18)

where Y is a vector of material balance residuals that need to be minimized in order to solve the equilibrium problem. In this study the optimization was done using the Newton–Raphson method as described by Westall.²³

In the above example the elements Y[4] to Y[6] are T_σ0,T_σβ and T_σd, the molar plane charges of the 0, β, and d planes (mol L⁻¹), which can be converted back to surface charge densities (C m⁻²) using a conversion factor

where s is the solid–liquid ratio (g L⁻¹) and a the specific surface area of the solid (m² g⁻¹).²⁴ The electrostatic balance can now be determined using the previously described methods of calculating the surface charges and potentials.

In the case of non-interacting surfaces, charge neutrality demands that at an infinite distance the electric charge is zero, meaning that all layer charges when added must be of equal charge and opposite sign to the diffuse layer charge

σ₀ + σ_β + σ_d = −σ_diff (20)

This assumption does not hold with interacting surfaces as the electric potential in the diffuse layer does not necessarily drop to zero with a non-infinite distance to the EDL of another surface. Instead, it is possible to use Gauss’ law to determine the integrated charge density in the Stern layer by calculating the derivative of the electric potential at the d-plane in the surface normal direction.

Imagining a thin cylindrical Gaussian surface with radius r (“pillbox”) on the charged surface we obtain the following equation for the plane charge through Gauss’ Law:

Because the D-field is perpendicular to the surface the equation for the flux can be calculated using Φ_D = πr²D⃑. The free enclosed charge can be calculated through the surface charge density on the analyzed surface Q_enc = πr²σ. Inserting these into eqn (20) leads to the following equation:

D⃑ = σ, with D⃑ = ε₀ε_rE⃑ (22)

where E⃑ = ∇ψ is obtained from numerical solution of the PBE.

For the current example this results in the equations

σ₀ = C₁(ψ₀ − ψ_β) (23)

σ_β = C₁(ψ_β − ψ₀) + C₂(ψ_β − ψ_d) (24)

σ_d = C₂(ψ_d − ψ_β) − ε₀ε_r · ∇ψ (25)

with the layer capacitances C_i defined for the TPM (Fig. 1). The surface charges and potentials get optimized iteratively in the same procedure as the rest of the chemical system.

Charge regulation between surfaces at small separations

The first calculations of interacting EDLs were done under the assumptions of either constant surface charge, allowing the potential to vary (CC) or constant potential while varying the charge (CP). While generally inaccurate to describe the equilibrium situation these two simplifications represent the upper and lower bounds of the actual charge-potential relationship in the case of interacting charged surfaces.¹⁶ Ninham and Parsegian¹⁴ were the first to describe a system where the electric potential between approaching surfaces is regulated by the equilibria of surface reactions, leading to a feedback between the surface charge and potential frequently called Charge Regulation (CR). A good approximation for full CR is the constant regulation approach in which a regulation parameter is used to obtain values between constant charge and constant potential solutions.¹⁶ The parameter can vary between 0 (CP) and 1 (CC), the values in between approximating charge regulation. It can be interpreted as a measure for how easy charge can be stored in the Stern layer(s) versus the diffuse layers. In cases where the Stern layer capacity dominates, the system exhibits constant charge properties. When the diffuse layer capacity dominates, the interaction is closer to constant potential.¹⁶ For decreasing surface separation the approximation error increases in comparison to an electrostatic SCM full regulation approach.¹¹

In previous studies CR is most commonly employed to model interaction forces between approaching surfaces in terms of DLVO theory.^11,16–18 However, the changes in surface speciation and diffuse layer potential are also significant for effects like anion exclusion in small pores²⁵ or trace element mobility mediated by surface sorption.²⁶ CR therefore needs to be included in the modeling of nanoporous or nano-structured media like clays or silica particles.^26–30

Full charge regulation in the current study is done by numerically solving the PBE. In the case of two or more interacting EDLs from multiple surfaces the calculation of the potential through the PBE is still valid, but the chemical system definition has to include every surface so that the chemical equilibria can be calculated simultaneously. The PBE gets solved on an entire domain that contains the relevant surfaces and Gauss’ law is applied by calculating the gradient of the electric potential perpendicular to the d-plane, and obtaining the respective values along the d-plane of the various surfaces. The resulting values are subsequently used for the next iteration of the Newton-Raphson solver until convergence (i.e., equilibrium) is reached. It is not necessary for the various surfaces to be described by the same chemical system as long as they are well defined within the chosen components, species, and equilibrium constant vectors as well as the stoichiometry matrix.

Finite-element solution for the Poisson–Boltzmann equation

Analytical solutions for the non-linear PBE are only known for simple geometries and electrolyte solution types. Exact determination of the diffuse layer charge for simple electrolytes is possible in the cases of a 1-D semi-infinite plane, and approximations for spheres and cylinders exist,³¹ but solutions for more complicated geometries and electrolyte types remain elusive. It is therefore common to solve the PBE, like other partial differential equations, numerically. A common mathematical method for this is the so-called Finite Element Method (FEM).³² A rigorous derivation of the mathematical fundamentals for the FEM is beyond the scope of this study, but a very brief overview for the steps needed to convert the PBE to a form that can be entered into a Python-based FE-solver will be provided. For a detailed description of the FEM and approaches to solving it in Python we refer to other works.^32–34

The very basic idea of the FEM is to solve the partial differential equation (PDE) only at certain nodes in space instead of covering the continuous domain and approximate the function between those nodes. The necessary step to achieve this is to rewrite the PDE as the variational formulation over an interval I = [0, L], also sometimes called the weak form. Taking as an example the 1D PBE this is achieved by multiplying both sides of the equation by a test function v(x), which is defined as vanishing at the endpoints of the interval and then integrating by parts:

The step from eqn (27) and (28) is only possible when the assumption that v(L) = v(0) = 0 is fulfilled and it simplifies the equation by reducing the order of the derivative of ψ(x) from two to one. Analogous to the test function v(x), the other function ψ(x) is often called the trial function. Both are part of a function space V₀. With the test function v(x), trial function ψ(x) and source function q(ψ(x)) defined, the weak form of the PDE can be used to solve the equation on a defined domain, with a properly defined function space, V₀, and boundary conditions.³²

In the case of this study the FE solver for the non-linear problem used was a Newton Solver implemented using FEniCSx.

Benchmark setup

The implementation of the numerical PBE solver was benchmarked against analytical solutions on a square domain, meshed using quadrilaterals and a first order Lagrange space. The side length of the domain was scaled to be a multiple of the Debye lengthas the ionic strength IS of the solution has a major impact on screening of the surface charge and therefore the extent of the diffuse layer. The domain was meshed with 40 by 40 nodes and the mesh density was increased towards the side of the domain to which the diffuse layer potential (ψ_d) was applied as a boundary condition. To achieve this, the side lengths of the cells were scaled by a geometric series with an exponent of 1.1. For the test cases with high d-plane potential (−100 mV) the mesh density was increased to 80 by 40 nodes. This was done to reduce the numerical error of the FE approximation in the region where rate of change of the potential is maximal.

The analytical solutions are valid for the assumption of a planar, semi-infinite solid surface. Analytical solutions were calculated for a symmetric electrolyte, 1 : 2 and 2 : 1 asymmetric electrolytes as well as an approximate analytical solution for 1 : 1–1 : 2 mixed electrolyte.^6,7,31 The concentrations of the different ionic species were 0.1 M, 0.01 M, or 0.001 M. In the asymmetric and mixed cases this was the concentration of the divalent ion and the monovalent ion had its concentration doubled. The d-plane potential applied to the surface as a Dirichlet boundary condition was chosen to be −100 mV, −10 mV or −1 mV.

An identical mesh to the one from the PBE solver benchmark was used for comparing the solution obtained by CRESCENDO against an equivalent system implemented in PHREEQC. The geometry is shown in Fig. 2(a). The surface and solution speciation of a quartz surface in equilibrium with a 0.01 M NaCl solution was modeled. In a second benchmark the sorption and speciation of Americium(III) on the same quartz surface was calculated. The speciation scheme from García et al.²⁰ was used, with Am as a chemical analogue to Eu and the assumption that the speciation and sorption behavior of Eu and Am are very similar.³⁵ The side length of the domain for this benchmark was fixed at 30 nm. This was roughly equal to 10 × λ_D as the ionic strength of solution was around 0.01 M for most pH values and Am-concentrations. Towards the very low and very high end of the pH scale the ionic strength increases due to the necessary addition of NaOH or HCl.

Fig. 2 Schematic images of the geometries considered for the benchmarks and the test case of charge regulation between perpendicular surfaces with isopotential lines (green). (a) Single surface in contact with solution as used for benchmarks 2 and 3. (b) Two parallel surfaces as used for benchmark 4. (c) two perpendicular surfaces as used for the test case of charge regulation.

As a final benchmark the SINFONIA code benchmarks from Gil-Díaz et al.¹⁹ were obtained from the GitHub repository linked in the paper to validate charge regulation. The model “Bench_1a_Silica” calculates the charge regulation of two approaching silica surfaces (Fig. 2(b)) for varying pH in a 0.01 M NaCl solution utilizing a simplified surface speciation scheme that includes only a single deprotonating surface site and no specific sorption of electrolyte ions from the solution. This same model was implemented in CRESCENDO, the d-planes of the two surfaces were defined as the left and right edges of the domain and for the distance variations the horizontal extent of the domain was changed from 28 nm down to 0.1 nm. The surface charge and d-plane potentials are compared for various pH values and distances between the surfaces. The calculation “Bench_3c_Goethite-Rutile” simulates two approaching surfaces of different chemistry, with [double bond splayed left]TiOH and [double bond splayed left]FeOH functional groups respectively, in a 0.001 M NaCl solution implemented using a 1-pK charging formalism with fractional charges. This is replicated in CRESCENDO through the use of a separate B-Matrix in the mass-balance equation to handle partial charges. The same method of varying the horizontal extent of the domain was used, modeling d-plane separation distances from 50 nm down to 0.1 nm at three different pH values. The values of ψ_d and σ₀ of the two surfaces are compared between the codes.

For the test case of CR between perpendicular surfaces (Fig. 2(c)) the mesh had the same geometric progression of cell side length applied to two edges (Fig. 3). The solution chemistry and speciation scheme were identical to the ones used in the previous benchmark. The number of nodes on the surfaces was chosen to be 40 and the surface speciation as well as layer charges and potentials were calculated at every node on the surface. The resulting mesh density is therefore maximal in the area of interest in the corner where the impact of charge regulation is the highest, decreasing towards the flat edges where interaction between surfaces is negligible and speciation is almost equivalent to an isolated flat surface. This also minimizes the numerical error of the finite element solution in the corner where the change in the slope of the potential is highest.

Fig. 3 The mesh used for the charge regulation benchmark. The d-planes of the surfaces were defined to be at the left and bottom edges of the domain, the quadrilateral mesh units are refined towards both surfaces as well as towards the corner.

The surface speciation, plane charges, and potentials for each segment of the surface were calculated individually, but simultaneously. The domain was defined to be symmetrical along the diagonal. When setting up the A matrix (similar to Table 1) each surface segment had its own surface species defined for proper calculation. This leads to a very large A matrix with 40 nodes requiring 40 × 10 surface species. The total amount of surface area in contact with solution for both surfaces combined was set to 40 m² L⁻¹ and the surface concentration of each segment was calculated to scale the resulting surface species' concentration to the varying sizes of these segments.

The code

All the tasks of solution and surface definition, mesh generation, finite-element solving of the PBE and Newton–Raphson optimization of the chemical speciation are handled by the code written in Python 3 (version 3.12.3). The mesh is generated using the Python api of gmsh (version 4.12.1). The non-linear PBE was solved using the open-source FEniCSx Python interface (version 0.1) for solving partial differential equations utilizing the FEM.

The inputs are entered in one Python file, requiring the definition of the solution and surface as well as the speciation reactions using the vector and matrix definitions described earlier. The setup and the mesh can be adapted to a variety of different applications and geometries.

The surface and solution speciation are optimized iteratively using the Newton–Raphson algorithm with the various layer charges and potentials also being recalculated and optimized at each step. Once convergence is reached the data is stored in a text file. The code for CRESCENDO, set up as it was used for calculation of the CR Benchmark, is available on GitLab via https://gitlab.kit.edu/kit/ine/Geochemistry/crescendo.

Results and discussion

Benchmarks

Benchmark 1: numerical vs. analytical PB solutions. The maximum absolute differences between the analytical and numerical solutions are summarized Table 2 and exemplary data is shown in Fig. 4. Scaling the domain size with λ_D causes the errors for different ionic strength values to be identical. Therefore, only varying electrolyte type and potential are shown.

Table 2 Maximum difference between the analytical and numerical solutions for various electrolyte types and d-plane potentials. Variation in ionic strength is not listed as scaling the domain with the Debye length resulted in identical errors for changes in ionic strength

Electrolyte type	d-plane potential (ψd)
	−1 mV	−10 mV	−100 mV
Symmetric (1 : 1)	3.316 × 10^−04 mV	3.305 × 10^−03 mV	2.524 × 10^−02 mV
2 : 1	5.016 × 10^−04 mV	3.173 × 10^−03 mV	1.809 × 10^−02 mV
1 : 2	3.327 × 10^−04 mV	3.412 × 10^−03 mV	2.736 × 10^−02 mV
1 : 1 : 2 : 1	3.307 × 10^−04 mV	3.205 × 10^−03 mV	1.874 × 10^−02 mV

---

Fig. 4 Potential profiles, difference, and relative error of some of the analytical and numerical solutions for ψ_d = −10 mV. The concentration of the lowest concentration ion was 10 mmol, leading to ionic strengths and λ_D values of 10 mmol and 3.04 nm for the symmetric electrolyte, 30 mmol and 1.76 nm for the 2 : 1 electrolyte, as well as 40 mmol and 1.52 nm for the mixed electrolyte. The top row shows the potential profiles with the blue line being the numerical solution and the dashed red line being the analytical solution. The middle row shows the difference Δ_{analytical–numerical} with the maximum of the absolute difference also shown. The third row shows the relative error, the red dashed lines are a visual aid for 2% error. The error between the numerical and analytical results are very low, the finite element solution is a very good approximation for the true values.

As shown in Fig. 4, the absolute difference reaches its maximum at the point in the potential profile where the change of the slope is still high but the mesh density is already decreasing, leading to a peak in numerical error. This effect is most noticeable when the analytical value is changing significantly in between finite element nodes. It is caused by the approximation between nodes, resulting in the calculated value oscillating between values slightly higher and slightly lower than the analytical value. The relative error for the range of interest close to the surface is less than one percent, increasing towards the far end of the domain due to the numerical noise where the analytical solutions approach 0. Overall, differences and relative errors are low for every tested configuration. The maximum absolute difference is always more than three orders of magnitude lower than the applied d-plane potential in every case and the relative error remains well below 2% for distances up to 5 × λ_D at which point the potential decayed to a value close to 0. These results show that the numerical solution of the PBE implemented in the code is a very good approximation for the actual values determined from analytical solutions.

Benchmark 2: surface and aqueous speciation – PHREEQC. The speciation scheme of García et al.²⁰ describes the silica surface using a two-site protolysis model with a more acidic site ([double bond splayed left]S_xOH, pK = −4.0, 1 site per nm²) and a regular site ([double bond splayed left]S_yOH, pK = −8.5, 3.7 sites per nm²). Overall surface concentration was defined to be 1 m² L⁻¹. This was implemented in both PHREEQC and CRESCENDO and the equilibrium speciation of the surface in contact with a 0.01 M NaCl solution was calculated.

The surface charge density and d-plane potential values were obtained from EDL calculations using a CD-MUSIC model in PHREEQC. A BSM was employed with a Stern layer capacitance of 2.0 F m⁻¹. CRESCENDO reproduces both the surface and aqueous speciation of the PHREEQC results (Fig. 5a) as well as the surface charge and d-plane potential values (Fig. 5b and c) very well. The ionic strength of the solution is increased towards low and high pH due to the necessary addition of HCl and NaOH to reach the pH values while maintaining charge neutrality of the bulk solution.

Fig. 5 Surface speciation (a), d-plane potential (b) and surface plane charge density (c) of a quartz surface in contact with a solution containing 0.01 M NaCl, calculated in PHREEQC (red squares) and CRESCENDO (solid lines). The ionic strength (d) increases towards both ends of the pH scale due to the addition of acid/base.

Benchmark 3: Am-adsorption – PHREEQC. For this benchmark a setup identical to the previous was used, but the solution now also contained 10⁻⁶ M of AmCl₃. The Am-species that were considered in solution and at the surface are Am³⁺, Am(OH)²⁺, Am(OH)₂⁺, Am(OH)₃, AmCl²⁺, AmCl₂⁺, [double bond splayed left]SOAm²⁺ and [double bond splayed left]SOAm(OH)₂. Bidentate sorption of Am on quartz is considered to be the actual binding mode,³⁶ but a monodentate sorption model was implemented here for simplicity and to avoid ambiguous handling of bidentate sorption in different speciation codes.³⁷ The implementation can be adjusted for multidentate sorption if necessary. The surface speciation of Am as well as σ₀ and ψ_d values were then compared between PHREEQC and CRESCENDO (Fig. 6a–c). The results of CRESCENDO again perfectly reproduce the values obtained from the PHREEQC implementation.

Fig. 6 Surface speciation (a), d-plane potential (b) and surface plane charge density (c) of a quartz surface in contact with a solution containing 0.01 M NaCl and 10⁻⁶ M AmCl₃, calculated in PHREEQC (red squares) and with CRESCENDO (solid lines). The relative fraction of the Am surface species (d) shows that at medium pH the positive surface species is dominating while at pH > 5.8 the neutral species is prevalent.

The values of ψ_d and σ₀ are significantly altered by the presence of Am at the surface (Fig. 6b and c). The d-plane potential no longer shows a monotonic decrease as it did for the silica surface without AmCl₃ in solution: after a slight decrease it increases towards a local maximum at pH 5.5–6 before it decreases again, indicating that the sorption of positively charged Am(III) species increases ψ_d at circumneutral pH. The graph of σ₀ is identical for low pH values until uptake of Am³⁺ begins at around pH 4. Towards higher pH values σ₀ is consistently more negative compared to the system without AmCl₃. The sorption of Am³⁺ in the Stern layer enhances the deprotonation of the hydrolyzed surface groups, leading to a more negative surface charge.

An Am sorption edge was calculated from the speciation data of both codes (Fig. 6d). Uptake begins at around pH 3.5. At pH < 5.8 the dominant surface species is the positively charged [double bond splayed left]SOAm²⁺. At pH ≈ 5.8 the fraction of both species is about equal and at pH > 7.5 all Am is present as the neutral species [double bond splayed left]SOAm(OH)₂. At pH > 7 almost all Am in the system is adsorbed to the surface. These observations explain the changes visible in the surface charge and d-plane potential graphs of the quartz surface with and without Am³⁺ present.

Benchmark 4: comparison of charge regulation codes. Surface charge density and d-plane potential values obtained from SINFONIA calculations are compared with the ones obtained from CRESCENDO. Fig. 7 shows the values for two interacting quartz surfaces. With decreasing distance between the surfaces, the absolute value of ψ_d increases towards a more negative potential while σ₀ becomes less negative, indicating suppressed surface deprotonation due to charge regulation effects. This effect is more noticeable at higher pH. CRESCENDO reproduces the values from SINFONIA for all pH values and interplanar distances.

Fig. 7 Comparison between the values for d-plane potential ψ_d (a) and surface charge density σ₀ (b) obtained from SINFONIA (red dots) and CRESCENDO (solid lines) for a simple quartz surface model in contact with a 0.01 M NaCl solution. The distance between the surfaces for which the charge regulation was calculated was varied between 28 nm and 0.1 nm, shown on the line labels. For low separations the d-plane potential becomes more negative while the surface charge density increases.

The effects of CR between surfaces of different chemistries are shown in Fig. 8. At high separation and low pH, the rutile surface (TiOH) is positively charged, with decreasing distance the surface charge density becomes more negative. The goethite surface (FeOH) shows an opposite trend, with very little variation at pH 3.8. ψ_dTiOH is lower than ψ_d FeOH at high separation. With decreasing separation between the two surfaces the d-plane potentials change, ψ_d TiOH increasing and ψ_d FeOH decreasing as the surfaces interact in close proximity.

Fig. 8 Comparison between d-plane potentials ψ_d (a) and (c) and surface charge density σ₀ (b) and (d) obtained from SINFONIA (red dots) and CRESCENDO (solid lines) for two approaching surfaces of different chemistry at various separations in contact with 0.001 M NaCl. One surface represents Rutile (top row) interacting with Goethite (bottom row), modeled at three different pH values.

Charge regulation in a corner of perpendicular flat surfaces

Two quartz surfaces were modeled in contact with a solution containing NaCl and AmCl₃.

The flat surfaces are perpendicular and interact with each other on a two-dimensional square domain with a side length of 30 nm. The surfaces meet at a right angle in the bottom left corner of this domain, which becomes the area of interest for CR. The mesh used for this model is shown in Fig. 3. The speciation, plane charges, and d-plane potential of the flat part of the surfaces far away from the corner are almost identical to a single surface in contact with the solution while the interacting diffuse layers in the corner cause significant deviation from that behavior (Fig. 9 and 10).

Fig. 9 Exemplary modeling result for charge regulation and constant potential at pH 6. The difference in potential across the domain (a) and (b) shows the impact of charge regulation in the bottom left corner (isopotential lines are shown in red, at intervals of 0.01 V). The electric field at the d-plane at y = 0 (c) and (d) shows strong influence of CR, with higher values when compared to the CP case which decreases down to 0 V m⁻¹. Potential profiles over the entire domain (e) and (f) show two-dimensional variation in potential. In the CP case the potential at y = 0 is flat whereas in the CR case it decreases towards the corner. The step size in y between each line is 0.75 nm.

Fig. 10 Effect of charge regulation at perpendicular quartz surfaces in contact with a solution containing 0.01 M NaCl and 10⁻⁶ M AmCl₃. The left column shows the amount of sorbed Am(III) (a), the surface charge density σ₀ (c), the β-plane charge density (e), and the d-plane potential (g) across the surface. The x axis represents the distance from the corner at x = 0. The right column shows the excess sorption of Am(III) in the corner (b), the excess surface charge in the corner (d), the excess β-plane charge in the corner (f) as well as the d-plane potential at the corner and a point far away from the corner, compared to the PHREEQC model of an isolated surface (h), versus pH.

Charge regulation is observable in the bottom left corner of the false-color potential maps, shown with isopotential lines (Fig. 9a and b). Comparison to a Constant Potential case shows that for most of the flat surface the potentials are identical while in the corner the potential deviates significantly. The derivative of ψ perpendicular to surface is the electric field used as input for Gauss’ law in eqn (25) to calculate the Stern-layer charge (Fig. 9c and d). pH and AmCl₃ concentration were varied in numerous simulations and the results are described in the following sections. Multiple effects of charge regulation can be observed.

Surface speciation, charges, and potentials vs. pH

The surface speciation of a solution containing 0.01 M NaCl and 10⁻⁶ M AmCl₃ was calculated for a range of pH 3 to 10. The results are shown in Fig. 10. Changes in the amount of Am³⁺ sorbed to the surface, σ_β, and ψ_d in relation to distance from the corner are immediately obvious (Fig. 10a, e and g), the variation in σ₀ is more subtle (Fig. 10c).

The total amount of Am at the surface increases with pH, leveling out on the flat part of the surface distant from the corner around pH 7 (in agreement with the sorption edge calculated in Fig. 5d). Charge regulation increases the uptake of Am in the corner compared to the flat surface at all pH values, with a peak of excess uptake between pH 5 and 6 (Fig. 10a and b). This value coincides with the pH at which the positively charged [double bond splayed left]SOAm²⁺ species is the dominant Am surface species (Fig. 5d). At pH > 6, when the dominant Am surface species is the neutral species [double bond splayed left]SOAm(OH)₂, the excess sorption decreases until the amount of Am sorbed in the corner and on the flat surface are almost equal.

The surface charge density σ₀ for different pH is shown in Fig. 10c and the difference in σ₀ between the corner and the flat part is shown in Fig. 10d. σ₀ of quartz in contact with this solution is negative for pH > 3, positive difference values therefore mean that σ₀ on the flat part is lower (more negative) than in the corner, negative difference values mean σ₀ in the corner is more negative. At low pH overall sorption is low and the difference is positive, decreasing towards higher pH. Between pH 5.5–7, where the excess sorption of Am in the corner is maximal, the difference is negative, indicating a more negative σ₀ in the corner. At higher pH values, where the excess sorption in the corner is negligible, the difference increases again.

The charge density in the β-plane σ_β increases with increasing pH, with the corner showing a more positive charge for all calculated pH values (Fig. 10e). The difference in σ_β between the corner and the flat part of the surface (Fig. 10f) increases with increasing pH up to a local maximum around the pH of maximum excess sorption. Towards higher pH the difference decreases slightly before increasing again at high pH > 9. Sorption of Am at pH < 7 is predominantly via the positively charged [double bond splayed left]SOAm²⁺ species, which involves deprotonation of the surface sites, leading to the negative σ₀ difference, and positive σ_β difference at the pH of maximum excess sorption. When the dominant Am species changes to the neutral [double bond splayed left]SOAm(OH)₂ species at higher pH σ_β difference decreases slightly because this species moves less charge into the β-plane compared to the positive species. The increase at high pH is due to sorption of Na⁺ in the β-plane.

ψ_d in the corner is more strongly negative for all pH values calculated here (Fig. 10g and h). The d-plane potential values of the flat part of the surface distant from the corner are identical to a PHREEQC model of the quartz surface without charge regulation. The difference between ψ_d in the corner and on the flat surface is relatively constant at around −9 mV for pH > 6, at which point the d-plane potential in the corner becomes even more negative compared to the flat part. At pH > 7 the difference is mostly constant again, around −17 mV. The pH value of this change correlates well with the excess Am sorption, the changes in surface charge density, and the change in dominant Am surface species. Higher sorption of positively-charged Am species in the corner shields the negative surface charge less than H⁺ ions on the flat part of the surface. At higher pH the overall neutral Am species cannot shield the charge more strongly than on the flat part of the surface so the difference remains constant between pH values.

Surface speciation, charges, and potentials at varying Am(III) concentrations

For pH 6 the surface speciation was calculated at varying AmCl₃ concentrations between 10⁻¹² M and 10⁻⁵ M (Fig. 11). pH 6 was chosen because the effects of charge regulation were strong for this value in the previous calculations. The values for concentrations of AmCl₃ < 10⁻⁸ M are almost indistinguishable but significant differences are observable for higher concentrations.

Fig. 11 Effect of charge regulation on perpendicular quartz surfaces in contact with a solution containing 0.01 M NaCl and varying amounts of AmCl₃ at pH 6. The changes in d-plane potential (a) and surface charge density σ₀ (b) are shown in the top row while the middle row shows the excess sorption of Am in the corner normalized to the flat part of the surface (c) and the β-plane charge density σ_β (d). The bottom row shows the relative amount of Am sorbed from solution, summed up over both surfaces, (e) and the speciation of Am across the surface (f). The x axis for all graphs except (e) is the distance from the corner at (x = 0).

ψ_d (Fig. 11a) shows more negative values in the corner for all concentrations. For higher concentrations the difference in d-plane potential between the corner and the flat part decreases and the absolute value of the d-plane potential decreases as well (becoming less negative).

The values of σ₀ are less negative in the corner for concentrations less than 10⁻⁷ M while at higher concentrations σ₀ in the corner is more negative than on the flat part of the surface (Fig. 11b). σ₀ for 10⁻⁵ M AmCl₃ is significantly more negative than for lower concentrations, the sorption of Am³⁺ to the surface enhances the deprotonation of the surface functional groups. At low concentrations the sorption of AmCl₃ in the corner does not shield the surface charge well, so H⁺ ions are also attracted into the corner which leads to an increase in σ₀. The trend in the β-plane charge density is opposite, with higher concentrations showing higher positive charge densities due to Am(III) sorption. In the corner the β-plane charge densities are more positive than on the flat part for all concentrations (Fig. 11d). In this plane the variation in charge density for concentrations less than 10⁻⁸ is again indistinguishable.

The relative amount of sorbed Am normalized to the concentration at the surface distant from the corner (Fig. 11c) shows that for lower concentrations of AmCl₃ the amount of Am in the corner is higher by factor of 3.5 compared to the flat part of the surface while for increasing concentrations this factor decreases. At 10⁻⁵ M AmCl₃ the distribution is almost flat, with only a slight increase towards the corner.

When investigating the amount of Am sorbed to the surfaces compared to the total amount of Am present in solution it is observable that about 99.9% of Am is sorbed to the surface for concentrations <10⁻⁶ mol L⁻¹, decreasing to 99.3% at 10⁻⁵ mol L⁻¹ of AmCl₃ (Fig. 11e). This is likely due to reduced attraction of Am³⁺ as the positive charge of sorbed Am³⁺ repels cations and the negative surface charge is shielded more strongly.

Comparison of the speciation of Am at the surface shows that at pH 6 the positively charged [double bond splayed left]SOAm²⁺ is the dominant surface species for all concentrations but the relative fraction of the neutral species [double bond splayed left]SOAm(OH)₂ increases with increasing AmCl₃ concentration (Fig. 11f) until the fraction of the two surface species is almost equal at 10⁻⁵ M. Charge regulation in the corner results in a higher fraction of the positive species compared to the flat part of the surface for all concentrations.

Comparison of computation time

Computation time for the chemical equilibrium simulations varies between the different programs. The time to calculate the solution and surface speciation of an equivalent simple system (Am(III) complexation on a silica surface like in Benchmark 3) on the computer used for this study was about 0.13 seconds for PHREEQC, 1.4 seconds for SINFONIA, and 22 seconds for version 1.0.0 of CRESCENDO.

PHREEQC, written in C++, is significantly faster than either of the Python-based codes. However, as their capabilities are significantly different this comparison is not simple. PHREEQC only calculates a 1D solution of a single surface in contact with the electrolyte using an analytical solution of the PBE and SINFONIA can calculate a 1D solution of two interacting parallel surfaces using an ODE solver for the PBE. Meanwhile CRESCENDO is able to calculate the chemistry of interacting planar surface segments, numerically solving the PBE on an arbitrary 2D domain with the FEM. Calculating the Jacobian matrix and solving the PBE on a 2D domain using the FEM every iteration are computationally intensive tasks and therefore slow, but other available software has so far not been able to calculate electrostatic interaction of surfaces either at all or in more complex geometries than parallel surfaces.

Furthermore, optimizations to improve the runtime of CRESCENDO calculations are being explored and will be implemented in future versions, with potential to reduce the calculation time by more than a factor of two.

Comparison to experimental studies

Few previous studies have used SCM to determine the chemistry of interacting surfaces while taking CR into account. The results of benchmarks included in the 1D CR code SINFONIA¹⁹ have been reproduced in this study and confirm the observations that CR effects cause a quartz surface in close proximity to another identical surface to show more positive surface charge and more negative d-plane potential than an isolated surface.

Some other studies investigating charging behavior of silica have been conducted using a simple chemical model for the silica surfaces. A CR model involving silica surfaces of varying roughness has recently been demonstrated by Ozcelik and Barisik.³⁸ In their model the dips of the rough surface show less negative surface charge densities and more negative potentials due to overlapping EDLs when compared to a flat surface in contact with an equivalent solution. Yang et al.³⁹ investigated charge regulation in silica nanopores of varying sizes in dependance of pH and ionic strength of the solution. Their results show that interacting EDLs in nanopores cause σ₀ of small pores to be less negative when compared to larger pores, with an asymptotic limit at the pore size at which opposing EDLs do not influence each other anymore. Sen and Barisik observed similar effects in silica nanochannels⁴⁰ and mesoporous silica.⁴¹ These results are consistent with the effects of charge regulation observed in the case of CR between perpendicular quartz surfaces shown in this work, with less negative surface charge density and more negative d-plane potential in the corner where the silica surfaces interact. These changes then lead to the altered sorption behavior of Americium(III) in the corner.

Some studies of trivalent Lanthanide sorption on mesoporous silica have also observed increased uptake of various Ln³⁺ ions in smaller pores, compared to larger pores with less diffuse layer interaction or non-porous surfaces. In the study by Ilgen et al.²⁹ overall Ln³⁺ uptake was higher for mesoporous SiO₂ with mean pore diameters of 4.4 nm than 7 nm. However, both were lower than the uptake on non-porous SiO₂, which the authors attributed to diffusion limited transport of the lanthanide ions through the long silica nanochannels. A later study by Fashina and Ilgen³⁰ then observed higher uptake by the porous SiO₂ compared to the non-porous fumed silica while using the same silica material but increasing equilibration time for the sorption experiments.

The main differentiating aspect of CRESCENDO when compared to previous 2D charge regulation approaches is the flexible layout of the chemical solver, allowing for comparatively simple implementation of complex solution and surface chemistries as demonstrated with the Am(III)-silica example. Previous implementations relied on hard-coded chemical reactions and charging equations, including only a very limited number of surface reactions and mostly disregarding solution chemistry.

Conclusions

Various benchmarks and test cases for CRESCENDO, a new Python-based surface complexation model code, have been described and their results presented. The code includes electrostatic calculations of up to a TPM, surface and solution speciation, and surface sorption under charge regulation conditions. The FE solution of the non-linear PBE has been validated against analytical solutions for a 1D semi-infinite planar geometry. The speciation code reproduced an equivalent model implemented in PHREEQC for solution composition, surface speciation and sorption. Benchmarks for 1D charge regulation of opposing planar surfaces implemented in SINFONIA were also reproduced in the new model. An example for 2D charge regulation has been presented for the case of two interacting perpendicular quartz surfaces and effects of charge regulation on surface speciation as well as interface charges and potentials have been extensively discussed. This work expands on previous charge regulation codes by enabling the modeling of two-dimensional systems including arbitrary configurations of flat surfaces on the nanoscale at the Poisson–Boltzmann level, whereas previous codes were limited to one-dimensional simplifications of opposing planar surfaces or simplified surface chemistries.

The consideration of charge regulation is vital for more accurate modeling of surface complexation, sorption, and interparticle forces on the nanoscale. Solving the PBE numerically allows for the consideration of complex geometries for which analytical solutions are not available, obtaining the electrostatic environment and surface speciation of nanoparticles or surfaces in nanoconfined settings. This can be of interest for a large variety of questions in many fields like material or environmental science. Being able to calculate accurate diffuse layer potentials for arbitrary electrolytes is another advantage of solving the full PBE as opposed to simplified analytical solutions which often only work for specific electrolyte types. This allows the code to be used for more complex systems with general solution compositions as encountered in biological applications, cementitious materials, and most environmental settings.

Next development steps will include application of the model to a variety of systems like nanopores and -particles with rounded surfaces or sorption on inhomogeneous surfaces, as well as problems like spillover effects on mineral edges. Further development of the model will allow it to be applied to problems like retention of radioactive isotopes in natural and engineered clay barriers in deep geologic repositories or trace element sorption in nanoporous media in general.

Author contributions

Conceptualization, F. H. and L. S., methodology, F. H. and L. S., software development, L. S., writing – original draft, L. S., writing – review & editing, L. S. F. H. and J. L., funding acquisition F. H.

Conflicts of interest

There are no conflicts to declare.

Data availability

The source code for CRESCENDO can be found within the GitLab repository at https://gitlab.kit.edu/kit/ine/Geochemistry/crescendo with https://doi.org/10.5281/zenodo.17750551. The version of the code employed for this study is version 1.0, which is set up to calculate the charge regulation of perpendicular surfaces test case.

Acknowledgements

The authors acknowledge the funding by Deutsche Forschungsgemeinschaft (Project number 545761212) within the French-German collaborative research project “Broadband dielectric and geochemical characterization of clay materials for underground storage (BARRIERS)”

References

1. Y. Zhang and W. Luo, BioResources, 2014, 9(2), 2484–2499.
2. Y. Liu, L. Zhang, Y. Wen, H. Zhai, Y. Yuan, C. Guo, L. Wang, F. Wu, C. Liu, J. Xiao, J. Liu, X. Yang, Y. Cai, J. Ji and Y. Liu, Sci. Total Environ., 2024, 948, 174856.
3. K. Fukushit and T. Sato, Environ. Sci. Technol., 2005, 39, 1250–1256.
4. S.-Y. Liang, W.-S. Lin, C.-P. Chen, C.-W. Liu and C. Fan, Appl. Sci., 2021, 11(13), 5879.
5. T. Hiemstra and J. Lützenkirchen, Rev. Mineral. Geochem., 2025, 91A, 13–84.
6. X. Liu, H. Li, R. Li and R. Tian, Surf. Sci., 2013, 607, 197–202.
7. Z. Chen and R. K. Singh, J. Colloid Interface Sci., 2002, 245, 301–306.
8. D. L. Parkhurst and C. A. J. Appelo, PHREEQC2 user's manual and program, US Geological Survey, Denver, Colorad, 2004.
9. J. P. Gustafsson, Visual MINTEQ 3.1, https://vminteq.com/, (accessed 27 Nov., 2025).
10. J. C. Westall, FITEQL: A computer program for determination of chemical equilibrium constants from experimental data, Department of Chemistry, Oregon State University, 1982.
11. G. Trefalt, S. H. Behrens and M. Borkovec, Langmuir, 2016, 32, 380–400.
12. T. Preocanin, A. Abdelmonem, G. Montavon and J. Luetzenkirchen, Clays, Clay Minerals and Ceramic Materials Based on Clay Minerals, 2016, ch. 3 DOI:10.5772/62082.
13. X. Meng and R. D. Letterman, Environ. Sci. Technol., 2002, 27, 970–975.
14. B. W. Ninham and V. A. Parsegian, J. Theor. Biol., 1971, 31, 405–428.
15. D. Chan, J. W. Perram, L. R. White and T. W. Healy, J. Chem. Soc., Faraday Trans. 1, 1975, 71, 1046–1057.
16. S. H. Behrens and M. Borkovec, J. Chem. Phys., 1999, 111, 382–385.
17. F. J. Ruiz-Cabello, G. Trefalt, P. Maroni and M. Borkovec, Langmuir, 2014, 30, 4551–4555.
18. C. Zhao, D. Ebeling, I. Siretanu, D. van den Ende and F. Mugele, Nanoscale, 2015, 7, 16298–16311.
19. T. Gil-Diaz, D. Jara-Heredia, F. Heberling, J. Lutzenkirchen, J. Link, T. Sowoidnich, H. M. Ludwig, M. Haist and T. Schafer, Adv. Colloid Interface Sci., 2021, 294, 102469.
20. D. García, J. Lützenkirchen, V. Petrov, M. Siebentritt, D. Schild, G. Lefèvre, T. Rabung, M. Altmaier, S. Kalmykov, L. Duro and H. Geckeis, Colloids Surf., A, 2019, 578, 123610.
21. T. Klačić, J. Katić, D. Kovačević, D. Namjesnik, A. Abdelmonem and T. Begović, Rev. Mineral. Geochem., 2025, 91A, 295–336.
22. D. A. Dzombak and F. M. M. Morel, Surface Complexation Modeling: Hydrous Ferric Oxide, Wiley, 1991.
23. J. Westall, in Particulates in Water, ed. M. C. Kavanaugh and J. Leckie, 1980, pp. 33–44 DOI:10.1021/ba-1980-0189.ch002.
24. F. Heberling, J. Lützenkirchen and T. Gil-Díaz, Encyclopedia of Solid-Liquid Interfaces, 2024, pp. 215–229 DOI:10.1016/b978-0-323-85669-0.00087-8.
25. C. Tournassat, I. C. Bourg, M. Holmboe, G. Sposito and C. I. Steefel, Clays Clay Miner., 2024, 64, 374–388.
26. C. Tournassat and C. I. Steefel, Rev. Mineral. Geochem., 2015, 80, 287–329.
27. X. Meng and R. D. Letterman, Environ. Sci. Technol., 2002, 27, 1924–1929.
28. C. Tournassat, I. C. Bourg, C. I. Steefel and F. Bergaya, in Natural and Engineered Clay Barriers, ed. C. Tournassat, I. C. Bourg, C. I. Steefel and F. Bergaya, Elsevier, Oxford, UK, 2015.
29. A. G. Ilgen, N. Kabengi, K. Leung, P. Ilani-Kashkouli, A. W. Knight and L. Loera, Environ. Sci.: Nano, 2021, 8, 432–443.
30. B. T. Fashina and A. G. Ilgen, Microporous Mesoporous Mater., 2026, 400, 113906.
31. H. Ohshima, in Interface Science and Technology, ed. J. Lützenkirchen, Elsevier, 2006, vol. 11, pp. 67–87.
32. M. G. Larson and F. Bengzon, The Finite Element Method: Theory, Implementation, and Applications, Springer, Berlin, Heidelberg, 2013.
33. H. P. Langtagen and K.-A. Mardal, Introduction to Numerical Methods for Variational Problems, Springer, Cham, 2019.
34. H. P. Langtagen and A. Logg, Solving PDEs in Python, Springer, Cham, 2016.
35. S.-G. Lee, K. Y. Lee, S. Y. Cho, Y. Y. Yoon and Y. Kim, Geosci. J., 2006, 10, 103–114.
36. S. Stumpf, T. Stumpf, J. Lutzenkirchen, C. Walther and T. Fanghanel, J. Colloid Interface Sci., 2008, 318, 5–14.
37. J. Lützenkirchen, R. Marsac, D. A. Kulik, T. E. Payne, Z. Xue, S. Orsetti and S. B. Haderlein, Appl. Geochem., 2015, 55, 128–137.
38. H. G. Ozcelik and M. Barisik, Phys. Chem. Chem. Phys., 2019, 21, 7576–7587.
39. J. Yang, H. Su, C. Lian, Y. Shang, H. Liu and J. Wu, Phys. Chem. Chem. Phys., 2020, 22, 15373–15380.
40. T. Sen and M. Barisik, Phys. Chem. Chem. Phys., 2018, 20, 16719–16728.
41. T. Sen and M. Barisik, Sci. Rep., 2019, 9, 137.

---