Influence of bedding structure on stress-induced elastic wave anisotropy in tight sandstones

Xiaying Li, Xinglin Lei, Qi Li,*, Dianguo Chen

a State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences, Wuhan,430071,China

b University of Chinese Academy of Sciences, Beijing,100049, China

c National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, 305-8567, Japan

d Chongqing Dizhiyuan Geology and Engineering Testing Co., Ltd, Chongqing, 401329, China

Keywords: Elastic wave anisotropy Stress-induced anisotropy Anisotropy reversal Bedding orientation Microstructure Tight sandstone

ABSTRACT To understand the evolution of stress-induced elastic wave anisotropy, three triaxial experiments were performed on sandstone specimens with bedding orientations parallel,perpendicular,and oblique to the maximum principal stress. P-wave velocities along 64 different directions on each specimen were monitored frequently to understand the anisotropy change at various stress levels by fitting Thomsen’s anisotropy equation. The results show that the elastic wave anisotropy is very sensitive to mechanical loading. Under hydrostatic loading, the magnitude of anisotropy is reduced in all three specimens.However, under deviatoric stress loading, the evolution of anisotropic characteristics (magnitude and orientation of the symmetry axis) is bedding orientation dependent. Anisotropy reversal occurs in specimens with bedding normal/oblique to the maximum principal stress. P-wave anisotropy ε′ is linearly related to volumetric strain Sv and dilatancy, indicating that stress-induced redistribution of microcracks has a significant effect on P-wave velocity anisotropy.The closure of initial cracks and pores aligned in the bedding direction contributes to the decrease of the anisotropy.However,opening of new cracks, aligned in the maximum principal direction, accounts for the increase of the anisotropy. The experimental results provide some insights into the microstructural behavior under loading and provide an experimental basis for seismic data interpretation and parameter selection in engineering applications.

1. Introduction

The anisotropy of elastic properties in most sedimentary rocks has long attracted the interest of geophysicists and geological engineers.Studies on elastic wave velocity anisotropy can be used to quantify the microstructural features of crustal rocks and explore the evolution characteristics of sedimentary sequences, which are particularly important for explanation of the seismic processing and imaging (Sayers,1999). For crustal rocks under in situ stress,the elastic wave anisotropy is generally characterized by two primary sources: intrinsic, alignment-based anisotropy and extrinsic,crack-dependent anisotropy(Barton and Quadros,2014;Allan et al.,2015).

In general, intrinsic anisotropy is a consequence of sedimentation and results from a variety of factors including mineral crystal elastic properties (Almqvist et al., 2010; Vasin et al., 2013), shape and volume of pores (Nishizawa,1982; Almqvist et al., 2011), and preferred orientation of cracks or penny-shaped pores (Anderson et al., 1974; Nishizawa and Kanagawa, 2010). Most investigations devoted to shales indicate that the alignment of clay minerals is the underlying cause for the elastic wave anisotropy (Johnston and Christensen, 1995; Sayers, 2005). However, some studies related to sandstones attributed to the presence of aligned cracks (Sayers and van Munster, 1991; Wu et al., 1991) and aligned mesoscale structures (Lei and Xue, 2009). It is widely accepted that a rock matrix embedded with randomly distributed cracks will exhibit isotropic mechanical behavior, whereas that containing aligned cracks will display an anisotropic behavior (Nishizawa and Kanagawa, 2010).

Extrinsic anisotropy is stress-sensitive and is mainly induced by the closure or opening of aligned, crack-like pores (Sarout and Guéguen, 2008a, b; Eslami et al., 2010; Allan et al., 2015; Koochak Zadeh et al., 2017; Grgic et al., 2019). Microcracks have been directly observed by Batzle et al.(1980)to close under compaction of increasing stress in rocks by a scanning electron microscope(SEM). Furthermore, cracks perpendicular to the maximum principal stress close, while those oriented parallel to the maximum principal stress remain open(Batzle et al.,1980). Recently,an SEM investigation of brittle creep experiments on shale samples performed by Geng et al.(2017)confirmed a time-dependent pressure solution, compaction localization, crack sealing/healing, and mineral rotation, which contributed to the evolution of anisotropy.

A large number of experimental studies have been concentrated on the evolution of elastic wave anisotropy in initially anisotropic rocks under the compression of confining pressure (Chan and Schmitt, 2014; Song et al., 2017) or deviatoric stress (Wu et al.,1991; Dewhurst and Siggins, 2006; Zhubayev et al., 2015; Grgic et al., 2019). Ultrasonic transmissions were performed on spherical rock samples in 132 independent directions at several levels of confining pressure (Pros et al., 2003). The change of threedimensional (3D) anisotropic elastic wave distributions in response to increasing confining pressure was recorded and compared for different rock samples. In most studied rocks, e.g.dunite and lherzolite, the anisotropic orthorhombic distribution caused by the presence of magmatic foliation and lineation at atmospheric pressure is mostly preserved even at high confining pressures,while for some rocks,e.g.pyroxenites and granulite,the distribution changes to almost isotropic or transversally isotropic.The experimental results of Li et al. (2018) performed on shale specimens confirmed that the confining pressure has a slight effect on the elastic wave anisotropy of the deformed shales, which is similar to the result of Kuila et al. (2011). However, the impact of confining pressure on the elastic anisotropy is generally ignored in field scales used in seismic data processing. By distinguishing the difference between intrinsic and extrinsic factors contributing to the elastic wave anisotropy, Allan et al. (2015) concluded that intrinsic anisotropy is the dominant source of anisotropy at all confining pressures.Extrinsic crack-based sources contribute to no more than 30% of the total anisotropy.

