Concentration dependent interfacial chemistry of the NaOH_(aq): gibbsite interface

Abstract

Caustic conditions are often employed for dissolution of a wide variety of minerals, where ion sorption, surface diffusion, and interfacial organization impact surface reactivity. In the case of gibbsite, γ-Al(OH)₃, the chemistry at the NaOH_(aq) interface is deeply intertwined with industrial processing of aluminum, including metal production and the disposition of Al-containing wastes. To date, little is known about the structure, speciation, and dynamic behavior of gibbsite interfaces (and that of many other minerals) with NaOH_(aq)—particularly as a function of ionic strength. Yet concentration-dependent interfacial organization and dynamics are a critical starting point to develop a fundamental understanding of the factors that influence dissolution. This work reports equilibrium molecular dynamics simulations of the γ-Al(OH)₃:NaOH_(aq) interface, revealing the sorption behavior and speciation of ions from 0.5–10 M [NaOH]. As inner-sphere complexes, Na⁺ primarily coordinates to the side of the gibbsite hexagonal cavities, while OH⁻ accepts hydrogen-bonding from the surface-OH groups. The mobility of inner-sphere Na⁺ and OH⁻ ions is significantly reduced due to a strong surface affinity in comparison to previous reports of NaCl, CaCl₂, or BaCl₂ electrolytes. At high [NaOH], contact ion pairing that is observed in the bulk solution is partially disrupted upon sorption to the gibbsite surface by the individual ion–surface interactions. The molecular-scale changes to surface speciation and competition between ion–surface vs. ion–ion interactions influence surface characterization of gibbsite and potential dissolution processes, providing a valuable baseline for starting conditions needed within future reactive molecular simulations.

---

1. Introduction

Gibbsite (γ-Al(OH)₃) is the most abundant aluminum hydroxide polymorph.^1–3 It consists of octahedral Al-atoms coordinated by three hydroxyl O-atoms each above and below the Al sheet^4,5 that form hexagonal cavities in the basal plane (001) (Fig. 1); its basal surface is terminated by hydroxyls coordinated to two Al³⁺ (Al₂OH), while the (100) and (110) edge surfaces consist of both AlOH and Al₂OH terminated by –OH and –OH₂ groups.^6,7 The dissolution of gibbsite from bauxite ores (which contain gibbsite) followed by its re-precipitation in highly concentrated NaOH_(aq) are the two critical steps of the Bayer process for aluminum and alumina production.^8–10 The caustic chemistry of gibbsite is relevant to other industrial settings, for example to the treatment and disposal of highly radioactive wastes such as those at the Hanford nuclear reservation. Therein, tank sludge, primarily composed of gibbsite and boehmite, is leached by highly concentrated NaOH solutions.¹¹

Fig. 1 Top: Illustration of hexagonal cavities in the basal plane (001) of gibbsite. Bottom: Simulation setup of an aqueous NaOH solution in contact with a gibbsite slab. The O atoms of gibbsite are shown in red, Al in pink, H atoms of gibbsite in white, O atoms of water in red, H atoms of water in white, Na⁺ in green, O atom of OH⁻ in magenta, H atom of OH⁻ in cyan. The shaded areas represent the inner-sphere (IS) and outer-sphere (OS) adsorption regions, respectively.

The nature of the sorption of Na⁺ and OH⁻ from highly concentrated NaOH_(aq) on gibbsite surfaces has not been comprehensively studied, with prior work focusing upon the development of a macroscopic or mesoscopic understanding of gibbsite dissolution by NaOH_(aq)^9,12–23 Yet the availability of these ions to participate within interfacial reactions could be intrinsically dependent upon their surface speciation and adsorption characteristics. Understanding the physical sorption of Na⁺ and OH⁻ from aqueous solution to the gibbsite basal surface is a critical first step toward characterizing the full surface chemistry behavior, as it provides relevant reference configurations for exploring the large phase space of potential dissolution reactions. Given that in alkaline conditions, the dissolution rate increases with [NaOH],¹³ there must be a catalytic effect of OH⁻ on the Al-O bond breaking that will need to be investigated by computation. Of course, this catalysis cannot be approached without understanding of the structure and dynamics of the interface prior to dissolution, i.e. the concentration-dependent interfacial organization, speciation, and diffusive behavior of Na⁺, OH⁻ and Na⁺·OH⁻. With this aim in mind, molecular dynamics (MD) simulations have been used in this work to study the speciation and characteristics of NaOH_(aq) at the (001) basal surface of gibbsite (from 0.5 M to 10 M). As inner-sphere (IS) complexes, Na⁺ primarily adsorbs to the side of the gibbsite hexagonal cavities, while OH⁻ accepts hydrogen-bonding from the surface hydroxyls. The ion mobility along the perpendicular direction is greatly inhibited due to the strong surface adsorption interaction and the mobility of Na⁺ and OH⁻ ions in the inner-sphere region is clearly reduced.

2. Computational methods

Simulation model

