Equilibrium Molecular Dynamics

Subscriber access provided by READING UNIV

Article

Elastic Characterization of S- and P- Wave Velocities in MarineLike Silica: The Role of Non-Equilibrium Molecular Dynamics Dolores Melgar, Marco Lauricella, Gareth S O'Brien, and Niall Joseph English J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b10552 • Publication Date (Web): 16 Jan 2018 Downloaded from http://pubs.acs.org on January 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Elastic Characterization of S- and P- Wave Velocities in Marine-like Silica: the Role of NonEquilibrium Molecular Dynamics Dolores Melgar*†, Marco Lauricella‡, Gareth S. O’Brien§, and Niall J. English*† †

School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland. ‡

Instituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche, Rome, Italy. §

Tullow Oil Limited, Central Park, Leopardstown, Dublin 18, Ireland.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 27

ABSTRACT

The α-quartz polymorph of SiO2 forms the basis of mineral sands stable down to 100 km depths below the surface, making it of central geoscientific relevance. The characterization of the nanoscale properties of these materials is of importance, especially for elastic properties governing phonon and sound propagation, and is of very high industrial relevance for oil exploration. Here, for the first time, we apply Non-Equilibrium Molecular-Dynamics simulation to analyze the propagation of an artificial velocity perturbation in silica systems, and, in so doing, determine Sand P- wave velocities, in a manner redolent in concept to seismic-based oil-exploration approaches. This propagation has been analyzed systematically by means of different metrics in terms of spatio-temporal system response; these produce consistent results, by and large. In particular, we find excellent quantitative agreement with experimental S- and P- wave velocities, in many cases.


2

Page 3 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION Silica is one of the most common materials in nature. For instance, silicates constitute 90% of the minerals in the Earth's mantle and crust. Furthermore, silicon dioxide is essential in many industrial applications and diverse fields, e.g., bio-medicine,1 optics,2

photonics,3 and

electronics.4 Moreover, one of the most important silica polymorphs, namely α-quartz, is a very common mineral and it is stable down to 100 km depth below the surface,5-6 rendering it of central geoscientific relevance. The characterization of the nano-scale physical properties of these materials is of importance and, indeed, of intrigue; several efforts have been reported in the literature in order to characterize, for instance, the structure and elastic properties of these systems, including shock-wave response, as well as the dislocation creep regimes and surface properties, among other interesting physicochemical characteristics.7-12 In view of these intriguing physical properties, particularly the elastic behavior and the shockwave response, acquiring a greater understanding of the mechanisms underpinning the elastic properties which govern phonon and sound propagation acts as a potent motivator, and is of very high industrial relevance for oil exploration. Not alone do silica-rich sandstone rocks form the majority of hydrocarbon-bearing reservoirs, but the elastic properties thereof determine their seismic response. Seismic exploration methods has acted as the cornerstone of subsurface investigation both in academia and in the hydrocarbon-extraction industry, where it can generate high-resolution subsurface images.13-14 Indeed, understanding SiO2 elastic properties, and their relationship with other material properties such as electrical conductivity, is of critical importance in assessing and de-risking hydrocarbon prospects, by means, for instance, of the Controlled-Source Electromagnetic (CSEM) technique. Since CSEM methods provide images with less resolution than ones originating from seismic investigations, a combination of both


3


Page 4 of 27

techniques would be a very interesting option for petroleum-exploration technology and it requires a detailed understanding of the relationship between elastic and electrical material properties. Regarding the elastic properties of geological relevant materials, one of the most interesting studied features is the velocity at which elastic waves propagate within the material, which depends on its density and stiffness. In a geological context, seismic waves are usually divided into two groups, namely body and surface waves, depending on the part of the planet the waves propagate (interior or surface, respectively). Regarding body waves, these are classified depending on the particularities of the motion of the particles as P- or S- waves. P-waves (or Primary waves) are longitudinal waves, e. g. the motion of the particles is parallel to the travel direction. They can propagate both in fluids and solids. They are the fastest body waves, and, therefore, the first ones to be detected by seismographs. S-waves (also known as Secondary or Shear waves) are characterized by a transversal motion of the particles with respect to the travel direction. As fluids do not support shear stresses, S-waves can only propagate in solids and they are always slower than the P-waves. Classical molecular dynamics (MD) is a very useful technique in order to obtain insights into the microscopic properties of silica materials. The interest in applying this methodology to SiO2 materials is reflected in the different forcefields developed in the last decades.15-20 Among the different MD applications, non-equilibrium MD (NEMD) has been used to calculate α-quartz properties, such as its thermal conductivity,21 and has also been applied widely to analyze shockwave effects on SiO2 materials.22-29 However, most of these studies related to shock waves focus on the microstructural response of the material under shock-loading conditions, and how shock compression can lead to a densification or phase transition, being the velocity of propagation of


4

