The concept of critical porosity was introduced to the rock physics community in the 1990s by Amos Nur and colleagues at Stanford University. Critical porosity has proved useful in reservoir characterization using seismic, borehole acoustic waves, and sonic logs, as well as in geomechanical applications such as prediction of rock strength and in-situ stress. This article reviews the concept and its relation to the phenomenon of mechanical percolation in statistical physics. Estimates of the critical porosity using laboratory measurements, computer simulations, and effective medium and effective field theories are discussed. Effective medium and effective field theories are not able to describe accurately the behavior of rocks in the vicinity of the critical porosity due to assumptions of idealized pore shapes and failure to explicitly account for pore connectivity. However, they do prove useful for understanding the effect of porosity, pore shape, and pore orientation distribution on elastic wave velocities, and suggest relations between porosity and elastic moduli useful for modeling of elastic waves in porous rocks.

Amos Nur (1938–2024) was a pioneering figure in the field of rock physics and geophysics. He founded the Stanford Rock Physics and Borehole Geophysics Project, significantly advancing the understanding of rock physics, seismic monitoring, and digital rock physics. His innovative work on ultrasonic velocity measurements in the laboratory and the interpretation in terms of rock properties and applied stress has been crucial in the characterization and monitoring of oil and gas reserves. Nur's interdisciplinary approach also bridged the gap between multiple disciplines. In addition to his well-known contributions to rock physics, he also united the diverse fields of geophysics and archaeology, exploring the impact of earthquakes on ancient civilizations, and his introduction of the concept of critical porosity established a bridge to percolation theory and statistical mechanics. He will be remembered for his creativity and his ability to inspire students and colleagues alike, leaving a long-lasting legacy in the geoscience community.

Crystalline rocks in the subsurface are of interest for geothermal energy extraction, nuclear waste storage, and, when weathered or fractured, as aquifers. Compliant discontinuities such as microcracks, cracks and fractures may nucleate and propagate due to changes in pore pressure, stress and temperature. These discontinuities may provide flow pathways for fluids and, if fracturing extends to surrounding rocks, may allow escape of fluids to neighbouring formations. Monitoring such rocks using sonic logs, passive seismic, borehole seismic and surface seismic requires understanding of the propagation of elastic waves in the presence of such discontinuities. These may have an anisotropic orientation distribution as in situ stress may be anisotropic. As crystalline rock may display intrinsic anisotropy due to foliation and the preferential orientation of anisotropic minerals, quantification of the relative importance of intrinsic and microcrack-induced anisotropy is important. This may be achieved based on the stress sensitivity of elastic wave velocities. A method that allows both the orientation distribution of microcracks and the stress dependence of their normal and shear compliance to be estimated independently of the elastic anisotropy of the background rock is presented. Results are given for anisotropic samples of gneiss from Bukov in the Czech Republic and granite from Grimsel in Switzerland based on the ultrasonic velocity measurements of Aminzadeh et al. The microcrack orientation distribution is approximately transversely isotropic for both samples with a preferred orientation of microcrack normals perpendicular to foliation. This preferred alignment is stronger in the sample of gneiss than in the granite sample, and the normal and shear compliance of the microcracks decreases with increasing compressive stress. This occurs because the contact between opposing faces of the discontinuities grows with increasing compressive stress, and this results in a decrease in elastic anisotropy with increasing compressive stress. At low stress, the ratio of microcrack normal compliance to shear compliance is approximately 0.25 for the granite sample and 0.7 for the sample of gneiss. The normal compliance Z(N) for both samples decreases faster with increasing compressive stress than the shear compliance Z(T), resulting in a decrease in Z(N)/Z(T) with increasing compressive stress.