The structural model of the gibbsite interface is illustrated in Fig. 1 and consists of a 4 × 6 × 2 supercell created from the monoclinic gibbsite unit cell of dimensions a = 8.68 Å, b = 5.08 Å, and c = 9.74 Å, with α = γ = 90°, and β = 94.54° and where the unit cell was cleaved along the (001) plane to form a gibbsite slab of dimensions 34.7 × 30.4 × 38.8 ų. The slab was put in contact with an aqueous NaOH solution, and five different NaOH concentrations were examined: 0.5, 1.0, 2.0, 5.0, and 10 M (molecular compositions and densities presented in Table S1, ESI†). Since a perfect cleavage of the basal surface was considered in the current work and that the maximum NaOH concentration is 10 M (pH = 15) the surface oxygens of the basal surface are all considered singly protonated. These are consistent with first-principles MD calculations of the pKa using the vertical energy gap method,²⁴ that predict a surface pKa value of 22 and are within 1–2 pK_a units of experiment.^24,25 As such a perfect basal surface cannot be proton active in the common pH range. Yet results from potentiometric and AFM measurements provide contradicting conclusions concerning the proton activity of basal Al₂OH groups.^6,26–30 Importantly, features that could give rise to a surface charge interpreted as proton activity²⁸ include ion pairs,^26,27 as well as the presence of defects (including steps and kinks) that reveal proton-active groups.^26,31 As such, it is important to examine and compare the concentrations of adsorbed ions as well as adsorbed contact ion pairs (CIP) and to understand their impact upon interfacial organization. The variation in ion speciation at the surface was compared to simulations of bulk NaOH solutions at the same concentrations.

Force field implementation and benchmarking

The ClayFF suite of interatomic potentials is essentially based on a nonbonded description of the metal–oxygen interactions associated with hydrated phases, which allows good transferability of the force field parameters while maintaining full flexibility of the simulation system. Furthermore, the use of a simple nonbonded interaction potential to describe metal–oxygen bonds and harmonic terms provides an efficient and accurate basis for simulating large systems.³² Here, the ClayFF potential was employed to describe the gibbsite mineral interactions, including an Al–O–H bending term³³ along with a harmonic O–H bond term and the flexible water SPC model.³⁴ ClayFF with the addition of the metal–O–H term has been fully validated for the gibbsite/water interface against AIMD,³³ and to our knowledge such validation has not been done for other force fields able to model gibbsite interfaces, such as INTERFACE FF.³⁵ As the metal–O–H term was initially associated with a Morse O–H potential but occurrences of unphysical bond breaking were observed at 10 M NaOH concentration, the original harmonic ClayFF O–H bond term was modified by selecting r₀(O–H) = 0.9572 Å instead of r₀(O–H) = 1.0 Å, resulting in a virtually identical surface H_hydroxyl–O_water RDF (related to hydrogen-bonding) compared to the Morse O–H bond term by Greathouse et al.³⁶ The Na⁺ Lennard-Jones parameters were taken from Joung et al.³⁷ and the OH⁻ parameters from Balbuena et al.³⁸ The Lorentz–Berthelot mixing rules were applied for the Lennard-Jones interactions of unlike particles. Combined with the flexible SPC model, the Na⁺ model predicts a Na⁺-O_water coordination number (CN) of 5.8 in water, within the margin of error of the experimental value obtained in dilute conditions.³⁹ The HO⁻⋯H_W coordination number of 5.8 predicted by the hydroxide model is an overestimation relative to the experimental value of 4.1,⁴⁰ however this is a common feature of 2-site point charge models.⁴¹ Although explicit polarization can improve the OH⁻ coordination number^42,43 it does not necessarily improve the general NaOH solution structure obtained from X-ray or neutron scattering⁴⁴ and is of course incompatible with the nonpolarizable ClayFF model employed for the gibbsite slab. The nonpolarizable model by Ufimtsev et al.⁴¹ is further computationally impractical because of the large number of sites involved. The electrostatic continuum correction was subsequently examined as a means to improve the OH⁻ CN. Test calculations indicated that although the H₂O⋯OH⁻ binding energies in the gas phase and HO⁻⋯H_water coordination numbers in solution were improved by scaling the charges by 75%, the Na⁺·OH⁻ CIP interaction energy and ion–surface interactions were significantly underestimated – a feature that could not be recovered by a modification of the Lennard-Jones terms. Given these observations, the original nonpolarizable OH⁻ potential was employed. Finally, the particle–particle particle–mesh (PPPM) method⁴⁵ was employed to compute the long-range electrostatic interaction with a tolerance of 10⁻⁵.

DFT calculations were performed to validate the ability of the force field to correctly describe the interactions of the Na⁺·OH⁻ CIP with the surface of gibbsite. A 3 × 4 gibbsite supercell (26.05 × 20.31 Ų) with the basal (001) surface exposed was constructed. Geometry optimization was performed on hydrated NaOH clusters of the form NaOH-(H₂O)_n (n = 0–7) placed on the (001) surface using the BAND module of AMS (Amsterdam Modeling Suite 2019) under periodic boundary conditions.⁴⁶ The PBE functional⁴⁷ with a dispersion correction through the Grimme approach^48,49 and a triple zeta basis set plus polarization (TZP) basis set us was employed. The DFT-optimized structures were then utilized to compute the binding interaction energies (ΔE) between the adsorbed NaOH(H₂O)_n hydrated complexes on the gibbsite surface. The (ΔE) is defined as:

ΔE = E_{gibbsite/NaOH–(H2O)n} − E_gibbsite − E_{NaOH–(H2O)n} (1)

where E_{gibbsite/NaOH–(H2O)n}, E_gibbsite, E_{NaOH–(H2O)n} are total electronic energies of the complexes between the gibbsite (001) surface and hydrated Na⁺·OH⁻ CIPs, of the gibbsite (001) surface, and of the hydrated Na⁺·OH⁻ CIPs, respectively. Fig. 2 shows the comparison of the DFT energies with the ones using the ClayFF force field (single point energies from the DFT-optimized structures), which presents good agreement with the DFT computation (RMS error of 3.83 kcal mol⁻¹). This consistency confirms that the classical force field used in this work can provide a reasonable description of the interfacial interactions between the aqueous Na⁺, OH⁻ ions and the gibbsite surface.

Fig. 2 Comparison of the DFT and ClayFF energies of formation for various hydrated NaOH clusters on the gibbsite basal surface and corresponding geometries (distances in Å) for select points. Geometry optimizations were done at the DFT level.

