Journal of Geophysical Research: Solid Earth RESEARCH ARTICLE 10.1002/2015JB012248 Key Points: • We mapped the mechanical properties of sandstones using indentation tests • Sandstones exhibit a lognormal distribution of their mechanical properties

Correspondence to: Y. Masson, [email protected]

Citation: Masson, Y., and S. R. Pride (2015), Mapping the mechanical properties of rocks using automated microindentation tests, J. Geophys. Res. Solid Earth, 120, doi:10.1002/2015JB012248.

Received 2 JUN 2015 Accepted 25 SEP 2015 Accepted article online 26 SEP 2015

Mapping the mechanical properties of rocks using automated microindentation tests Yder Masson1 and Steven R. Pride2 1 Institut de Physique du Globe de Paris, Paris, France, 2 Earth Sciences Division, Lawrence Berkeley National Laboratory,

Berkeley, California, USA

Abstract A microindentation scanner is constructed that measures the spatial fluctuation in the elastic properties of natural rocks. This novel instrument performs automated indentation tests on the surface of a rock slab and outputs 2-D maps of the indentation modulus at submillimeter resolution. Maps obtained for clean, well-consolidated, sandstone are presented and demonstrate the capabilities of the instrument. We observe that the elastic structure of sandstones correlates well with their visual appearance. Further, we show that the probability distribution of the indentation modulus fluctuations across the slab surfaces can be modeled using a lognormal probability density function. To illustrate possible use of the data obtained with the microindentation scanner, we use roughly 10 cm × 10 cm scans with millimeter resolution over four sandstone planar slabs to numerically compute the overall drained elastic moduli for each sandstone sample. We show that such numerically computed moduli are well modeled using the multicomponent form of the Hashin-Shtrikman lower bound that employs the observed lognormal probability distribution for the mesoscopic-scale moduli (the geometric mean works almost the same). We also compute the seismic attenuation versus frequency associated with wave-induced fluid flow between the heterogeneities in the scanned sandstones and observe relatively small values for the inverse quality factor (Q−1 <10−2 ) in the seismic frequency band 102 Hz < f < 104 Hz. The numerically computed frequency dependence in the attenuation varies from one type of sandstone to another, and we observe significant anisotropy in the attenuation associated with waves propagating in different directions.

1. Introduction Geological materials have heterogeneities in their poroelastic properties on different scales, ranging from pore sizes to the reservoir scale. The spatial distribution of these heterogeneities at scales smaller than core sizes, at millimeter or slightly finer resolution, is important as such heterogeneity influences the mechanical and transport properties being measured on the core. This information can be used in upscaling techniques (e.g., homogenization or effective-medium theories) to better model fluid flow [e.g., Farmer, 2002] and seismic wave propagation [e.g., Pride et al., 2004] in porous rocks.

©2015. American Geophysical Union. All Rights Reserved.

MASSON AND PRIDE

The elastic and viscoelastic properties of rocks can be obtained using either direct measurements (e.g., by applying stress on a sample and measuring the strain response) or indirect measurements (e.g., by propagating a wave through a sample and deducing the elastic moduli from secondary observables such as the wave speed). Both approaches are complimentary and can be used to obtain the viscoelastic properties of entire core samples [e.g., Spencer, 1981; Batzle et al., 2006; Adam et al., 2009; Birch and Bancroft, 1938; Winkler and Nur, 1982; Murphy, 1982; Lucet et al., 1991; Cadoret, 1993]. These macroscopic measurements at the core scale provide constraints on the seismic velocities, the anisotropy, and the attenuation that can enhance the quality of seismic images and their interpretation. However, these measurements do not necessarily identify the key physical mechanisms at work at smaller scales that are responsible for the macroscopic observation. For example, mapping the fluctuations in the elastic moduli down to a scale of a few grain diameters is important in determining whether wave-induced flow associated with the presence of heterogeneities [e.g., Pride and Berryman, 2003a, 2003b; Pride et al., 2004; White et al., 1975; White, 1975; Norris, 1993; Gurevich and Lopatnikov, 1995] is a mechanism that can explain the observed attenuation and dispersion in the seismic band (10 Hz to 104 Hz). Further, the fluctuation of elastic properties in actual rocks has not been widely measured over scales between millimeters and tens of centimeters. We call these scales “mesoscopic” because they are large THE MECHANICAL PROPERTIES OF ROCKS

1

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

enough for a porous-continuum description to be used but are smaller than seismic wavelengths. How the overall elastic moduli of a large sample are influenced by the mesoscopic-scale elastic structure it contains is largely unknown. Transmitted [e.g., Murphy et al., 1984] and refracted [e.g., Winkler and D’Angelo, 2006] ultrasonic waves have been used to produce high-resolution velocity images of rocks that provide information about their elastic structure. Microindentation and nanoindentation techniques can be used to directly measure the elastic modulus of a rock sample at a given point [e.g., Ibrahim, 2008; Wang et al., 2009; Leite and Ferland, 2001]. Zhu et al. [2007] show that nanoindentation mapping of the Young’s modulus is feasible in cement pastes and natural rocks; however, the dimensions of their experimental maps do not exceed the size of a few grains. We found no studies using indentation techniques to measure the spatial fluctuations in the elastic properties of rocks at scales ranging from a few grain diameters up to a few tens of centimeters. To access this information, we have developed an automated microindentation device that maps the indentation modulus of rocks at roughly 1 mm resolution. If a sandstone being probed is dry, each point measurement of this device is of the local (millimeter scale) drained modulus associated with the elastic response of the local framework within a volume on the order of 10 mm3 . Previously unanswered questions can be addressed with this device. What is the percent fluctuation of the drained elastic moduli over mesoscales from millimeters to tens of centimeters? In a relatively homogeneous sandstone like Berea, are the mesoscale fluctuations 10% or 1000% or somewhere in between? Do the fluctuations of the drained elastic properties over mesoscales follow a consistent statistics from one type of sandstone to the next? We will show here, for the first time, that the drained elastic moduli obey nearly perfect lognormal statistics. How are the overall drained moduli of say a 10 cm × 10 cm sample of rock related to the statistical fluctuations of the drained moduli over the smaller mesoscopic scales? In other words, which average of the mesoscale elastic structure best characterize the overall drained moduli? Answers to these questions and more will be provided in this paper. In what follows, we first present the instrument we have developed for this article that performs multiple automated indentation tests across flat rock surfaces; both the design of the instrument and the processing of the collected data are discussed. We then present maps of the indentation modulus obtained from clean, well-consolidated sandstone samples that demonstrate the capabilities of the instrument. Last, we show that the data collected can help to better understand and to better model wave propagation through porous rocks. Specifically, we show that the drained moduli of four well-consolidated sandstones are well modeled using the multicomponent form of the lower Hashin-Shtrikman bound and that despite the elastic fluctuations being on the order of 100% of the mean elastic moduli, rather low-attenuation levels associated with wave-induced flow between the mesoscale patches are computed in the seismic frequency band.

