Investigating the chemical-mechanical of silica requires a force field which both correctly estimates the mechanical properties of silica and the water/electrolyte interaction with the surface. Very few MD force fields exhibit both capabilities. These include one (Mahadevan and Garofalini, 2007) which is an adaptation of a dissociative water force field (Guillot and Guissani, 2001) and the other, a bond-order based force field (ReaxFF) (Van Duin et al., 2003; Fogarty et al., 2010; Pitman and Van Duin, 2012; Yeon and van Duin, 2015). ReaxFF allows for the dissociation and creation of chemical bonds and redistribution of charges through environmentally dependent charge equilibrium methods. Previous studies have used the Yeon (Yeon and van Duin, 2015) and the Pitman (Pitman and Van Duin, 2012) parametrizations of ReaxFF to investigate silica surface energies (Rimsza et al., 2017a; Yu et al., 2018) and the Fogarty (Fogarty et al., 2010) parameterization to calculate mechanical properties (Yu et al., 2016). For this work, the Yeon version of ReaxFF was selected because it is a reparametrization of the Fogarty version that improves the energetics of water-surface interactions (Rimsza J. M. et al., 2016), which are critical to chemical-mechanical fracture. Additionally, our previous investigations of silica fracture toughness (K_IC) process zone size, and radius of curvature (Rimsza et al., 2017b) have been consistent with experiment. Some concerns have been raised surrounding the over estimation of bulk modulus of α-quartz, reported as 602 GPa by Yu et al. (2018). Elastic constants for the Yeon ReaxFF parametrization were calculated as E = 76.1 ± 0.3 GPa, K = 55.2 ± 5.2 GPa, μ = 21.8 ± 0.8 GPa, and ν = 0.32 ± 0.02 for the elastic modulus, bulk modulus, shear modulus, and Poisson's ratio, respectively. Reported error is from variation in strain rates. These values do not compare exactly with experimental results (E = 73.1 GPa, K = 36.7 GPa, μ = 31.2 GPa, and ν = 0.17) (Varshneya, 2013) but are within an order of magnitude. The elastic constant that is of primary concern is E, since the loading methodology (discussed below) is focused on mode I or pure tension. The overestimation of G in the Yeon parametrization (55.2 ± 5.2 GPa from simulation and 36.7 GPa from experiment) would be a larger concern for studies which focused on the analysis of silica in shear and compression. All simulations were performed using the LAMMPS MD code with the USER-REAXC package with the Yeon and van Duin Si/O/H ReaxFF parameterization (Aktulga et al., 2012; Yeon and van Duin, 2015) as the force field (included in Supporting Information).

To investigate crack propagation in silica, thin, quasi-2D simulation cells (140 × 140 × 28 Å) were created through a melt-and-quench procedure. Three replicas of the silica glass were generated to enable statistical averaging. The initial cristobalite configurations contained 38,400 atoms in a 2:1 oxygen to silicon ratio. These systems were heated to 4,000 K, then held for 100 ps to melt the crystal structure under periodic boundary conditions. They were then cooled from 4,000 to 300 K at a rate of 5 K/ps, which is reported to accurately replicate bond distance and angles in glassy silica (Lane, 2015). Densities were controlled during cooling at the measured density of silica, 2.2 g/cm³ (Mozzi and Warren, 1969), through a canonical (NVT) ensemble which controls the number (N) of atoms, volume (V), and temperature (T) of the simulation with a damping coefficient of 100 time steps and a time-step of 0.5 fs. After the systems had reached 300 K, the density fluctuated in simulations that alternated isothermal-isobaric dynamics with energy minimization to achieve low energy configurations. Final silica structures had a density of 2.187 ± 0.001 g/cm³.

