A novel inversion method for the laboratory determination of Thomsen’s δ anisotropy parameter on cylindrical rock specimens from ultrasonic data has been recently reported in the literature. We further assessed this method through a direct comparison of the results of the traditional method (involving a single off-axis P-wave velocity measurement at 45°) and the new method (involving 65 P-wave velocity measurements at several angles to the symmetry axis). We prepared and characterized two vertical shale specimens from the same preserved vertical core to assess their similarity in terms of structure, mineralogy, porosity, and density. The shale was assumed to be transversely isotropic in view of the observed (horizontal) bedding. We subjected both specimens to the same brine saturation and effective stress state. Using the two methods, we obtained similar results for Thomsen’s α (vertical P-wave) and ε (P-wave anisotropy) parameters. However, a significant discrepancy was observed for Thomsen’s δ parameter: We obtained results of 0.13 using the new method and 0.39 using the traditional method. As a result of the overdetermined nature of the P-wave velocity measurements used in the new method, we believe that the corresponding δ value is more reliable. Also, the value derived with the new testing method seems to match more closely the reported field data.


Shales are extremely variable in terms of physical, mechanical, and chemical properties. The relationships between these properties and petrophysical measurements performed downhole or in the laboratory are generally not straightforward. The current understanding of the elastic properties and associated transverse isotropy (TI) of shales at realistic in situ conditions is particularly insufficient in view of the need for such information in field applications such as exploration or monitoring (e.g., 4D seismic monitoring), especially at small offset angles (normal moveout). This is due to a combination of several factors including: (1) the fact that among the scarce amount of laboratory data available in the literature, only a few experiments were performed under realistic in situ conditions, especially in terms of specimen preservation, pore-fluid saturation, and pore-pressure control; (2) limited availability of preserved samples; and (3) limited availability of laboratory measurements caused partly by the time-consuming nature of testing resulting from extremely low shale permeability (nanodarcy range) and partly by the lack of well-preserved samples.

This letter focuses on data from laboratory ultrasonic measurements acquired using a dense array of transducers attached to a single cylindrical shale specimen (Sarout et al., 2014). These data are interpreted in terms of Thomsen’s anisotropy parameters (Thomsen, 1986) using a novel inversion algorithm recently published by Nadri et al. (2012) for rocks with tilted TI. A simplified version of this algorithm is used for the specimens studied here for which vertical TI (VTI) is assumed. Special emphasis is placed on the δ parameter that controls the geometry of the propagating wavefront and plays a crucial role in seismic processing and therefore interpretation of field data (Brevik et al., 2007). The center frequency of the ultrasonic waves propagated in the specimens is of the order of 0.5 MHz, so that the wavelength is of the order of 7 mm for the measured velocities of 3.4–4.2km/s. This wavelength is significantly larger than the lamination thickness (see Figure 1). Dispersion or scattering effects are expected to be minimal.

The novel aspect of this experiment and associated data-inversion method (Nadri et al., 2012; Sarout et al., 2014), as compared with similar work, is that the determination of two Thomsen’s anisotropy parameters, ε (P-wave anisotropy) and δ, is now based on the measurement of P-wave velocity on a single shale specimen along numerous raypaths having variable inclination to the bedding plane (Figure 1). In contrast, the traditional method usually reported (Podio et al., 1968; Jones and Wang, 1981; Lo et al., 1986; Vernik and Nur, 1992; Hornby, 1998; Dewhurst and Siggins, 2006; Sarout and Guéguen, 2008; Delle Piane et al., 2011; Dewhurst et al., 2011) relies on only three P-wave velocities, involving a single off-axis propagation path inclined at 45° to the bedding plane (Figure 1).

The use of a single off-axis velocity measurement to determine δ (Figure 2) raises the question of whether this laboratory method is sufficiently accurate for direct comparison with or interpretation of field data. Although in the field P-wave velocities can be estimated along many different raypaths, in the laboratory, only a single off-bedding velocity is usually reported. The claimed large difference in the magnitude of Thomsen’s δ parameter determined in the laboratory and observed in the field (Leaney, 1994; Brevik et al., 2007) triggered the need for a better assessment of this specific anisotropy parameter at the laboratory scale. By increasing the number of independent measurements at numerous angles to bedding in the laboratory (Figure 2), the new method is expected to better probe the wavefront geometry and therefore provide a more accurate estimate of the δ parameter.

This novel approach is expected to be particularly beneficial in situations for which the minimum P-wave velocity does not coincide with the symmetry axis of the VTI shale. For instance, Tsvankin (2005) shows that in such a medium at vertical incidence (along the symmetry axis where the ray angle ϕ and phase angle θ are zero), the convexity of the P wavefront is controlled by the sign of the δ parameter; i.e., d2VP/dθ2|θ=0=2δVP(θ=0°). Negative convexity implies that the minimum P-wave velocity does not coincide with the (vertical) symmetry axis, but it is slightly offset. In the worst-case scenario, at small offset, a VTI medium with negative convexity could be mistaken with a horizontal TI medium.