Page 5 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the elastic waves a sidelight. For isotropic systems, such velocities are directly related to the bulk and shear modulus. Examples can be found in the literature were Ab Initio Molecular Dynamics simulations were performed in order to calculate the shear and bulk modulus of a given material and, then, simply compute the isotropically averaged sound velocities.30-31 For non-isotropic systems as α-quartz, the relationship between the elastic properties of the material and the seismic waves velocities becomes more complicated: the Christoffel equation32 allows the calculation of wave velocities in different directions from the elastic stiffness tensor components. Several examples of the application of this equation can be found in the literature, where equilibrium MD techniques were employed.33-35 Here, for the first time, we apply Non-Equilibrium Molecular-Dynamics simulation to analyze the propagation of an artificial velocity perturbation in silica, α-quartz systems, and, in so doing, determine S- and P- wave velocities, in a manner similar in concept to seismic-based oilexploration approaches. The perturbation is directly introduced in the system by altering the velocities of the atoms of a particular region of the material. To the best of our knowledge, this methodology was never used in NEMD simulations for determining seismic wave velocities at an atomistic scale. A similar approach was used some years ago to study the sound wave acceleration in granular materials,36 but in that case, a discrete element method was applied, in which a set of spherical particles were meant to emulate the material grains and no atomistic resolution was considered. The perturbation propagation will be analyzed by means of different system-response metrics (in a spatio-temporal sense), and the corresponding results will be compared with the velocities calculated for trigonal systems by Farnell37 and the measured velocities by McSkimin;38 we find excellent quantitative agreement with experimental S- and P- wave velocities in many cases. The


5


Page 6 of 27

agreement with the experimental results for α-quartz validates our methodology. The importance of these results is the development of a new methodology, consisting in explicitly introducing a perturbation in a given region of the system and then following its propagation, which is applied here for the first time. This new approach opens the door for future applications on different materials and conditions. Seismic velocities are not particularly hard to measure experimentally, but the NEMD method we present in this paper permits the complete control on the features of the geological sample in terms of composition, porosity, lattice heterogeneities, stress, temperature, and other petrophysical characteristics of the sample, as well as the particularities of the seismic wave itself, as its magnitude and direction, because of the full control that MD simulations allow on the design of the sample represented in the simulation box and also on the creation of the perturbation. Therefore, relationships between the seismic waves velocities and different characteristics of the sample can be established in a more systematic, direct and simpler way.

COMPUTATIONAL DETAILS We created an orthorhombic bulk model of α-quartz from the unit cell reported by Levien et al.7 (see Figure 1).The x and z laboratory axes were aligned onto a and c crystallographic ones. The y-laboratory axis is perpendicular to x and z, and forms a right-handed system. This orientation of the crystal axis is in accordance with the IRE Standards on Piezoelectric Crystals 39 and, therefore, consistent with the work of McSkimin.38 The selected unit cell was replicated 30 times in each direction, containing 81000 SiO2 units. The box volume was equilibrated under three-dimensional periodic boundary conditions in the NPT (isothermal-isobaric) ensemble. Then, the box lengths were selected so as to fit to the experimental density of α-quartz at 300 K


6

Page 7 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 1 atm (ρ=2.6485 g/cm3),40 but respecting the ratios obtained by NPT equilibration. The resulting box lengths were: Lx=147.43 Å, Ly=127.68 Å, and Lz=162.11 Å. Owing to its relatively high quality and accuracy, the silica forcefield proposed by Vujić and Lyubartsev was used to describe both bonded and non-bonded interactions,41 with the smooth particle-mesh Ewald method for long-range electrostatics.42 The bonded interactions are described by means of harmonic potentials for both Si-O bonds and O-Si-O angles. The non-bonded interactions are represented by assigning to the atoms an electrostatic charge and the corresponding Van der Waals parameters. All of the simulations were performed using the DL-POLY Classic package.43 Classically-propagated molecular dynamics were carried out for 10 ns (with a time step of 1 fs) to equilibrate the system in the NVT ensemble (canonical ensemble) at 300 K using a NoséHoover thermostat with a relaxation time of 0.5 ps,44 followed by a 1 ns run under the microcanonical (NVE) ensemble, with excellent energy conservation and essentially no temperature drift; the final configuration was used as the reference system to be perturbed in NEMD.

Figure 1. Top figures represent the unit cell of α-quartz reported by Levien et al. from different perspectives, inspired by the wave travel directions [100], [120], and [001] in the direct-lattice


7


Page 8 of 27

vector basis, in a ball-and-stick representation (yellow for Si atoms and red for O atoms). Bottom figures are the surfaces resulting from the replication of the top figures. The compasses in the left of the figures indicate the lattice vectors orientation. In order to analyze how a velocity perturbation is propagated along different directions in a quartz crystal, we chose three travel directions, parallel to the x-, y-, and z- axes (or [100], [120], and [001] in the direct-lattice basis vectors). The propagation along the [010] direction was not studied because, in the case of α-quartz, it is symmetrically equivalent to [100]. This perturbation was created by altering the initial velocities of the atoms belonging to a particular designated orthorhombic slab from x, y or z = 0 Å to x, y or z = 2 Å, depending on the travel direction of the wave (x-, y- or z- axis, respectively). This perturbation was performed in two different ways: (i) to gauge P-wave velocity (vP), the change in the starting velocities is along the travel direction, whilst (ii) for estimating S-wave velocity (vS), a velocities component perpendicular to the travel direction is perturbed. For instance, when considering zˆ as the travel direction, to simulate a Pwave, vz is altered, while for mimicking an S-wave, both vx and vy modifications are explored. The reader is addressed to Figure 2 for an illustrative description. The modification of a single component of the selected atoms' velocity (perpendicular or parallel to the travel direction, depending on the type of wave) avoids the possibility of mixing longitudinal and transversal modes. As three-dimensional periodic boundary conditions are applied to the system, reflection on the boundaries of the sample and, then, possible coupling of the longitudinal and transversal modes, cannot occur. In both P- and S- wave cases, the respective laboratory-axis velocity components of the selected atoms’ initial velocities were set up to 20 Å/ps, with equal negativedirection momentum distributed evenly to existing atomic velocities over all other atoms in the system, so as to conserve total linear momentum. This perturbation amplitude was chosen


