Figure 3. a) Rupture speed (V R) of the crack tip as a function of initial tensile stress. Above a critical speedV oscillate = 0.15CR, whereC R is the Rayleigh wave speed, the crack tip oscillates. Above a critical speed V branch = 0.26CR, crack branching occurs. b) Change in potential energy of the simulation domain from the relaxed initial system and the relaxed state after the crack has propagated through the domain, as a function of rupture speed. ɣ s is the surface energy.
The perfect quartz lattice is an energy minimum for our system. One way to characterize the damage caused by crack propagation is to compare the potential energy of the relaxed initial system with the potential energy of the relaxed state after the crack has propagated through the system. To obtain high crack speeds, the system has to be loaded substantially above the fracture toughness of the material, and this leads to thermal energy being released during crack propagation. This energy is absorbed by the thermostat during and after crack propagation and we make sure to measure the potential energy at the same temperature for the initial state and the cracked state. The potential energy is averaged over the last 50 ps of each simulation, and this averaging window starts at least 5 ps after the crack has reached the far end of the sample. For a perfectly straight crack, this energy difference should be equal to 2 s, where ɣ s is the surface energy of a flat crack surface and A is the surface area of the new crack. In our simulations, the surface energy for straight cracks was measured by dividing the change in potential energy over the crack surface area, and found to be equal to 0.93 J/m2. For comparison, the surface energy of the (\(0\ 0\ 0\ 1\)) plane of quartz in air has been measured experimentally to be 0.92 J/m2 (Parks, 1984). Figure 3b shows how the dissipated energy varies with rupture speed. The results show that the energy dissipated in the system increases with increasing rupture speed due to fracture oscillation and microbranching that damages the material. The contribution of potential energy of new surfaces and damaged quartz structures varies for the oxygen and silicon atoms (Figures S10a and S11a). For the silicon atoms, the new surfaces provide a higher potential energy than the damaged structures. Conversely, for the oxygen atoms, the energy in the new surfaces is lower than in the damaged structure. This can be explained by the energy from the steric repulsion and the Coulomb interactions between oxygen pairs (see Supplementary Information).
4 Discussion
4.1 Morphology of the fracture surface
Fracture surfaces have been studied extensively in brittle materials, like polymethyl methacrylate (e.g., Fineberg et al., 1991; Sharon and Fineberg, 1996), Homalite (Ravi-Chandlar and Knauss, 1984), silica glass (Sharon and Fineberg, 1998) or neo-Hookean brittle gels (Livne et al., 2005). These studies showed that, above a threshold rupture speed, the crack surface changed from smooth to a jagged structure that undergoes instabilities characterized by either oscillations or microbranching. Our results show that the morphology of a crack surface in a model α-quartz varies with rupture speed. At low speeds, the crack is propagating straight (Figure 2a). When the rupture speed increases, the crack undergoes path instabilities first by oscillations, where the amplitude of the oscillations increases as the crack speed increases (Figures 2c, d and f) and then, at higher rupture speeds by branching (Figure 2e).
Because polymethyl methacrylate is an amorphous material, Fineberg et al. (1991) and Sharon and Fineberg (1996) argued that the source of the observed oscillations in the experiments could not be linked to the material structure. Instead, the occurrence of the oscillations was attributed to a sudden increase in rupture speed. In contrast to an amorphous material, we observe that fracture propagation is controlled by the crystalline structure of α-quartz to some degree. For example, the atom layers control the transition from mirror to mist surface as the crack jumps between different layers. At the onset of oscillations, the oscillations have a small amplitude and are irregular; they do not seem to follow the crystalline structure (e.g., Figure 2c). However, when the oscillations reach a more regular amplitude, they follow the crystalline structure of quartz (e.g., Figures 1c and 2d). During microbranching, certain areas with oscillations are controlled by the crystalline structure while, in other areas, zones of damage show a slightly rounded shape (Figure 1c), which cannot be directly linked to the crystalline structure of α-quartz.
4.2 Rupture velocity and crack instabilities
The threshold speed above which path instabilities are observed has been measured in several materials and found to be in the range 0.3-0.74C R (e.g., Zhou et al., 1996; Buehler and Gao, 2006; Bleyer and Molinari, 2017), and experiments on single-crystal quartz observed branching when the crack approached 0.45C R (Leong et al., 2018). In our simulations, we have identified two threshold speeds:V oscillate=0.15C R is the speed of the first path instability that produces oscillations, andV branch=0.26C R is the speed at which the first microbranching is observed. These two speeds are slightly lower than experimental measurements performed in quartz on larger samples and significantly lower than the theoretical limit ofC R. In the next paragraph, we argue that this difference is controlled by the functional form of the cohesive force between silicon and oxygen atoms.
Buehler and Gao (2006) modeled an elastic 2D material and tuned the cohesive force between interacting atoms to study how this parameter controls rupture speed. The cohesive force is dependent on the breaking of atomic bonds, which occurs at a critical atom separation (r break). This critical atom separation can be expressed as the ratio between the interatomic distance where the potential energy slope is steepest, and the minimum potential energy. It is one of the key parameters governing the crack instability speed. Whenr break=1.135, a dynamic instability can occur for rupture velocities below 0.2C R, while whenr break increases to 1.25, a dynamic instability occurs for rupture velocities above 0.73C R. In the interatomic potential used in our simulations,r break for the two-body interactions between silicon and oxygen atoms is equal to 1.144 (Figure S1). This value falls in the lower range of investigated r breakparameters (Buehler and Gao, 2006) and is consistent with the threshold speeds for crack propagation instabilities in the range 0.15-0.32C R observed in our simulations. It is therefore likely that the shape of the two-body interaction between silicon and oxygen, and in particular the critical separationr break governs the rupture speed. In future adjustments of the parameters of silica potentials, one should therefore consider this quantity explicitly.
4.3 Production of damage during dynamic rupture
Damage zones in faults contain a high density of fractures at all scales, which may be created by several processes occurring during or after the fault formation, such as Andersonian fracturing, early fault tip migration, fault tip linkage (e.g., Johri et al., 2014; Mitchell and Faulkner, 2009), and dynamic ruptures (Rudnicki, 1980; Wilson et al., 2003; Paul et al., 2007). Models of dynamic rupture calculate the stress field around a fault tip during propagation, which allows quantifying the damage along the length of the rupture (e.g., Madariaga, 1976; Andrews, 1976). Our simulations indicate that a propagating tensile crack creates nanoscale damage. This damage is generated in the material surrounding the propagating crack and aborted branches, or in areas where oscillations shifts direction during propagation (e.g., Figure 2d).
Another type of damage produced by dynamic rupture is pulverized rocks, which have been shattered in situ without significant shear strain. These rocks consist of very fine grains (e.g., Dor et al., 2006). If our damaged crystal were exposed to further strain through shear deformation, this would lead to the formation of nanosized fragments (Figure 1c). In the simulations where rupture propagation produces microbranches, we have measured the height and length of nascent fragments. Independently of initial tensile stress, the height of the fragments would be around 20 nm, while their length would vary in the range 40-110 nm. For cracked grains in a gouge, the lower limit for fragment diameter is reported to be 30 nm; however, the fragments can also reach sizes up to 100 µm (Keulen et al., 2007). For grain sizes below 1.2 \(\mathbf{\pm}\) 0.3 µm, Keulen et al. (2007) proposed that these grains are produced by another mechanism than grinding. Sammis and Ben-Zion (2007) explored mechanisms that might produce fragment below the grinding limit. These authors showed that nanometer-sized fragments rarely form by grain crushing in simple shear under compressive loading, but could be produced during tensile loading at high strain rates. An estimate of the smallest fragment produced by tensile stress can be found by\(\mathbf{K}_{\mathbf{\text{IC}}}\mathbf{=\sigma}\sqrt{\mathbf{\text{πa}}}\), where K IC is the fracture toughness and 1 MPa m1/2 for quartz (e.g., Iwasaki and Torikai, 1993),\(\mathbf{\sigma}\) is the tensile stress and a is the radius of the flaw. In our simulation, the tensile stress is in the range 1-2 GPa, which will give a =80-320 nm. If the flaws are almost equal in size and distributed uniformly in the rock, a will indicate the dimensions of the fragments, otherwise even smaller fragments could be formed (Sammis and Ben-Zion, 2007). From our simulations, even the longest fragments (110 nm) are much smaller than the grinding limit for quartz and in the lower range of estimated fragment sizes due to tensile loading condition. Therefore, dynamic rupture at molecular scale could provide an explanation for the formation of rock fragments with dimensions smaller than the grinding limit. This is supported by Wilson et al. (2005) who argued that the smallest fragments in gouges might be explained by damage from a crack tip during dynamic rupture.
5 Conclusion
We perform molecular scale simulations of dynamic rupture propagation in α-quartz loaded under tensile stress. We identify critical rupture speeds for crack oscillations and crack branching consistent with results from Gao and Buehler (2006). When the rupture speed is below 460 m/s, the crack propagates straight, while at speeds above 460 m/s, it starts oscillating. This rupture speed corresponds to 0.15C R. For microbranching to occur, the rupture speed must exceed 800 m/s, which corresponds to 0.26C R. During a crack branching event, two branches propagate simultaneously until one branch dominates. These two branches may provide an initial structure that could lead to quartz nanosized fragments with dimensions below the grinding limit. These fragments in fault damage zones have a high surface area and thus may react relatively fast in the presence of aqueous fluids, leading to fast fault healing and strength recovery.