First a slit crack has been introduced into the silica on a half-plane by deleting bonds between atoms on opposite sides of the half-plane, then the crack was mechanically loaded by displacement of far-field atoms in a region outside a cylinder centered on the crack tip (Figure 1). This displacement field was taken from the linear elastic fracture mechanics (LEFM) continuum solution for a semi-infinite slit crack in mode I (tension) loading that is spatially varying and parameterized by K_I (Lawn and Wilshaw, 1993, Equation 2.11). Prior to loading, the in-plane (x and y) boundary conditions were changed to non-periodic since periodicity is lost during loading of the system (Figure 1). Via the prescribed displacement field, a pre-existing slit crack on a halfplane was first opened to a K_I of 0.05 MPa√m. Alternating minimizations and low temperature NVT relaxations were then performed to reach an equilibrium structure. Subsequently, the systems were loaded by incrementally increasing K_I by 0.071 MPa√m. While holding the periphery fixed to enforce the prescribed mode I displacement, the interior region of atoms was integrated at 300 K with a 0.1 fs time step for 5 ps with NVT dynamics at each mechanical loading step. The effect of different relaxation times between each loading step was previously tested for values of 2.5, 5, and 10 ps (Rimsza et al., 2018a). It was found that changing the relaxation time did not significantly alter the eK_IC values. Therefore, it appears that variation in the glassy silica structure has a greater impact on the mechanical response than the relaxation time used between loading steps. Through 100 loading steps the system was fractured with a final loading state of K_I = 0.65 MPa√m. During loading, four different chemical environments (vacuum, water, 1M NaCl, and 1M NaOH) were introduced in separate simulations to investigate the role of different aqueous solution compositions on fracture properties.

We use effective fracture toughness (eK_IC) in this paper, which is different than the K_IC used in classical fracture mechanics. Classically, the K_IC is the point at which a pre-existing fracture converts from subcritical to critical behavior. Due to our small system size, a single fracture event, in which a few bonds break at the fracture tip, results in fracture velocities on the order of 10 m/s. These high fracture velocities are indicative of a critical behavior, but several subsequent events would need to occur to produce an experimentally measurable change in fracture length. Therefore, we are defining the loading at which the first fracture event occurs as eK_IC to highlight that it is (1) defined differently than the K_IC in classical mechanics and (2) is likely to be an underestimation of the experimentally-derived K_IC values, which do not have the atomistic level accuracy exhibited in these simulations. Also, as we will see, eK_IC is impacted by the chemical environment unlike the strict definition of K_IC defined as an intrinsic material property in fracture mechanics text books (Lawn and Wilshaw, 1993).

For the vacuum systems, no water was introduced, and the internal surfaces stay unhydroxylated during the simulation. When aqueous solutions (water, 1M NaCl, 1M NaOH) are introduced it is necessary to add water to the fracture volume as the crack opens, so that water is available to wet the surface as the system is loaded. A grand canonical Monte Carlo (GCMC) method was used to add water into the crack during the simulation. The GCMC input used a chemical potential of −250 kcal/mol, corresponding to a target density of liquid water, and 10,000 attempts were made to insert water molecules at the beginning of every loading step. To avoid inserting water directly into the glass matrix, a modified GCMC method was used that limited insertion to within 3 Å of an existing water molecule. To maintain stable NaCl and NaOH concentrations, one NaCl or NaOH was added following the addition of water. Each NaCl or NaOH was added as a molecule at the widest point of the fracture, rather than inserted near the crack tip, to emulate the diffusion of ions from the bulk fluid. A reflective wall at the crack opening was used to keep the solution inside the fracture without adding external pressure.

After insertion of the water and NaCl/NaOH molecules, 5 ps of NVT dynamics at 300 K was performed to allow for equilibration of the system. The relaxation time for water-silica interactions was investigated by performing several fracture simulations with shorter (2.5 ps) and longer (10 ps) periods between loading steps. Previous studies have found that a range of relaxation times resulted in overlapping silica energy curves and minimal differences in the results (Rimsza et al., 2018a), indicating that variation in glass structure dominate differences in fracture due to simulated relaxation times.

Due to the high variation in the glassy silica structure, we used 12 replicas of the pre-cracked systems in each of the four loading conditions (48 simulations in total). These 12 replicas were constructed by rotating the three original silica systems by 0, 90,180, and 270° in plane and then introducing slit cracks on the -x half-plane. The results from these 12 systems are reported as the average and standard error unless otherwise noted.

As part of the analysis of the systems, the radius of curvature of the crack tip was calculated by fitting the locations of the surface atoms to a second order polynomial (y = f(x) = ax² + bx + c) and then obtaining the crack tip radius (ρ) from the formula (Howard, 1988):

The radius of curvature has been previously calculated for silica fracture in MD simulations using this method (Rimsza et al., 2017b).