Simulation and equilibration protocol

All simulations were performed using LAMMPS⁵⁰ with a time step of 1 fs. The volume of the aqueous phase was first equilibrated using the NPT ensemble according to the equations of Shinoda et al.⁵¹ at 300 K and 1 atm along the z-axis, resulting in a slit pore of length ∼70 Å along this dimension, ensuring the absence of confinement. Within this stage of the equilibration, the integrity of the mineral slab was preserved by constraining the velocity of its Al- and O-atoms to zero, which still allowed fluctuations of the volume of the aqueous phase. In the subsequent NVT run, two 1 ns temperature cycles from 300 K to 500 K and back were performed. At the final temperature of 300 K the mineral slab constraint was released; during this phase, the radial distribution functions (RDFs) for all O–H and Na–O pairs were monitored and the system was considered equilibrated when they reached convergence (after 30–40 ns). A final 1 ns equilibration of the density was done in the NPT ensemble at 300 K. All the properties of the system were computed based on a subsequent 40 ns production phase in the NVT ensemble.

Umbrella sampling

Umbrella sampling was employed to investigate the free energy profiles of ion adsorption to the basal gibbsite surface. Note that the inclusion of multiple ions would requires tens of ns of ensemble averaging over 75 umbrella sampling windows. Separate MD simulations were performed using the gibbsite slab, water, and a single ion present (either Na⁺ or OH⁻). The collective variable was defined as the z-coordinate of the ion relative to the plane formed by the outermost mineral O-atoms (set to z = 0). One hundred ps umbrella sampling windows were separated by 0.1 Å with k = 125 kcal mol Å⁻² for 1 Å < z < 6 Å and 0.2 Å with k = 30 kcal mol Å⁻² for 6 Å < z < 11 Å. To prevent the Na⁺ from moving laterally across the hexagonal surface cavity, the distance between the ion and the line orthogonal to the surface and passing through the site was constrained to zero. However, the force constant of the restraint k = 1 kcal mol Å⁻², was much weaker than the restraint of the biasing collective variable, ensuring that this secondary, lateral restraint had a negligible influence on the potential of mean force (PMF). As the gibbsite surface is neutral, the electrostatic interaction between every ion and the surface is zero and the interaction at z = 10 Å is close to zero.

3. Results and discussion

Adsorption characteristics

We first consider the general adsorption characteristics of NaOH_(aq) moving from the bulk solution to the surface of gibbsite, as shown in Fig. 3via the number density profiles for Na⁺, OH⁻ and the contact ion pairs (CIP) between Na⁺ and OH⁻, abbreviated Na⁺·OH⁻. Contact ion pairs are defined as those ion pairs with a direct ion–ion interaction and no intervening solvent molecules between them. In the trajectory analysis, one CIP is counted every time the Na⁺-OH⁻ distance is less than 2.85 Å (distance corresponding to the minimum after the first Na-O(OH⁻) RDF peak). For clarification, a single ion can have multiple direct contacts with an ion of opposite charge, i.e. several CIPs can be associated with a single ion. The coordinates of the Na⁺·OH⁻ CIP are chosen as the center of the Na-O axis.

Fig. 3 Concentration profiles (averaged over both interfaces within the simulation) of Na⁺, OH⁻, Na⁺·OH⁻ CIPs (a–f) the number density profiles of water O_W (solid) and H_W (dashed). The vertical dashed lines indicate the boundaries between inner sphere (IS) and outer-sphere (OS) regions (at 3.3 Å) and between the OS and the rest of the solution (at 5.3 Å). The value of the z-coordinate of the outermost layer of O-atoms on both surfaces was defined as 0, and the density profiles averaged over the two surfaces of the mineral slab.

The density profiles further allow the definition of inner-sphere (IS) sorbed ions as those whose distance is z < 3.3 Å, meaning their direct interaction with surface-atoms, while outer-sphere (OS) sorbed ions reside within 3.3 < z < 5.3 Å (having intervening H₂O between the ion and the surface atoms). Even if a third, less structured ion layer exists (visible for both ions at 10 M in Fig. 3e), the z > 5.3 Å portion of the solution is considered “bulk” for the sake of convenience. The density profiles clearly show the IS and OS species as two distinct layers from z = 0 Å, respectively. The distribution of Na⁺ within the IS layer consists of three peaks at 1.4 Å, 1.8 Å and 2.3 Å that correspond to different adsorption sites (vide infra). The distribution of OH⁻ within the IS layer consists of one peak at 2.6 Å. The Na⁺·OH⁻ CIPs are mostly located within the IS layer at a distance of 2.2 Å. The O_W and H_W number density profiles are shown in Fig. S1 (ESI†). The most prominent effect therein is the decrease of the height of the first H_W peak (z < 2 Å) associated with H_WO–H_W⋯O_gibbsite) hydrogen bonding, presumably caused by hindrance due to the increased number of Na⁺ coordinated to the surface O-atoms.