Deviatoric stress has been proven to have a significant impact on the anisotropy of elastic wave velocity,even in an isotropic medium(Sarout and Guéguen,2008a;Eslami et al.,2010).The elastic wave anisotropy will be significantly modified by the applied deviatoric stress. However, this change is dependent on the maximum principal stress orientation with respect to the bedding structure(Piane et al.,2011;Bonnelye et al.,2017;Wang and Li,2017;Li et al.,2018).Previous attempts to evaluate stress-induced anisotropy with increasing deviatoric stress were primarily devoted to two extreme conditions, in which the maximum principal stress orientation is either perpendicular or parallel to the bedding structure.When the maximum principal stress is bedding parallel, the elastic wave anisotropy is enhanced as the deviatoric stress increases.However,when the maximum principal stress is normal to bedding, the elastic wave anisotropy decreases; for instance stress-induced additional isotropy planes turn horizontal transversely isotropic material into seismically orthorhombic medium (Sviridov et al.,2019). Recently, a few studies were reported on anisotropy reversal behavior in shales, which were deformed under constant stress (Li et al., 2018), constant strain (Bonnelye et al., 2017) and creep loading (Geng et al., 2017) conditions. The anisotropy reversal, which indicates that the velocity perpendicular to the bedding will exceed that parallel to the bedding,occurs only in the case of deformed specimens with bedding-perpendicular loading of the deviatoric stress. Anisotropy reversal is extremely typical in the case of creep-loading experiments, in which the magnitude of the elastic wave anisotropy is twice the initial value (Geng et al.,2017).

In many sites, especially in folded formations, the maximum principal stress is neither perpendicular to nor parallel with the bedding planes. Thus, the question arises as to whether such anisotropy reversal will occur for the case in which the maximum principal stress is inclined to the bedding plane.Indeed,studies on the evolution of elastic wave anisotropy on a deformed specimen with an inclined bedding structure have been very limited.In order to understand the influence of bedding orientation on anisotropy evolution under both isotropic and anisotropic stress fields, we carried out triaxial deformation and velocity measurement experiments on tight sandstones with bedding orientations parallel with,perpendicular to, and oblique to the maximum principal stress,respectively.

2. Specimens and methods

2.1. Specimen description

The experimental sandstone specimens used in the present study were sampled from a quarry located in Xiangshi Town,Longchang County in Sichuan Province, southwestern China(Fig.1a). Fig.1a shows that the quarry is within the Sichuan Basin where the sandstone is a stretch of the Lower Shaximiao Formation in the Middle Jurassic series. The Lower Shaximiao Formation spreads widely throughout Sichuan Province and is composed of gray-green or light-green sandstone with interlayers of purple-red argillaceous sandstone or sandy mudstone. The thickness of the Lower Shaximiao Formation ranges from 500 m to 2000 m in the Sichuan Basin. The sandstone is mainly composed of medium-tofine-grained quartz or feldspar sandstone, with the quartz content ranging from 40%to 70%(Ye and Lv,2010).Fig.1b and c show the photographs of geologic section of the Lower Shaximiao Formation and an outcrop block in the quarry,respectively.As shown in Fig.1b and c, the sandstone formation is characterized by welldeveloped bedding structures.

The mineral composition of the sandstone was determined using an X-ray fluorescence spectrometer(D8 ADVANCE).The volume fraction of each mineral phase is listed in Table 1. The Lower Shaximiao sandstone is mainly composed of SiO2, with a total volume fraction of 70.38%.The porosity of the sandstone specimens varies from 9% to 12% and the permeability is approximately 0.2 mD.

In order to investigate the influence of bedding orientation on elastic wave anisotropy, we prepared three cylindrical sandstone specimens, with length of 125 mm and diameter of 50 mm, as shown in Fig.1d.Three specimens,labeled SS-00,SS-90,and SS-35,were cored parallel with,perpendicular to,and at a 35°angle with respect to the bedding plane, respectively.

2.2. Methods

2.2.1. Anisotropy model

In order to quantify the anisotropic property of rocks, several simplified velocity equations and models were developed to fit the raw velocity data. The simplest method is to calculate a P-wave anisotropy parameter AP,which is the fractional difference between the P-wave velocity parallel and perpendicular to the bedding plane (Stanchits et al., 2006; Piane et al., 2011; Allan et al., 2015):

Fig.1. (a) Map showing the sampling location in the Sichuan Basin, China. (b) Photograph of an outcrop of the lower Shaximiao formation with clear sedimentary direction. (c)Photograph of a sandstone block characterized by bedding structure.(d)Photograph of a test specimen with a diameter of 50 mm and a length of 125 mm.The GIS data are from the National Geomatics Center of China (affiliated with National Administration of Surveying, Mapping and Geoinformation of China).

where VPHand VPVis the P-wave velocity parallel and perpendicular to the bedding plane, respectively.

In general, the VPHis higher than VPV. In some cases (e.g. Pros et al., 2003), the APis calculated by the fractional difference of Pwave maximum and minimum velocities.

An elliptical model will also be used to describe the velocity distribution within the stressed samples as follows (Yanagidani et al.,1985):

where V⊥and V‖are P-wave velocities perpendicular and parallel to the axial direction of sample (maximum principal stress direction),respectively;? is the ray path angle of the measured velocity with respect to the axial direction; V（?） is the P-wave velocity along the ray path with an angle ?.Eq.(2)defines an ellipse for the reciprocal of P-wave velocity,which is named P-wave slowness.Pwave velocity fitting will also be used to describe the velocity field,but it has been found that the form of P-wave slowness shows a better fitting result than the form of the velocity itself on the elliptical dependence on propagation direction (Stanchits et al.,2003).