The strong sensitivity of velocity to stress observed in many sandstones originates from the response of stress-sensitive discontinuities such as grain contacts and microcracks to a change in effective stress. If the change in stress is anisotropic, then the change in elastic wave velocities will also be anisotropic. Characterization of stress-induced elastic anisotropy in sandstones may enable estimation of the in situ three dimensional stress tensor with important application in solving problems occurring during drilling, such as borehole instability, and during production, such as sanding and reservoir compaction. Other applications include designing hydraulic fracture stimulations and quantifying production-induced stresses which may lead to rock failure. Current methods for estimating stress anisotropy from acoustic anisotropy rely on third-order elasticity, which ignores rock microstructure and gives elastic moduli that vary linearly with strain. Elastic stiffnesses in sandstones vary non-linearly with stress. Using P- and S-wave velocities measured on Gulf of Mexico sandstones, this non-linearity is found to be consistent with a micromechanical model in which the discontinuities are represented by stress-dependent normal and shear compliances. Stress-induced anisotropy increases with increasing stress anisotropy at small stress but then decreases at larger stresses as the discontinuities close and their compliance decreases. When the ratio of normal-to-shear compliance of the discontinuities is unity, the stress-induced anisotropy is elliptical, but for values different from unity, the stress-induced anisotropy becomes anelliptic. Although vertical stress can be obtained by integrating the formation's bulk density from the surface to the depth of interest, and minimum horizontal stress can be estimated using leak-off tests or hydraulic fracture data, maximum horizontal stress is more difficult to estimate. Maximum horizontal stress is overpredicted based on third-order elasticity using measured shear moduli, with estimates of pore pressure, vertical stress and minimum horizontal stress as input. The non-linear response of grain contacts and microcracks to stress must be considered to improve such estimates.

Unconsolidated sandstones are attractive targets for underground storage of carbon due to their high porosity and permeability. Monitoring of injection and movement of CO 2 in such formations using elastic waves requires an understanding of the acoustic properties of the sandstone. Current approaches often use the so‐called soft‐sand model in which a Hertz–Mindlin model of the acoustic properties at high porosity is mixed with the acoustic properties of the mineral phase to predict the acoustic properties over the entire porosity range. Using well‐log data from two unconsolidated sand formations of interest for CO 2 storage, we discuss the limitations of this model and provide an alternative approach in which the mechanical properties of grain contacts are obtained by inversion, and the properties of infill material lying within the pore space are estimated. The formations considered are the Paluxy Formation in Kemper County, Mississippi, and the Frio Formation near Houston, Texas. The ratio of the normal to shear compliance of the grain contacts is found to be significantly less than unity for both formations. This implies that the grain contacts are more compliant in shear than in compression. However, the grain contact compliance is higher and the ratio of the normal to shear compliance is lower for the Frio example than for the Paluxy case, and this may lead to sliding at grain contacts with low shear compliance and transport of grains during fluid flow, particularly if CO 2 acts to weaken any cement that may be present at the grain contacts. Such transport was suggested by Al Hosni et al. in explaining why the magnitude of the time‐lapse effect due to the injection of CO 2 at the Frio CO 2 injection site is greater than predicted using conventional rock physics models. A simple model of the mechanical properties of infill material lying within the pore space suggests that the bulk and shear moduli of infill material in the Paluxy case are significantly higher than the Frio case, consistent with the lower grain contact compliance in the Paluxy case.

Organic-rich shales contain large amounts of oil and gas and are anisotropic because of fine-scale layering and the partial alignment of organic matter and anisotropic clay minerals with the bedding. An example is the Wolfcamp Shale in the Permian Basin. Elastic anisotropy needs to be accounted for in the characterization of such formations using seismic data and plays a role in hydraulic fracturing and in the evaluation of stress changes and geomechanical effects resulting from production. Using extensive well log data acquired in the Midland Basin, the eastern sub-basin of the Permian Basin, we estimate and compare the elastic anisotropy in the Middle and Upper Wolfcamp Shale by combining data from a vertical pilot well with two lateral wells, one (6SM) drilled in the Middle Wolfcamp and one (6SU) drilled in the Upper Wolfcamp. The data used were acquired at the Hydraulic Fracture Test Site 1, located in the eastern part of the Midland Basin. Thomsen's anisotropy parameter gamma$\gamma $ calculated from the fast and slow shear sonic is higher on average for the 6SM lateral than for 6SU, consistent with there being less carbonate content in 6SM than in 6SU. However, the anisotropy parameter gamma$\gamma $ in some regions with higher carbonate content in well 6SU is higher than in well 6SM. This may indicate the influence of natural fractures. The primary set of steeply dipping fractures observed in the lateral wells at Hydraulic Fracture Test Site 1 acts to increase gamma$\gamma $ if the ratio of the normal-to-shear fracture compliance is less than about 0.5. Sub-horizontal fractures may also increase gamma$\gamma $ and could affect the vertical extent of hydraulic fractures. Relations between elastic moduli C33 and C55 in the Upper and Lower Wolfcamp in a vertical pilot well allow C33 to be predicted in a lateral well using measurements of C55 in that well. Comparison of Thomsen's anisotropy parameters gamma$\gamma $ and epsilon$\varepsilon $, with gamma$\gamma $ calculated from the measured values of C55 and C66 and epsilon$\varepsilon $ calculated from the measured values of C11 and predicted values of C33, show that epsilon$\varepsilon $ is mostly greater than gamma$\gamma $.