We report here a direct comparison of the novel and traditional methods of estimating Thomsen parameters on the same shale lithology, with special emphasis on the estimation of Thomsen’s δ parameter. To this end, two cylindrical shale specimens are sampled a few centimeters apart from a preserved parent core oriented with their axes perpendicular to the bedding planes. The sample size is as follows: length L80mm and diameter D=38mm. Synthetic brine (3.5% NaCl) was used as a pore fluid throughout the test. Following a preliminary resaturation and equilibration stage, each specimen is subjected to a confining pressure of 50 MPa and pore pressure of 32 MPa; this isotropic loading path was conducted over a period of approximately three months to ensure saturation and pore-pressure equilibration within the low-permeability samples. Both specimens are fitted with several ultrasonic transducers directly attached to their surfaces. Specimen no. 1 is equipped with 14 P-wave transducers, whereas specimen no. 2 is fitted with the classical configuration usually reported in the literature and involving six P-wave transducers and four S-wave transducers, with a single off-axis raypath (Figure 1) allowing the measurement of five independent velocities:

  1. VP(90°): P-wave propagating parallel to the bedding,

  2. VP(0°): P-wave propagating normal to the bedding,

  3. VP(45°): P-wave propagating at 45° to the bedding,

  4. VSH(0°): S-wave propagating normal to the bedding with polarization parallel to the bedding,

  5. VSH(90°): S-wave propagating parallel to the bedding with polarization parallel to the bedding.

The estimation of Thomsen parameters for a transversely isotropic medium from the data for specimen no. 2 is based on the following equations:  
where ρ is the bulk density of the tested specimen and the indicated angles stand for the direction of propagation with respect to the symmetry axis of the TI shale (orthogonal to the bedding plane). The subscript H for the S-waves indicates the (horizontal) direction of polarization.

The stiffness tensor components Cij require phase-velocity measurements as inputs. Phase and group (ray) velocities along the TI symmetry coordinates (at 0° or 90°) coincide. However, in view of the transducer size (6 mm in diameter) compared with the propagation distances, the off-axis velocities are expected to be group (ray) velocities at a ray angle of 45° (Dellinger and Vernik, 1994). Following a method reported by Dewhurst and Siggins (2006), a correction is implemented to recover the corresponding phase velocity and phase angle, which allows for the calculation of the required phase velocity at the true phase angle of 45°.

For Specimen no. 1, a more involved inversion procedure is required, the details of which have been published previously (see Nadri et al., 2012). However, the steps of the inversion workflow are summarized in Figure 3, and the important features of the method are as follows:

  1. The inverse problem is only weakly sensitive to the magnitude of the slow S-wave velocity β=VSH(0°). However, the inversion algorithm inverts for it as one of the sought Thomsen’s parameters, and as such requires a prior input value for it (for details, see Nadri et al., 2012). The inverted value β=VSH(0°) is not expected to be reliable because no actual S-wave data from this experiment are used.

  2. Ray tracing is used so that Thomsen’s parameters are inverted for directly from the measured group (ray) velocities, so that no group-to-phase correction is required a posteriori.

  3. An estimate of the uncertainties associated with the inverted parameters is derived from the stochastic inversion procedure (see Table 1).

In the particular application of this inversion method to specimen no. 2, VTI is assumed a priori because the specimen is cored orthogonal to the shale’s bedding plane. Because no P-wave velocity measurement is available in the direction perpendicular to the bedding plane for specimen no. 1, this velocity is calculated after Thomsen’s parameters have been inverted for. An additional outcome of the new inversion method is that based on the number of input measurements available versus the number of parameters to be inverted for, an uncertainty in the determination of the anisotropy parameters is computed (for details, see Nadri et al. [2012] and the results in Table 1).


The mineralogy of the tested specimens quantified by powder x-ray diffraction comprises kaolinite (38%), mixed layer illite-smectite (23%), quartz (14%), and calcite (14%). Other minor constituents include chlorite, aragonite, and mica. The bulk density of the specimens is 2530kg/m3, and the porosity, estimated from the difference in weight between the saturated and dry specimen, is 8.9%.

As shown in Figure 1, the two neighboring samples show a similar internal structure when visualized via x-ray tomography with weak sedimentary bedding traces at a nominally orthogonal angle to the core plug axis and the occurrence of dispersed dense (bright in the figure) phases. Moreover, the density distributions along their axes are similar as shown by the CT number profile (i.e., the scaled x-ray attenuation value, which is a function of the density, atomic number, and thickness of the sample being analyzed) showing a range of values between 2350 and 2550 for both samples. On the basis of these observations, it is therefore assumed that there is reasonable homogeneity in the bulk properties of the two plugs.