From 0.5 M to 2 M [NaOH] the surface coverages of Na⁺ and OH⁻ continuously increase, and saturation is reached 5 M with a coverage of ∼7–8 ions nm⁻². Interesting comparisons can be made to a similar study of NaCl_(aq) sorption by Ho et al.⁷ within the same concentration range and using similar quality force fields. Although the current work reports a surface coverage of Na⁺ and OH⁻ at 1 M that is greater than 2 per nm², in the case of Na⁺ and Cl⁻ from NaCl_(aq) surface coverage values of 0.4 and 0.7 ions nm⁻² were observed, respectively. At 2 M NaCl_(aq), 1.1 Na⁺ nm⁻² were observed to be sorbed to the basal surface of gibbsite as compared to our 4 Na⁺ nm⁻² from 2 M NaOH_(aq) (note that the Cl⁻ surface coverage was not mentioned by the authors at this concentration). Excluding the possible effects of the box size, the force field, and the equilibration procedure, which are likely minor, these data indicate a likely smaller affinity of Cl⁻ for the surface than OH⁻, which results in less Na⁺ adsorbed due to charge compensation. Ho et al. found even smaller adsorption values for Ca²⁺ (in CaCl_2(aq)) and Ba²⁺ (in BaCl_2(aq)), where 0.8–0.9 ions nm⁻² were sorbed from 2 M electrolyte solutions. The surface coverage of Na⁺ in the IS layer is systematically larger than that of the OS layer, in agreement with observations regarding NaCl_(aq) on the gibbsite basal surface. However, the IS/OS ratio for Na⁺ decreases with increasing concentration.

Very different trends are observed for OH⁻ sorption, as similar amounts of anions are found in the IS and the OS layers irrespective of [NaOH]. The lower degree of adsorption of OH⁻ in the IS layer compared to Na⁺ indicates a cumulatively weaker interaction of OH⁻ with the hydroxyls caused by hydrogen-bonding and coulombic forces than the coulombic interaction of Na⁺ with surface O-atoms. At higher concentrations, the balance between Na⁺ and OH⁻ in the inner sphere tends to equalize as a result of saturation being reached, even if a slight imbalance in favor of Na⁺ remains. Interestingly, the number of OH⁻ in the OS slightly overcompensates this imbalance at 5 M and 10 M, though an approximately neutral net charge is obtained when the third layer is accounted for (not shown here). Finally, H₂O surface coverage is very modestly affected by ion sorption, with a decrease in surface coverage of ∼1 H₂O nm⁻² (or a 10% decrease) between 0.5 M and 10 M in the first water layer, and virtually no change in the second layer (Fig. 4). Note that the coordinate perpendicular to the surface was allowed to relax during the equilibration, and thus the small decrease of the number of H₂O at the interface compared to the large increase in the number of sorbed ions between 0.5 M and 10 M may be explained by the overall contraction of the H₂O⋯H₂O distances in the volume unoccupied by the ions.

Fig. 4 Surface adsorption coverage for inner-sphere (IS) and outer-sphere (OS) Na⁺ and OH⁻ ions as well as 1st water layer (W_L1) and 2nd layer (W_L2) as a function of ion concentration. The value indicated for the water surface coverage is relative to its value at 0.5 M; the absolute values at 0.5 M are 9.32 and 8.75 per nm² for W_L1 and W_L2, respectively.

As indicated from Fig. 5a, between 0.5 and 2 M NaOH the overall Na⁺·OH⁻ contact ion pairing is slightly increased at the gibbsite surface relative to bulk NaOH_(aq) simulations, suggesting that the surface facilitates ion pairing at moderate concentrations. In contrast, under higher [NaOH] contact ion pairing is strongly reduced at the solid/liquid interface where the IS and OS layers are saturated (Fig. 4, 5 and 10 M). Across the concentration range, there is competition amongst the inter-particle interactions associated with adsorption of individual ions vs. contact ion pairing. Fig. 5b shows the number of CIPs between an ion in the IS and an ion of opposite charge either in the IS (IS–IS) or in the OS (IS–OS), divided by the total number of Na⁺ in the solution. This increases up to 2 M NaOH, then decreases with higher concentrations, indicating that ion saturation at the surface causes an energetic penalty for CIP interactions of sorbed species.

Fig. 5 Ratio between the number of contact ion pairs and the number of Na⁺ (or OH⁻) as a function of concentration, where N(Na⁺·OH⁻) is the total number of contacts between Na⁺ and OH⁻. (a) The number of CIPs normalized by the number of Na⁺ in the bulk NaOH_(aq) model and gibbsite/NaOH_(aq). Note that the “interface” values reflect only the CIPs and Na⁺ in the IS and OS. (b) The number of CIPs in different inner-sphere and outer-sphere regions, normalized by the total number of Na⁺. In this case the X–Y notation reflects contact ion pairing between an ion in the X layer and an ion in the Y layer with X, Y indicating IS, OS, or B – where IS = inner-sphere; OS = outer-sphere; B = rest of the solution (z > 5.3 Å).

Ion adsorption sites

The planar (xy) density distributions were then examined to elucidate the adsorption sites of Na⁺ and OH⁻. Fig. 6 depicts the preferential adsorption sites along the xy plane for inner-sphere adsorbed Na⁺ and OH⁻. Two adsorption sites are observed for Na⁺ (Fig. 6a): #1 above the center of a hexagonal cavity formed by the Al-atoms and coordinating with three surface O-atoms; #2 above the Al-atoms and coordinated by three surface O-atoms. On average, 0.76 Na⁺ and 0.60 OH⁻ per hexagonal cavity are adsorbed as IS complexes at saturation (5 M). For OH⁻, the surface adsorption site is primarily located near the surface H-atoms belonging to the surface hydroxyls that point toward the solution (Fig. 6b) through a hydrogen bond (HB) interaction. Fig. 6c shows the planar density distributions of H₂O in the 5 M NaOH_(aq). It is clear that most H₂O are distributed on the edge of the hexagonal cavities on the surface (001), just like the inner-sphere bound OH⁻. The anticipated effect this may have upon an estimated surface charge is discussed in the ESI.†

Fig. 6 Planar density distributions of the inner-sphere (IS) layer for the 5 M NaOH solution. (a) Na⁺ (b) OH⁻ (c) O_W. (d) Planar density distributions of O_W of the first H₂O layer for the pure water system. In (c) and (d), the balls in beige represent the Al atoms.

