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
2Aɣ 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.