A relatively accurate anisotropy notation was raised by Thomsen(1986) to quantify the anisotropic properties of rock samples. The above-mentioned studies on quantitative estimation of elastic wave anisotropy are mostly based on the Thomsen anisotropy equation (Piane et al., 2011; Bonnelye et al., 2017; Li et al., 2018).Under the assumption of transverse isotropy (TI), P- and S-wave velocity along arbitrary directions can be expressed by five parameters (α, β, δ, ε, and γ), which can be conveniently converted from the elastic constants of the stiffness matrix.P-wave velocities along arbitrary direction can be approximated as follows (Mavko et al., 2009):

where φ is the angle of wave propagation orientation with respect to the symmetry axis; VP（φ） is the pseudo-longitudinal wave velocity; α is the P-wave velocities propagating along the symmetry axis;ε is a measure of P-wave anisotropy;δ is the anellipticity of the wave front, controlling the wave front shape complexity. For a special case when ε = δ, the anellipticity vanishes and the wave front becomes an ellipse,while for the other TI media with ε≠δ,the wave front will deviate from being ellipse.Note that the symmetry axis defined in the paper is the symmetry axis of the velocity curve and along the orientation of the minimum velocity.

Table 1 Mineral composition of the Lower Shaximiao specimens (as measured by X-ray fluorescence spectrometer).

Berryman (2008) extended the validity of Thomsen’s weak anisotropy formulation to stronger anisotropy and a wider range of angles without significantly affecting the equations’simplicity.The expression of P-wave velocity is as follows:

The parameter φmis a key value needed in the extended formula and can be easily deduced from the measured P- and SV-wave velocities:

where cijis element of the elastic stiffness tensor, VSVis the velocity of pseudo-shear wave.

Compared with the Thomsen’s equation,the most improvement of this extended formula is to give the angular location of the peak(or trough)of the quasi-SV-wave closer to the correct location,thus a more accurate quasi-SV-wave can be obtained.While for P-waves,the extended Berryman’s formula will sometimes yield better results at small angle than the Thomsen’s approximation and sometimes worse,depending on the actual value of φm.

For calculation of the anisotropic parameters, P-wave and Swave velocities should be measured along several different directions both in the Thomsen’s anisotropy equation and the extended Berryman’s formula. For instance, most previous studies relied on the measurement of VP（90°）, VP（0°）, VP（45°）, VSH（90°）and VSH（0°）, on three samples at 0°, 45°, and 90°to the bedding plane(Johnston and Christensen,1995;Zhubayev et al.,2015),or on a single specimen with transducers attached at these three angles(Dewhurst and Siggins, 2006; Sarout and Guéguen, 2008a; Kuila et al., 2011). However, in the laboratory, P-waves can be easily generated and recorded than S-waves which come at later time.The measurement of S-wave velocity is more difficult because it requires special S-type transmitters and receivers. Moreover, both P-and S-waves will be recorded by the S-type PZT receivers.Thus,the arrivals of S-waves are unclear and usually hidden in the coda of P-waveforms or even noise, increasing difficulty in determination of S-wave arrivals and large uncertainty of S-wave velocity.Due to the less accurate measurement of S-wave velocity,a relatively large error exists in the determination of the parameter δ. Considering the uncertainties in arriving time(±0.1 μs for P-waves and±0.4 μs for S-waves), the error of parameter δ estimated by Kuila et al.(2011) can be as high as ±33%.

In order to address this problem,a new approach was developed for estimating the Thomsen’s parameters (ε and δ), as well as the symmetry axis orientation (Li et al., 2018). In consideration of the following facts: (1) the laboratory measured P-velocities are ray velocity and the calculation of phase velocity (used in Thomsen’s equation) from ray velocity is not straightforward; and (2) the inaccurate measurement of S-wave velocity. For practical applications,we represented directly the measured ray velocities using an equation similar to the Thomsen’s anisotropy equation as follows:

where α′, δ′, ε′, and φ′have the same physical meaning with the Thomsen’s equation. In facts, as shown in Section 3, our experimental data can be represented very well. Therefore, these parameters are very close to the parameters of Thomsen’s anisotropy equation. The advantage of this method is to avoid measuring Swave velocities and only P-wave velocities are needed to be measured along several different directions on a single sample.Both the parameters(ε′and δ′)and the orientation of the symmetry axis can be directly inverted from 64 independent P-wave velocities using Eq.(7).A detailed description can be found in Li et al.(2018).By taking advantage of an array of PZTs (piezoelectric ceramic transducers) glued to a single sample, the error caused by the difference between two samples of the same lithology (or at least assumed to be in the same physical state)could be minimized.The overdetermined P-wave velocity measurements could significantly improve the accuracy and reliability of the parameter δ, as confirmed by Sarout et al. (2015).

In this paper, the same method from Li et al. (2018) will be adopted to investigate the influence of bedding orientation on the elastic wave anisotropy in the sandstone samples. Since the microfabric parameters of sandstone samples, such as different types of microcracking, micropores, mineral compositions and qualities of grain boundaries, are quite different from the shale samples studied in Li et al.(2018),similar investigation should also be performed in the sandstones. In addition, a bedding-inclined sandstone sample is added in this paper for understanding the influence of bedding orientation on stress-induced elastic wave anisotropy.