High‐porosity sandstones are important for hydrocarbon production, underground CO 2 storage, extraction of geothermal energy and freshwater aquifers. Porosity of sandstones may be estimated using elastic wave velocities, but these depend also on fluid saturation, clay content, pore shape and contacts between sand grains. An understanding of how elastic properties of sandstones depend on these factors is important for characterizing their storage potential and for geomechanical issues, such as sanding, borehole stability, reservoir compaction and fracturing. Ultrasonic velocity measurements in clay‐bearing sandstones indicate that much of the clay in shaly sandstones is non‐load‐bearing. This enables a simple approach for modelling the elastic properties of shaly sandstones that includes the effect of pore concavity and agrees with ultrasonic P‐ and S‐velocities measured in the laboratory. Despite this agreement, some clay may reside within the contacts and may act to inhibit the development of quartz cement, thus reducing porosity loss and helping to preserve storage volume. This appears to be the case for the Lower Mt. Simon Sandstone, a target formation for underground storage of CO 2 in the Illinois Basin, for which the bulk moduli agree with the predicted bulk moduli, but the shear moduli are lower than predicted. This appears to result from an increase in a shear compliance of the grain contacts that may enable sliding along the grain contacts and increase the tendency to shear failure.

Geothermal resources have potential for providing cost-effective and sustainable energy. Monitoring of production-induced changes in geothermal reservoirs using seismic waves requires understanding of the elastic properties of the rock and how they change due to injection of fluids and opening and closing of natural and hydraulic fractures. P- and S-wave velocities measured in a granitic geothermal reservoir using sonic logging are systematically lower than those predicted using the composition of the rock. Cracks may occur in granitic rocks from tectonic stresses and from the thermal expansion mismatch between differently oriented anisotropic crystals. An isotropic orientation distribution of microcracks causes a significant reduction in both the P- and S-velocities, consistent with the observed sonic P- and S-velocities. Vertical fractures cause a difference in the velocity of vertically propagating shear waves polarized parallel and perpendicular to the fractures. An assumption that the lower measured velocities are caused by the presence of vertical fractures is inconsistent with the sonic data. This is because vertical fractures cause a decrease in slow S-wave velocity that greatly exceeds the decrease in P-wave velocity, in contrast to the observed data. The growth of vertical fractures in the geothermal reservoir may be monitored using the difference in velocity of the fast and slow shear waves, while the change in P-velocity in a crossplot of measured P- and slow S-velocities is useful for estimating the ratio of the normal-to-shear compliance of the fractures.

Explicit formulae are presented for anisotropic parameters for the orthorhombic case of a polar anisotropic body (e.g. a shale) penetrated by a single set of vertical, rotationally variant fractures (e.g. joints).

The elastic frame moduli for the Lower Mt. Simon Sandstone are estimated by accounting for grain boundary compliance and concavity of pores in the vicinity of grain contacts. This sandstone is the target formation for underground storage of CO2 in the Illinois Basin – Decatur Project, a 1-million-ton carbon capture and storage project in Decatur, Illinois. Knowledge of the elastic frame moduli allows the feasibility of time-lapse seismic for monitoring the injection and storage of CO2 injection to be assessed.

ABSTRACTThe Marcellus Shale is one of the largest shale gas resources in the world and is anisotropic due to fine layering and the partial alignment of anisotropic clay minerals and organic matter with the bedding. This anisotropy can be approximated as Vertically Transversely Isotropic (transversely isotropic with a vertical axis of rotational symmetry) with five independent density‐normalized elastic stiffnesses A11, A13, A33, A55 and A66 in the two‐index notation and with axis of rotational symmetry along x3. Compressional and dipole shear‐wave data acquired in a horizontal well allows estimation of A11, A55 and A66, while the same data in a vertical pilot well allows estimation of A33 and A55. The ratio of vertically propagating P‐ to S‐velocity is and has a quadratic dependence on in the Upper Marcellus with minimum occurring for the largest volume fraction of kerogen. This relation allows estimation of A33 along a lateral well using measured values of A55 in the well. Comparison of Thomsen's anisotropy parameter ε is found to be mostly greater than anisotropy parameter γ. Estimating A13 as the average of A33 – 2A55 and A11 – 2A66 as proposed recently by Yan and Vernik allows Thomsen's anisotropy parameter δ and a parameter K0 that relates the horizontal and vertical effective stresses to be estimated. The results are expected to help in estimating horizontal stress needed for design of hydraulic fractures, and in interpreting seismic amplitude variation with offset required for reservoir characterization.