Ion hydration and interfacial hydrogen-bonding

RDFs were used to characterize the structure of the ion hydration shells where IS complexation significantly decreases hydration numbers in the first coordination shell (Fig. 7). For example, at 2 M [NaOH] the ion-O_W coordination number (CN) in the IS region is reduced by 1.8 and 2.0 for OH⁻ and Na⁺, respectively, compared to the bulk. We note that the bulk CN(Na⁺–O_W) of 5.1 is smaller than the value mentioned in the Methods section as a result of significant ion-pairing with OH⁻ that obviously does not occur at high dilution.

Fig. 7 RDFs between ions and water molecule O-atoms (O_W) within the typical 2 M NaOH solution when the ions exist as (a) inner-sphere and (b) outer-sphere complexes, or (c) when they are in the bulk. The corresponding configuration snapshots are presented in (d), where the immediately surrounding H₂O are indicated with ion-water distances reported for reference.

As in the bulk, the gibbsite basal surface features intralayer hydrogen-bonding between hydroxyl groups in the basal plane and out of the basal plane,⁵² the latter representing 58% of the OH groups (Fig. S2, ESI†). Hydroxyl groups also form HBs with adsorbed H₂O (both as HB donors and acceptors) and OH⁻ (only donated). The evolution of the average number of HBs as a function of concentration is represented in Fig. 8c. The corresponding HB autocorrelation function⁵³ is defined as

where h(t) = 1 if the tagged water pair maintains a HB from time 0 to time t; otherwise, h(t) = 0. Fig. 8a and b present the S_HB(t) curves for the electrolyte-surface HBs and OH⁻-surface HBs as a function of NaOH concentration, while the HB lifetimes τ^HB_S determined from a triexponential fit of S_HB(t) are presented in Fig. 8d (Tables S2, S3 and Fig. S3, ESI†). As shown in Fig. 8a, the higher [NaOH] concentration, the slower the decay of S_HB(t), which indicates a strengthening of the surface⋯water HBs. The number of HBs (N(HB)) per HB type donated or accepted by the surface hydroxyls is shown in Fig. 8c. As far as hydroxyl groups are concerned, the number of HBs donated to neighboring hydroxyl groups is 0.45–0.49 depending on [NaOH]—in agreement with the fact that 42% of the hydroxyls are parallel to the plane (Fig. S2, ESI†). Most of the rest of the hydroxyls donate their HB to H₂O or to OH⁻, and a minority does not donate HBs as the total N(HB) does not reach 1.00 even at 10 M (where N(HB)_total = 0.95, Fig. 8c). As [NaOH] increases, more surface -OH groups donate HBs to the hydroxide anion, with an increase in N(HB) from ∼0 (0.5 M) up to 0.12 (5 and 10 M, Fig. 8c). The number of HBs donated to H₂O is however nearly constant (0.46–0.49), is agreement with the very minor release of H₂O when [NaOH] increases (Fig. 2). The seemingly contradictory decrease in the number of HBs accepted by the surface -OH from the H₂O from 0.28 to 0.15 is likely due to the change in water orientation caused by OS sorption. The increase in the number of surface-hydroxide anion HBs is accompanied by a strong increase in their lifetime from ∼2 to ∼20 ps (Fig. 8d), a feature that is related to a strengthening of the network of interactions that include hydrogen bonding, ion–dipole, and ion–ion interactions. In contrast, the increase in the lifetime of HBs between the surface and sorbed H₂O is much more modest, from 0.5 to 2 ps. This indicates that water diffusion is not considerably slowed down by the creation of the network of interactions that involve sorbed OH⁻ and Na⁺.

Fig. 8 The S_HB(t) curves of the water-surface HB (a) and the hydroxyl-surface HBs (b) for different NaOH concentrations. Inserts show the corresponding configurations. (c) Average number of hydrogen bonds N(HB) between different atom pairs, based on the H···O distance. The cutoff (ranging from 2.27 to 2.40 Å) corresponds to the first minimum of the RDF of the corresponding pair after the first intermolecular peak. Within the notation “g” denotes the gibbsite (001) surface; “w” denotes H₂O. (d) Average hydrogen bond lifetime τ^HB_S.

It is interesting to consider whether the observed modifications to interfacial water structure would be detectable by experiment, for example in atomic force microscopy (AFM). Nakouzi et al.⁵⁴ recently utilized fast force mapping AFM, complemented by classical molecular dynamics rare event simulations, to observe the molecular structure of the interface of boehmite (γ-AlO(OH)) with pH 11 water, a surface related to the α-Al(OH)₃:water studied here. That work found that the molecular structure of the interface measured by AFM was strongly sensitive to the local interfacial water density, which likely reflects a templating effect of the mineral surface structure upon that of interfacial water. If true, it suggests that surface charge at even more elevated pH would not affect the average interfacial water structure until the surface charge density is sufficiently high to affect a comprehensive change in the mineral–water interface structure (e.g., several charged sites per nm²). This in turn might create an unstable surface, leading to dissolution. Regardless, since interfacial water density is not significantly modified by ion sorption at saturation, and given that an ideal gibbsite basal surface site is not proton active,⁵⁵ the aforementioned data on boehmite indicates that the basal gibbsite interface structure measured by AFM spectroscopy would not vary significantly at high pH as long as the density of step edges remains small.

Adsorption energetics