The coordination environment of the Na⁺ ion was analyzed to identify how the salt solution was altering the silica surface. The coordination was calculated by using interatomic distances with cut-offs: Si-O (2.3 Å), H-O (1.25 Å), Na-O (2.20 Å), O-O (1.5 Å), and H-H (1.0Å). Cut-off values were selected based on previous investigations of Na⁺ adsorption on silica surfaces from NaOH solutions using the ReaxFF force fields (Rimsza et al., 2018b). See Rimsza et al. (2018b) for the pair distribution functions used for cut-off selection and further discussion.

During fracture the energy balance is controlled by the energetic cost associated with creating new surfaces, namely the surface energy, as described through the Griffith criterion (Lawn and Wilshaw, 1993). Therefore, a detailed understanding of the surface chemistry as the fluid infiltrates the fracture and reacts with freshly fractured surface is critical to identify the mechanisms that are responsible for environmentally-assisted fracture. Previous work has been done on the interactions between electrolyte solutions and silica surfaces (Ho et al., 2011; DelloStritto et al., 2016; Pfeiffer-Laplaud and Gaigeot, 2016); however, these have not typically included the application of external loading. A stress field can alter the structure of the silica by introducing voids, stretching bonds and ring structures (Muralidharan et al., 2005), and generally increasing the reactivity of the silica matrix (Rimsza J. M. et al., 2016).

Here we identify four different surface species. The most dominant surface species when the surface is contact with an aqueous solution is a silanol group, with a dangling bond terminated by one hydrogen (Si-OH). Several other surface species are present including protonated silanols (Si-OH₂), deprotonated silanols (Si-O⁻), and protonated siloxane oxygen atoms (Si-OH-Si). In this work the crack surface is dominated by either Si-OH or Si-O⁻, with concentrations of Si-OH₂ and Si-OH-Si being <0.06 #/nm².

The concentration of silanol groups on silica surfaces varies, but is reported as 4.0–5.0 #/nm² (Zhuravlev, 2000) and represents the total number of surface sites including neutral, deprotonated, and protonated silanol groups. We analyze all the available surface sites for protonation to determine how the silica surface has relaxed in the presence of the different solution compositions. Fewer surface sites would suggest that more silanol groups have recombined, forming Si-O-Si from Si-OH and Si-O⁻ sites. In our simulated silica, the initial concentration of available sites for formation of Si-OH groups is ~3.05–3.5 #/nm² (Figure 2A). The slight differences in surface site concentrate at 0.15 MPa√m is due to initial relaxation of the system before loading is applied. The silica surfaces exposed to salt solutions exhibit more sites (Si-OH, Si-O⁻) than those in pure water, suggesting that these surfaces have higher concentrations of broken Si-O bonds.

By separating the types of sites (Figures 2B,C), it becomes apparent that the surface species formed in 1M NaCl solution are more likely to be Si-OH in contrast to those formed in the 1M NaOH solution. The concentration of surface sites (Si-O⁻, Si-OH, Si-OH₂) and level of relaxation has been explicitly tied to the K_IC of the material (Lawn and Wilshaw, 1993) through the surface energy and elastic modulus.

The surface structure is therefore influenced by solution composition. These variations can be simulated when a reactive force field is used (compared with traditional MD force fields that do not allow for molecular disassociation). In these simulations the bonding states of atoms can change and the introduction of a NaOH or NaCl molecule into solution does not guarantee that the molecules will remain intact as the simulation evolves. Investigations of NaOH solution wetting a silica surface with a reactive force field (Rimsza et al., 2018b) noted the development of many different aqueous species including OH⁻, H₃O⁺, H₂, and O₂. Here, the solution composition is tracked over time to look at the development of these species.

In other work, the steric effects of solution infiltration in fractures focused on what, if any, impact occurred with changing crack tip accessibility (Jones et al., 2018), identifying that the size of the molecules was less important than the reactivity in altering the fracture properties of the material. Salt solutions are known to alter the silica surface due to adsorption (Dove and Craven, 2005). In the present context, the addition of OH⁻ creates a basic solution that increases the dissolution rates of glass (Du and Rimsza, 2017), while the addition of Cl⁻ has a less significant effect on the pH of the solution. Studies of the dissolution rate of glasses have shown the rate increases in both acidic and basic conditions, but is slower in acidic conditions (Alexander et al., 1954). In Figure 3, the OH⁻ and H₃O⁺ concentrations influencing pH in the simulations are reported as molar concentrations based on the number of water molecules in the fracture volume. Clearly, the most rapid change occurs during the initial loading, before the crack propagates, where OH⁻ decreases rapidly and H₃O⁺ increases slightly. The formation of a steady-state concentration occurs at ~0.25–0.3 MPa√m, due to a balance of the rate of water infiltration and addition of NaCl or NaOH molecules as the fracture is loaded.

