The chemistry of glass fracture has been under scientific debate for a long time (Fairbarn and Tate, 1859; Griffith, 1921). However, a quantitative description of the underlying mechanisms of crack propagation, critical to the further development of crack-resistant glasses, is presently unavailable (Wondraczek et al., 2011; Freiman, 2012). This is largely due to the fundamental nature of cracking as it occurs in the broad group of brittle oxides, and especially in glassy silicates: crack propagation reflects a complex and, as of yet, not fully understood interplay of chemical and physical interactions (Ciccotti, 2009). In particular, below a critical value of applied stresses (the subcritical regime), chemical reactions and the effects of mechanical loading are so intertwined that the separation of the governing parameters becomes an extremely challenging task. The subcritical crack growth and its dependence on environmental conditions have profound consequences for the glass industry and the suitability of glass as a structural material. Here, parametric knowledge is needed not only for the focused development of new material compositions and toughening procedures but also for enabling better time-to-failure estimations (Freiman et al., 2009). Beyond materials development, this need encompasses other related areas such as the understanding of geophysical fracture processes (Atkinson, 1982; Brantut et al., 2013; Zhang and Zhao, 2014).

Generally speaking, the kinetics of crack propagation below Mode I critical stress concentration K_Ic can be divided into three regimes: in region I, even though the applied stresses are much smaller than the critical value, the crack propagates through a directional corrosion reaction on stressed bonds. For the case of molecular water as the reaction partner, a theoretical framework has been developed (Michalske and Freiman, 1982; Michalske and Bunker, 1984, 1987) (Figure 1). In this region, the crack propagation velocity depends strongly on the applied stresses, temperature, and environment. In region II, the crack growth velocity becomes increasingly limited by the transport of the reactive molecules to the crack tip, resulting in a constant crack velocity with increasing stress intensity. Finally, in region III, the applied stresses are high enough to break the stressed atomic bonds without environmental help, so the crack velocity is again strongly dependent on the applied stresses and temperature, but not on the environment.

Generally, the rate of crack growth v is modeled as ν = λ⋅k (Michalske and Bunker, 1984; Schoeck, 1990), where λ is a geometrical parameter depending on the matrix through which the crack is growing, and k is the rate of success for the crack to advance by λ units. The reaction rate is then framed within the transition state theory, and written in the form of an Arrhenian expression, k = k₀ exp(−ΔF/RT), where k₀ is the attempt frequency based on the oscillations of the system around an equilibrium position and ΔF is the molar energy difference between the initial and final states, in this case representing the energy barrier to break the bond. While working on the delayed fracture of glasses, Wiederhorn and coworkers (Wiederhorn, 1967; Wiederhorn et al., 1974a,b) pioneered a semiempirical equation to describe the phenomenon of subcritical crack growth in the presence of water: (1)v=v0a(H2O)nexp−E∗−bKIRT.

Here, v₀ is the pre-exponential term, a(H₂O) is the water activity, n is the number of water molecules which interact with the Si–O–Si bond during its cleavage, E* is the activation energy, K_I is the stress concentration factor at the crack tip, and b is a function of the activation volume and the crack tip radius (Wiederhorn et al., 1974a,b). The mathematical formulation for the activation energy in Eq. 1 is analogous to the enthalpy of gasses under constant pressure, only the pressure term is substituted by the applied stresses (Wiederhorn, 1967; Wiederhorn and Bolz, 1970). This expression has found widespread use in the literature, influencing a great part of the developments on the understanding of subcritical crack growth since its inception.

For the moment, we note that the primary feature of Eq. 1 is the reduction of E* by the strain energy U_strain = bK_I. This can be understood as an offset contribution in the potential well of the considered bond, increasing the probability for breakage (case I in Figure 2).

The derivation of Eq. 1 assumes that the applied mechanical field offsets the potential energy at the considered bond, thus reducing the activation energy of the breakage event and increasing breakage probability. However, this appears as only one of two extremes: the second is an increase in the overall excitation of the system (case II in Figure 2) without directly affecting the activation barrier. This interpretation is a natural result of the democratic fiber bundle model (DFBM), in which the material is represented as a series of springs connecting two parallel bars and oriented perpendicular to the direction of crack propagation. These springs are considered to be ideally elastic, ideally brittle with random failure threshold (modeling inherent structural disorder), and the load on each spring is also subjected to thermal noise (Kun et al., 2000; Roux, 2000; Ciliberto et al., 2001; Scorretti et al., 2001).