The mechanical properties of sandstones, including porosity, density and elastic moduli, can be estimated non-destructively through elastic wave-velocity measurements. Here, the variation of elastic wave velocity with porosity in sandstones is modelled using Maxwell's effective field theory, extended to the elasticity of heterogeneous media by Sevostianov and coworkers. Comparing measured and predicted elastic wave velocities shows that on deposition, pores in sandstones are less stiff than spherical pores, but that their stiffness increases as porosity decreases. This suggests that concavity of pores in sandstone decreases with decreasing porosity. This interpretation is confirmed by the simple model of Sevostianov and Giraud in which concave pores are represented as superspherical pores, defined by a shape parameter that allows the effect of pore concavity on elastic wave velocities to be investigated. Inversion of measured velocities for this parameter indicates that pore concavity decreases with decreasing porosity. Moreover, values of the shape parameter obtained by inverting measured P-velocities alone are found to give a good prediction of both P- and S-wave velocities, confirming the applicability of the model.

The SEG Advanced Modeling (SEAM) Barrett model was designed to model typical land basins found in the North American midcontinent that host unconventional reservoirs, such as fractured shale reservoirs. This model was used recently in several studies to assess whether shale bodies could be resolved using azimuthal 3D P-P reflection seismic data. In one study, it was claimed that near-surface complexity prevents the identification of the shale bodies using azimuthal analysis and concluded that velocity variation with azimuth (VVAz) and amplitude variation with azimuth (AVAz) are not worth running in the Permian Basin. However. another study by different authors applied a different seismic processing sequence to successfully resolve the reservoir geobodies and showed promising AVAz and VVAz results. We have focused on the SEAM Barrett model itself. Despite some advantages. the limitations of the Barrett model prevent general conclusions to be drawn about the usefulness of VVAz and AVAz to characterize fractured unconventional reservoirs.

Measurements of elastic wave velocities enable non-destructive estimation of the mechanical properties, elastic moduli and density of snow and firn. The variation of elastic moduli with porosity in dry snow and firn is modeled using a differential effective medium scheme modified to account for the critical porosity above which the bulk and shear moduli of the ice frame vanish. A comparison of predicted and measured elastic moduli indicates that the shear modulus of ice in snow is lower than that computed from single crystal elastic stiffnesses of ice. This may indicate that the bonds between snow particles are more deformable under shear than under compression. A partial alignment of ice crystals also may contribute. Good agreement between elastic stiffnesses of the ice frame obtained from elastic wave velocity measurements and the predictions of the theory is observed. The approach is simple and compact, and does not require the use of empirical fits to the data. Owing to its simplicity, this model may prove useful in a variety of potential applications such as construction on snow, interpretation of seismic measurements to monitor and locate avalanches and estimation of density within compacting snow deposited on glaciers and ice sheets.

The elastic properties of clay in shales can be derived from a simple model consisting of domains of partially aligned anisotropic clay platelets embedded in a softer interplatelet medium representing clay-bound water and interparticle contacts. Approximating clay platelets as ellipsoidal inclusions and the interplatelet medium as isotropic, with effective bulk and shear moduli, allows the effect of the spatial distribution of inclusions upon the elastic properties of the clay matrix to be investigated. The pair distribution function governs the respective locations of clay platelets and is similarly assumed to be ellipsoidal. Due to compaction, the spacing of clay platelets perpendicular to their long axis is typically less than their spacing parallel to it. Hence, the aspect ratio of the pair distribution function is less than that of the clay platelets. This paper examines the impact of the pair distribution function upon the elastic stiffnesses and anisotropy of clay domains. Assuming a lower aspect ratio for the pair distribution function than for the clay platelets is found to reduce in-plane elastic stiffnesses C11, C12 and C66, but to increase out-of-plane stiffnesses C33, C13 and C55. Since in-plane elastic stiffnesses decrease whereas out-of-plane elastic moduli increase as the aspect ratio of the pair distribution function decreases below that of the clay platelets, the anisotropy parameters ϵ, γ and η decrease. However, the change in Thomsen’s δ parameter is positive.