Fig. 2. Schematic of experimental apparatus for triaxial compression and velocity measurements.

2.2.2. Data acquisition and experimental procedure

Fig.3. (a)Specimen SS-90 with 16 PZTs attached to the surface in two parallel lines along axial direction.(b)Configuration of 16 PZTs and four strain gauges on the unfolded surface of the specimen.(c)64 different ray paths obtained during the velocity measurements.PZTs S1 through S8,denoted in red,are treated as sources.PZTs R1 through R8, denoted in blue, are acted as receivers.

Fig.4. (a)Waveforms recorded by receivers R1 through R8 when a pulse ultrasonic sourced from S5.The trigger time is set at 0 μs and the arrival time is indicated by the red line.(b) Fitting results for three specimens SS-00, SS-90 and SS-35, under an unstressed condition. The red arrow indicates the orientation of the symmetry axis.

Fig.2 shows a schematic of the experimental apparatus used in this study. The apparatus consists of loading controlling system,waveform recording system, and stress/strain recording system.Details of the same apparatus for velocity measurement can be found in previous studies (Lei and Xue, 2009; Li et al., 2018). In order to obtain P-wave velocities along different directions as possible,16 P-wave transducers were glued directly onto the surface of the test specimens in two parallel lines along axial direction by using an instant adhesive with a negligible thickness,as shown in Fig. 3a. The location of each PZT is shown in Fig. 3b. The PZTs,with a diameter of 5 mm and a thickness of 2 mm,have a resonant frequency of 1 MHz.Eight PZTs on the left side,denoted as S1 to S8,were treated as sources by connecting them to a pulse generator through a high-speed switch box.Eight PZTs on the right-hand side,denoted as R1 to R8,were used as receivers by feeding the recorded signals to a high-speed waveform recording system, which had a dynamic range of 16 bits and a sampling rate of up to 100 MHz.During each velocity survey, PZTs S1 to S8 were sequentially and consecutively pulsed with an ultrasonic pulse with a voltage of up to 400 V, and the remaining transducers were used to record the ultrasonic waveforms. In order to understand the microstructural characteristics within the specimens,signals recorded by PZTs R1 to R8 will be used for investigation of the anisotropy evolution.Fig.4a shows an example of waveforms recorded by receivers R1 to R8 when an ultrasonic pulse is excited from S5.The arrival time of each waveform, marked as red bars, was determined using the method of Akaike Information Criteria (AIC). The trigger time, when the ultrasonic pulse was excited from S5,was set to 0 μs.Therefore,the P-wave velocity along a selected ray path can be determined from the length of the ray path and the travel time, which is the difference between the arrival time and the trigger time. In total, 64 Pwave velocities along different ray paths could be obtained during each velocity measurement.

By fitting the measured velocity data with the anisotropy equation(Eq.(7)),the P-wave velocity α′along the symmetry axis,anisotropy parameters ε′and δ′,and orientation φ′of the symmetry axis could be extracted from the overdetermined set of independent P-wave velocities. Fig. 4b shows the fitting curves of Eq. (7)and the experimental data points for three specimens SS-00, SS-90,and SS-35 under unstressed condition.Since the symmetry axis orientations of the three specimens are different, for the sake of uniformity,the orientation of each ray path for different specimens is indicated by the angle θ between the x1axis and the ray path.Therefore,each experimental data point(θ,Vp),indicating a P-wave velocity along a given ray path, can be projected in a polar coordinate system. The blue line indicates the fitting curve of the Eq. (7). The red arrow indicates the orientation of the symmetry axis. Under the unstressed condition, all symmetry axes are along the minimal velocity orientation and normal to the bedding plane,whereas the maximum velocity orientation is along the bedding plane, which is consistent with the results of most studies (e.g.Sarout and Guéguen, 2008a; Piane et al., 2011; Li et al., 2018).

A triaxial compression experiment was performed on each specimen in order to investigate the evolution of the anisotropy parameters(α′,δ′,ε′and φ′)under a changing stress field.Before the experiment, each specimen was dried in an oven with a temperature of 60°C. Similar to 16 PZT transducers, four cross-type strain gauges, used for local axial and circumferential strain measurements,were also attached to the central surface of each specimen.In order to isolate the confining oil, a silicone sealant with a thickness of 5-6 mm was painted on the surface of both the sample and transducers.In the present study,three specimens were loaded under both isotropic and anisotropic stress fields. Specimens were first hydrostatically loaded by increasing the confining pressure from 0 MPa to 22.5 MPa. The confining pressure was then maintained constant for several minutes before axial stress loading.Deviatoric stress(σ1-σ3)was increased from 0 MPa to 100 MPa at a rate of 2.5 MPa/min(Under this axial loading rate,the creep effect is very small and will not be taken into account).During the stress loading process, velocity surveys were performed as the confining pressure increased by 2 MPa or the deviatoric stress increased by 10 MPa. In addition, the axial and circumferential strains were monitored continuously in order to characterize the evolution of the mechanical behavior of the experimental specimens under the changing stress fields.

Fig.5. (a)Average axial(Sa)and circumferential(Sc)strains,(b)P-wave velocities along different propagation directions as a function of the confining pressure in the experiment of specimen SS-00.

Fig.6. (a)Average axial(Sa),circumferential(Sc),and volumetric(Sv)strains,(b)P-wave velocities along different propagation directions as a function of the deviatoric stress in the experiment of specimen SS-00.

3. Results

3.1. Stress dependence of strain and P-wave velocity