As shown in the Fig. 9, the free energy profiles for ion sorption (obtained from potential of mean force simulations) feature two separated free energy minima, corresponding to the IS and OS adsorption positions in the corresponding density profiles (Fig. 3). In Fig. 9a, the Na⁺ on site #2 (side of the hexagonal cavity) has a smaller energy minimum in IS, which agrees with the density profiles. For Na⁺ and OH⁻ the first PMF minimum (corresponding to IS adsorption) is respectively at 1.9 Å and at 2.7 Å, consistent with the observation that the first adsorption layer of Na⁺ is closer to the gibbsite surface (001) than OH⁻. In the PMF profiles, the first free energy minimum of the Na⁺ on site #2 is approximately 2.3 kcal mol⁻¹ lower in energy than that of the OH⁻ (−2.0 kcal mol⁻¹), demonstrating that Na⁺ has a stronger adsorption stability. The PMF in Fig. 9a reveals that Na⁺ adsorption to the side of the hexagonal cavity consists of three sites that correspond to three minima observed at z-values between 1.9 and 2.4 Å depending on whether Na⁺ is directly above the Al-atom or slightly offset (Table 1).

Fig. 9 PMF of Na⁺ (a) and OH⁻ (b) adsorbing on the gibbsite (001) surface. Two surface sites for inner-sphere Na⁺ and one for OH⁻ are shown above the image. The blue and red vertical dashed lines, respectively, indicate the separation between the inner-sphere layer and the outer-sphere layer and the separation between outer-sphere layer and the bulk phase. Given the absence of long-range electrostatics and the fact that the density profiles at 0.5 M and for z > 6.5 Å are almost flat, the interaction energy is considered to be equal to zero at z = 11 Å.

Table 1 Equilibrium constants (K) for the interconversion of Na⁺ and OH⁻ between the inner-sphere (IS) and outer-sphere (OS) states and between the OS and bulk state. Given that adsorption to site #2 is more stable than site #1 by ∼2 kcal mol⁻¹ in Fig. 8, the corresponding PMF was selected for the calculation of K for Na⁺

		Na^+	OH^−
IS → OS	K = [OS]/[IS]	0.21	0.34
OS → bulk	K = [dissolved]/[OS]	0.88	0.83

---

However, the difference in stability and the energy barriers between those sites is less than kT (0.6 kcal mol⁻¹) and thus inter-site hops are not rare events and thermodynamic distinction between these sites is not meaningful. Less importantly, for the PMF to site #1, near z = 2.4 Å (the third IS minimum) some incursions near the center of the cavity occur due to the weak lateral restraint.

The adsorption equilibrium constants (K) for the NaOH near the gibbsite surface, were calculated according to:

with “#” referring to a transition state. We take the leftmost defined reaction coordinate (∼1.3–1.6 Å) as the left boundary for IS and the right boundary is at 9.1 Å for the bulk based upon the transition from the OS region and the bulk shown in Fig. 9. These equilibrium constants are in apparent contradiction with the results from Fig. 3, e.g., the IS → OS transformation for OH⁻ corresponds to an equilibrium constant of 0.34, yet the ratio between the OS concentration and the IS concentration at 10 M is 1.22. This is due to the fact that at high concentration, both the adsorption equilibrium and the contact ion pairing equilibrium need to be taken into account. However, the reported value for the adsorption of Na⁺ on gibbsite from dilute NaNO_3(aq) is K = 0.4, which is an intermediate value between that of the IS → bulk transformation (0.21 × 0.88 = 0.18) and the OS → bulk transformation (0.88).⁵⁶

Dynamic behavior of adsorbed ions

In order to elucidate the dynamic behavior of the Na⁺ and OH⁻ ions within the interfacial region, their in-plane mean square displacement along the xy plane (MSDxy) and the out-of-plane mean square displacement along the z-direction (MSDz) have been computed. The MSDs are defined as:

MSD_xy = 〈Δx(t)²〉 + Δy(t)²〉 ∝ t^k (4)

The parameter k is an indicator of the diffusion mode. For Fickian diffusion k = 1, while k = 0.5 represents a single-file diffusion mode. The intermediate states k > 1 and k < 1 correspond to superdiffusion and subdiffusion^57–59 regimes, respectively.

The simulated MSD_xy and MSD_z for both ions at 1 M NaOH within the interfacial IS and OS regions are shown in Fig. 10 and those for other concentrations in Fig. S5 (ESI†). In general, Na⁺ and OH⁻ ions show identical concentration dependent trends in diffusion behavior. In the IS, the slopes of all the logarithmic MSD curve are nearly 0.5 in the period of 1–100 ps. Beyond that, the diffusion of ions exhibits Fickian diffusion (k ≈ 1). The ion mobility it evaluated using the self-diffusion coefficient D_s with the Einstein relation.⁶⁰ The in-plane two-dimensional diffusion coefficient D_s can be calculated from the slope of the in-plane mean squared displacements (MSDxy) as,

where t is the simulation time. We used the MSD_xy during the 700–1000 ps to calculate the self-diffusion coefficients of 1 M system. The calculation results of self-diffusion coefficients D_s are presented in Fig. 10 and the diffusion coefficients for all concentration conditions are listed in Table S4 (ESI†). Unsurprisingly, the self-diffusion coefficients of the Na⁺ and OH⁻ ions in IS and OS are smaller than these of the bulk occurring in the order IS < OS < bulk as a result of the strong surface confining interaction. Further, the in-plane self-diffusion coefficients are larger than the out-of-plane ones. A drastic decrease of the self-diffusion coefficients for Na⁺ (relative to OH⁻) is observed for those ions in the inner-sphere layer, and thus Na⁺ mobility is more inhibited than OH⁻. In addition, the self-diffusion coefficients of Na⁺ and OH⁻ ions in the OS are larger than that in IS, which is consistent with the fact that the surface interaction becomes weaker in the OS region. The obviously reduced mobility of Na⁺ and OH⁻ ions in the IS region is in good agreement with their stronger adsorption interaction (Fig. 9). As shown in Table S4, (ESI†) the self-diffusion coefficients of the Na⁺ and OH⁻ ions in the IS and OS layers become smaller with increasing concentrations. This can be attributed to the enhanced steric hindrance effect under higher concentration and the observed increase in the HB lifetimes. Ion mobility necessitates the ability of other sorbed ions and H₂O to rearrange on the surface, and as the concentration increases ion mobility will require more displacements of the nearby strongly adsorbed ions and the subsequent breakage of intermolecular interactions.

