Figure 1 . a) Sketch of the simulation setup. The slab is made of α-quartz and contains a notch that guides tensile fracture propagation. The domain is a three-dimensional block dimensions lx=300 nm, ly=12 nm, and lz=60 nm. The length of the initial notch is 60 nm, and the aperture is 16 Å on the left side. The upper and lower parts of the block (red areas) are fixed. The block is exposed to a compressive normal stress along the x -direction (\(\sigma_{x}\)), followed by a tensile strain along the z -direction (\(\sigma_{z}\)).b) Orientation of the block with a x -, y - andz -cut. Both the crystallographic axes (a1, a2, a3 and c) and the coordinates axes (x , y , and z ) used in our simulations are shown.c) A sketch of a fracture, drawn from a simulation snapshot. Solid lines represent observed damage while the dashed line represent potential damage at the atomic scale. Box 1 (red) illustrates a propagation path influenced by the crystallographic structure (see also Figure 2d), whereas box 2 (red) illustrates an area where the propagation path is less influenced by the crystallographic structure (see also Figure 2e).
First, a compressive normal stress is imposed along thex -direction using an isothermal-isobaric thermostat (NPT). Then, an extensional displacement controlled by the canonical ensemble (NVT) is imposed along the z -direction, perpendicularly to the notch. The imposed tensile strain leads to an initial tensile stress in the range 1000-2000 MPa, and the stress is measured just before the crack starts to propagate. These levels of tensile stress agrees with expectations for dynamic rupture propagation in the crust (e.g., Reches and Dewers, 2005). For most simulations, the α-quartz block has dimensions 300x12x60 nm, the imposed compressive stress in thex -direction is 500 MPa, the extensional displacement along thez -axis is imposed over a period of 30 ps, and the temperature is set to 300 K.
The simulations are run for 500 ps with a time step of 1 fs, during which a tensile (mode I) dynamic crack propagates. During this time, depending on the initial tensile stress imposed on the system, one of three outcomes are observed: 1) the crack does not propagate from the initial notch; 2) the crack starts propagating and stops inside the simulation domain; or 3) the crack propagates throughout the simulation domain along the x -direction. If the crack does not propagate, the simulation is discarded because the loading was insufficient for dynamic crack propagation. If the crack propagates, but does not cross the whole domain, we increase the duration of the simulation such that the crack reaches the end of the simulation domain. Thus, we obtain a set of simulations where a crack has propagated throughout the entire system. We sample the position of the crack tip in thex -direction through time and calculate the average rupture speed. To assess the sensitivity of the rupture speed on the details of the simulation setup, we varied the loading plane, compressive stress, temperature, and loading rate. The effect of these variables on the rupture speed is described in the Supplementary Information, but overall the rupture speed is not particularly sensitive to these parameters and we observe the same qualitative behavior.
2.2 Interatomic interaction potential for α-quartz
We perform the simulations using the molecular dynamics package LAMMPS (Plimpton, 1995). We model an α-quartz with the interatomic potential initially defined by Vashishta et al. (1990) using the updated parametrization from Broughton et al. (1997). The 1997-version of the potential is modified from the Nakano et al. (1994) three-body potential, which was modified from the original version of the potential (Vashishta et al., 1990). The first versions of the potential were designed to reproduce the properties of molten, crystalline, and amorphous silica. Broughton et al. (1997) adjusted three parameters to optimize the material properties towards those of crystalline α-quartz. Tables S1 and S2 display the parameters used in the simulations.
The 1997-version of the Vashishta potential reproduces the anisotropy of α-quartz, including the elastic parameters and the compressibility of the lattice lengths. To ensure that the simulated material has properties similar to α-quartz within the pressure and temperature range of our simulations, we have calculated the elastic parameters at T=1 K and P=0.1 MPa (Table S3), the density in a temperature range 1-520 K at 0.1 MPa and in a pressure range 0-10 GPa at 298 K, as well as the compressibility parameters of the lattice along different crystalline directions (a1 and c) at 298 K in a pressure range 0-3 GPa (see Supplementary Information). Overall, when the applied pressure is below 3 GPa, the compressibility along the x - andz -axis, corresponding respectively to the crystallographic a1-axis and c-axis, agrees with experimental studies (Jorgensen, 1978; d’Amour et al., 1979; Hazen et al., 1989; Glinnemann et al., 1992). The density of our simulated material is consistent with the density of α-quartz in the temperature range 1-520 K. At 298 K, the density (ρ ) is slightly lower, but has a similar dρ /dP -slope in the range 0-3 GPa, which is within the range of our simulations, while it deviates at higher pressures (Figure S3). Based on the elastic parameters (Table S3) and the density, the Rayleigh wave speed is 3048 m/s when applying an extensional displacement on the (0 0 0 1) plane and 3097 m/s when applying an extensional displacement on the (\(0\ \overset{\overline{}}{1}\ 1\ 0\)) plane, using the procedure of Vinh and Ogden (2005). For comparison, Leong et al. (2018) measured the Rayleigh wave speed to be 3841 m/s when applying an extensional displacement on the (\(5\ \overset{\overline{}}{10}\ 5\ 8\)) plane and 3694 m/s when applying an extensional displacement on the (\(1\ \overset{\overline{}}{2}\ 1\ 3\)) plane. The difference between our measured speeds and the speeds reported by Leong et al. (2018), apart from different loading planes, is that the modelled material is slightly softer, as indicated by the elastic parameters (Table S3).
3 Results
For all the simulations, we describe the fracture surface roughness after the crack has propagated through the simulation domain by classifying the fracture geometry as straight cracks, oscillating cracks, and cracks exhibiting microbranching (Figure 2). When oscillations occur, we measure the amplitude of these oscillations (Figure 2f). For the straight cracks (Figure 2a), the propagation produces a combination of mirror and mist surface roughness. In some areas, the crack propagates along a single atomic plane, producing a mirror surface (Figure 2a, box 1). In other areas, the crack propagation creates a mist surface, by shifting between atomic layers (Figure 2a, box 2). Wherever the crack changes course or turns abruptly (Figure 2d), or where the crack branches (Figure 2e), damage is created in the surrounding quartz structure. The damage is characterized by a lower atomic density, and for a critical atomic distance of 1.9 Å, the coordination of silicon and oxygen atoms decreases from 4 and 2 in α-quartz to 3 and 1 in damaged quartz.