2. The Microindentation Scanner 2.1. Overview The microindentation scanner consists of an indentation probe mounted on an industrial x -y-z motion system. It performs multiple automated indentation tests at different locations on the surface of a flat rock sample. Pictures of the instrument are shown in Figure 1. The heavy-duty motion system has been chosen to be stiff enough to minimize frame deformation when performing the indentations. It is also sufficiently fast to perform high-resolution mapping in a reasonable amount of time of, say, a few hours for a rock slab several centimeters in dimension. The indentation probe of the microindentation scanner is schematized in Figure 2 and consists of an indenter mounted on a rigid rod. A simple circuit diagram of the automated device is shown in Figure 3. The force transmitted to the sample through the rod is monitored using a piezoelectric load cell, and the indentation depth (i.e., the distance between the sample’s surface and the indenter’s tip) is monitored using three displacement sensors (linear variable differential transformers or LVDTs). The piezoelectric transducer used to monitor the applied load has been chosen because it offers a faster response than the more classical load cells based on strain gages. The LVDTs have been selected for their high reliability and accuracy that offer submicron resolution. As shown in Figure 1b, three LVDTs positioned at equal angular intervals around the main rod are used to monitor the indentation depth. This reduces the noise level and prevents measurement error due to an eventual angle between the indenter and the sample’s surface due to surface roughness. MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

2

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 1. (a) A photograph of the indentation probe poised over the Rainbow sandstone slab. (b) A zoomed view showing the indenter’s tip and the bottom end of three linear variable differential transformers (LVDTs) used for monitoring the indentation depth.

The indentation probe is mounted on a linear ball bearing that is attached to the vertical axis of the motorized frame as schematized in Figure 2. A spring is placed between the block supporting the indentation tip and the block mounted on the z axis. This spring is compressed when the indentation tip is in contact with the surface of the sample. The advantage of this setup is that the servomotor driving the vertical axis can be used both to raise and to lower the indentation probe, and to control the force applied on the rock sample’s surface by the indenter. When the indenter is in contact with the sample’s surface, the applied force is directly proportional to the length of the compressed spring. The force applied to the rock sample can be controlled with a very high precision as it is proportional to the spring’s rigidity. The motion system can theoretically achieve positioning with an accuracy of about a micron. Finally, all the operations performed during the microindentation scanning are automated and controlled by a computer. This setup was custom built by the first author just for this study. 2.2. Indentation Tests The microindentation tests are performed using a loading, relaxation, unloading sequence. For each test, the indentation probe is first lowered until the indenter establishes contact with the sample’s surface. Then, the MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

3

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 2. Mechanical principle of the microindenter. The indenting probe is mounted on a rod fixed to the z axis rail of the motorized frame. (a) When the indenter tip is away from the surface of the rock sample, the probe is at rest and no force or displacement is recorded. (b) The indenting probe is then lowered until the indenter tip first touches the surface of the rock sample. (c) The applied load is increased by lowering the indenting probe until it reaches its maximum value. These operations are then performed in reverse.

force applied on the sample’s surface by the indenter is progressively increased until it reaches its maximum value (this is done by slowly lowering the z axis so the spring becomes compressed and pushes on the indentation tip). Once the applied force reaches the maximum value, the indentation probe is maintained in place for a certain period allowing for an eventual viscoelastic relaxation. Finally, the force applied is released slowly until contact between the sample’s surface and the indenter is broken (this is done by slowly raising the z axis so the spring releases the pressure on the indentation tip). Generally, during the loading stage, most materials exhibit a nonlinear elasto-plastic behavior. The relaxation stage can eliminate the viscosity impact and thereby ensure that the unloading data used for computational purposes are almost elastic. During each test, the indentation depth and applied force are recorded continuously. Once the test is done, the results may be plotted as a load-displacement curve similar to the one shown in Figure 6. These curves are numerically analyzed in real time by the instrument to compute an indentation modulus. In the next section, we detail the numerical processing employed to obtain an indentation modulus from the load-displacement curve. 2.3. Data Analysis On the assumption that the indentation data arise from a purely elastic contact (Hertzian contact), numerous models have been proposed for indentation data analysis [e.g., Sneddon, 1965; Doerner and Nix, 1986; Oliver and Pharr, 1992]. The form most often used is the one presented by Oliver and Pharr [1992] and is known as the Oliver and Pharr method. In this method, load-displacement curves are analyzed using the analytical model

Figure 3. Simplified circuit diagram of the microindenter. The arrows symbolize the direction of the data flow.

MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

4

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 4. Adapted from Oliver and Pharr [1992]: A schematic representation of a section through an indentation showing various quantities used in the analysis.

S=

𝜕P 2 √ = 𝛽 √ Er A, 𝜕h 𝜋

(1)

where the contact stiffness S = 𝜕P∕𝜕h is the slope of the tangent to the unloading part of the loaddisplacement curve at the maximum load (as illustrated in Figure 5), h is depth of the indenter tip below the surface, A is the contact area, and the reduced indentation modulus Er is an averaged modulus that accounts for the effects of using a nonrigid indenter. In this study, the dimensionless parameter 𝛽 is taken as unity because only axisymmetric indenters are considered. We refer the reader to Oliver and Pharr [2004] for a discussion on appropriate values of the 𝛽 parameter if nonaxisymmetric indenters are employed. We have 2 1 1 − 𝜈 2 1 − 𝜈i + = Er E Ei

(2)

where 𝜈 and E are the Poisson’s ratio and Young’s modulus of the rock sample, and 𝜈i and Ei are the Poisson’s ratio and the Young’s modulus of the indenter. Equation (1) was first introduced by Sneddon [1948] for contact between a rigid indenter of defined shape and a homogeneous isotropic elastic half-space. Later, Pharr et al. [1992] showed that this equation is robust and

Figure 5. Adapted from Oliver and Pharr [1992]: A schematic representation of load versus indenter displacement data for an indentation experiment. The quantities shown are hmax , the indenter displacement at peak load; hs , the projected indentation depth; hf , the final depth of the contact impression after unloading; and S, the initial unloading stiffness. A graphical interpretation of the contact depth is presented as well.

MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

5

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