Fig. 10 MSDs in log–log scale for ions in the interface region (IS and OS) in the 1 M system. For both Na⁺ and OH⁻, the MSDs results include the in-plane MSD along the xy plane and the out-of-plane MSD along the z direction (the analogous self-diffusion coefficients D_s are presented in the same corresponding color). Dashed lines specifying the plot slopes are provided for comparison.

4. Conclusions

MD simulations have been employed to investigate the concentration dependent interfacial organization and ion adsorption behavior of NaOH_(aq) on the (001) basal surface of gibbsite. Both Na⁺ and OH⁻ are observed to adsorb primarily as inner-sphere species, with Na⁺ residing either within or on the edge of a hexagonal cavity of the surface O-atoms. For OH⁻, the surface adsorption sites are located above the gibbsite hydroxyls from which they accept hydrogen bonds. The adsorption capacity of Na⁺ and OH⁻ ions on the gibbsite surface is higher than those previously reported for aqueous NaCl, CaCl₂, or BaCl₂ electrolytes. This behavior could be attributed to the enhanced affinity of OH⁻ ions with the gibbsite surface (and consequently Na⁺ due to contact ion pairing). Up to 2 M NaOH, the surface is found to enhance contact ion pairing within the IS region, however it strongly disrupts the formation of interfacially bound CIPs above 2 M NaOH. The interfacial diffusion properties of the ions indicate that ion mobility is greatly inhibited in the inner-sphere region, while the motion of ions in the outer-sphere region is more unrestrained due to the weaker surface interaction.

Note that although dissolution may occur primarily at edges (including step edges), at very high pH (typically pH > 13) the (100) surface (which is a predominant edge surface according to Hiemstra et al.⁶) only consist of OH groups,⁵⁵ as in the basal (001) surface considered in this work. Such edges may be quantitatively modulated by the O–Al coordination, surface charge, among other aspects, yet the concentration dependent behavior of ion sorption on the (001) including the modulation of contact ion paired species at the surface, in combination with the changes to surface diffusion, provide valuable new insight toward ongoing and future studies of gibbsite dissolution under NaOH_(aq) conditions relevant to Al industrial processing and metal production.⁶ More generally, this work expands a fundamental understanding of ion-sorption under highly caustic conditions, which is relevant to a wide variety of processes that leverage NaOH (or other hydroxide electrolytes) for dissolution.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported by the Interfacial Dynamics in Radioactive Environments and Materials (IDREAM), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences. Molecular simulations used resources from the Center for Institutional Research Computing at Washington State University.

References

