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.