Here, we consider first the energetic equilibrium for cracks growing at rates below the limit of dynamic stability (Fineberg and Marder, 1999), (2)G=1−μ2EAI(v)KI2=2γ, where G is the strain release rate, μ is Poisson’s ratio, E is the elastic modulus, A_I is a material and crack speed dependent variable, and γ is the surface energy. Equation 2 can be reasonably approximated by Freund (1990), (3)G=(1−v∕cR)KI2E∕1−μ2=2γ, (4)∴v=cR1−2Eγ1−μ2KI2, with c_R being the Rayleigh wave speed. However, Eq. 4 is only valid for K_I > K_Ic, so to extend the range for which the equation is applicable, it has to be rewritten. Here, we consider the term 1−2Eγ∕(1−μ2)KI2 as equivalent to the first two terms of a Taylor series expansion, leading to this original equation: (5)v=cRexp−2Eγ(1−μ2)KI2.

Equation 5 fulfills the requirement of convergence to Eq. 4 for K_I > K_Ic, and additionally it produces a non-zero crack growth velocity at K_I = K_Ic and v = 0 for K_I = 0. From the transition state theory, we can write that (6)v=k0exp−ΔFRT=v0exp−2γ(NA∕ρA)RT, where the energy difference between the initial and final states ΔF is approximated by the ratio between the energy needed to create two surfaces and the number of chemical bonds that make up this surface, so N_A is Avogadro’s constant and surface bond density (ρ_A). This approach is consistent with the approximation of ideally brittle fracture and the lattice trapping understanding of fracture (Hsieh and Thomson, 1973; Lawn et al., 1980; Michalske and Bunker, 1984; Schultz et al., 1994; Zhu et al., 2004). Putting Eqs 5 and 6 together yields: (7)v=v′0exp−2γRTρA∕NA+(1−μ2)KI2∕E or alternatively, (8)v=v′0exp−E∗RT+CKI2, where v′0 is the pre-exponential term, E* = 2γN_A/ρ_A and C = (1 − μ²)N_A/Eρ_A. As noted above, Eq. 8 shares the Arrhenian scaling with the approach of activation offset, and by including both thermal as well as structural noise into its variance, there is a clear parallel with the average time-to-failure predictions derived through DFBM (Kun et al., 2000; Roux, 2000; Ciliberto et al., 2001; Scorretti et al., 2001). The main difference to Eq. 1 is that instead of decreasing the magnitude of the energy barrier to bond failure, the strain energy U_strain = CK_I² acts concertedly with the thermal energy RT, increasing the excitation level of the bond in question. This collaborative action can be intuitively thought as stemming from the fact that elastic deformation and thermal expansion have its origin on how the acting stress and temperature interact with the interatomic bond potential, leading to the same outcome: a shift on the interatomic equilibrium distance.

While the empirical usefulness of Eq. 1 is widely accepted in glass science, here, we focus on testing the parallel validity of Eq. 8. This is challenging because literature provides a relatively broad scatter in the values of activation energy for Si–O bond rupture, ranging from approximately 50–150 kJ mol⁻¹ (Lindsay et al., 1994; Walsh et al., 2000; Del Bene et al., 2003; Zhu et al., 2005). Such a variation may result in prediction errors of crack velocity of several orders of magnitude. Therefore, we used a highly accurate multilevel CCSD(T):MP2//DFT approach for generating data on the homogeneous dissociation of the Si–O bond. Specifically, we evaluated the energy of 624 kJ mol⁻¹ for the central Si–O bond of a dimeric hydridosilsesquioxane molecule. This model is used to reflect the energetics of straining and breakage of a single silica bond at the crack tip because it is small enough to be computationally accessible, but provides a level of accuracy of <2 kJ mol⁻¹, unmatched by alternative methods on larger systems (Feller et al., 2011). Since the energy was used to create two “surfaces,” we will be using the value of 312 kJ mol⁻¹ for comparison to experimental data on region III cracking. Calculations were also performed for Si–O bond cleavage in the presence of one, two, and three H₂O molecules (mimicking region I cracking), yielding values of 164.2, 166.6, and 158.3 kJ mol⁻¹, respectively. No simulation was performed for more than three water molecules interacting with the DHS molecule; experimental data suggest that the hydration number of the Si–O–Si bond in amorphous silica lies between 2.4 and 1.8 (Fournier and Marshall, 1983). These values will be used in the following for reference.

Applying Eq. 8 to published experimental data on subcritical crack growth of silica (Suratwala and Steele, 2003) and a soda–lime silicate glass (Wiederhorn, 1967; Wiederhorn and Bolz, 1970), we first calculate the constant C. Theoretical calculations yield values of 1.13 × 10⁻⁶ and 1.05 × 10⁻⁶ m⁴ (N mol)⁻¹, respectively, for the two different glasses. In parallel, rearranging Eq. 8 allows for the experimental determination of C, i.e., using the slope of a plot of T versus K_I², (9)T(KI)=−CRKI2+E∗Rln(v0∕v).