Fig.5 shows the measured local strains and P-wave velocities of SS-00 during isotropic compression. After the confining pressure was increased to 22.5 MPa, both the axial and circumferential strains increased with a progressively decreasing strain rate, indicating a nonlinear compression of specimen volume in both the axial and circumferential directions. However, the circumferential strain is larger than the axial strain.The difference between the two strains also exhibit progressive increases with increasing confining pressure. Even if specimen SS-00 was compressed under an isotropic stress field, specimen SS-00 underwent an anisotropic mechanical compaction with larger deformation in direction normal to the bedding plane.

For simplicity,11 P-wave velocities along different directions are chosen as shown in Fig.5b.A significant difference(up to 0.3 km/s)is observed among these P-wave velocities, indicating a weak velocity anisotropic characteristic of specimen SS-00.As expected,the P-wave velocity propagating normal to the bedding (θ= 0°) is the lowest;the larger the angle of the ray path,the higher the P-wave velocity. In addition, all P-wave velocities are pressure-sensitive and increase by approximately 32% during the compression of the confining pressure.As the confining pressure increased from 0 MPa to 22.5 MPa, the P-wave velocity along the ray path of θ = -64.6°increased by approximately 0.8 km/s(from 2.99 to 3.77 km/s).The P-wave velocity along the ray path of θ=0°increased from 2.69 to 3.56 km/s, showing an increase of 0.87 km/s, which was slightly larger than those of other directions.Even if the P-wave velocity at θ=0°exhibited a larger increase,the absolute value of the P-wave velocity was still the minimum,and the relative magnitude among other velocities was unchanged under the compression of the confining pressure.

Fig. 6 shows the local strains and P-wave velocities of SS-00 specimen under a deviatoric stress up to 100 MPa. Axial strain shows a linearly increasing trend with increased deviatoric stress.Specimen SS-00 was compacted in the axial direction and dilated in the circumferential direction. Volumetric strain (Sv= Sa+ 2Sc)increased at first and then decreased. The transition position is defined as the yield point.

The evolution of the P-wave velocities due to the increasing deviatoric stress exhibits dependence on wave propagation direction.As shown in Fig.6b,the P-wave velocities along the directions with larger angles (e.g.θ = -64.6°, 61°, and ±56.4°) are increased by up to approximately 0.3 km/s as the deviatoric stress increases from 0 to 57.5 MPa. In contrast, the P-wave velocities along the directions at smaller angles (e.g. θ = 0°and ±16.8°) exhibit very small changes of less than 0.05 km/s approximately.As soon as the deviatoric stress exceeds the yield point, velocities with larger angles increase slowly or even remain unchanged, whereas those with smaller angles begin to decrease with increasing deviatoric stress. Due to the different changes in response to the increasing deviatoric stress, the difference among different velocities is significantly enhanced. Therefore, we can qualitatively infer that the velocity anisotropy of specimen SS-00 is enlarged under the compression of the deviatoric stress(see Section 3.2 for details).

Fig. 7. (a) Average axial (Sa) and circumferential (Sc) strains, (b) P-wave velocities along different propagation directions as a function of confining pressure in the experiment of specimen SS-90.

Fig. 8. (a) Average axial (Sa), circumferential (Sc) and volumetric (Sv) strains, (b) P-wave velocities along different propagation directions as a function of deviatoric stress in the experiment of specimen SS-90.

Fig. 7 shows the experimental results for local strains and Pwave velocities of specimen SS-90 under isotropic compression.As the confining pressure increased from 0 MPa to 22.5 MPa, the responses of axial and circumferential strains to the increasing pressure were similar to those for specimen SS-00, but showed a relatively small difference between the axial and circumferential strains. These two strains were approximately equal when the confining pressure was less than 10 MPa. As expected, P-wave velocities along paths of larger angles had lower values. Velocities along angles of -64.6°and 0°were increased by up to approximately 0.9 km/s and 0.7 km/s, respectively, when the confining pressure increased from 0 MPa to 22.5 MPa.

Fig. 8 shows the experimental results of specimen SS-90 under deviatoric compression. The strain-stress curves were very similar to those of specimen SS-00. Before the yield point, the specimen displayed a linearly elastic compression. P-wave velocities along larger angles (e.g. θ = ±64.6°and ±56.4°) increased more significantly than those with smaller angles (e.g. θ = ±16.8°and 0°). As specimen SS-90 entered the yield region, the P-wave velocities along larger angles continued to increase with increasing deviatoric stress, but at a lower rate. In contrast, the P-wave velocities along smaller angles began to gradually decrease with the increasing deviatoric stress.Note that P-wave the velocities along larger angles(θ = ±42.1°, ±56.4°, and ±64.6°) exceeded the ones along smaller angles(θ=±16.8°and 0°)as the deviatoric stress reached 50 MPa.In other words,P-wave velocity anisotropy reversal was induced by application of the bedding-normal deviatoric stress, and the orientation of the maximum velocity was no longer parallel with the bedding plane, which was quite different from that in the unstressed condition.In addition,as shown in Fig.8b,the gap among these P-wave velocities is gradually narrowed with the increasing deviatoric stress before the yield point, but widened as soon as specimen SS-90 enters into the yield region.