The addition of the NaOH molecules results in a relatively steady OH⁻ concentration of 1M (Figure 3A). The concentration of OH⁻ drops more dramatically to <0.25 M for both the 1M NaCl and pure water solutions. The addition of NaCl molecules results in a higher concentration of both OH⁻ and H₃O⁺ in solution compared to pure water; however, both solutions become slightly acidic, with higher concentrations of H₃O⁺ than OH⁻. The concentration of H₃O⁺ increases with decreasing pH: 1M NaOH < water < 1M NaCl, as expected. The deviation from neutral pH indicates that in both salt solutions, silica dissolution should be higher than in pure water, increasing the chemical impact on chemical-mechanical subcritical fracture growth.

For the salt solutions, the silica surface interacts not just with water, but also with dissolved Na⁺ and either Cl⁻ (1M NaCl) or OH⁻ (1M NaOH) ions. Both salt solutions change the concentrations of surface species (Si-O⁻, Si-OH, Si-OH₂) when compared with water; this can be partially attributed to adsorption of Na⁺ to maintain local charge balance. In a previous study, we used ReaxFF (Rimsza et al., 2018b) to investigate Na⁺ adsorption environments with OH⁻ counter ions and identified several adsorption structures with Na⁺ ions coordinated by a mixture of water, silanols, and siloxane oxygen atoms. In this work, eight distinct Na⁺ structures are found in the NaCl and NaOH solutions. Different Na⁺ structures are identified by the coordinating oxygen species, which include the oxygen atoms in water molecules, silanol groups, and siloxane groups. The most common species is a free Na⁺ ion coordinated exclusively by water molecules (Figure 4A). These cations have at least one water molecule between them and the silica surface, and account for ~75% of the Na⁺ ions in the 1M NaCl simulations compared to only 61.5% in the 1M NaOH simulations (Table 1). The next three structures are inner-sphere monodentate (Figure 4B), bidentante (Figure 4C), and tridentate (not pictured), which include Na⁺ coordinated by water molecules and by one, two, or three Si-OH groups, respectively. These types of configurations (including the tridentate) were also observed on flat surfaces at 1M NaOH concentrations (Rimsza et al., 2018b). In contrast, the coordination of Na⁺ by Si-O⁻ species did not occur on flat surfaces at the same NaOH concentration, where Si-O⁻ was more commonly coordinated with either neighboring silanol groups or water (Rimsza et al., 2018b). However, in the fracture simulations, coordination of Na⁺ by Si-O⁻ does occur and is more common than the tridentate species (Table 1). We hypothesize that the confining conditions at the crack tip results in increased Si-O⁻ interactions with Na⁺ (Figure 4D) compared with flat surfaces. The final three configurations include Na⁺ that is at least partially coordinated by siloxane oxygen species. The first is a buried species (Figure 4F), which is a Na⁺ coordinated only by siloxane oxygen and indicates the Na⁺ has been entirely incorporated into the silica structure. Na⁺ diffusing into bulk silica occurs readily due to the open structure of the silica and Na⁺ is a common component of industrial glasses (Smit et al., 1978; Varshneya, 2013). The next is an adsorbed species (Figure 4E), defined as a Na⁺ ion coordinated by siloxane and water oxygen, but no silanol. Finally, there is a mixed coordination structure at the surface with Na⁺ ions coordinated by water, silanol, and siloxane oxygen atoms. Overall, the higher free and lower adsorbed Na⁺ ion concentrations in the 1M NaCl simulations may be leading to the higher Si-OH surface concentration and the lower deprotonation that is seen on the surface.