1. R. Demichelis, Y. Noel, P. Ugliengo, C. M. Zicovich-Wilson and R. Dovesi, J. Phys. Chem. C, 2011, 115, 13107–13134.
2. S. Nagendran, G. Periyasamy and P. V. Kamath, Z. Anorg. Allg. Chem., 2015, 641, 2396–2403.
3. M. F. Peintinger, M. J. Kratz and T. Bredow, J. Mater. Chem. A, 2014, 2, 13143–13158.
4. J. T. Kloprogge, L. V. Duong, B. J. Wood and R. L. Frost, J. Colloid Interface Sci., 2006, 296, 572–576.
5. S.-L. Wang and C. T. Johnston, Am. Miner., 2000, 85, 739–744.
6. T. Hiemstra, H. Yong and W. H. Van Riemsdijk, Langmuir, 1999, 15, 5942–5955.
7. T. A. Ho, J. A. Greathouse, A. S. Lee and L. J. Criscenti, Langmuir, 2018, 34, 5926–5934.
8. A. R. Hind, S. K. Bhargava and S. C. Grocott, Colloids Surf., A, 1999, 146, 359–374.
9. J. A. M. Pereira, M. Schwaab, E. Dell’Oro, J. C. Pinto, J. L. F. Monteiro and C. A. Henriques, Hydrometallurgy, 2009, 96, 6–13.
10. H. B. Yang, X. L. Pan, H. Y. Yu, G. F. Tu and J. M. Sun, Nonferrous Met. Soc., 2015, 25, 4151–4159.
11. J. G. Reynolds, J. K. McCoskey and D. L. Herting, Ind. Eng. Chem. Res., 2016, 55, 5465–5473.
12. R. F. Scotford and J. R. Glastonbury, J. Chem. Eng., 1971, 49, 611–616.
13. R. F. Scotford and J. R. Glastonbury, Can. J. Chem. Eng., 1972, 50, 754–758.
14. A. Bauer and G. Berger, Appl. Geochem., 1998, 13, 905–916.
15. P. Benezeth, D. A. Palmer and D. J. Wesolowski, Geochim. Cosmochim. Acta, 2008, 72, 2429–2453.
16. J. A. Chermak, Miner, 1992, 40, 650–658.
17. J. A. Chermak, Miner, 1993, 41, 365–372.
18. T. R. Graham, M. Dembowski, E. Martinez-Baez, X. Zhang, N. R. Jaegers, J. Hu, M. S. Gruszkiewicz, H.-W. Wang, A. G. Stack, M. E. Bowden, C. H. Delegard, G. K. Schenter, A. E. Clark, S. B. Clark, A. R. Felmy, K. M. Rosso and C. I. Pearce, Inorg. Chem., 2018, 57, 11864–11873.
19. H. Grenman, T. Salmi, D. Y. Murzin and J. Addai-Mensah, Eng. Chem. Res., 2010, 49, 2600–2607.
20. C. D. Peskleway, G. S. Henderson and F. J. Wicks, Am. Mineral., 2003, 88, 18–26.
21. N. Brown, J. Cryst. Grow., 1972, 12, 39–45.
22. Z. Shen, S. N. Kerisit, A. G. Stack and K. M. Rosso, J. Phys. Chem. Lett., 2018, 9, 1809–1814.
23. R. Schliemann and S. V. Churakov, Geochim. Cosmochim. Acta, 2021, 293, 438–460.
24. M. Sulpizi and M. Sprik, Phys. Chem. Chem. Phys., 2008, 10, 5238–5249.
25. T. De Meyer, B. Ensing, S. M. J. Rogge, K. De Clerck, E. J. Meijer and V. Van Speybroeck, ChemPhysChem, 2016, 17, 3447–3459.
26. J. Rosenqvist, P. Persson and S. Sjöberg, Langmuir, 2002, 18, 4598–4604.
27. M.-C. Jodin, F. Gaboriaud and B. Humbert, J. Colloid Interface Sci., 2005, 287, 581.
28. Y. Gan and G. V. Franks, 6.
29. S. Mitchell, Theses and Dissertations.
30. B. R. Bickmore, C. J. Tadanier, K. M. Rosso, W. D. Monn and D. L. Eggett, Geochim. Cosmochim. Acta, 2004, 68, 2025–2042.
31. A. Klaassen, F. Liu, D. van den Ende, F. Mugele and I. Siretanu, Nanoscale, 2017, 9, 4721–4729.
32. R. T. Cygan, J.-J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255–1266.
33. M. Pouvreau, J. A. Greathouse, R. T. Cygan and A. G. Kalinichev, J. Phys. Chem. C, 2017, 121, 14757–14771.
34. O. Teleman, B. Jönsson and S. Engström, Mol. Phys., 1987, 60, 193–203.
35. H. Heinz, T.-J. Lin, R. Kishore Mishra and F. S. Emami, Langmuir, 2013, 29, 1754–1765.
36. J. A. Greathouse, J. S. Durkin, J. P. Larentzos and R. T. Cygan, J. Chem. Phys., 2009, 130, 134713.
37. I. S. Joung and T. E. Cheatham, J. Phys. Chem. B, 2008, 112, 9020–9041.
38. P. B. Balbuena, K. P. Johnston and P. J. Rossky, J. Phys. Chem., 1996, 100, 2706–2715.
39. M. Galib, M. D. Baer, L. B. Skinner, C. J. Mundy, T. Huthwelker, G. K. Schenter, C. J. Benmore, N. Govind and J. L. Fulton, J. Chem. Phys., 2017, 146, 084504.
40. S. E. McLain, S. Imberti, A. K. Soper, A. Botti, F. Bruni and M. A. Ricci, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 094201.
41. I. S. Ufimtsev, A. G. Kalinichev, T. J. Martinez and R. J. Kirkpatrick, Chem. Phys. Lett., 2007, 442, 128–133.
42. R. Vácha, T. Megyes, I. Bakó, L. Pusztai and P. Jungwirth, J. Phys. Chem. A, 2009, 113, 4022–4027.
43. J. S. Hub, M. G. Wolf, C. Caleman, P. J. van Maaren, G. Groenhof and D. van der Spoel, Chem. Sci., 2014, 5, 1745.
44. D. Semrouni, H.-W. Wang, S. B. Clark, C. I. Pearce, K. Page, G. Schenter, D. J. Wesolowski, A. G. Stack and A. E. Clark, Phys. Chem. Chem. Phys., 2019, 21, 6828–6838.
45. J. W. Eastwood, R. W. Hockney and D. N. Lawrence, Comput. Phys. Commun., 1980, 19, 215–261.
46. G. Velde, F. M. Bickelhaupt, E. J. Baerends, C. F. Guerra, S. J. A. V. Gisbergen, J. G. Snijders and T. Ziegler, J. Comput. Chem., 2001, 22, 931–967.
47. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
48. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
49. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
50. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
51. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
52. R. Demichelis, B. Civalleri, Y. Noel, A. Meyer and R. Dovesi, Chem. Phys. Lett., 2008, 465, 220–225.
53. A. Chandra, J. Phys. Chem. B, 2003, 107, 3899–3906.
54. E. Nakouzi, A. G. Stack, S. Kerisit, B. A. Legg, C. J. Mundy, G. K. Schenter, J. Chun and J. J. De Yoreo, J. Phys. Chem. C, 2021, 125, 1282–1291.
55. X. Liu, J. Cheng, M. Sprik, X. Lu and R. Wang, Geochim. Cosmochim. Acta, 2013, 120, 487–495.
56. R. Weerasooriya, H. J. Tobschall, H. K. D. K. Wijesekara, E. K. I. A. U. K. Arachchige and K. A. S. Pathirathne, Adsorp. Gibbsite, Chemosp., 2003, 51, 1001–1013.
57. A. Striolo, Nano Lett., 2006, 6, 633–639.
58. J. P. Bouchaud, A. Georges, J. Koplik, A. Provata and S. Redner, Phys. Rev. Lett., 1990, 64, 2503–2506.
59. H. F. Ye, H. W. Zhang, Y. G. Zheng and Z. Q. Zhang, Nanofluid, 2011, 10, 1359–1364.
60. E. H. P. Liu and B. Berne, J. Phys. Chem. B, 2004, 108, 6595.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp01997c

‡ These authors contributed equally to this work.

---