In the case of specimen SS-35, which was cored along a 35°angle to the bedding structure, the mechanical behavior was similar to those of specimens SS-00 and SS-90 under both isotropic and anisotropic stress conditions (see Figs. 9 and 10). During the loading of the confining pressure,a sudden increase in axial strain was induced by a wrong step of axial stress.The missed load was in the elastic range and thus did not damage the specimen because the axial strain was quickly reduced as soon as the axial load was released.The magnitude of the velocity increase was dependent on the propagation orientation as the confining pressure increased.Pwave velocities with negative angles (with small angles to the bedding plane) are always larger than those with positive angles(with large angles to the bedding plane). A clear gap among different P-wave velocities can also be found in specimen SS-35.During the increase of the confining pressure, P-wave velocities with angles of 31.1°and 42.1°exhibited the largest improvement and were increased by up to 1.5 km/s,whereas the ones with angles of -42.1°and -56.4°were increased by 1.3 km/s. Therefore, this gap was gradually narrowed as the confining pressure increased.

In addition, the response of the P-wave velocities to the deviatoric stress is also propagation orientation dependent. As the deviatoric stress increased from 0 MPa to 80 MPa,P-wave velocities with angles of θ = -42.1°, -56.4°, and -64.6°were increased slowly. However, P-wave velocities with angles of θ = 42.1°, 56.4°,and 61°increased more significantly, whereas the ones with smaller angles tended to decrease as soon as the specimen SS-35 entered into the yield region. The gap among different velocities tends to diverge under the applied deviatoric stress, especially in the yield region.

Fig.9. (a)Average axial(Sa)and circumferential(Sc)strains,(b)P-wave velocities along different propagation directions as a function of the confining pressure in the experiment of specimen SS-35.

3.2. Stress dependence of anisotropy

In order to quantitatively assess the anisotropic characteristic of each specimen, the measured P-wave velocities, as represented in Section 3.1, were used to invert four parameters (α′, δ′, ε′and ψ′)using a simple least-squares fitting method.Fig.11 shows the fitting curves projected in polar coordinates under the confining pressure.The points denote experimental data measured along the selected ray paths, as shown in Fig. 3c, which indicates that the curves gradually expand outward with slowly decreasing anellipticity under the compression due to the confining pressure.Fig.12 shows the fitting curves of three specimens under the applied deviatoric stress.The shapes of the velocity structures of the three specimens were significantly modified by the applied deviatoric stress. The anisotropy of specimen SS-00 was enhanced with larger anellipticity.The shapes of specimens SS-90 and SS-35 tended to become similar to that of specimen SS-00 at the end of the deviatoric stress loading.The symmetry axes(directions of minimum velocity)of specimens SS-90 and SS-35 gradually rotated to the direction of x1-axis under the compression of deviatoric stress.For details,the first quadrants of these curves in polar coordinates were magnified and shown in a way of curves of P-wave velocity with respect to the angle of the ray path. The anisotropic characteristics of the three specimens were significantly modified by the deviatoric stress. Anisotropy reversal could be easily identified in specimens SS-90 and SS-35.

Fig.10. (a) Average axial (Sa), circumferential (Sc) and volumetric (Sv) strains, (b) P-wave velocities along different propagation directions as a function of deviatoric stress in the experiment of specimen SS-35.

Fig.11. Fitting curves of the Eq. (7) projected in polar coordinates under compression of the confining pressure for specimens (a) SS-00, (b) SS-90, and (c) SS-35. The solid lines correspond to velocity structures under different confining pressures,which are inverted from Thomsen’s anisotropic equation.The points denote the experimental data obtained by the velocity measurements. The arrows indicate the orientation of the symmetry axis.

Fig.12. Fitting curves of the Eq. (7) projected in polar coordinates under the compression due to the deviatoric stress for specimens (a) SS-00, (b) SS-90, and (c) SS-35. The first quadrants of these curves are magnified on the right-hand side and are shown in a way of curves of P-wave velocity with respect to the angle of the ray path. The color bar represents the deviatoric stress.The solid lines correspond to the velocity structures inverted from the anisotropic equation.The points denote the experimental data obtained by velocity measurements. The arrows indicate the orientation of the symmetry axis.

Fig. 13. Fitting results for (a, e) the P-wave velocity, α′, along the symmetry axis, (b, f) Thomsen’s anisotropy, ε′, (c, g) the anellipticity, δ′, and (d, h) the orientation ψ′ of the symmetry axis with respect to the confining pressure and deviatoric stress, respectively.

Fig.13 shows the evolution of the anisotropic parameters as a function of the confining pressure and deviatoric stress. Under isotropic compression, the P-wave velocity, α′, along the symmetry axis for each of the three specimens increased greatly.Both the P-wave anisotropy, ε′, and the anellipticity, δ′, decreased significantly with increasing confining pressure.The angles of the symmetry axes of three specimens are almost unchanged under the isotropic compression, thus the symmetry axis keeps normal to the bedding plane for each specimen as shown in Fig. 11. In other words, the orientation of the maximum velocity is always along the bedding plane. Under anisotropic compression, the three P-wave velocities, α′, along the symmetry axis, first increased and then decreased with increasing deviatoric stress.However,the change in α′for specimen SS-90 was up to 0.4 km/s,which was much larger than the changes for the other two specimens. The P-wave anisotropy, ε′, and the anellipticity, δ′, of specimens SS-00 and SS-35 increased monotonically with increasing deviatoric stress,whereas those of specimen SS-90 first decreased and then increased. Both the minimums of ε′and δ′appeared when the deviatoric stress was approximately 37.5 MPa.In terms of the symmetry axis,the angles of the symmetry axes of specimens SS-90 and SS-35 changed from the bedding orientation to nearly 0°, whereas the angle of specimen SS-00 remained unchanged at approximately 0°.