applies to tips with a wide range of shapes. Sneddon’s contact solution predicts that the unloading data for an elastic contact for many simple axisymmetric indenter geometries (sphere, cone, flat punch, and paraboloids of revolution) follow the power law )m ( P(h) = 𝛼 h − hf , (3) where P is the applied load, h is the elastic displacement of the indenter, hf is the final elastic displacement of the indenter as defined in Figure 5, and 𝛼 and m are fitting constants. Oliver and Pharr [1992] use this model to fit the unloading part of the load-displacement curve and to obtain the indentation stiffness )m−1 ( 𝜕P(h) || Smax = = 𝛼 m hmax − hf . (4) 𝜕h ||h=hmax To estimate the projected contact area A, the tip shape of the indenter must be known accurately so A can be expressed as a function of the indentation depth hc A = f (hc ), (5) where the function is the usual expression for the surface area of either a cone, a spherical lens, or a cylinder (flat punch) of height hc . Oliver and Pharr derive the following relationship for the contact depth hc from Sneddon’s solutions (see Figure 5 for definition of these depths) P hc = hmax − hs = hmax − 𝜖 max , (6) Smax where 𝜖 = 0.72, 0.75, and 1, for cone-, sphere-, and flat-punch geometries, respectively.

In summary, the Oliver and Pharr method consists of first fitting a power law function to the unloading segment of the load-displacement curve. Taking the slope of this function at maximum load yields the contact stiffness Er . This contact stiffness in addition to the appropriate value of 𝜖 is used in order to determine the actual contact depth hc so that it is finally possible to derive the indentation modulus Es that we now define as 1 1 − 𝜈2 . = Es E

(7)

Figures 4 and 5 show a schematic sketch of such an analysis. Following the Oliver and Pharr method, the microindenter performs the following processing for each indentation test: 1. Record load and displacement versus time during the indentation cycle. 2. Preprocess the data: filter, remove outliers, and extract the unloading part of the load-displacement curve (an example of a load-displacement curve is pictured in Figure 6). 3. Perform a least squares fit (power law) of the unloading part of the load-displacement curve to determine the parameters 𝛼 and m in equation (3). 4. Compute the slope Smax of the tangent to the unloading part of the load-displacement curve at the maximum displacement point using equation (4). 5. Use equation (6) to obtain the indentation depth hc and then equation (5) to obtain the contact area A. 6. Compute Er using equation (1). All the steps above are performed in real time by the instrument and only the values of the reduced modulus Er are written in the output files. The value of the indenter’s modulus Ei needed to evaluate Es is obtained independently when calibrating the instrument as described in the next section. 2.4. Measurement Error and Calibration The calibration of the indenter consists of determining the ratio (1 − 𝜈i2 )∕Ei . This can be achieved by using equation (2) where the reduced modulus Er is measured for a material with known Young’s modulus E = Eref and Poisson’s ratio 𝜈 = 𝜈ref . Though the ratio (1 − 𝜈i2 )∕Ei can be estimated using a single test on a given material, to reduce error and check for consistency in the results, we perform numerous estimates using different materials with Young’s moduli ranging from 3 GPa up to 70 GPa (e.g., glass, plexiglass, copper, and aluminum). An estimate of the measurement error can be obtained after calibration by mapping the indentation modulus Es for homogeneous samples with known Young’s modulus Eref and Poisson’s ratio 𝜈 = 𝜈ref . For example, we calculate the standard deviation from the reference value using √ √ N √1 ∑( )2 E(xi ) − Eref , 𝜎=√ (8) N i=1 MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

6

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

where N is the number of measurements in the map, xi is the spatial position of the ith measurement, and E(xi ) is the measured value at location xi . Using this approach, we estimate the standard deviation of the measurement error to be less than 2.5% of the reference value for a range of materials having widely varying Young’s moduli.

Figure 6. An example of an actual force-displacement curve obtained using a well-consolidated sandstone sample.

When the indenter tip is maximally compressed into a point on the rock surface, the radius of the contact area is amax . As shown by the Sneddon [1948] elastic solution of the indentation problem, the volume of rock that dominantly controls the elastic unloading that is being measured lies within a few amax of the contact area in both the horizontal and vertical directions. Crudely, one could take (4amax )3 as a reasonable estimate of the rock volume being sampled during the unloading portion of the indentation test.

3. Experiments on Four Sandstone Samples In Figures 7–10, we present maps of the indentation modulus Es obtained with the microindentation scanner for four different sandstone samples. Before we performed these measurements, in order to determine what is best adapted for these types of rocks, indenters with different shapes and made of different materials were tested: conical (diamond), spherical (carbide), and flat(steel). In the final setup used for all the maps presented in the present paper, a flat indenter was used because it induces the least amount of damage to a sandstone sample’s surface. It also offers the advantage of not having to determine the contact area from the load-displacement curve (hc above), which reduces uncertainty in the estimated indentation modulus. The maximum load applied for the indentation tests has been chosen based on two criterion: (1) it has to be large enough so that the indenter’s tip establishes full contact with the sample’s surface; and (2) it has to be as small as possible to minimize damage in the sample. Based on these two constraints, we performed multiple tests and established that indentation loads ranging from 50N to 150N may be considered optimal for the sandstone samples we worked with. This roughly corresponds to indentation depths of tens of microns for the samples considered in this section (see Figure 6). We also made sure that when repeating indentation tests at a given location with the aforementioned maximum indentation load, the measured value for the indentation modulus did not vary significantly. This implies that additional damage induced by repeated indentation tests is minimal. All the samples presented in this section have been cut through bulk rocks using a diamond saw, and no specific treatment or polishing has been applied to the surfaces mapped. All the maps were performed within the central part of the samples, staying roughly 1 cm away from the edges of the slab to avoid potential side effects. The total time needed to complete each map did not exceed 48 h for a maximum of 40,000 tests (or roughly 4 s per test). In the analysis that produces each point measurement of indentation modulus, the rocks are taken to be locally homogeneous over scales of twice the diameter of maximum contact 4amax . For the cylindrically shaped indenter used in this work, we have 2amax = 1 mm. Therefore, the maps presented in this section show the fluctuations of a local effective modulus associated with a local rock volume of roughly 8 mm3 and should be understood as such. We repeat that the indenter is only pushed into each such volume a distance of roughly 0.01 mm. Further, the grain diameters in our sandstones are on the order of 0.1 to 0.3 mm, so the cylindrical indenter is in contact with between 10 and 100 sand grains. MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

7

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 7. (a) Photograph of the Berea sandstone slab of dimensions 8.5 cm × 7.5 cm and 5 cm thick. (b) Map showing the spatial fluctuation in the indentation modulus Es . The spatial sampling interval is 1 mm in both the x and y directions. (c) Probability density histogram showing the indentation modulus distribution as well as the corresponding distribution for K and G, assuming a Poisson’s ratio 𝜈 = 1∕8. Because Young’s modulus of the sample E is almost equal to the indentation modulus Es , we have E = (1 − 𝜈 2 )Es ≈ 0.98Es . The solid blue line is the lognormal probability density function in equation (9) where the parameters 𝜇 , 𝜎p , and x0 have been adjusted to fit the observed histogram.

The first sample, photographed in Figure 7a, has been cut through a Berea sandstone drill core. This Berea is a fine-grained, well-sorted grayish sandstone with closely spaced planar bedding. The second sample, photographed in Figure 8a, has a yellowish-beige color, coarser grains of varying sizes, and exhibits some curved bedding patterns. In addition, this sample has two horizontal fractures that span the entire mapped area and that are roughly perpendicular to the sample surface. To prevent the sample from breaking during the indentation tests, it was first cased into a plaster shell and sliced afterward. The third sample, photographed in Figure 9a, is a Rainbow sandstone paving tile purchased at a local tile dealer. It is a very fine-grained and well-sorted sandstone with bright colors and a quasi-planar visual bedding. The fourth sample, photographed in Figure 10a, is a beige Massillon sandstone cut through a drill core. It has coarse grains of various sizes and exhibits cross bedding. The physical dimensions of each sandstone slab are given in the captions to these figures. The indentation maps obtained are plotted in Figures 7b, 8b, 9b, and 10b. In these maps, at most 0.5% of the measurements were considered as outliers and replaced by the average values of measurements performed at neighboring locations. To the first order, we observe that the fluctuation in the indentation modulus matches the visual characteristics of the sample. However, the indentation maps appear to be noisier at smaller scales compared to the photographs. Part of this high-frequency noise can be attributed to measurement error. However, given that our estimate of the uncertainty in the measured values is no more than a few percent error, we think that most of the elastic-property fluctuation observed at small scales is related to the rock’s structure. Our interpretation of the fluctuation observed in these samples is that the large-scale fluctuations that correlate with the rock’s bedding are due to changes in the mix of components forming the grain skeleton MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

8

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 8. (a) Photograph of the fractured sandstone sample of dimensions 6.5 cm × 13.5 cm and 3.5 cm thick. (b) Map showing the spatial fluctuation in the indentation modulus Es . The spatial sampling interval is 1 mm in both the x and y directions. (c) Probability density histogram showing the indentation modulus distribution as well as the corresponding distribution for K and G, assuming a Poisson’s ratio 𝜈 = 1∕8. The solid blue line is the lognormal probability density function in equation (9) where the parameters 𝜇 , 𝜎p , and x0 have been adjusted to fit the observed histogram.

(grain size, mineralogy, etc.), while the small-scale fluctuations are more likely due to variations in the way these components are locally and randomly packed and cemented together. In Figure 7b, we observe that the elastic structure of the sample correlates very well with the rock’s bedding. Within layers, the indentation values are relatively homogeneous and the fluctuation around the mean value is likely due to the random nature of grain packing. The oblique stripes are artifacts and correspond to visible saw marks at the surface of the sample. These may be due to topography or damage induced by the sawing. The computed Young’s modulus is E ≈ 16.6 GPa which is consistent with available values obtained using standard methods that lie in the range (12.4 GPa–16.7 GPa) at 0 MPa confining pressure according to American Society for Testing and Materials Standard D7012-14 [Anonymous, 2004]. In Figure 8b, we observe that the two horizontal fractures at heights 3 cm and 11 cm appear as thick lines with low indentation modulus. Notice that as the indenter’s tip is significantly larger than the thickness of the fractures, the measured values for the indentation modulus correspond to the effective properties of the sample in the vicinity of the fractures and not to the actual properties of the fractures themselves. It is interesting to note that the thickness of the weak area around the fractures is significantly larger than the indenter’s tip diameter (i.e., the material remains weak for an additional 1 or 2 mm on either side of the fracture surface). This suggests that as the fracture broke initially, it weakened or broke grain contacts in a small 1 or 2 mm neighborhood perpendicular to the fracture plane. It is also interesting to observe that the fracture at 3 cm height separates layers with contrasting indentation modulus. We emphasize that the indentation modulus values measured on and around these fractures should depend strongly on the angle between the cut plane and the fracture plan. We would expect much lower moduli in cases where the indenter pushes in a direction MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

9

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 9. (a) Photograph of the Rainbow sandstone sample of dimensions 20 cm × 20 cm and 1.5 cm thick. (b) Map showing the spatial fluctuation in the indentation modulus Es . The spatial sampling interval is 1 mm in both the x and y directions. (c) Probability density histogram showing the indentation modulus distribution as well as the corresponding distribution for K and G, assuming a Poisson’s ratio 𝜈 = 1∕8. The solid blue line is the lognormal probability density function in equation (9) where the parameters 𝜇 , 𝜎p , and x0 have been adjusted to fit the observed histogram.

with a significant component that is perpendicular to the fracture plane. In this sample, the fractures are more or less perpendicular to the sample surface and thus have a relatively strong measured indentation modulus. This has implications when we numerically compute attenuation for this sample. Last, when looking at the bottom layer in Figure 8b, we observe centimeter-scale patches of stronger material; these visually correlate with areas exhibiting some degree of enhanced quartz cementation. In Figure 9b, most of the fluctuation in the indentation modulus occurs at small scales. This is consistent with the fact that the grain size is very homogeneous throughout the whole sample, and it is mostly, but not entirely, the variation in pigmentation (iron content) that gives the sample its layered visual appearance. In Figure 10b also, most of the fluctuation in the indentation modulus is observed at small scales which correlates with the heterogeneous nature of the Massillon sandstone at grain scales. Large-scale fluctuations are nonetheless visible and reflect the rock’s bedding. The computed Young’s modulus is E ≈ 11GPa which is a little smaller than the value E ≈ 15 GPa [Lumley et al., 1992] obtained by measuring compressional and shear ultrasonic velocities. The probability density histograms in Figures 7c, 8c, 9c, and 10c share a bell shape with significant positive skew. The solid blue lines in Figures 7c, 8c, 9c, and 10c are lognormal probability density functions (PDFs) of the form ( ( )2 ) (ln x − x0 ) − 𝜇 1 , f (x) = (9) exp − √ 2𝜎p2 (x − x )𝜎 2𝜋 0

p

where the parameters 𝜎p , 𝜇 , and x0 have been adjusted to match the data. We observe that the lognormal PDF fits the histogram estimates almost perfectly. This finding honors the late Albert Tarantola who promoted the use of lognormal PDFs to describe the fluctuation in the elastic properties of rocks [e.g., Mosegaard and MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

10

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 10. (a) Photograph of the Massillon sandstone core sample of circular diameter 10 cm and 1 cm thick. (b) Map showing the spatial fluctuation in the indentation modulus Es . The spatial sampling interval is 1 mm in both the horizontal and the vertical directions. (c) Probability density histogram showing the indentation modulus distribution as well as the corresponding distribution for K and G, assuming a Poisson’s ratio 𝜈 = 1∕8. The solid blue line is the lognormal probability density function in equation (9) where the parameters 𝜇 , 𝜎p , and x0 have been adjusted to fit the observed histogram.

Tarantola, 2002; Tarantola, 2005]. He did so primarily because the lognormal distribution requires the property being measured to be positive which is the case for most rock properties. The central limit theorem further tells us that lognormal distributions arise when what is being measured is the result of the product of random variables. If, for example, the indentation modulus is well modeled as the geometric mean of the underlying grain-contact strengths, a lognormal PDF for the indentation modulus should be expected. We will not speculate further on why lognormal distributions are being measured but it is a key finding of this study.

4. Numerical Computation of Each Sample’s Complex Elastic Properties In this section, we use the 2-D maps of the indentation modulus to numerically simulate the 2-D response of each sample to time-varying displacements applied normal to the sample boundaries. The local stress and strain rate are computed throughout each sample’s 2-D plane using the explicit finite-difference algorithm of Masson et al. [2006] that solves Biot’s poroelasticity equations on a staggered grid. The spatially averaged stress tensor and strain-rate tensor are then determined as a function of time. A temporal Fourier transform is taken next to give both the complex frequency-dependent stress tensor 𝜎ij (𝜔) and strain tensor 𝜖ij (𝜔) for each sample. The details of this modeling algorithm are given by Masson and Pride [2007]. From these simulations, we compute the five complex and frequency-dependent elastic moduli M11 (𝜔), M22 (𝜔), M12 (𝜔), M21 (𝜔), and G(𝜔) using the relation ⎡ 𝜎xx ⎤ ⎡ M11 M12 ⎤ ⎡ 𝜖xx ⎤ ⎥⎢ 𝜖 ⎥. ⎢𝜎 ⎥=⎢M M ⎢ yy ⎥ ⎢ 21 22 ⎥ ⎢ yy ⎥ G ⎦ ⎣ 𝜖xy ⎦ ⎣ 𝜎xy ⎦ ⎣

MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

(10)

11

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Table 1. Hydro-Physical Properties Assumed When Computing the Seismic Attenuations Presented in Figure 11a Solid Grain Material Bulk modulus Ks (GPa) 36.0 Density 𝜌s (kg/m3 ) 2650 Skeletal Framework of Grains Poisson ratio 𝜈 Massillon

Berea

Rainbow

Fractured

0.125

0.125

0.125

0.125

Berea

Rainbow

Fractured

0.2

0.3

0.3

Porosity 𝜙 Massillon 0.22 Permeability k

(m2 )

Massillon 0.6 × 10−12

Berea

Rainbow

Fractured

0.4 × 10−12

1 × 10−12

1 × 10−12

Fluid Bulk modulus Kf (GPa) 2.25 Density 𝜌f (kg/m3 ) 1000 Viscosity 𝜂 (N s m−2 ) 10−3

a K is the bulk modulus of the solid grain material, 𝜌 is the density of the solid grain material, 𝜈 is the Poisson ratio of s s the skeletal framework of grains, 𝜙 is the porosity of the skeletal framework of grains, k is the permeability of the skeletal framework of grains, Kf is the bulk modulus of the fluid, 𝜌f is the density of the fluid, and 𝜂 is the viscosity of the fluid.

As will be seen, only small amounts of anisotropy (less than 1%) are observed in the overall elastic moduli; however, the attenuation exhibits more anisotropy. In order to perform the 2-D numerical experiments, we must translate the maps of indentation modulus into local estimates of the poroelastic moduli which, through the Biot and Willis [1957] relations, require maps of the local drained bulk moduli K and shear moduli G over the grid. To do so, we assume that at each grid point, the drained Poisson ratio is 𝜈 = 1∕8 which means that K = G at all points in our simulation domain. In an experimental study, Castagna et al. [1985] show that the drained bulk modulus roughly equals the shear modulus for all 30 sandstones in their study. In this manner, we simply map the indentation modulus in Figures 7–10 that were taken on dry samples into equivalent maps of drained bulk and shear frame moduli (that are identical due to the 𝜈 = 1∕8 assumption). We assume the solid grains are made of quartz (solid-mineral bulk modulus Ks = 36 GPa and density 𝜌s = 2650 Kg/m3 ) for all samples. The local porosity and permeability were not directly measured; we instead took their values to be uniform across each sample and equal to the overall porosity and permeability level as measured on lab cores by others. The overall elastic moduli and levels of peak attenuation being computed are independent of permeability, which only influences the timing of fluid-pressure equilibration and thus the frequency dependence of dispersion and attenuation. In the future, we can imagine making maps of porosity and permeability using existing equipment such as an X-ray scanner for the porosity and a micropermeameter for the permeability. The values of the different parameters assumed when computing the complex frequency-dependent moduli of equation (10) are given in Table 1. 4.1. Numerical Estimates of the Drained Elastic Moduli of the Sandstone Samples In this section, we use the above procedure to numerically compute the overall drained elastic moduli for each sample and compare these computational results to a range of simple averages that use the observed probability distributions of the local elastic moduli in an estimate of the overall elastic moduli. MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

12

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

To obtain the drained Mij of equation (10), we use a pore fluid that has fluid properties that are 106 smaller than that of a liquid. Specifically, we use a fluid bulk modulus of Kf = 103 Pa, a fluid viscosity of 𝜂 = 10−9 Pa s, and fluid density of 𝜌f = 10−3 kg/m3 . Due to the tiny fluid bulk modulus, the fluid pressure does not noticeably change when the rock is deformed which means that no induced flow occurs and that the imaginary part of the measured moduli Mij is negligibly small. Nonetheless, we evaluate the real part of each complex Mij at the low end of the frequency spectrum (at 100 Hz). To obtain an isotropic drained bulk modulus K2D for the entire sample from the computed drained Mij , we imagine an experiment where 𝜖xx = 𝜖yy and 𝜖xy = 0 (pure compression) and use the definition (1,2) K2D =

𝜎xx + 𝜎yy

(11)

2(𝜖xx + 𝜖yy )

along with equation (10) to obtain (1,2) K2D =

M(1) + M(2) + M(1) + M(2) 11 22 12 21 4

.

(12)

The superscript (1, 2) represents that two numerical computations need to be performed to determine the Mij : (1) 𝜖xx ≠ 0 and 𝜖yy = 0 to determine M(1) and M(1) ; and (2) 𝜖xx = 0 and 𝜖yy ≠ 0 to determine M(2) and 11 12 22 (2) (1,2) (2) (1) (2) M21 . Similarly, to obtain an isotropic shear modulus G2D from knowledge of M(1) , M , M , and M for the 11 22 12 21 sample, we imagine an experiment where 𝜖yy = −𝜖xx and 𝜖xy = 0 (pure shear) and use the definition G(1,2) = 2D

𝜎xx − 𝜎yy

(13)

2(𝜖xx − 𝜖yy )

along with equation (10) to obtain G(1,2) = 2D

+ M(2) − M(1) − M(2) M(1) 11 22 12 21 4

.

The other way (3) to compute G2D is to apply 𝜖xy ≠ 0 along with 𝜖xx = 𝜖yy = 0 to obtain G(3) = 2D

𝜎xy 𝜖xy

.

(14)

(15)

Due to numerical error, we have that G(1,2) differs from G(3) by a couple of percent and so we elect to use for 2D 2D our estimate of the shear modulus G(1,2) + G(3) 2D G2D = 2D . (16) 2 Finally, to obtain estimates of the desired 3-D moduli, we use Hooke’s law for an isotropic 3-D material and assume that our 2-D numerical experiments correspond to plane strain 𝜖zz = 0 in a 3-D sample to obtain (1,2) K 3D = K2D − G3D ∕3

(17)

G3D = G2D .

(18)

The values for all the above moduli computed for each sandstone sample in this study are given in Table 2. Note that M11 and M22 are less than 1% different and that each rock’s elasticity is approximately isotropic despite the visual layering in the samples. Also note that M12 and M21 differ by a couple of percent. Thermodynamic arguments (Maxwell relations) require that M12 = M21 . The difference we observe seems related to error in the numerical experiments involved with how the spatially fluctuating shear moduli are used to update the local stress components on the staggered grid. We we were able to reduce the difference between M12 and M21 by reducing the grid size below the resolution of the indentation maps. We also observe that K3D is not very different from the corresponding G3D . This may not seem surprising since we assumed each local grid point to have K = G. However, in a heterogeneous sample, even if K = G locally, that by itself does not guarantee that a pure shear experiment on a heterogeneous sample will yield the same overall modulus as a pure compression experiment. Details of how the heterogeneity is put together can make the two macromoduli different. MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

13

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Table 2. The Numerically Computed Modulia (1)

(1)

(2)

M11

M12

M21

(2)

(1,2)

(1,2)

(3)

M22

K2D

G2D

G2D

G2D

(1,2)

K3D

K 3D

(GPa)

(GPa)

(GPa)

(GPa)

(GPa)

(GPa)

(GPa)

(GPa)

(GPa)

Massillon

11.38

1.68

1.63

11.32

6.50

4.85

4.78

4.82

4.88

(GPa) 4.89

Berea

17.26

2.49

2.38

17.08

9.80

7.37

7.23

7.30

7.34

7.37

Fractured

5.05

0.72

0.73

4.91

2.85

2.13

2.09

2.11

2.14

2.15

Rainbow

40.05

5.88

5.46

39.97

22.84

17.17

16.77

16.97

17.12

17.18

a The superscript denotes the experiment used to compute the moduli. Numerical experiment(1) is applying 𝜖 ≠ 0 xx with 𝜖yy = 0 and 𝜖xy = 0. Experiment (2) is applying 𝜖yy ≠ 0 with 𝜖xx = 0 and 𝜖xy = 0. Experiment (3) is applying (1,2) (1) (2) (1) (2) (1,2) (1) (2) (1) (2) 𝜖xy ≠ 0 with 𝜖xx = 0 and 𝜖yy = 0. We have K2D = (M11 + M22 + M12 + M21 )∕4, G2D = (M11 + M22 − M12 − M21 )∕4, (1,2) (3) (1,2) (1,2) (1,2) (1,2) G2D = (G2D + G2D )∕2, K3D = K2D − G2D ∕3, and K 3D = K2D − G2D ∕3.

Given the measured lognormal probability distribution function fM (M) for a local elastic modulus M, where we consider either M = K or M = G but have fK (K) = fG (G) due to the assumption that 𝜈 = 1∕8, we can use a range of standard averages and bounds to see which average does the best job predicting the numerically computed elastic moduli K 3D and G3D of each sandstone sample. Note that fM (M) is normalized so that ∫ fM (M) dM = 1. The arithmetic average is defined ⟨M⟩a =

The harmonic average is defined

∫

fM (M)M dM.

(19)

fM (M) 1 = dM. ∫ ⟨M⟩h M

(20)

The geometric average is defined ln⟨M⟩g =

∫

fM (M) ln M dM.

(21)

∑N If we write the integral over the finite range of M as a sum ∫ fM (M) ln M dM = ΔM i=1 fM (Mi ) ln Mi where Mi+1 = Mi + ΔM and M1 is the smallest modulus in the distribution and MN the largest, then upon taking the exponential we have N ∏ f (M )ΔM ⟨M⟩g = Mi M i . (22) i=1

Thus, the geometric mean involves the product of the mesoscale moduli in the sample as raised to the appropriate exponent that reflects the volumetric probability of each modulus in the sample. The lower (HS−) and upper (HS+) Hashin-Shtrikman bounds in d dimensions for a multicomponent mixture of isotropic solids are given by Torquato [2002] [ ]−1 fM (M) ⟨M⟩HS− = dM − M∗min (23) ∫ M∗min + M and

[ ⟨M⟩HS+ =

fM (M) dM ∫ M∗max + M

]−1

− M∗max ,

(24)

where for M = K, G

MASSON AND PRIDE

∗ = Kmin

2(d − 1) Gmin , d

(25)

∗ Kmax =

2(d − 1) Gmax , d

(26)

G∗min = Hmin ,

(27)

G∗max = Hmax .

(28)

THE MECHANICAL PROPERTIES OF ROCKS

14

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Table 3. The Numerically Computed Value of K 3D as Given in Table 2 for Each Sample’s Overall Bulk Modulus Compared to Various Possible Averagesa Bulk Modulus K 3D (GPa) Hashin

Hashin

Shtrikman

Shtrikman

Hashin

Simple

Harmonic

Geometric

Lower

Upper

Shtrikman

Computed

Mean

Mean

Mean

Bound

Bound

Mean

Massillon

4.89

5.21

4.70

4.93

4.87

5.08

4.98

Berea

7.37

7.74

7.13

7.42

7.33

7.58

7.46

Fractured

2.15

2.24

2.07

2.15

2.12

2.20

2.16

Rainbow

17.18

18.53

16.46

17.38

17.11

17.97

17.54

a The estimate in bold italics is the one that came closest to the computed value.

Here Gmin is the smallest shear modulus in the distribution and Gmax is the largest. Similarly, Hmin, max are the smallest and largest values in the distribution of the elastic modulus H defined as [ H=G

] dK∕2 + (d + 1)(d − 2)G∕d . K + 2G

(29)

We are interested in 3-D moduli and take d = 3. The results of the above averages are given in Tables 3 and 4. We observe that the lower Hashin-Shtrikman bound does the best overall job in matching the calculated values of K 3D and G3D ; however, the geometric mean also does an excellent job in matching the calculated values. To within the roughly 2% error of the numerical computations, one can say that the geometric mean and the lower Hashin-Shtrikman bound do an equally good job predicting the calculated bulk and shear moduli. Note that although we have K = G locally throughout our samples, the Hashin-Shtrikman bounds for the overall K and G differ from each other due to the different local response in pure compression and pure shear experiments involving heterogeneous composites. 4.2. Numerical Estimates of the Seismic Attenuation of the Sandstone Samples When seismic waves are traveling through saturated porous materials, they locally compress or dilate the propagating medium and induce a heterogeneous fluid-pressure response that correlates with the incompressibility structure present. These wave-induced fluid-pressure gradients equilibrate via pore pressure diffusion and result in viscous loss [e.g., Pride and Berryman, 2003a, 2003b; Pride et al., 2004; White et al., 1975; White, 1975; Norris, 1993; Gurevich and Lopatnikov, 1995]. When heterogeneities are all of the same size, attenuation, as measured by the inverse quality factor Q−1 , is at a maximum when the fluid pressure has just enough time in a wave period to diffuse across the heterogeneous patches present. Peak attenuation thus occurs at a frequency fc that scales as k∕L2 where k is the local permeability and L is the characteristic size Table 4. The Numerically Computed Value of G3D = G2D Compared to the Various Possible Averagesa Shear Modulus G3D (GPa)

Computed

Hashin

Hashin

Shtrikman

Shtrikman

Hashin

Simple

Harmonic

Geometric

Lower

Upper

Shtrikman

Mean

Mean

Mean

Bound

Bound

Mean

Massillon

4.82

5.21

4.70

4.93

4.84

5.05

4.94

Berea

7.30

7.74

7.13

7.42

7.29

7.54

7.41

Fractured

2.11

2.24

2.07

2.15

2.11

2.18

2.15

Rainbow

16.97

18.53

16.46

17.38

16.97

17.82

17.40

a The estimate in bold italics is the one that came closest to the computed value.

MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

15

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Figure 11. Attenuation versus frequency computed following the method presented in Masson and Pride [2007] for the four sandstones in this study and using the indentation modulus map in Figures 7–10 as input.

of the patches. Further, the maximum attenuation Q−1 (fc ) is proportional to the square of the amplitude ΔM of the fluctuation in the incompressibility M with respect to its average value ⟨M⟩; i.e., Q−1 ∝ (ΔM)2 [Masson and Pride, 2007]. Last, the shape of the heterogeneous patches influences the anisotropy in the attenuation [e.g., Masson and Pride, 2014]. In the general case of a material containing an arbitrary amount of heterogeneous patches having various sizes and shapes, the attenuation and dispersion versus frequency can be obtained using numerical methods [e.g., Masson and Pride, 2007; Quintal et al., 2011]. In this section, again following the method introduced in Masson and Pride [2007], we compute the attenuation due to wave-induced flow for the four samples presented in the previous section but with water in the pores. Due to the fluid-pressure diffusion inside the sample, the response to the imposed stress is not instantaneous. Therefore, as stated earlier, the computed elastic moduli are complex and frequency dependent. For a given elastic modulus M(𝜔) measured on a sample, the attenuation can be obtained as a function of frequency 𝜔 by taking the ratio (a proof is given in the first appendix of Masson and Pride [2014]) Q−1 (𝜔) = M

Im{M(𝜔)} . Re{M(𝜔)}

(30)

In this manner, we obtain Q−1 for each of the five modes in 2-D corresponding to M11 , M22 , M12 , M21 , and G. M In Figure 11, we present the attenuation curves versus frequency computed using the numerical simulations of induced mesoscopic flow. Each panel in Figure 11 corresponds to one of the samples in Figures 7–10. The different lines correspond to the inverse quality factors associated with the moduli M11 , M12 , M22 , and G. For all samples, the attenuation levels are relatively low within the seismic frequency band (say Q−1 ≤ 2 × 10−3 ) which suggests that an elastic solid might be appropriate to represent these media. Apart from the Massillon sample, we observe no significant frequency dependence in the attenuation levels across the seismic band which suggests that constant Q−1 might be appropriate when modeling seismic waves through consolidated sandstones. The significant increase in attenuation with increasing frequency for the Massillon sandstone is due to the fact that most of the fluctuation in the elastic properties occurs at small scales. The decrease MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

16

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

in the attenuation toward low frequencies is due to the finite size of the samples (i.e., we have no heterogeneities with size greater than the sample’s size). For larger samples, we would expect the attenuation to stay more or less constant toward low frequencies. Only the attenuation curves associated with the Berea sample seems to suggest that higher attenuation levels might be expected at lower frequencies. When comparing the quality factors associated with the moduli M11 and M22 , we see that the Berea sandstone, the fractured sandstone, and the Rainbow sandstone are significantly anisotropic. We see more attenuation for waves traveling in the direction perpendicular to the rock’s bedding compared to waves traveling in the direction parallel to the bedding. For all samples, the attenuation associated with the shear modulus is almost nonexistent in the seismic frequency band.

5. Conclusion We introduce a microindentation scanner to measure the spatial distribution over mesoscopic scales (mm to 10 cm) of an indentation elastic modulus within natural rocks. We show that mapping the elastic properties of rocks with a fine resolution of millimeters or less is feasible using automated microindentation tests. The indentation maps obtained for clean, well-consolidated sandstones share some common features. For all samples scanned, we observe a spatial correlation between the fluctuation in the elastic properties and the rock’s bedding. At smaller scales within the sample, the fluctuations likely reflect spatial variations in the grain-packing arrangements and degree of cementation between grains. At larger scales, the fluctuations are associated with variations in the properties of the materials composing the rock. The indentation modulus probability distributions are well characterized using lognormal probability density functions. We use the indentation modulus maps measured over slabs of sandstone to numerically estimate the overall drained elastic moduli of each slab. Using the observed lognormal distribution measured on each slab, the multicomponent lower Hashin-Shtrikman bound is shown to do an excellent job of predicting the numerically computed moduli (however, to within the error of the numerical simulations, the geometric mean does an equally good job). We show that the indentation scans can be used along with the numerical experiments to compute estimates of seismic attenuation versus frequency. We observe that wave-induced flow in the well-consolidated sandstones we measured does not produce significant levels of attenuation. A limitation of 2-D scans of the elastic moduli of rock slabs is that fractures perpendicular to the slab will not have their large compliance normal to the fracture plane be represented properly in the scans. Acknowledgments This work was funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under award DE-ACO2-05CH11231. Y. Masson has recently been supported through the European Union’s Seventh Framework Program(FP-7-IDEAS-ERC), ERC Advanced Grant (WAVETOMO). People interested to obtain the data presented in this manuscript may contact Yder Masson via email at [email protected].

MASSON AND PRIDE

References Adam, L., M. Batzle, K. T. Lewallen, and K. van Wijk (2009), Seismic wave attenuation in carbonates, J. Geophys. Res., 114, B06208, doi:10.1029/2008JB005890. Anonymous (2004), Test methods for compressive strength and elastic moduli of intact rock core specimens under varying states of stress and temperatures, ASTM Int., doi:10.1520/d7012-14. Batzle, M. L., D.-H. Han, and R. Hofmann (2006), Fluid mobility and frequency-dependent seismic velocity—Direct measurements, Geophysics, 71(1), N1–N9. Biot, M. A., and D. G. Willis (1957), The elastic coefficients of the theory of consolidation, J. Appl. Mech., 24, 594–601. Birch, F., and D. Bancroft (1938), The effect of pressure on the rigidity of rocks. II, J. Geol., 46, 113–141. Cadoret, T. (1993), Effet de la saturation eau/gaz sur les propriétés acoustiques des roches. Étude aux fréquences sonores et ultrasonores, PhD thesis, Université de Paris VII, Paris, France. Castagna, J. P., M. L. Batzle, and R. L. Eastwood (1985), Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks, Geophysics, 50, 571–581. Doerner, M., and W. Nix (1986), A method for interpreting the data from depth-sensing indentation instruments, J. Mater. Res., 1, 601–609. Farmer, C. (2002), Upscaling: A review, Int. J. Numer. Methods Fluids, 40, 63–78. Gurevich, B., and S. L. Lopatnikov (1995), Velocity and attenuation of elastic waves in finely layered porous rocks, Geophys. J. Int., 121, 933–947. Ibrahim, N. (2008), Caractérisation des propriétés mécaniques des géomatériaux par technique de micro indentation, PhD thesis, Lille 1, Villeneuve-d’Ascq, France. Leite, M. H., and F. Ferland (2001), Determination of unconfined compressive strength and Young’s modulus of porous materials by indentation tests, Eng. Geol., 59, 267–280. Lucet, N., P. N. J. Rasolofosaon, and B. Zinszner (1991), Sonic properties of rocks under confining pressure using the resonant bar technique, J. Acoust. Soc. Am., 89, 980–990. Lumley, D. E., D. Bevc, J. Ji, and S. Talas (1992), Lab measurement of elastic velocities in dry and saturated Massillon sandstone, Stanford Exploration Project Report, 73, Stanford Univ., Stanford, Calif. Masson, Y., and S. R. Pride (2014), On the correlation between material structure and seismic attenuation anisotropy in porous media, J. Geophys. Res. Solid Earth, 119, 2848–2870, doi:10.1002/2013JB010798. Masson, Y. J., and S. R. Pride (2007), Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopic-scale heterogeneity, J. Geophys. Res., 112, B03204, doi:10.1029/2006JB004592. Masson, Y. J., S. R. Pride, and K. T. Nihei (2006), Finite-difference modeling of biot’s poroelastic equations at seismic frequencies, J. Geophys. Res., 111, B10305, doi:10.1029/2006JB004366. Mosegaard, K., and A. Tarantola (2002), Probabilistic approach to inverse problems, Int. Geophys. Ser., 81, 237–265.

THE MECHANICAL PROPERTIES OF ROCKS

17

Journal of Geophysical Research: Solid Earth

10.1002/2015JB012248

Murphy, W. F. (1982), Effects of partial water saturation on attenuation in Massillon sandstone and Vycor porous glass, J. Acoust. Soc. Am., 71(6), 1458–1468. Murphy, W. F., J. N. Roberts, D. Yale, and K. W. Winkler (1984), Centimeter scale heterogeneities and microstratifaction in sedimentary rocks, Geophys. Res. Lett., 11, 697–700. Norris, A. N. (1993), Low-frequency dispersion and attenuation in partially saturated rocks, J. Acoust. Soc. Am., 94, 359–370. Oliver, W. C., and G. M. Pharr (1992), An improved technique for determining hardness and elastic-modulus using load and displacement sensing indentation experiments, J. Mater. Res., 7, 1564–1583. Oliver, W. C., and G. M. Pharr (2004), Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology, J. Mater. Res., 19(01), 3–20. Pharr, G. M., W. C. Oliver, and F. R. Brotzen (1992), On the generality of the relationship among contact stiffness, contact area, and elastic-modulus during indentation, J. Mater. Res., 7(3), 613–617. Pride, S. R., and J. G. Berryman (2003a), Linear dynamics of double-porosity and dual-permeability materials. I. Governing equations and acoustic attenuation, Phys. Rev. E, 68, 036–603. Pride, S. R., and J. G. Berryman (2003b), Linear dynamics of double-porosity and dual-permeability materials. II. Fluid transport equations, Phys. Rev. E, 68, 036–604. Pride, S. R., J. G. Berryman, and J. M. Harris (2004), Seismic attenuation due to wave-induced flow, J. Geophys. Res., 109, B01201, doi:10.1029/2003JB002639. Quintal, B., H. Steebr, M. Frehner, and S. M. Schmalholz (2011), Quasi-static finite element modeling of seismic attenuation and dispersion due to wave-induced fluid flow in poroelastic media, J. Geophys. Res., 116, B01201, doi:10.1029/2010JB007475. Sneddon, I. (1965), The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile, Int. J. Eng. Sci., 3(1), 47–57. Sneddon, I. N. (1948), Boussinesq’s problem for a rigid cone, Math. Proc. Cambridge Philos. Soc., 44(4), 492–507. Spencer, J. W. (1981), Stress relaxations at low-frequencies in fluid-saturated rocks—Attenuation and modulus dispersion, J. Geophys. Res., 86, 1803–1812. Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Pa. Torquato, S. (2002), Random Heterogeneous Materials, 701 pp., Springer, New York. Wang, W., S. Corn, and P. Ienny (2009), Characterization of mechanical properties of geomaterials by microindentation test. Calibration and applications, Eur. J. Environ. Civil Eng., 13(3), 305–328. White, J. E. (1975), Computed seismic speeds and attenuation in rocks with partial gas saturation, Geophysics, 40, 224–232. White, J. E., N. G. Mikhaylova, and F. M. Lyakhovitsky (1975), Low-frequency seismic waves in fluid-saturated layered rocks, Izv. Acad. Sci. USSR, Phys. Solid Earth, 11, 654–659. Winkler, K. W., and R. D’Angelo (2006), Ultrasonic borehole velocity imaging, Geophysics, 71, F25–F30. Winkler, K. W., and A. Nur (1982), Seismic attenuation: Effects of pore fluids and frictional-sliding, Geophysics, 47(1), 1–15. Zhu, W., J. J. Hughes, N. Bicanic, and C. J. Pearce (2007), Nanoindentation mapping of mechanical properties of cement paste and natural rocks, Mater. Charact., 58, 1189–1198.

MASSON AND PRIDE

THE MECHANICAL PROPERTIES OF ROCKS

18