8

Page 9 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


because it corresponds to the highest initial velocity observed after the equilibration for the atoms to be perturbed. Therefore, the amplitude of the perturbation is large enough to be detected over the thermal noise, but small enough to produce a linear elastic wave. As the choice of the perturbation magnitude could seem arbitrary, we followed the same procedure described in the following sections for two different perturbation amplitudes (10 Å/ps and 40 Å/ps), in order to check possible dependencies of our results on this particular feature of the perturbation. For the three traveling directions, P- and S- waves' velocities were calculated and then compared with the ones obtained for the 20 Å/ps perturbation. The difference on the seismic velocities was in all the cases less than 5% (the average change is less than 2%). In this context, it is worth to highlight that our results have an associated uncertainty around 3% on average. Therefore, we concluded that, at least for this order of magnitude, there is not a strong dependence of our results on the amplitude of the initial perturbation.

Figure 2. Schematic representation of the change in velocities performed in order to simulate the S-wave (a) and P-wave (b) and how the simulation box is divided into slabs in order to analyze their propagation, considering the z axis as the travel direction. This division in slabs is applied only at the measuring stage, while the sample constitutes an infinite and perfect crystal during


9


Page 10 of 27

the NEMD simulations. As three dimensional periodic boundary conditions were used, the schematic representation must be interpreted as it was replicated along x, y and z. This means, for instance, that slab #1 and #10 will behave as in direct contact. Moreover, it is worth noting that this is a schematic representation to understand how the velocities (and not the positions) are directly altered. However, it seems obvious that forcing the particles to have a given velocity in a particular direction will have associated a particle displacement in that direction just after the perturbation is created. After the perturbation is created, an NVE simulation was then performed on each perturbed system, with the perturbation propagating along the two directions. The 1 fs-step NEMD simulations were extended until the collision of each wave with its periodic ‘mirror-image’ wave propagating away from the slab along the propagation axis but in opposite direction; the NEMDsimulation time was 1 and 2 ps for the P- and S- wave cases, respectively. It is worth noting that three dimensional periodic boundary conditions were applied to all the simulations presented in this paper. In Figure 2, this implies, for instance, that slab #1 is in contact with slab #10. Therefore, when a wave is traveling in the z direction and reaches slab #10, after a given time it will reach slab #1 again. The analysis of each perturbation’s propagation called for division of the simulation box into 2 Å-width slabs along the travel direction (see Figure 2), followed by spatio-temporal evolution of various slab-based system-response metrics: temperature, difference in (number-) density from equilibrium, the forces, and mass-velocity along the initial perturbation direction. The reader must take into account that this slab-division of the system is only a way of measuring different properties at a given distance of the source of the perturbation and that the sample within the studied waves travel is still an infinite and perfect α-quartz crystal. In the S-wave case, as the


10

Page 11 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


initial perturbation occurs in a direction perpendicular to the travel direction, no marked difference in density from equilibrium is expected along the travel direction, due to the periodic boundary conditions. With the information of the spatio-temporal evolution of the different metrics, the time-step at which the maximum of each metric was reached in each slab (charting the wave’s ‘time of arrival’) was plotted versus the position of the corresponding slab. This result in a set of straight lines (time-distance curves) whose slope is related to the wave-propagation velocity (vide infra). Unfortunately, atomic thermal ‘noise’, inter alia, can hinder crisp and unambiguous identification of these maxima. For instance, in the case of slab temperature, after the local maximum due to the perturbation, there is a ‘drift’ that can shift this local maximum somewhat in time. Therefore, in order to be sure that the considered maximum is the one corresponding to the wave’s ‘time of impact’ at each slab, two different methods were applied. The simplest considers the metric’s absolute maximum, independently on times of occurrence, whilst the other defines a time window and uses it to find the correct local maximum: for the slab adjacent to the initial perturbation, this local maximum should be found within the time window defined. Then, for subsequent slabs, the local maximum should be located within a time window centered in the timestep in which the previous slab’s local maximum was found. To be sure that the definition of this window does not affect our results, we used several values, i.e., 100, 200, 300, 400, 500, and 600 fs. Therefore, for a given perturbation, and for each metric, we may derive a set of seven values for the wave-propagation velocity: the one calculated with no time window and the ones from the six time windows studied. For this set of values, a weighted mean was calculated, with its corresponding standard deviation. The weight of each estimate was computed taking into account the goodness-of-fit parameters of each plot. Another possible approach would be to perform several NVE simulations from a different initial configuration


11


Page 12 of 27