Fig.14. (a) Relationship between the P-wave anisotropy, ε′, and the volumetric strain when the confining pressure ranges from 0 MPa to 22.5 MPa.(b)Relationship between the P-wave anisotropy, ε′, and the volumetric strain when the deviatoric stress ranges from 0 MPa to the yield stress.(c)Relationship between the P-wave anisotropy,ε′,and the dilatancy when the specimens enter the yield region.

4. Discussion

The present study highlighted the sensitivity of the P-wave velocity and anisotropy to the changing stress fields(both isotropic and anisotropic) in three specimens with different bedding orientations. Significant responses of P-wave velocity and anisotropy to the applied stress were shown in three samples. The stress sensitivity of P-wave velocity anisotropy suggested a possibility of the in situ stress field measurement from inversion of velocity in formations and rock masses. For example, crack closure induced anisotropy can be deduced from velocity measurements under hydrostatic pressure through a simple transformation, which is independent of any idealized crack geometry and limitation of low crack concentration (Mavko et al.,1995).

As mentioned earlier,such velocity response under loading was generally attributed to the change of pore space structure, e.g. the difference in closing of the preexisting cracks or the opening and expansion of new cracks aligned in a statistical sense.Fig.14 shows the relationship between the P-wave anisotropy, ε′, and the volumetric strain, Sv, or the dilatancy Δ, during different loading processes.Significant linear or polynomial function relationships were observed in three different loading stages. Therefore, the state of primary pores and induced microcracks was identified as the extrinsic sources governing the stress sensitivity of velocity and associated anisotropy.

According to the response of the P-wave velocity, as well as its anisotropy, the evolution of the internal microstructure to the changing stress field can be characterized as shown in Fig. 15. In unstressed state,three specimens demonstrate significant velocity anisotropic characteristics, as shown in Fig. 4. The maximum velocity was along the bedding plane,and the minimum velocity was normal to the bedding, indicating that most primary pores and microcracks were distributed along the bedding direction (see Fig.15a).It is because the minimum velocity generally occurs in the direction normal to the plane of cracks. Stronger anisotropy indicates a higher level of alignment along the bedding direction.

During the isotropic loading, the nonlinear increase in the volumetric strain,Sv,indicated that the mechanical behaviors of the specimens were consistent with the closure of the initial flaws including cracks and penny-shaped pores as the confining pressure ranged from 0 MPa to 22.5 MPa. The P-wave anisotropy, ε′, decreases linearly with the increasing volumetric strain,Sv,in all three specimens (see Fig. 14a). Since a larger increase of velocities occurred in the direction normal or nearly normal to the bedding plane, as compared to the velocities parallel or subparallel to the bedding, the preferential closure of cracks with planes parallel to the bedding could be induced.Therefore,most of the closed cracks and pores were inferred to be aligned in the bedding direction, as shown in Stage I (see Fig. 15). The greater the reduction of the anisotropy parameter, ε′, is, the more the aligned microcracks are closed.

Conversely, under anisotropic mechanical loading, the three specimens demonstrate a similar mechanical behavior. It can be inferred that the evolution of the P-wave velocities along different directions is a result of the combination of anisotropic compression in the axial direction before yield stress and dilatancy from microcracking activity in the circumferential direction after the yield stress. This is also indicated by acoustic emission (AE) monitoring (Lei et al., 2011; Li et al., 2016). Before the yield point, the linear increase in volumetric strain indicates that these specimens undergo linear elastic compression and thus become denser as most of the initial flaws are compressed and become closed(Stage II in Fig.15).P-wave velocities along directions with larger angles(e.g.θ = ±64.6°and ±56.4°, angle between the propagation direction and the bedding plane) increased more than that along other directions (e.g. θ = ±16.8°and 0°). A polynomial relationship between ε′and Svis more suitable for describing this linear elastic compaction process.In this stage,it can be inferred that the elastic compression of mineral grains and stiff pores likely account for most of the compression of the volumetric strain.

Fig.15. Sketch of the evolutions of the internal microstructure for three samples under the loading of both confining pressure and deviatoric stress.The nonlinear compaction refers to as the stage when the confining pressure ranges from 0 MPa to 22.5 MPa. The linear compaction is the stage when the deviatoric stress ranges from 0 to the yield point. The dilatancy is the stage when the samples enter into the yield stage.

Fig. 16. (a) P-wave velocities, PH and PV, parallel and perpendicular with the compression direction, respectively, with respect to the deviatoric stress for three samples. The values of PH and PV are calculated from the fitting curves (Fig. 13). (b)Crack densities,ΓH and ΓV,with planes parallel and perpendicular to the compression direction, respectively, with respect to the deviatoric stress for three samples.

When the axial stress exceeds the yield point, the volumetric strain deviates from a straight line θ(extension of Stage II in Fig.15)and changes from compaction to dilatancy.The corresponding yield point is in good agreement with crack initiation in the stress-strain curves. The dilatancy, Δ, is defined as the difference between the straight line θ and the volumetric strain,Sv.In this stage,velocities along larger angles increased at a decreasing rate, whereas those with smaller angles gradually began to decrease with the increasing deviatoric stress. In addition, as shown in Fig. 14c, ε′increases linearly with the increasing dilatancy of the rock volume for all three specimens, indicating that the change in velocity anisotropy results primarily from the volume dilatancy during the yield stage.