In addition to the coordination environment, the location of the Na⁺ ions can vary inside the fracture geometry. Several hypotheses on the role of pH on the fracture of silica suggest that the accessibility of the crack tip by OH⁻ is the primary means of enhancing fracture in basic solutions (White et al., 1987); therefore, understanding how the confined geometry of the crack tip alters the fracture properties is critical. Snapshots of the fracture geometry and the infiltrating fluid (Figure 5) provide some geometric insight into the crack tip accessibility by solution species. For silica fracture in water, the crack tip is primarily filled with complete water molecules along with a few H₃O⁺; the OH⁻ ions present do not migrate to the crack tip, in contrast to those in NaOH and NaCl solutions. In the NaCl solution, the fracture volume is filled with both Na⁺ and Cl⁻ ions, along with some OH⁻ and H₃O⁺ ions. Therefore, it appears that the presence of Na⁺ with a Cl⁻ counter ion may allow for easier migration toward the crack tip. In the breakdown of the different coordination environments (Supporting Information, Table S2) there is some limited coordination of the Na⁺ ion by Cl⁻, with an average of 0.52 ± 0.02 Cl⁻ ions coordinating each free Na⁺ ion in the first solvation shell. In the NaOH system, the Na⁺ ions have limited diffusion into the fracture volume, so that most of the aqueous species near the crack tip are either surface coordinated OH⁻ or free H₃O⁺ molecules. NaCl and NaOH solutions result in different fracture properties (discussed below), which is unexpected considering that both are promoting the diffusion of OH⁻ to the crack tip and potentially encouraging the hydroxide–mediated reaction (Equation 2). It is hypothesized that the hydration shells of the Na⁺ and Cl⁻ ions are affecting the OH⁻/H₃O⁺/H₂O-silica interactions, altering the energy barriers or the frequency of interactions (White et al., 1987) and the presence of different counter ions at the crack tip alters the activity of the OH⁻.

Siloxane bond strength in the glass structure varies with the number of oxygen bonds connecting two silica tetrahedra. Previous MD (Mahadevan and Garofalini, 2008) and DFT (Criscenti et al., 2006) simulations have been used to investigate the relative Si-O bond strength of different linkages to elucidate glass dissolution mechanisms as a function of glass composition. Here, the Q_n (n = 0–4) descriptor is used to distinguish silicon species, where n is the number of siloxane oxygen atoms attached to a silicon. A perfectly coordinated silicon atom is bonded to exactly four siloxane oxygen atoms and, hence, is a Q₄ species. Surfaces, defects, and network modifiers (e.g., Na⁺, Li⁺, Ca²⁺) can disrupt the glass structure, introducing silicon with lower concentrations of siloxane bonds and therefore lower Q_n species. Occasionally Q₅ species can also be formed, typically in low concentrations (<3%); this is a known defect in silica glass simulated with ReaxFF that has been previously reported for different silica systems (Rimsza J. et al., 2016; Rimsza and Du, 2018). Q₅ species have been identified by both MD simulations (Rimsza J. et al., 2016) and nuclear magnetic resonance (NMR) characterization (Zotov and Keppler, 1998), but at lower concentrations. By looking at the change in Q_n species in the structure, we can compare the relative strength of the siloxane bonds that link Q_n species in the glass.

The change in Q_n species from the first to the last loading step (0.15 MPa√m to 0.65 MPa√m) is provided in Figure 6 and generally, we expect higher Q_n species to transition to lower Q_n species during dissolution (Q₄ → Q₃ → Q₂ → Q₁ → Q₀). Regardless of environment, our simulations show that there is a decrease in Q₄ species. This is balanced by an increase in Q₃, Q₂, and Q₁ species. Interestingly, a water environment results in a decrease in the Q₀ species during loading, suggesting that the few Q₀ species that are formed reprecipitate in pure water. Previous computational results showed that during dissolution without mechanical stress, the Q₁ → Q₀ conversion has the lowest activation energy (Mahadevan and Garofalini, 2008). Therefore, when dissolution occurs we would expect a decrease in Q₁ species and an increase in Q₀. In contrast, the Q₄ → Q₃ and Q₃ → Q₂ transitions have higher activation energies and, hence, these species should be maintained in the system under dissolution-only conditions. As Figure 6 demonstrates, Q₃ and Q₂ concentrations are also maintained in these chemo-mechanical fracture simulations. When water is present, Q₃ and Q₂ appear to be more stable than in dry conditions, consistent with predictions from dissolution-only studies.