In Figure 3, the relation above is applied to subcritical crack growth data of soda–lime silicate glass measured between 2 and 90°C (Wiederhorn and Bolz, 1970). The value of C does not remain constant through the complete range of temperature–stress intensity data, and even as it stabilizes there is a disparity of 4 orders of magnitude between the experimental fit and the theoretical calculation. This effectively means that the actual “conversion factor” of strain to thermal energy is much smaller than the theoretical one, which might be due to many different effects which are not covered by our approach: first, for lower crack growth velocities, the applied stress intensity factor is much closer to the threshold value K_th [the static fatigue limit, below which the crack ceases to grow (Fett et al., 2005)]. This means that the stresses around the crack tip become increasingly shielded from the applied stresses due to diffusive exchange between mobile ions in the glass and the hydronium ions in the environment, leading to an increased slope of the T × K_I² graph. Second, in parallel with the DFBM approach (Scorretti et al., 2001), the parameter C can be interpreted as a measure of structural disorder, related to how the applied stresses are preferentially concentrated on the bonds located closer to the crack tip. However, the stress concentration is assumed to follow a Gaussian distribution, and that might not reflect the actual distribution, leading to the disagreement between the experimental and theoretical values. It is also worth noting that Eq. 8 is derived under the assumption of static crack growth with constant K_I, while the experimental tests were conducted under constant load. Therefore, following the observed disparity, we used the experimentally fitted values to construct the Arrhenius plots shown in Figure 4.

The authors are unaware of published crack growth rates of silica glass under vacuum; therefore, the comparison has to be made using crack growth in water vapor. Silica glass exhibits accelerated water diffusion under stress (Tomozawa et al., 1991), which causes the crack surfaces to expand, creating compressive stresses and shielding the crack tip from the applied stresses (Fett et al., 2005; Wiederhorn et al., 2011). This explains the apparently anomalous results reported for silica (Suratwala and Steele, 2003), and silica aerogels (Despetis et al., 2004), where the measured crack propagation speed decreases with increasing temperature. Wiederhorn et al. (2013b) have developed a formalism that enables the direct calculation of the shielding stresses as a function of both applied stresses and temperature. The estimated activation energies for region I are 114 ± 8 kJ mol⁻¹ for the soda–lime silicate glass and 156 ± 5 kJ mol⁻¹ for silica glass; for region III the activation energy for the soda–lime silicate is 314 ± 9 kJ mol⁻¹.

As commented before, on a first approximation of ideal brittleness (Célarié et al., 2003a,b; Guin and Wiederhorn, 2004, 2006; Wiederhorn et al., 2013a,b), the only contribution to the strain release rate G during fracture is the creation of the fracture surface. Therefore, with proper knowledge of the Si–O–Si surface bond density, the energy required to break each Si–O bond can be estimated by the surface energy. As its value depends only on the type of bond that is being broken, the calculated values for glassy and crystalline silica should be similar. The surface energy of silica glass under N₂ atmosphere and at 300 K is reported (Wiederhorn, 1969) to be approximately 4.4 J m⁻²; while MD simulations (Leed and Pantano, 2003) yield an area density of Si–O–Si bonds of 7.2 bonds nm⁻², resulting on an activation energy of 360 kJ mol⁻¹. This is remarkably consistent with monocrystalline α-quartz, with values of 310 and 340 kJ mol⁻¹ [for cracks propagating on the a<112¯0> plane, on directions normal to the planes z<011¯1> and r<101¯1>, respectively] calculated using fracture (Atkinson, 1979; Timms et al., 2010) and bond density data (Bloss and Gibbs, 1963) under inert atmosphere.

The activation energy values calculated from cleavage in monocrystals and the simulation data obtained by CCSD(T):MP2//DFT agree to a remarkable extent with the results from Eq. 8, which leads us to believe that our proposed equation provides accurate predictions of crack propagation rates.

As noted initially, the principal difference between the two approaches rests on how the applied mechanical strain affects the chemical bonds directly at the crack tip. The offset approach [as seen in Eq. 1, but also revised elsewhere (Henderson et al., 1970; Bartenev, 1973)] considers that the applied stresses decrease the energy barrier for bond breakage. Some simulation studies (Gagnon et al., 2001) support this interpretation, but they also give a very different mechanistic description of crack propagation, following the coalescence of voids forming ahead of the crack tip; a model that seems to be at odds with the most recent experimental investigation (Wiederhorn et al., 2013a). A possible route for resolving this apparent discrepancy would be to combine these considerations with the aspect of concerted activation, in which the applied stresses play a role analogous to temperature, increasing the overall excitation level of the system without altering the energy profile of the chemical bonds (Stillinger, 1995; Sciortino, 2005). Thereby, it is important to acknowledge the empirical nature of both Eqs 1 and 8, if applied individually. We hypothesize that in reality, the physical reactions taking place may be best understood on the basis of simultaneous reactions of activation offset and concerted activation.

We considered parametrization of models for crack propagation in brittle glasses. The action of strain energy can be modeled in two empirical ways, i.e., acting either as an offset to the activation barrier for bond rupture, or concertedly with the thermal energy of the system. Both approaches provide accurate reproductions of experimental data. However, non-empirical treatment is expected to require simultaneous implementation.

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

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.