Preferred orientation of newly formed microcracks leads to the anisotropic velocity field in response to the change of axial loading.In a matrix containing aligned ellipsoidal cracks,major reduction of velocity generally occurs in the direction normal to the plane of the cracks (Anderson et al.,1974). Therefore,it is likely to characterize the stress-induced crack distribution from elastic velocity measurement.Followed by the reveal of long,straight and narrow crack induced by the dilatancy (Tapponnier and Brace, 1976), a lot of models,which are used to invert the crack characteristics from the measured P- and S-wave velocities, were developed based on the assumption of flat spheroidal cracks.Not only the crack density and orientation distribution (Soga et al.,1978; Sayers et al.,1990), but also the anisotropic permeability(Gibson Jr.and Toks?z,1990)can be estimated from the velocity inversion. Due to the lack of the Swave measurements,the model proposed by Soga et al.(1978)was used to invert the crack densities.This model has been applied for investigation of the anisotropic crack orientation with increasing stress on Basalt and Granite by Stanchits et al. (2006). The results are in good qualitative agreement with the observed acoustic emission activity and the distribution of AE source types(Stanchits et al.,2006).The equations used to invert anisotropic crack density parameter are as follows:

where ΓHand ΓVdenote the crack densities with planes parallel with and perpendicular to the compression direction,respectively;P0is the P-wave velocity of the crack-free rock; PHand PVare Pwave velocities parallel(θ=90°)with and perpendicular(θ=0°)to the compression direction, respectively. Since the maximum velocity of crack-free is unknown,the assumption that P0exceeds by 5% the maximum velocities (Stanchits et al., 2006) will also be adopted in this paper.The coefficients a1= 1.452 and a2= 0.192 were obtained from Anderson et al. (1974), which had also been used for Granite (Soga et al.,1978) and Basalt samples (Stanchits et al., 2006).

The values of PHand PVare calculated from the fitting curves(Fig.13)and shown in Fig.16a with respect to the deviatoric stress.In the three samples, the axial velocity PVincreased more significantly than the horizontal velocity PHduring the axial loading.For the SS-90 sample, the anisotropy reversal was evidently observed because PVexceeded PHas the deviatoric stress was larger than 35 MPa.The result of crack density with crack planes parallel with and perpendicular to the axial direction is shown in Fig. 16b. It shows that ΓHincreases monotonically but ΓVdecreases significantly with the increasing deviatoric stress, revealing that the density of dilatancy-induced cracks with axially parallel planes increases and further the density of cracks with axially normal planes diminishes, as shown in the Stage Ш (Fig.15). Such a scenario is similar to that for crystalized rocks (e.g. Lei et al., 1992).Therefore,it can be concluded that the axial stress closes the cracks perpendicular to the axial stress and leaves the cracks oriented parallel with the axial stress open.

5. Conclusions

The present study was devoted to the investigation of the influence of bedding orientation on the evolution of elastic anisotropy on three Lower Shaximiao sandstone specimens,which were subjected to both increasing confining pressure and deviatoric stress in the laboratory. The elastic anisotropy was determined by means of ultrasonic measurements and was quantified by Thomsen’s anisotropic parameters. In order to minimize the errors caused by specimen heterogeneity and measurement uncertainty, velocity measurements along 64 directions were performed on a single specimen in order to calculate these parameters.The measured ray velocities under both isotropic and anisotropic stress conditions are represented well by the proposed equation, even for the case of a relatively large anisotropy (ε > 0.2 in conditions of low confining pressure and high deviatoric stress). As an alteration of the exact Thomsen’s anisotropic model,the simplified model for ray velocity is of high accuracy to describe the anisotropic properties of the sandstone and is applicable in geophysical investigations such as velocity tomography,in which ray velocity is used.

The experimental results reveal that the Lower Shaximiao sandstone behaves anisotropically, which is very sensitive to mechanical loading (both isotropic and anisotropic). The increasing confining pressure contributes to a relatively larger increase in Pwave velocities in directions normal or nearly normal to the bedding plane as compared to those parallel or subparallel with the bedding and therefore accounts for the significant decrease in both ε′and δ′. However, the orientation of the symmetric axis remains normal to the bedding plane under the compression due to the confining pressure. The linear relationship between ε′and Svindicates that the mechanical behavior of the specimens is consistent with the closure of most microcracks aligned in the bedding direction.

Under deviatoric compaction, the evolution of the anisotropy parameters is dependent of bedding orientation. The anisotropic characteristics of three experimental specimens are significantly modified by the applied deviatoric stress,especially for specimens with bedding orientations normal and oblique to the maximum principal stress. The orientation of the symmetry axis changes to the direction normal to the maximum principal stress, indicating that the P-wave anisotropy characteristics of these two specimens are reversed. The stress-induced redistribution of crack directions has a significant effect on the velocity distribution, and then the velocity anisotropic characteristic. The opening of new cracks or expansion of old cracks, which are aligned in the maximum principal direction, accounts for the dilatancy of the volumetric strain.

The experimental results establish a linear connection between the velocity anisotropy and the dilatancy, and thus reveal the evolution of the microstructures of sandstone specimens with different bedding orientations under both isotropic and anisotropic compressions. However, further experiments, including microscopic structure measurement, are urgently needed in order to well characterize the stress-induced change of the microstructures.

Data availability

The data of measured P-wave velocities are online available for Figs.5-10.

Declaration of competing interest

The authors wish to confirm that there are no known conflictsof interest associated with this publication and there has beenno significant financial support for this work that could have influenced its outcome.

Acknowledgments

The research was partially supported by the National Natural Science Foundation of China (Grant Nos. 41902297, 41872210),the Natural Science Foundation of Hubei Province (Grant No.2018CFB292) and Open Research Fund of State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences (Grant No.Z017006).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2020.06.003.