This study uses well log analysis along with empirical and theoretical equations to characterize properties of the unconventional resource in the Bakken Formation of North Dakota. Cross-plot analysis of petrophysical and elastic properties indicates that density (ρ), P- and S-impedance, Young’s modulus (E), shear modulus (μ), and μρ are good total organic carbon (TOC) content indicators. Shales with high TOC content (greater than 10 wt%) have lower velocities and densities (<2.31 g/cm3) than lower TOC cases. High and low TOC rocks exhibit Vp/Vs values between 1.65 and 1.75. Both the calibration of Greenberg and Castagna’s trends and Vernik’s predictive methods can be helpful for shear-wave velocity and impedance prediction in the Bakken Formation. A total porosity range between 6% and 10% corresponds to Bakken shales with TOC content ≥ 10 wt% with P-impedance (AI) and S-impedance (SI) cutoff values of 25,000 g/cm3·ft/s and 15,000 g/cm3·ft/s, respectively. Elastic moduli decrease as TOC content and hydrocarbon saturation increase. The high TOC shales correlate with E, μ, and μρ values between 11–24 GPa, 2.5–8.75 GPa, and 10–18 GPa·g/cm3, respectively. We present empirical equations, calibrated for the Bakken Formation members, to be used for quick-look calculations and comparisons. The methodology presented here can help inform analysis for selecting promising target areas in unconventional resource plays.

A fracture may be characterized by its normal and shear compliances that relate the discontinuity in displacement of the fracture faces to an applied traction. To assess the feasibility of using elastic waves for monitoring natural and hydraulic fractures in anisotropic shale reservoirs, an understanding of how these compliances relate to the properties of the background rock and any fracture infill material is required. Three models of vertical fractures in a VTI (transversely isotropic medium with vertical axis of rotational symmetry) background suitable for feasibility studies are discussed. The first two models are suitable for open vertical fractures and represent the fractures as aligned low aspect ratio elliptic cylindrical inclusions with long axis of the ellipse either horizontal or vertical. In these models the horizontal and vertical shear compliance differ in a way that depends on the anisotropy of the background medium and the orientation of the fractures. The third model is suitable for estimating the properties of hydraulic fractures containing proppant. In this model, proppant is represented as a random packing of identical elastic spheres subject to a uniaxial compression due to the normal stress acting on the fracture. The results are expected to help improve the reliability of fracture models used for assessing the feasibility of elastic waves for characterizing natural and hydraulic fractures.

Seismic shear wave velocity (S-velocity) shows a decrease towards the base of ice sheets in Antarctica and Greenland that is not accompanied by a corresponding decrease in compressional velocity (P-velocity). This decrease has been interpreted as arising from liquid water below the melting point (pre-melt water) at grain boundaries, but the lack of a corresponding decrease in P-velocity has not been explained. Representing grain boundaries as displacement discontinuities allows the change in P- and S-velocities to be written as functions of the normal and shear compliance of the grain boundaries. This allows the normal-to-shear compliance ratio of the grain boundaries to be constrained, and seismic anisotropy resulting from a partial orientation of grain boundaries to be estimated. This approach demonstrates that the observed reduction in S-velocity with no significant decrease in P-velocity near the base of ice sheets in Antarctica and Greenland can be explained by pre-melt water at small aperture grain boundaries. Such water may enable sliding along the grain boundaries and so may enhance creep of ice near the base of ice sheets. If stress state is anisotropic the aperture of water-containing grain boundaries may vary with azimuth, with the most open grain boundaries oriented with strikes perpendicular to least compressive stress. Microcracks and fractures may be treated also as displacement discontinuities and, together with oriented grain boundaries, may contribute to shear wave splitting as observed in West Antarctica in a fast-moving ice stream.

Unconsolidated sands provide zones of high porosity and permeability important for freshwater aquifers, hydrocarbon production, and CO2 sequestration. An understanding of the acoustics of unconsolidated sands enables the characterization of such formations using ultrasonics, borehole acoustics, and seismic methods. Inversion of ultrasonic compressional and shear velocities measured for unloading as a function of confining pressure for room-dry unconsolidated sands allows information on the mechanical properties of the grain contacts to be obtained using an approach based on the divergence theorem. This allows the effective compliance of sand to be written as the sum of the compliance of the pores and of the grain contacts, and it does not assume that the grains are identical spheres, in contrast to previous approaches. Grain contacts are found to be more compliant under shear than under normal compression, and the ratio of the normal-to-shear compliance decreases with decreasing confining pressure, implying that the shear compliance increases faster with decreasing confining pressure than the normal compliance. This is of importance in understanding the role of shear in the failure of unconsolidated sands, such as occurs in shallow water flow, sanding, and failure around injectors, where the change in stress is a function of the normal-to-shear compliance ratio of the grain contacts.