The results of the two experiments are shown in Figure 4 as velocities versus angle to the symmetry axis. It shows that

  1. measured group velocities for the available ray angles in specimen no. 1 (gray dots) and in specimen no. 2 (gray squares),

  2. computed group velocities at the same ray angles in specimen no. 1 using the inverted Thomsen parameters (black dots),

  3. computed phase velocities at the corresponding phase angles in specimen no. 1 (open circles) and in specimen no. 2 (open squares).

Table 1 summarizes these results for quantitative comparison purposes at fewer important angles and in terms of Thomsen anisotropy parameters.

The results in Table 1 show a reasonably good agreement between the two methods in terms of P-wave velocity orthogonal to the bedding plane (α) and in terms of P-wave anisotropy (ε). The difference in velocity between the specimens at any given angle is lower than the accepted uncertainty for P-wave velocity measurements in the laboratory (5%; see, for instance, Delle Piane et al., 2011). Note that α has been directly measured on specimen no. 2, whereas it has been obtained by inversion for specimen no. 1. The fact that similar values of α have been obtained for both specimens suggests that the inversion for α in specimen no. 1 is reliable.

Note that the slow S-wave velocity recovered for specimen no. 1 through the new inversion method (β=VSH(0°)=2.00km/s) is not expected to be reliable, and it is indeed relatively different from the corresponding velocity measured directly on specimen no. 2 (β=VSH(0°)=2.48km/s).

A significant discrepancy between the two methods is observed for the δ parameter. It is believed that the δ value found for specimen no. 1 is a closer representation of the actual shale property due to the overdetermined nature of the P-wave velocity measurements available for this particular specimen. In Figure 5, the experimental data for specimen no. 1 (gray dots) and the computed group velocity at all ray angles for δ=0.13 (thick black line) are shown again (the same as in Figure 4). The black dot represents the computed group velocity for a ray angle of 45°, and the error bar associated with it represents a simulated experimental uncertainty of ±5%. The thin black lines represent the computed evolution of the group velocity with ray angle for various values of δ ranging between δ=0.2 and δ=0.9 (all the other parameters are taken from specimen no. 1 in Table 1).

If only the group velocity at 45° were measured for specimen no. 1, the measured value would lie somewhere along the uncertainty interval, that is, between 3.44 and 3.80km/s. Using this group velocity alone and the traditional computation method yields a value of δ ranging from 0.1 to 0.5. In contrast, using the additional independent velocity measurements at various ray angles (gray dots) and the new inversion method yields δ=0.13±0.05(0.08<δ<0.18).

Further to the above analysis, it is of interest to compare our experimental results of δ evaluation with those available in the literature. Figure 6 shows a comparison of the Thomsen δ parameter estimated for various shales in the field and in the laboratory. In particular, the plot shows a range of values obtained via the analysis of

  1. multioffset borehole seismic data for a shale formation reported in Leaney (1994),

  2. a set of wireline logs acquired in vertical and deviated wells at various deviation angles through a shale formation reported in Brevik et al. (2007),

  3. a comprehensive database compiled from many articles, conference proceedings, theses, presentations, and other public-domain sources that includes core measurements and borehole seismic measurements (Horne, 2013),

  4. the results of the present laboratory study for specimens no. 1 and no. 2.

Note that for the literature data, the mean value (dashed lines) and distribution (vertical bars) are shown. For the data set compiled by Horne (2013), a gray box shows the interval of ±1 standard deviation. It is obvious that δ for shales can span across a rather large range of values (positive and negative) depending on the specific shale lithology considered. In addition, note that the experimental value derived from specimen no. 1 with 65 raypaths falls well within the range reported by the field investigation of Brevik et al. (2007) and to a lesser extent within the range reported by Leaney (1994).

On the other hand, the value derived for specimen no. 2 with only four raypaths is coincidentally very close to the mean estimated from the data set compiled by Horne (2013).


Using a traditional and the new method, similar results for Thomsen’s α (vertical P-wave) and ε (P-wave anisotropy) parameters were obtained on the two shale specimens. However, a significant discrepancy is observed for Thomsen’s δ parameter: 0.13 using the new method in specimen no. 1 and 0.39 using the traditional method on specimen no. 2. As a result of the overdetermined nature of the P-wave velocity measurements used in the new method, we believe that the corresponding δ value is more reliable. Also, the value derived with this new method seems to match more closely the reported field data. The samples used for these experiments were brought to failure after each experiment to extract their strength properties (triaxial tests) so that repeating the velocity data acquisition using the other ultrasonic array was not possible. In addition, the arrangement of ultrasonic transducers for the new method does not allow us to carry out the traditional analysis (no raypath at 45° and 90° to the bedding). Because of the heterogeneous nature of rocks, we are currently preparing a manuscript in which we report similar results for synthetic TI samples made of phenolic resin.


Financial support from the sponsors of the SHARC research project is gratefully acknowledged (BG-Group, Chevron, ConocoPhillips, ExxonMobil, Sinopec, Statoil, and Total).

Freely available online through the SEG open-access option.