The impact of solution composition is most apparent when examining the Q₁ species data. Generally, there is an increase in Q₁ from vacuum to water to 1M NaOH to 1M NaCl. The only simulations that exhibit silica monomers released from the surface are those with a 1M NaCl solution, and the Q₀ is still in very low concentrations. In contrast, bulk glass dissolution is expected to occur most readily for the 1M NaOH solution (due to the higher pH). In the 1M NaCl solution the ~33% increase seen in Figure 6 is an average increase of total Q₀ speciation of ~0.5 silicon atoms, which is large only with respect to the low initial concentration in the replica systems. Despite the low change, several previous investigations have not seen spontaneous dissolution (Dove, 1999; Kubicki et al., 2012; Zapol et al., 2013; Rimsza and Du, 2017, 2018).

It is well-established that the strength of the silica decreases in the presence of an aqueous environment, which is typically attributed to dissolution inside the fracture volume and around the crack tip. The radius of curvature of the crack is a structural parameter that is connected to the amount of local inelasticity in a brittle material, with a blunter tip being associated with more inelasticity. It has also been investigated (Rimsza et al., 2017b) for its role in identifying the level of dissolution in the fracture volume, and how this impacts the chemical-mechanical fracture. Previous work (Rimsza et al., 2018a) noted that, in the presence of water, the radius of curvature decreased, resulting in a narrow crack tip, suggesting a higher stress concentration and a lower threshold for fracture propagation. Interestingly, the simulations with a 1M NaCl solution exhibit fractures with a larger radius of curvature than those with a 1M NaOH solution or pure water (Figure 7), and are the simulations that exhibit the highest increase in Q₁ species, and the only formation of simulations in which Q₀ species. Based on these results, we hypothesize that the presence of a highly reactive fluid like the NaCl solution under confinement may cause rapid increases in dissolution and in the radius of a crack tip, even within the picosecond timeframe in these MD simulations.

Due to the metastable nature of glass, the silica network within the bulk material is somewhat inelastic and is modified when an external stress is applied. The energy associated with the structural changes in the bulk material is described by the energy dissipation (G_diss) term in the well-known strain energy release rate (G_IC) equation:

which accounts for the energetic cost of fracture propagation in terms of the surface energy (2γ_s) associated with new surfaces and a non-conservative process to create them (i.e., irreversible changes in bulk silica structure quantified by G_diss). The G_diss term is extracted from the atomistic simulations through the change in internal energy of the silica (dU) and the added surface area (dA) (Rimsza et al., 2017b).

A derivation for Equation (5) is included in Appendix I. Here, surface area was calculated through construction of a surface mesh along the crack face with a spherical probe molecule with a radius of 4.5 Å as implemented in the Ovito Visualization Tool (Stukowski, 2009). In previous investigations, the energy dissipation has been associated with the amount of inelasticity occurring around the crack tip (Rimsza et al., 2017b, 2018a). Using this description, the following trend in inelastic energy generation was calculated for the four fracture environments considered: vacuum > 1M NaOH solution > 1M NaCl solution > pure water (Figure 8). The higher G_diss value for the 1M NaOH solution compared with the other aqueous solutions suggests that there is more unrecoverable inelastic strain in the region surrounding the crack tip.

Identification of eK_IC values for the four silica systems was based on step-wise changes in the potential energy of the silica during loading, as outlined in previous studies of silica fracture (Rimsza et al., 2017b). A sudden change in energy is accompanied by bond breakage, indicating a crack propagation event. The loading at which these crack propagation events occur are averaged over multiple simulations to estimate the eK_IC for the system. The results identify that the silica systems fractured in vacuum exhibit the highest eK_IC value, the fewest fracture events, and the highest G_IC value (Table 2). This is consistent with experimental observations that silica is stronger in a non-reactive or dry environment (Rimsza et al., 2018a). Experiments have shown that the addition of any electrolytes that perturb the pH of the solution make the system fracture at lower loads than in neutral water (Dove, 1995). From the simulation data reported here, the eK_IC values overlap, but it appears that the silica in 1M NaCl solution has a slightly higher eK_IC than that in water or a 1M NaOH solution. It is hypothesized that a broader crack tip radius (Figure 7) reduces the drive for fracture propagation at these time scales.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.