and, then, average the wave's time of arrival over the set of trajectories. This average would remove the thermal noise and, thus, could avoid mistaking the maximum of a metric indicating the wave's impact with a post-impact maximum due to the aforementioned 'drift', but the need of running several simulations for each perturbation will considerably increase the computational cost of the seismic velocities calculation. Therefore, we chose the time-window approach in order to optimize our resources.

RESULTS AND DISCUSSION Figure 3 provides a representative example of time evolution of slab-based mass velocity; for clarity, only the values around the maximum for the 20 first slabs are shown. Nevertheless, salient and important features of wave propagation emerge: these velocities form a peak for each slab, indicating that the perturbation actually reached the analyzed slab. From Figure 3, it can be observed that, as the wave-front’s distance from the originally-perturbed slab increases, these peaks become lower in height and wider, as well as less symmetric and smooth, indicating that attenuation occurs in the wave’s wake as it propagates, as indeed expected.


12

Page 13 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Mass-velocity evolution with time for the first 40 slabs. Only the values around the maximum for each slab are represented for clarity. Also, the maximum for the initially perturbed slab (slab #1) is not shown, because it is much higher than the other maximum values (19.98 Å/ps). Once the timestep at which the P- and S-waves reaches every slab has been computed, they were plotted against the distance from the origin of the perturbation. For all the metrics studied and all the propagation directions two straight lines with the same slope were found (see Figure 4 for an example). This double rectilinear pattern has two important interpretations. On one hand, the fact that the timestep at which the metrics’ maximum value takes place for each slab forms a straight line when represented as a function of the distance from the perturbation origin means that the perturbation created in slab #1 is propagated at constant velocity through the quartz crystal, this velocity being related with the slope of the abovementioned straight line. On the other hand, the fact that two straight lines with almost the same slopes were observed means that the perturbation created propagates in two directions with a very similar velocity. When the initial perturbation is created, the velocity component in the propagation direction for the selected atoms is changed to a given value. This means that these atoms will have a tendency at the very beginning of the simulation to displace in that direction, even if their initial positions were not directly altered. The crystal structure represented by the bonded interactions in the forcefield will oppose to such move and would try to bring back those atoms to their equilibrium position, causing a displacement in the opposite direction, and, therefore, another wave with similar characteristics but traveling in the opposite direction. From Figure 1, as the initial perturbation was created in slab #1, the first wave will propagate from slab #1 to slab #2 and so


13


Page 14 of 27

one, while the second wave will almost instantly propagate from slab #1 to slab #10, due to the periodic boundary conditions, in the opposite direction of the first wave. Tables 1-3 compares the seismic velocities calculated for both travel directions (v+ and v-) along the propagation axis to show their similarity. Tables 1-3 include the averaged values for the velocities over the different time windows considered. If some of the v+ and v- values and their standard deviations appear to differ for some cases, it is worth highlighting that each calculated velocity itself bears an intrinsic uncertainty from least-squares fitting. This uncertainty is not included in Tables 1-3, but was taken into account to affirm that both waves propagating in opposite directions travel with the same velocity.

Figure 4. Timestep at which the mass-velocity maximum value is reached for each slab along the z-axis. Both waves are represented: the one in the z-axis positive direction (z > zpert, blue line) and the one in the negative direction (z < zpert, red line).

Table 1. Summary of the velocities calculated for P- and S-waves propagated along the x- axis in both positive (v+) and negative (v-) directions averaged over all the time windows considered


14

Page 15 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and the corresponding standard deviation. The last column collects the results averaged over all the time windows for both directions. The reader is addressed to the text for more details about the time windows used to perform this analysis. The metrics analyzed are temperature (T), mass velocity (mv), the difference in density with respect to the equilibrium (Δρ), and the sum of forces in the selected direction (ΣF). Type of wave

Metric Analyzed

v+ (km/s)

v- (km/s)

v+,- (km/s)

P

T

7.73±0.24

7.80±0.00

7.79±0.07

mv

7.65±0.01

7.62±0.02

7.63±0.02

Δρ

8.37±0.73

8.67±1.68

8.49±1.25

ΣF

7.87±0.54

7.14±0.65

7.43±0.69

T

4.34±0.03

4.44±0.09

4.40±0.09

mv

4.35±0.06

4.43±0.05

4.39±0.07

ΣF

4.49±0.40

4.40±0.39

4.46±0.39

T

4.66±0.01

4.68±0.01

4.67±0.02

mv

4.67±0.01

4.68±0.01

4.67±0.01

ΣF

4.64±0.12

4.63±0.13

4.63±0.12

Sy

Sz

Table 2. Summary of the velocities calculated for P- and S-waves propagated along the y- axis in both positive (v+) and negative (v-) directions averaged over all the time windows considered and the corresponding standard deviation. The last column collects the results averaged over all the time windows for both directions. The reader is addressed to the text for more details about the time windows used to perform this analysis. The metrics analyzed are temperature (T), mass velocity (mv), the difference in density with respect to the equilibrium (Δρ), and the sum of forces in the selected direction (ΣF).


15


Page 16 of 27

Type of wave

Metric Analyzed

v+ (km/s)

v- (km/s)

v+,- (km/s)

P

T

7.83±0.12

8.24±0.11

8.12±0.22

mv

7.91±0.18

8.27±0.17

8.10±0.25

Δρ

8.26±0.62

7.83±0.54

7.98±0.59

ΣF

9.40±0.55

8.04±1.31

8.65±1.22

T

4.24±0.08

4.43±0.08

4.35±0.12

mv

4.54±0.06

4.60±0.12

4.57±0.10

ΣF

4.52±0.60

4.31±0.42

4.40±0.50

T

4.57±0.03

4.53±0.06

4.55±0.04

mv

4.58±0.02

4.71±0.08

4.64±0.09

ΣF

4.60±0.28

4.66±0.05

4.64±0.16

Sx

Sz

Table 3. Summary of the velocities calculated for P- and S-waves propagated along the z- axis in both positive (v+) and negative (v-) directions averaged over all the time windows considered and the corresponding standard deviation. The reader is addressed to the text for more details about the time windows used to perform this analysis. The last column collects the results averaged over all the time windows for both directions. The metrics analyzed are temperature (T), mass velocity (mv), the difference in density with respect to the equilibrium (Δρ), and the sum of forces in the selected direction (ΣF). Type of wave

Metric Analyzed

v+ (km/s)

v- (km/s)

v+,- (km/s)

P

T

8.05±0.12

7.97±0.00

7.98±0.05

mv

8.11±0.00

7.98±0.00

8.03±0.07

Δρ

8.04±0.69

7.56±0.23

7.81±0.56

ΣF

7.86±0.16

8.09±0.08

7.96±0.17


16

Page 17 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Sx

Sy

T

4.93±0.04

5.01±0.05

4.97±0.06

mv

5.00±0.05

5.02±0.05

5.01±0.05

ΣF

4.72±0.28

5.06±0.27

4.93±0.32

T

4.94±0.04

4.98±0.04

4.96±0.04

mv

5.03±0.03

5.03±0.03

5.03±0.02

ΣF

4.21±0.62

4.48±0.22

4.37±0.44

In order to test that these results are not dependent on the particular features of our simulation box, the same procedure was applied to a set of boxes with different size. The same unit cell reported by Levien et al. was replicated a different number of times in the three directions of the same, resulting in a set of simulations boxes with the same structure but with different box lengths combinations (from 25 Å to 325 Å). The seismic velocities calculated from these simulations are all very similar. Therefore, we concluded that our results are not dependent on the system size. The simulation box with lengths Lx=147.43 Å, Ly=127.68 Å, and Lz=162.11 Å was chosen for a more detailed study because it represents a point of compromise between the quality of the statistics and the computational cost of the NEMD simulations. To check the consistency of the vP and vS values calculated from the analysis of the different metrics, a set of single-factor ANOVA tests were performed.45 The objective on these tests was to find out if the data sets from different metric analysis were representing the same mean value. The tests were performed for a significance value of α = 0.01. For the waves propagating along the x-axis, all of the metrics should be considered for calculating vP, although only temperature and mass velocity are reliable metrics for calculating both vS. For the waves propagating along the y-axis, the data from the analysis of the temperature, the mass velocity, and the difference in


17


Page 18 of 27

density for the P-wave are consistent, whilst for the vSx, only the temperature and the sum of forces should be taken into account. In contrast, for the vSz, all the metrics produce reliable data. Lastly, for the waves propagating along the z-direction, all of the metrics should be included for the calculation of vP and vSx. However, in the case of vSy, the three metrics analyzed seem to produce different mean values, in spite of their closeness. Therefore, in this case we will take the metric whose least-squares fittings produce better goodness-of-fit parameters, i.e., the mass velocity. Bearing the ANOVA-test conclusions in mind, averaged vP and vS values and the corresponding standard deviations were computed for the three travel directions considered. These results are presented in Table 4.

Table 4. Averaged values for P- and S- waves' velocities taking into account the metrics that the ANOVA tests with a significance value of 0.01 show as consistent for each direction of propagation. Travel direction

Type of wave

v (km/s)

x

P

7.70±0.28

Sy

4.40±0.07

Sz

4.67±0.01

P

8.07±0.36

Sx

4.36±0.25

Sz

4.61±0.10

P

8.01±0.11

Sx

4.99±0.10

Sy

5.03±0.05

y

z


18

Page 19 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results from Table 4 were compared with the ones from the works by Farnell37 and McSkimin.38 The phase velocity of a plane wave propagating in an anisotropic but homogeneous solid, such as α-quartz, varies depending on the travel direction. From Hooke's law, which relates the stress and the strain tensors of a material, the relationships between both tensors with the particles displacement, and considering a plane wave as an assumed solution for the equations of motion, non-trivial solutions are found for three phase-velocity values: one corresponding to a compressional wave and two associated with shear waves (referred to as fast and slow shear wave, FS and SS, respectively, from this point on). In the work of Farnell,37 the value of these velocities was calculated from the experimental values of the stiffness constants and taking into account the particularities of the α-quartz crystal symmetry, whilst in the case of the work of McSkimin,38 the approach was to measure the velocities in order to calculate the stiffness constants. Both sets of data are in excellent agreement, and they considered x-, y-, and z-cuts of the α-quartz crystal, with the same orientation than the NEMD simulations presented in this paper. These results are collected in Table 5.

Table 5. Velocities for the P-, FS-, and SS- waves reproduced from the works of Farnell37 and McSkimin.38

Direction of propagation

vP (km/s)

vFS (km/s)

vSS (km/s)

Farnell

McSkimin

Farnell

McSkimin

Farnell

McSkimin

x

5.73

5.75

5.08

5.11

3.32

3.30

y

5.99

6.01

4.32

4.32

3.88

3.92

z

6.35

6.32

4.66

4.69

4.66

4.69


19


Page 20 of 27

Given the symmetry of the α-quartz crystal, the x- and z- axes are pure-mode axes for longitudinal waves. The x-axis is also a pure-mode axis for shear waves, whilst the z-axis is a degenerate axis for transverse waves, in which both shear waves will travel with the same velocity. This means that for the waves propagating along the x-axis and for the compressional wave along the z-axis, our NEMD results can be directly compared with the ones from Farnell 37 and McSkimin.38 For the other cases, one has to be more careful, as the angle between the rod and the direction of the energy propagation is not zero. But, as this angle is always less than 30º, and considering the length scale of our NEMD experiments, as well as the way of measuring the wave time of arrival, given that we are not examining solely the particular direction in one dimension, but rather a maximum value of a series of metrics in a set of volumetric slabs, the deviation of the direction of the energy propagation from a particular axis should not strongly affect our results. Therefore, the NEMD results from Table 4 and the experimental results from Table 5 can be compared, bearing in mind that the Sy- and Sx- waves are the SS-waves for the perturbations travelling along the x and y directions, respectively. Seismic velocities are expected to increase under stress. To check the physical reliability of our methodology, NEMD simulations at different levels of stress were performed. The introduction of different levels of stress was done by means of changes in the volume of the simulation box. We proceed as follows: from the equilibrated configuration whose volume produces the right density for α-quartz at 300 K and 1 atm, the volume was reduced by 5% keeping the ratio between box lengths constant. Then, a 1 ns NVE run was performed in order to check that the temperature drift was not relevant. This final configuration was used to generate the velocity perturbation in the same way as explained in the Computational Details section, and also to produce the starting point for the next simulation, whose volume is reduced again by 5%. The


20

Page 21 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


process was repeated to produce 5 equilibrated simulation boxes with volumes 95%, 90%, 85%, 80%, and 75% the original box, whose volume was optimized to reproduce the experimental density of α-quartz at 300 K and 1 atm. The seismic velocities calculated for the three modes traveling in the three directions of the space at different degrees of stress follow the same pattern: the closer the pressure is from 0, the lower the velocities. Pressure is related to the structural stress of the system. Therefore, our results show that the less stressed the structure, the slower the seismic waves. This behavior of the calculated seismic velocities is in perfect agreement with the knowledge on seismic waves and, thus, it strengthens the physical meaning of our methodology.

CONCLUSIONS A set of classically-propagated NEMD simulations were performed in the NVE ensemble, in which α-quartz SiO2 bulk was perturbed artificially to study how this perturbation evolves spatio-temporally within the system. This perturbation consists of altering the parallel and orthogonal components to three directions of propagations for atomic velocities in a slab, depending on the type of wave (P- or S-) velocity to be calculated. The directions of propagation considered were those of the laboratory axes. By means of dividing the box into slabs along the travel direction and calculating the maximum of different metrics in each slab, we gauge spatiotemporal evolution of temperature, the difference in density with respect to equilibrium, massvelocity and forces acting on the particles (in the orthogonal or parallel directions, depending on the type of wave to be studied). For all of the properties measured, both the P- and S- waves travel along both directions of the propagation axis with constant velocity, as evidenced by the straight-line motif in Figure 3 (and repeated over all metrics analyzed). Moreover, single-factor


21


Page 22 of 27

ANOVA tests were performed in order to consider only consistent velocities extracted from the analysis of different properties for the average velocity values. An absence of system-size- and magnitude-perturbation-dependence of these results was also confirmed. Moreover, we have also proved that the velocities calculated with this new methodology react under stress as expected. Simulations performed at different levels of stress show that the more stressed the structure, the higher the wave velocities, for both longitudinal and transversal modes in the three travel directions. This is a very important point that strengths the physical meaning of our method to predict seismic velocities. Comparing results from Table 4 and 5, several conclusions can be extracted: 

The NEMD-calculated velocities are over-estimated, but by less than 35% in all cases. In the case of experiments (both for calculating the velocities or the stiffness constants) the crystal sample is finite and not-perfect. In our NEMD simulations, we used (infinitely) periodic boundary conditions and a perfect α-quartz unit cell. These two features can justify an over-estimation of 20-35% of the seismic velocities for a reasonably good forcefield. When only 1% of the SiO2 units were randomly removed, the velocities decrease for all the modes and travel directions. This decrease is about 6% on average, but in some cases it represents a reduction of 11%. Therefore, it seems reasonable to assume that the imperfection, as well as the boundary and size effects, of the experimental sample could justify the over-estimation in our NEMD-calculated velocities.



For the three directions studied, the relation between vS and vP is between 54% and 63%, which matches perfectly the well-known estimation in rock mechanics of vS as


22

Page 23 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 3 vP. This feature supports the previously mentioned idea that longitudinal and

transversal modes were not being mixed. 

Our methodology reproduces the existence of a fast and a slow shear wave for the xand y- directions, as well as the degeneracy in the S-wave velocity for the z-direction.



The experimental results establish that the P-wave is faster in the z- than in the ydirection, as well as it is faster in the y- direction than in the x- direction. Our results show a slower P-wave along the x- axis compared to the other ones. In the case of the y- and z- directions, our results are very similar, but, due to the higher standard deviation for the vP in the y- direction, our results are largely consistent with the mentioned experimental observation.

The agreement between the experimental observations and the values extracted from our simulations is sufficiently quantitative to consider our method robust. The application of this method to other geologically relevant systems opens up a new avenue for research, underscoring the role of molecular simulation in acting as a rich prototyping environment for geological science and oil/seabed exploration. Quartz is an important constituent for other rock types with embedded pore fluids. The good results for vP and vS calculated from NEMD results encourage future investigations on attenuation and petrophysical relationships.

Supporting Information. The following files are available free of charge. Specific details on the time-distance curves parameters, including slope, calculated velocity, and goodness-of-fit parameters, for all the metrics, directions of propagation, and types of waves,


23


Page 24 of 27

together with the information of the corresponding ANOVA tests (PDF)

Corresponding Authors *E-mail: [email protected] Phone: +353 1 7161812 *E-mail: [email protected] Phone: +353 1 7161646

Author Contributions The NEMD simulations, as well as their analysis, were performed by DM. Discussions about the results obtained were carried out with the participation of all authors. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

ACKNOWLEDGMENTS The authors thank Christian Burnham, Zdenek Futera and Mohammad Reza Ghaani for helpful discussions. Regarding funding sources, DM thanks the Irish Research Council and Tullow Oil plc for an Enterprise postdoctoral fellowship (Project ID: EPSPD/2016/26), together with NE to Science Foundation Ireland (Grant Code: SFI 15/ERC-I3142) for funding of computing resources.


24

Page 25 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES 1. Feng, J.; Li, J.; Wu, H.; Chen, Z. Metabolic responses of HeLa cells to silica nanoparticles by NMR-based metabolomic analyses. Metabolomics 2013, 9, 874-886. 2. Shi, F.; Tian, Y.; Peng, X.; Dai, Y. Combined technique of elastic magnetorheological finishing and HF etching for high-efficiency improving of the laser-induced damage threshold of fused silica optics. Appl. Opt. 2014, 53, 598-604. 3. Jiang, L.-Y.; Wu, H.; Jia, W.; Li, X.-Y. Polarization-independent negative refraction effect in SiO2-GaAs annular photonic crystals. J. Appl. Phys. 2012, 111, 023508. 4. Campbell, M. G.; Liu, Q.; Sanders, A.; Evans, J. S.; Smalyukh, I. I. Preparation of nanocomposite plasmonic films made from cellulose nanocrystals or mesoporous silica decorated with unidirectionally aligned gold nanorods. Materials 2014, 7, 3021-3033. 5. Palache C.; Berman H.; Frondel C. Dana's System of Mineralogy. GFF, New York, US, 1962. 6. Cohen, L. H.; Klement, W. High-low quartz inversion: Determination to 35 kilobars. J. Geophys. Res. 1967, 72, 4245-4251. 7. Levien, L.; Prewitt, C. T.; Weidner, D. J. Structure and elastic properties of quartz at pressure. Am. Mineral. 1980, 65, 920-930. 8. Hirth, G.; Tullis, J. Dislocation creep regimes in quartz aggregates. J. Struct. Geol. 1992, 14, 145-159. Iglauer, S.; Salamah, A.; Sarmadivaleh, M.; Liu, K.; Phan, C. Contamination of silica 9. surfaces: Impact on water–CO2–quartz and glass contact angle measurements. Int. J. Greenhouse Gas Control 2014, 22, 325-328. 10. Bragg, W.; Gibbs, R. E. The structure of alpha and beta quartz. Proc. R. Soc. A 1925, 109, 405-427. 11. Howm, R. A. Silica: physical behavior, geochemistry and materials applications. Mineral. Mag., 1996; 60, 390-391. 12. Wackerle, J. Shock‐wave compression of quartz. J. Appl. Phys. 1962, 33, 922-937. 13. Yilmaz, Ö. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data. Society of exploration geophysicists, Tulsa, US, 2001. Aki, K.; Richards, P. Quantitative Seismology, University Science Books, New York, 14. US, 2002. 15. Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D. Development of an empirical force field for silica. Application to the quartz−water interface. J. Phys. Chem. B 2006, 110, 2782-2792. 16. Vessal, B.; Amini, M.; Leslie, M.; Catlow, C. R. A. Potentials for molecular dynamics simulation of silicate glasses. Mol. Sim. 1990, 5, 1-7. 17. van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Force fields for silicas and aluminophosphates based on ab initio calculations. Phys. Rev. Lett. 1990, 64, 1955-1958. 18. Munetoh, S.; Motooka, T.; Moriguchi, K.; Shintani, A., Interatomic potential for Si–O systems using Tersoff parameterization. Comput. Mater. Sci. 2007, 39 (2), 334-339. 19. Carré, A.; Horbach, J.; Ispas, S.; Kob, W., New fitting scheme to obtain effective potential from Car-Parrinello molecular-dynamics simulations: Application to silica. Europhys. Lett. 2008, 82, 17001. 20. Cowen, B. J.; El-Genk, M. S. On force fields for molecular dynamics simulations of crystalline silica. Comput. Mater. Sci. 2015, 107, 88-101.


25


Page 26 of 27

21. Yoon, Y.-G.; Car, R.; Srolovitz, D. J.; Scandolo, S. Thermal conductivity of crystalline quartz from classical simulations. Phys. Rev. B 2004, 70, 012302. 22. Wang, J.; Rajendran, A. M.; Dongare, A. M. Atomic scale modeling of shock response of fused silica and α-quartz. J. Mater. Sci. 2015, 50, 8128-8141. 23. Farrow, M. R.; Probert, M. I. J. Atomistic molecular dynamics simulations of shock compressed quartz. J. Chem. Phys. 2011, 135, 044508. 24. Barmes, F.; Soulard, L.; Mareschal, M. Molecular dynamics of shock-wave induced structural changes in silica glasses. Phys. Rev. B 2006, 73, 224108. 25. Renou, R.; Soulard, L.; Lescoute, E.; Dereure, C.; Loison, D.; Guin, J.-P. Silica glass structural properties under elastic shock compression: experiments and molecular simulations. J. Phys. Chem. C 2017, 121, 13324-13334. 26. Kubota, A.; Caturla, M.-J.; Davila, L.; Stolken, J.; Sadigh, B.; Quong, A.; Rubenchik, A. M.; Feit, M. D. Structural modifications in fused silica due to laser-damage-induced shock compression, Proc. SPIE 2002, 4679, 108-116. 27. Dávila, L. P.; Caturla, M.-J.; Kubota, A.; Sadigh, B.; Díaz de la Rubia, T.; Shackelford, J. F.; Risbud, S. H.; Garofalini, S. H. Transformations in the medium-range order of fused silica under high pressure. Phys. Rev. Lett. 2003, 91, 205501. 28. Su, R.; Xiang, M.; Chen, J.; Jiang, S.; Wei, H. Molecular dynamics simulation of shock induced ejection on fused silica surface. J. Appl. Phys. 2014, 115, 193508. 29. Izvekov, S.; Rice, B. M. Mechanism of densification in silica glass under pressure as revealed by a bottom-up pairwise effective interaction model. J. Chem. Phys. 2012, 136, 134508. 30. Shigeaki, O. Elastic properties of CaSiO3 perovskite from ab initio molecular dynamics. Entropy 2013, 15, 4300-4309. 31. Martorell, B.; Wood, I. G.; Brodholt, J.; Vočadlo, L. The elastic properties of hcpFe1−xSix at Earth's inner-core conditions. Earth. Planet. Sci. Lett. 2016, 451, 89-96. 32. McSkimin, H. J.; Jr., P. A.; Thurston, R. N. Elastic moduli of quartz versus hydrostatic pressure at 25° and − 195.8°C. J. Appl. Phys. 1965, 36, 1624-1632. 33. Bedrov, D.; Ayyagari, C.; Smith, G. D.; Sewell, T. D.; Menikoff, R.; Zaug, J. M. Molecular dynamics simulations of HMX crystal polymorphs using a flexible molecule force field. J. Comput. Aided Mater. Des. 2001, 8, 77-85. Sewell, T. D.; Menikoff, R.; Bedrov, D.; Smith, G. D. A molecular dynamics simulation 34. study of elastic properties of HMX. J. Chem. Phys. 2003, 119, 7417-7426. 35. Zhang, Y.; Zhao, D.; Matsui, M., Anisotropy of akimotoite: A molecular dynamics study. Phys. Earth Planet. Inter. 2005, 151, 309-319. Mouraille, O.; Mulder, W. A.; Luding, S. Sound wave acceleration in granular materials. 36. J. Stat. Mech. Theory Exp. 2006, 2006, P07023. Farnell, G. W. Elastic waves in trigonal crystals. Can. J. Phys. 1961, 39, 65-80. 37. 38. McSkimin, H. J. Measurement of the 25°C zero‐field elastic moduli of quartz by high‐frequency plane‐wave propagation. J. Acoust. Soc. Am. 1962, 34, 1271-1274. 39. Brainerd, J.; Jensen, A.; Cumming, L.; Batcher, R.; Begun, S.; Black, H.; Grown, G.; Burrows, C.; Busignies, H.; Cady, W. Standards on piezoelectric crystals. Proc. IRE 1949, 37, 1378-1395. 40. Sosman, R. B., The Properties of Silica; an Introduction to the Properties of Substances in the Solid Non-Conducting State. The Chemical Catalog Company, inc.: New York, US, 1927.


26

Page 27 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


41. Bojan, V.; Alexander, P. L. Transferable force-field for modelling of CO 2 , N 2 , O 2 and Ar in all silica and Na + exchanged zeolites. Modell. Simul. Mater. Sci. Eng. 2016, 24, 045002. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A 42. smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577-8593. 43. Smith, W.; Forester, T. R. DL_POLY_2.0: A general-purpose parallel molecular dynamics simulation package. J. Mol. Graph. 1996, 14, 136-141. 44. Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635-2643. 45. Peck, R.; Devore, J. L., Statistics: The exploration & analysis of data. Duxbury Press, Belmont, US, 2011.

TOC GRAPHIC


27

Equilibrium Molecular Dynamics

Recommend Documents