Characterizing the elastic signatures of overpressure of shale caused by the smectite-to-illite transition relies on a good understanding of this mechanism and is also necessary for pore-pressure prediction. Methods of pore-pressure prediction in shales that have undergone smectite-to-illite transition are mostly based on empirical fitting without a quantitative interpretation based on a micromechanism analysis. With upscaled wireline-logging data, two trends of smectite-to-illite transition are categorized by using the crossplot of sonic traveltime and density. Trend I associated with a fluid-expansion scenario exhibits a decrease of sonic velocity with little change in the bulk density, whereas trend II induced by a fluid-loss scenario contains an increase of density with little change in the sonic velocity. The fluid expansion typically gives rise to high-magnitude overpressure and tends to happen when the overlying formations have more shaly contents and low permeability. The fluid loss case tends to have relatively deeper overpressure onsets, and its overlying formations tend to have more sandy contents with relatively high permeability. We develop a modeling framework to capture the elastic and pore-pressure evolution characteristics in shale during the smectite-to-illite transition. With proper bulk volume models, the velocity, density, and pore pressure increase of shale can be computed in the fluid expansion, fluid loss, and a mixture of these two scenarios. After calibration with logging data, rock-physics modeling can quantitatively interpret the rock-property evolution characteristics within the smectite-to-illite transition zone.
Overpressure exists widely from active plate boundaries to passive basins and has been extensively studied in academia and the petroleum industry. Overpressure tends to happen once the pore-fluid pathways to the surface (seabottom for offshore and water table for onshore) are cut off so that the interstitial fluid cannot escape to restore the hydrostatic pressure. It is necessary to investigate pore pressure because excess pore pressure may cause geohazards (Swarbrick and Osborne, 1998; Chatterjee et al., 2011) and can reflect the status of pore-fluid existence, accumulation, and migration (Bruce, 1984; Hao et al., 2007).
Low permeability is the main reason for the retardation of vertical or lateral fluid movement in sedimentary rocks. Mudrocks including shales constitute approximately 70% of sedimentary rocks, and they can have an extremely low permeability (from micro- to nanodarcy) due to their fine-grained () particles. Overpressure can occur in shales in the process of mechanical compaction, volume expansion of pore fluids from physical or chemical changes, or both (Swarbrick and Osborne, 1998; Dutta, 2002; Mukerji et al., 2002; Singha and Chatterjee, 2014; Dasgupta et al., 2016). Common mechanisms include disequilibrium compaction or undercompaction, kerogen-to-gas transition, and smectite-to-illite transition or clay diagenesis. Therefore, to perform reliable pore-pressure prediction, we need to understand the mechanisms of overpressure in a region and their effects on the physical properties of sedimentary rocks and thus geophysical responses.
The signatures of overpressure with various mechanisms on velocity and density profiles are sketched in Figure 1. The dashed lines describe the loading curves of velocity and density increase similarly with depth at the normal compaction condition when pore pressure is identical to the hydrostatic pressure; the solid curves represent the behaviors of velocity and density varying with depth when disequilibrium compaction (Figure 1a) or unloading (Figure 1b) happen. Velocity and density vary with depth similarly for overpressure caused by disequilibrium compaction (Bowers, 2001; Lahann et al., 2001; Katahara, 2003; Hoesni, 2004). This is because disequilibrium compaction develops during the compaction loading process, which is controlled by effective stress (the difference between overburden stress or confining pressure and pore pressure). If fine-grained deposits are buried at a rapid sedimentation rate, pore fluid does not have adequate time to be expelled from the pore space due to the poor permeability. When fluid retention helps to inhibit compaction, the porosity also stops reducing, and the bulk density and velocity increase is retarded.
The velocity trend deviates from the density trend in Figure 1b, and the onset depth of overpressure is caused by unloading (Bowers, 2001; Ramdhan and Goulty, 2011). The unloading mechanism refers to geologic reasons for a decrease of effective stress, such as the smectite-to-illite transition, kerogen-to-gas transition, erosion, uplift, and aquathermal expansion (Bowers, 1995; Carcione and Gangi, 2000; Lahann et al., 2001; Katahara, 2006; Hao et al., 2007; Lahann and Swarbrick, 2011; Ramdhan and Goulty, 2011; Couzens-Schultz et al., 2013; Tingay et al., 2013; Yu and Hilterman, 2014; Dutta, 2016). Overpressure caused by smectite-to-illite or kerogen-to-gas transitions without tectonic activity is attributed to a secondary mechanism. When the volume of pore fluids expands, an internal increase of pore pressure can be quickly built up within a narrow range of depths. Because compaction is irreversible (decreasing the effective stress will not increase the porosity), bulk properties such as porosity and density change little when the effective stress is reduced. By contrast, seismic velocity and resistivity are more sensitive to a decrease of effective stress than bulk properties (Bowers, 2001), in that the effective stress affects the degree of grain contact, opening or closure of cracks/fractures, and the frame modulus. The main objectives of the study are (1) to relate the evolution characteristics of shale elastic properties to the concomitant excess pore pressure during the smectite-to-illite transition and (2) to quantitatively interpret or predict the rock-physics properties and pore-pressure evolution of shales during the smectite-to-illite transition.
To interpret the geophysical properties variation and predict pore pressure in shales containing clay diagenesis, a few methods are used. These methods include increasing the coefficient of Eaton (1975), shale compaction model by Dutta (1986), elastic unloading method by Bowers (1995), and load-transfer method by Lahann et al. (2001). Dutta’s (1986) approach is theoretical, whereas methods of Eaton (1975), Bowers (1995), and Lahann et al. (2001) are empirical. Empirical methods can work with local calibration (Gutierrez et al., 2006; Lahann and Swarbrick, 2011; Dutta et al., 2015). The elastic unloading model assumes a purely elastic process during the smectite-to-illite transition, which means that the shale’s velocity and density will recover along the unloading path when the effective stress increases. However, this method only caters to a single compaction trend and requires the derivation of the historical maximum effective stress. By contrast, the load-transfer model assumes that the smectite-to-illite transition shifts the shale’s properties in a plastic behavior, so that the compaction relation between porosity and effective stress has changed, i.e., a reloading will not track back the unloading curve but along a new compaction trend. Even though the effective stress decreases, the density can increase as a result of chemical compaction or clay diagenesis. Nevertheless, the load-transfer model does not explain the gap between different compaction trends and the relationship established between porosity and velocity. Little work has focused on the fundamental micromechanism of the smectite-to-illite transition and the quantitative description of how shale’s elastic stiffness and pore pressure change during the transition.
To understand the shale-property evolution within the smectite-to-illite transition zone, we select 37 wells from hundreds of wells from offshore Louisiana that have relatively good data quality, where the smectite-to-illite transition is reported as a primary mechanism of overpressure (Verm et al., 1998; Yu and Hilterman, 2014). This is followed by an identification of two trends of smectite-to-illite transition (fluid expansion and fluid loss/escape) based on the wireline-logging measurements and a detailed discussion of how geologic settings control the occurrence of both. Their signatures on the crossplots of the P-wave velocity and density can be explained by the hydration status of the interlayer space and distribution of transferred water from the interlayer bound water. During the transition, if the interlayer water in smectite transfers to free water and the free bulk water is preserved in the pore system, fluid expansion causes a velocity decrease, little density change, and a greater degree of overpressure. If the released water escapes, the shale has a density increase with little velocity change and with a relatively small pore pressure increase.
To quantitatively assess the shale elastic properties variation, we propose a framework to model the elasticity of the illite-smectite mineral with various amounts of interlayer water and elasticity of shale within the smectite-to-illite transition zone. We present three bulk volume models corresponding to fluid expansion, fluid loss, and the hybrid of these two. The modeled velocities and densities accord with the trends of the well-logging data. The issue of pore pressure increase within the transition zone where smectite alters to illite is also resolved through the volumes and compressibility of the pore fluid and pore space. After calibrating with two wells containing fluid expansion and fluid loss, respectively, the modeling result can be used to interpret how shale properties change within the transition zone based on the corresponding micromechanisms.
The structure of this paper is organized as follows: First, we review the structure and elasticity of clay minerals. Then, we diagnose the elastic characteristics of the smectite-to-illite transition with logging data and categorize two end-member trends for this transition. Next, we elaborate our framework of rock-physics modeling and apply the established rock-physics modeling to interpret the variations of shale P-wave velocity, density, and pore pressure during this transition. Finally, further discussion and conclusions are presented.
CLAY-MINERAL STRUCTURE AND ELASTICITY
The structure of clay mineral
Clay minerals have a sheet-like or layered structure, consisting of two basic units: a silicon tetrahedral sheet (T) and an alumina octahedral sheet (O). Clay minerals composed of a three-layered structure (T-O-T) belong to the smectite group (Colten-Bradley, 1987; Moyano et al., 2012), which includes smectite and illite. In illite, the distance between the first silica layer of one T-O-T sheet and the next silica layer of a neighboring T-O-T sheet (basal spacing) is approximately 10 Å (Figure 2). The basal spacing of smectite depends on hydration and can vary between 9.3 and 20 Å (Ougier-Simonin et al., 2016), resulting from smectite adsorbing various amounts of polar water molecules present in between the clay platelets. This type of water is called an interlayer or bound water, which can be eliminated between 110°C and 200°C in a laboratory (Sarout and Guéguen, 2008). The diameter of a water molecule is close to 3 Å, and basal spacing with one, two, or three layers of interlayer water is approximately 12, 15, or 17.8 Å, correspondingly (Ebrahimi et al., 2012; Sayers and den Boer, 2016).
The elastic stiffness of clay mineral
The elastic stiffness of clay minerals is difficult to measure directly because of the fine-grained nature of clay minerals and the uncertainty of the clay hydration status. Additionally, the components of clay, the dry clay sheet (T-O-T), and the interlayer water, are also difficult to determine. Although we almost never observe “dry clay” underground and clay always contains bound water and/or interlayer water, some authors estimate the properties of the dry clay mineral (e.g., density, bulk, and shear moduli) through laboratory measurement (Han et al., 1986; Wang et al., 2001) or calibration with rock-physics modeling (Holt and Fjær, 2003). For the interlayer water, its largest difference with free water (bulk water) is that the motion of interlayer water is highly restricted because of the electrostatic force between the hydrated ions and the negatively charged clay surface. To characterize this, Holt and Fjær (2003) and Kolstø and Holt (2013) set a nonzero shear modulus for interlayer water when they modeled the elastic properties of shale. Holt and Kolstø (2017) also confirm that by using a discrete particle flow model, the shear stiffness of the clay interlayer region increases from zero to a few gigapascal when clay surfaces are charged.
Another solution to obtain the elastic stiffness of clay is by using a theoretical approach. The elastic stiffness of Na-montmorillonite as a function of hydrated status has been calculated by Ebrahimi et al. (2012). They obtain the full elastic tensor of Na-montmorillonite with zero to 258 water molecules added into the interlayer space by using molecular dynamics with a force field for structural simulation of clay minerals. In this quasistatic method, the stiffness exhibits linear elasticity at strains up to 0.2%. The basal spacing ranges from 9.3 to 19.6 Å with the corresponding numbers of water molecules added. The rotational symmetry axis of the clay mineral is assumed to be perpendicular to the dry clay platelet to make the clay elastic properties easy to use. We approximate the general elastic tensor of Na-montmorillonite with 21 elastic constants to a transverse isotropic (TI) medium with five elastic constants through an algorithm by Dellinger (2005). Figure 3 plots the best-fitting elastic constants of TI symmetry: , , , , and , which are similar to the results of Sayers and den Boer (2016).
The values of , which range from 150 to 300 GPa, appear to be higher than some other authors’ work on muscovite’s property (Alexandrov and Ryzhova, 1961; Vaughan and Guggenheim, 1986; McNeil and Grimsditch, 1993), which is believed to be close to the elasticity of illite. This might be due to the effects of boundary conditions on the horizontal stiffness in the molecular simulation. However, the trend is insightful, indicating that during clay dehydration (the number of water molecules in the interlayer space decreases from 258 to zero), clay becomes stiffer, which is in agreement with the experimental work by Ghorbani et al. (2009). If we treat interlayer bound water as a component of Na-montmorillonite, it has shear rigidity (a few gigapascal) because and increase during dehydration. The interlayer shear rigidity value indicated by the dehydration/hydration process is close to the values when clay surfaces are charged (Kolstø and Holt, 2013), which resembles the condition of the smectite-to-illite reaction.
ELASTIC CHARACTERISTICS OF THE SMECTITE-TO-ILLITE TRANSITION
In the northern Gulf of Mexico (GOM), overpressure exists extensively and can be caused by different mechanisms, such as disequilibrium compaction and the smectite-to-illite transition. Disequilibrium compaction takes place in young sediments with a fast sedimentation rate, and the smectite-to-illite transition (Lahann et al., 2001; Katahara, 2003; Lahann and Swarbrick, 2011; Yu and Hilterman, 2014) happens in smectite-rich shale at a temperature greater than 70°C.
We used wireline-logging data from 320 wells located in the shelf of the GOM, offshore Louisiana. The locations of the wells are mapped in Figure 4. The characteristics of the data are similar to the data in Yu and Hilterman (2014). The data include P-wave velocity, resistivity, and density averaged on each 60 m (200 ft) depth intervals. The 60 m interval data filter out high-frequency fluctuations in rock properties, which may obscure the low-frequency pattern. An interval is defined as sandstone with a shale volume less than 50% or as shale with a volume larger than 50%. The onset depth of overpressure in each well is recorded based on the logging run. Mud weights and temperature are read directly from the logging run and interpolated at each 60 m interval. No high-quality pressure measurement is available, so the pore pressure is approximated from mud weights. Because mud weights are typically (0.5 lb/gal) higher than the pore-pressure gradients (Dutta and Khazanehdari, 2006), we focus on the trend of pore-pressure development and subtract this value from mud weights to obtain the pore pressure. The hydrocarbon intervals and intervals with questionable logging quality were omitted. With the upscaled data of water-saturated shale formations, it is convenient for us to study the compaction trend. Based on the work of Verm et al. (1998), in Figure 4b, the smectite-to-illite transition is responsible for the overpressure in most wells in the north of the study area, whereas the disequilibrium compaction dominates the overpressure in the other wells (Yu and Hilterman, 2014).
Based on Figure 1b, the depth trends differ between velocity and density for unloading, which motivates us to select 37 wells with relatively good-quality data containing unloading overpressure. Specifically, this unloading is caused by the smectite-to-illite transition rather than by hydrocarbon generation. First, it is because there is little source rock reported in the Miocene, Oligocene, and Pliocene formations, which are generally shallower than 5000 m and thermally immature within this area (Galloway, 2008). Second, lacking an increase in wireline-log resistivity, together with the small amount of total organic carbon (generally less than 1%) estimated from DT-log resistivity (Passey et al., 1990), indicates that the pore space charging with hydrocarbon may not have occurred (Lahann et al., 2001). In the unloading zone, the low formation temperature (generally less than 120°C) with relatively high porosity (larger than 15% in average) does not indicate that aquathermal expansion would make a tremendous contribution to overpressure.
Figure 5 plots density in grams per cubic centimeter against the sonic traveltime in microseconds per foot of shale from the 37 wells with color-coded depth. Two empirical compaction trends of smectite and illite are added in the dark and light blue as references (Lahann et al., 2001; Dutta, 2002; Gutierrez et al., 2006). At depths greater than 3000 m, most of the shale data fall near the smectite line, whereas at deep depths less than 3000 m, most of the data follow the illite line. So, smectite-rich shale converts to illite-rich shale with increasing depth, and most of the transitions take place between 2500 and 3500 m. However, it is challenging to distinguish how the transition shifts the rock properties from the smectite line to the illite line.
Two trends of smectite-to-illite transition
Two end-member trends can be deduced for the smectite-to-illite transition zone through a mind experiment. We can simplify shale and focus on the clay-mineral transformation. Before the transition, due to the fine grain and large surface area of clay minerals, other minerals can be suspended within clay aggregates, such as quartz and feldspar. When the environment is conducive, smectite-rich shale transforms to illite-rich shale, which involves releasing the interlayer bound water as free water (Powers, 1967; Swarbrick and Osborne, 1998; Lahann et al., 2001). Based on the possible distribution of the released water, there are two end-member cases: fluid expansion and fluid loss.
Because the immovable interlayer water transfers to movable water, the volume of the pore liquid increases and the volume of clay decreases. The neighboring clay sheets tend to squeeze the interlayer space, thereby compressing the transferred water and increasing the pore pressure. If the smectite-to-illite transition proceeds in a closed system, the released water is preserved in the pore space to keep the system approximately at a constant bulk density, which is termed fluid expansion (Lahann, 2002; Mukerji et al., 2002). So, the effective stress decreases dramatically as the pore pressure increases internally, resulting in unloading overpressure. If the mineral composition of the shale does not change much with depth (Hower et al., 1976), there will be a rightward trend on the crossplot of sonic slowness (DT) and density as shown in Figure 6a because as the soft liquid component increases, the modulus decreases, density changes little, and sonic velocity decreases. This characteristic is shown in trend I.
An alternative to the closed-system reaction is the fluid loss. The transition takes place in a system with its boundary allowing the released fluid to escape entirely and thus generating a smaller pore pressure increase. The rotation and sliding of the surrounding grains are quite likely to occur so as to achieve a new balance between the porosity and effective stress; thus, the rock’s density would increase under the effect of continued compaction. As shown in Figure 6b, data in the transition zone trend are upward because the light and soft components’ volume fractions decrease and the modulus and density increase and sonic velocity changes little (the modulus and density increase and cancel each other). This characteristic is shown in trend II.
Sonic slowness-density crossplot analysis for single well
Sonic slowness (DT)-density crossplot analysis has been performed on the 37 wells selected from the offshore Louisiana, and some of them have the characteristics described in the mind experiment. Figure 7 displays three wells with trend I characteristic. Within their transition zones, the primary characteristics are that acoustic slowness increases and density changes little with increasing depth, similar to the increase in pore pressure caused by fluid expansion. Their onset depths of overpressure are, respectively, at 2600, 2524, and 3287 m based on the logging run. The data are averaged for each 60 m depth section (color-coded) along the well. Generally, at depths shallower than 3000 m, data of each well define the compaction trend, which is consistent with the smectite loading curves (the dark-blue line). Solid black arrows indicate the approximate paths of DT and density varying with increasing depth within the smectite-to-illite transition zones, spanning from 2800 to 3500 m approximately with formation temperatures between 80°C and 110°C falling within the clay diagenesis window. After completion of the smectite-to-illite transition, a reloading will occur if no extra water is released and pore fluid can discharge. For depths greater than the transition zones, the data trend for each well is consistent with the illite loading curves (the light-blue line), i.e., acoustic velocity and density increase. Therefore, the decrease in velocity and little change in density within the transition zone and the new loading curve support that the trend I smectite-to-illite transition is close to a process of inelastic unloading. It is different from the hydrocarbon gas generation case in Tingay et al. (2013) because the gas-generation case can be explained by an elastic unloading process that reloading will track back the unloading path and continue to move along the original loading curve. Moreover, a small decrease in density values in Figure 7c might suggest an increase in porosity by high-magnitude overpressure if no fluid is expelled and the shale composition does not change with the depth. Dense data near the smectite loading curves above the transition zone might result from the overpressure induced by the disequilibrium compaction, which has the relatively high values of mud density (up to or 15 pounds per gallon of mud weight).
Figure 8 shows the crossplots of shale density and sonic slowness (DT) for another three wells with a trend II pattern. Within their transition zones, it is the primary characteristic that density increases along with either little change (Figure 8a and 8c) or a slight reduction (Figure 8b) in velocity, indicating that transferred interlayer water cannot be preserved and partially or entirely escapes out of the boundary. The onset depths of overpressure are, respectively, 2546, 2908, and 3104 m based on the logging run. Depths are color coded in Figure 8a and 8c, and mud density in gram per cubic centimeter is color coded in Figure 8b. Likewise, data at depths shallower than 3000 m follow the smectite loading curves (the dark-blue line), and data at depths greater than 3500 m follow the illite loading curves (the light-blue line). The solid black arrows indicate how DT and density vary with depth when the smectite-to-illite transition occurs, spanning from 2500 to 3300 m approximately. The formation temperatures of the transition zones are also between 75°C and 110°C, which are within the clay diagenesis window. Also, the slight reduction in velocity (Figure 8b) might be caused by the fluid expansion to a mild degree, which also reflects that P-wave velocity is sensitive to the effective stress change due to excess pore pressure. Data near the smectite loading curves are not as dense as in the trend I wells, which have the relatively small values of mud density (up to ) and might be less affected by the disequilibrium compaction, i.e., vertical fluid transfers are easier in trend II than in trend I.
Comparing the six wells displayed with other wells bearing the trend I or trend II smectite-to-illite transition, the trend I transition zones start with a density near in average, which is larger than the density of the onset of the trend II transition zone near . We speculate that data near smectite-rich shale loading curve (formations overlying the transition zone) of trend I have smaller total porosity and lower permeability than trend II. So, disequilibrium compaction tends to accompany the trend I transition because of the relatively low permeability of the overlying formations. It also suggests that the smectite-to-illite transition can be inhibited by the disequilibrium compaction (Colten-Bradley, 1987) due to the relatively small porosity of the onset for the trend I transition as a result of a more thorough degree of compaction. Some features of the two trends are shown in Table 1. Another result of the smaller porosity for trend I is that a small increase of the pore-fluid volume can lead to a larger increase in the pore pressure.
We checked the logging curves of the six wells and found that the trend I transition tends to take place if the overlying strata have a relatively large amount of clay and thus relatively low permeability to impede vertical fluid transfer. For instance, we plot well A (trend I) and well F (trend II) in Figure 9, from left to right: gamma ray (API), spontaneous potential (mV), the difference between neutron porosity and density porosity, and volume fraction of sand from petrophysical interpretation. All curves are at a sampling rate of 0.3 m (1 ft). The left three properties indicating high clay content present a high gamma ray, high spontaneous potential, and a big difference between the neutron and density porosity (Katahara, 2006). The distance between the two wells is approximately 50 km. Above the transition zones of the two wells (above the dashed line), the average volume fraction of sand (approximately 12%) in well A is much less than that in well F (approximately 55%), and the clay content in well A is higher than well F, e.g., neutron-density difference at 0.20 in average for well A and at 0.11 for well F. Within the transition zone (below the dashed line and above the dotted line), although there are some intercalated sandstone formations, the mineral compositions do not change much within well A and do not change much within well F.
The trend I or fluid expansion tends to takes place if smectite-rich shales have relatively high clay content and the shales can be first overpressured by disequilibrium compaction. The red circle in Figure 6a represents a speculated onset of overpressure (caused by disequilibrium compaction), which is shallower than the onset of the smectite-to-illite transition. Nevertheless, trend II or fluid loss tends to occur if smectite-rich shales are more compacted, and the released fluid from the interlayer water may escape to a large extent due to good pore connectivity. The overpressure onset (the red circle in Figure 6b) is deeper than the onset of smectite-to-illite transition.
The locations of wells of two trends of smectite-to-illite transition are not randomly distributed and follow a particular pattern as shown in Figure 10. Most wells with the trend I transition fall within the white triangle-like area, whereas wells with the trend II transition fall outside of the area. Further analysis can focus on checking clay content or permeability of layers above the onset of the smectite-to-illite transition and looking for seals.
Although the two trends are categorized and interpreted by fluid expansion or fluid loss, there are still lots of wells to be attributed to a combined result of these two trends, which represent the geologic reality (Lahann et al., 2001; Kataraha, 2014). Bruce (1984) also suggests that water expelled from smectite during the process of smectite-to-illite transition might migrate out of the host rock early or might be trapped partially or entirely and released gradually during the geologic time. In other words, if the smectite-to-illite transition occurs and the system allows a part of the released water to escape and another part of the released water to be preserved, this is insufficient to maintain the constant effective stress. Therefore, the interlayer hydration status and transferred water distribution might be first-order factors controlling the visible change in sonic measurements and seismic responses to the changes in rock properties and the concomitant overpressure.
Another goal of this paper is to link the micromechanism of the smectite-to-illite transformation at the grain scale to the corresponding physical property evolution of shale. To quantitatively assess the shale elastic property variation, we propose a framework to model the elasticity of the illite-smectite mineral with a various amount of interlayer water and elasticity of shale within the smectite-to-illite transition zone.
The smectite-to-illite transition involves the transfer of large amounts of interlayer water in between the clay T-O-T platelet to normal free water (Powers, 1967; Swarbrick and Osborne, 1998; Lahann et al., 2001). We presume that the evolution of elastic properties of shale during the smectite-to-illite transition can be approximated by a procedure of smectite clay dehydration (Burst, 1969). Figure 11 illustrates our framework for modeling shale elastic properties during smectite-to-illite transformation. The rock consists of a solid phase, interlayer water, and free water. First, we model the elastic properties. Second, we present bulk volume models. Third, we solve the pore pressure increase. The detailed steps are elaborated as follows:
Step 1: Modeling the elastic properties
Modeling the elasticity of the clay mineral with interlayer water: The elastic properties of the dry clay platelet and smectite are assumed to have intrinsic TI symmetry, including five elastic constants: , , , , and . This is because of (1) theoretically, the presence of strong covalent bonds within the clay hexagonal platelet layer (T-O-T) and the relatively weak electrostatic bonds between platelet layers (Sayers, 1994) and (2) experimentally, the observation of increasing confining pressure increasing the anisotropy of the shale samples (Tosaya, 1982), which is contributed from the intrinsic anisotropy of clay instead of the closure of aligned pores.The interlayer space consisting of water molecules and cations (e.g., Na+ and Ca2+) is between the clay platelets in Figure 2. The elastic properties of clay are influenced by the hydration status, which can be estimated through a layered medium model. With the increasing illite/smectite (I/S) ratio, the average interlayer thickness of clay will decrease during the smectite-to-illite transition because illite is nonexpandable with a basal layer spacing approximately 10 Å and smectite is expandable with a basal layer spacing ranging from 9.3 to 20 Å. Due to such tiny thicknesses of clay platelets compared with the seismic or ultrasonic wavelength, we use the Backus average (Backus, 1962) to characterize the smectite mineral’s stiffness. The anisotropic Backus averaging theory gives a transversely isotropic medium described by five effective stiffnesses (Backus, 1962; Mavko et al., 2009):where indicates averages of the enclosed properties weighted by their volumetric proportions. As an example, , where is the volume fraction of interlayer water in clay and and are the elastic constants of the interlayer water and dry clay platelets. The yielded TI characteristics depend mainly on the contrast between the elastic properties of the interlayer water and dry clay. To model the clay mineral properties during the smectite-to-illite transition, we next determine the properties of interlayer bound water and dry clay successively.(5)
Though the stiffness of interlayer water has not been measured directly and reported in the literature, the interlayer water property can be inverted through the theoretically calculated stiffness of Na-montmorillonite by Ebrahimi et al. (2012). They solve the quasistatic elastic stiffness with 0 up to 258 water molecules and the corresponding basal layer spacing varying from 9.3 to 19.6 Å. We can calculate the elastic stiffness of a two-layer medium comprised of dry clay (Na-montmorillonite without interlayer cation) and interlayer water by testing a series of values for bulk (2.25–15 GPa) and shear moduli (0−22.5 GPa) for interlayer water. Two inversion strategies are implemented. One is through seeking the best-fitting values of bulk and shear moduli of the interlayer water at each basal spacing (Figure 12) with error analysis (Figure 13). In the first inversion, horizontal stiffnesses and are underestimated (less than 15%), whereas the modeling vertical stiffnesses of and have tiny errors less than 1% and 5%, respectively. The second inversion is to invert the bulk and shear moduli of the interlayer water as invariants regardless of the basal spacing or hydration status (Figure 14). A notable difference between Figures 14 and 12a is the (blue triangles) at a basal layer spacing of less than 11.5 Å because the setting of the constant bulk and shear moduli of the interlayer water does not demonstrate a sharp decrease in from the dry state to the monolayer hydration. For other hydration states, two inversions present similar results. For inversion details, see Appendix A. The first inversion demonstrates that the elastic moduli of the interlayer water as a function of interlayer thickness are not monotonic but fluctuate approximately 10 and 2 GPa for bulk and shear moduli. Moreover, the second inversion obtains the best-fitting results of the bulk and shear moduli for interlayer water are 9.7 and 1.8 GPa.
To determine the elasticity of dry clay, we notice the chemical composition difference between Na-montmorillonite () and illite . In fact, during the smectite-to-illite transition, the reaction degree and illite/smectite (I/S) ratio increase progressively (Charpentier et al., 2003). For instance, there are several intermediate products by time sequences (Wilson et al., 2016): Initially a random I/S mixed-layer structure (smectite-dominant), then a regular 1:1 I/S interstratification, subsequently an ordered structure as IIISIIIS (illite-dominant), and finally pure illite.
Because the mixed illite/smectite can be found over the whole depth column (Hower et al., 1976; Bruce, 1984; Lahann et al., 2001; Wilson et al., 2016), we adopt the dry sheet silicates of the 1:1 ratio of illite-smectite () as the dry clay unit in the modeling. We neglect the elasticity variations of the dry clay during the smectite-to-illite transition because, first, the thickness variations in the T-O-T platelets are small, and, second, the interlayer hydration status is assumed to dominate the clay property evolution. Based on first-principle calculations and density function theory, Militzer et al. (2011) derive the elastic tensor for the dry sheet silicates of the 1:1 ratio of illite-smectite. The approximated TI elastic constants for , , , , and are 83.6, 49, 22.6, 31.5, and 19.3 GPa, respectively. These elastic constants are combined with the bulk and shear moduli of the interlayer water to simulate the effective stiffness of clay as a TI medium.
There is no agreement regarding the density of interlayer water, which depends on temperature, pressure, clay chemical composition, types of hydrated cations, and the number of water molecules/interlayer thickness. For instance, Bathija (2008) estimates the density of interlayer water at 1 GPa uniaxial stress, which is approximately for smectite with no less than one-layer interlayer water (32 water molecules) by using molecular simulation. Ebrahimi et al. (2012) estimate the density of the first layer (corresponding to a basal layer spacing of approximately 12 Å) at approximately at zero pressure. According to Powers (1967) and Colten-Bradley (1987), the interlayer water density would exceed , if the basal spacing is less than 20 Å. We apply 1.15 and for the densities of the interlayer water and dry clay in our modeling.
It is noteworthy that we use a different clay component model compared with Sayers and den Boer (2016). They use thickness fractions of the stiff layer (T-O-T structure) and the soft layer (interlayer) to express smectite clay as two isotropic layers. In our modeling, we use the thickness fraction of dry clay (smectite without water) and interlayer water, whose thickness is the difference between hydrated smectite and nonhydrated smectite. Dry clay is treated as a TI medium, whereas the interlayer water is treated to be isotropic. Although the chemical expression of clay and illite/smectite ratio can change and the exchanged cations may also affect the interlayer water properties, the varying soft component fraction will dominate the clay property evolution during the smectite-to-illite transition. In this intuitive way, the elastic properties of clay undergoing the smectite-to-illite transformation can be modeled through the interlayer water fraction because of the similar thickness between dry smectite and illite.
Modeling the elasticity of shale solid phase: Besides clay minerals, other minerals, such as quartz and feldspar, exist to form the solid phase of shale. However, due to the fine-grained (mud or clay size) nature and major components of the solid phase, clay minerals can serve as the load-bearing materials of a shale’s matrix. Quartz and feldspar can be considered as in a suspension status, surrounded by clay minerals. Therefore, we propose to use the self-consistent approximation in the anisotropic domain (Hornby et al., 1994) to generate the solid phase properties.
- iii)Modeling the elasticity of porous medium: We use the differential effective medium (DEM) method in the anisotropic domain (Nishizawa, 1982) to simulate a porous medium’s elastic properties because DEM treats each constituent asymmetrically. It can set smectite as a host material and add aligned pores filled with free water as inclusions. In the DEM algorithm, a small amount of inclusions of one phase is incrementally added to a background host medium in an iterative behavior, and the process is continued until the desired proportion of the constituents is reached. The change in effective elastic stiffness due to an increase of can be represented as (Nishizawa, 1982; Hornby et al., 1994)where denotes the stiffness of free water filled, represents the stiffness of the porous medium containing the solid phase and the pore-fluid water, and is a fourth-rank tensor given by the strain Green’s function integrated over the inclusion shape (Mura, 1987). The mechanically compacted pores are the inclusions having an ellipsoidal shape with an aspect ratio of approximately 0.08−0.25, which are typical ranges for interparticle pores of shales (Moyano et al., 2012; Zhao et al., 2016). The bulk modulus and density of water are set at 2.55 GPa and , respectively.(6)
Notice that our modeled clay mineral elasticity is at the grain scale, which has a symmetry perpendicular to the layered bedding. The clay aggregates vary in orientation but align locally, and a considerable portion of clay particles are randomly oriented. Overall, the degree of alignment presents a broad range of dispersion and depends on the compaction and diagenesis (Vernik and Liu, 1997). Ho et al. (1999) and Charpentier et al. (2003) observe that the alignment of clay minerals is progressively enhanced within a narrow depth window corresponding to the smectite-to-illite zone, which results from the effect of fluid loss and redistribution of effective stress. The more the released water escapes, the more aligned the clay aggregates become. Because fluid expansion does not involve much bulk volume change, the clay mineral alignment degree seldom changes. Although fluid loss leads to continued compaction, the clay minerals become more aligned. Hence, an orientation distribution function (ODF) is used to relate the clay alignment degree with the total pore space, which is the sum of pores filled with free water and interlayer space filled with interlayer water. Also, the ODF considers the randomization effect of quartz content on the clay-alignment degree. For the details of ODF and how to average the elasticity of clay aggregates, see Appendix B.
Step 2: Bulk volume models
It is necessary for us to use proper bulk volume models accounting for the scenarios of fluid expansion, fluid loss, and the combination of them. Bulk volume models are shown in Figure 11, before the smectite-to-illite transition, the shale consists of solid matrix and mechanical-compacted pores filled with free water (the light-blue area) denoted with “fw.” The solid phase is comprised of quartz and clay, including dry clay and interlayer water (the dark-blue area) denoted with “iw.” When the environment is conducive, smectite-to-illite transition occurs, and interlayer water begins to transfer to free water. The remaining interlayer water (the dark-blue area) in the interlayer space is denoted by “iwr,” reflecting the degree of transformation. The distribution of released interlayer water plays a key role in deciding whether the system is open or closed. The released interlayer water can either drain away (the white area denoted with “lw”) or be preserved in the pore space (the light-blue area denoted with “rw”).
Fluid expansion: If the released water does not drain and is preserved in the pores corresponding to the case of “fluid expansion” shown in Figure 11, the solid part dissolves and leaves a void filled with transferred free water (Katahara, 2007). The transferred interlayer water is replaced by soft pores with an equal volume to model the softening effects during fluid expansion (velocity decrease). The soft pores are filled with free water and have an ellipsoidal shape with an aspect ratio ranging from 0.005 to 0.02 due to the flat geometry of the interlayer space. In addition to the mechanically compacted porosity with an aspect ratio ranging from 0.2 to 0.5, there are two types of pores in the modeling for fluid expansion. This is in agreement with the rock-physics modeling work by Ruiz and Azizov (2011) and Moyano et al. (2012): Pores in shale can be considered to have a bimodal distribution of pore aspect ratios when effectively modeling the elasticity of shale.
Fluid loss: If the released interlayer water can fully escape as the “fluid loss” case shown in Figure 11, the total porosity (the sum of interlayer water porosity and free water porosity) decreases and the bulk density increases. After subtracting the volume of the interlayer water loss, no internal pore pressure is generated and the velocity changes trivially because of the increases in the bulk density and stiffness.
Combined case of fluid expansion and fluid loss: For a combined case, we modify the bulk volume model to accommodate that one part of transferred interlayer water escapes and another part is preserved.
Step 3: Estimating pore pressure increase
There are two major contributors to shale’s pore pressure increase: the gain in pore water volume during the smectite-to-illite transition and the continual burial or increased confining pressure for an undrained rock. We will solve each factor respectively and add them together.
The fluid expansion in the smectite-to-illite transition is conceptually modeled through a process of opening soft pores by the excess free water transferred from the interlayer bound water and accommodating the excess free water. The pore pressure is larger than the one before the smectite-to-illite transition so long as not all of the released water escapes out of the system or pure fluid loss occurs. Imagine that after the interlayer water transfers to free water in a closed system under constant confining pressure, the volume fraction of the solid framework decreases and the free water porosity enlarges. The original space taken by the solid interlayer water is replaced with liquid free water, and it continues to resist the confining pressure. Hence, an increasing part of the confining pressure acts on the pore fluid and the pore pressure increases. Meanwhile, if there is a density difference between the interlayer and the free water due to the water molecular packing, the excess fluid volume is different from the volume of the corresponding interlayer water. We can solve the pore pressure increase if the pore space accommodates the pore fluid, at an assumed constant density and compressibility for each component.
Figure 15 schematically illustrates how to solve pore pressure increase during the smectite-to-illite transition based on the contrast of volumes between the pore fluid and the pore space. In the whole process, we assume that the confining pressure does not change. Before the transition, we have a pore-fluid volume denoted by , which represents the initial free water porosity. The initial free water porosity equals the pore-space volume at the initial pore pressure , . In this paper, we use the italic letter “” in lowercase to refer to pore, the capital letter “” in italic to stand for pressure, and the letter “p” in lowercase to stand for the P-wave; the italic letter “” represents volume, and the capital letter “V” without italic denotes velocity.
Modeling results and analysis
The elastic velocity, density, and pore pressure increase () are modeled during the smectite-to-illite transition for a series of possible situations, and the results are displayed in Figure 16. The dark-blue and light-blue lines are the empirical loading curves for smectite-rich shale and illite-rich shale (Lahann et al., 2001; Dutta, 2002; Gutierrez et al., 2006). Some modeling parameters are set as follows: The initial compacted porosity is 17% with a pore aspect ratio of 0.1, the soft pore has an aspect ratio of 0.008, the initial quartz volume fraction is 5% in the whole rock, and the grain density of illite is .
Measured well-log data can be interpreted using our modeling result, which is composed of the scattered points linked by dotted lines. Each point has its particular status of the interlayer hydration and the distribution of released water. The dark-blue point at the upper-left corner without the pore pressure increase is the initial point, and the clay has an equivalent basal layer spacing of 11.5 Å. The value represents the smectite takes 75% of the mixed illite/smectite because the thickness of smectite with one-layer water molecules and illite is 12 and 10 Å. With a sampling of 0.4 Å, the interlayer thickness gradually decreases to 10.3 Å, corresponding to a volume fraction of smectite in clay at 15%, to simulate the smectite-to-illite transition. Additionally, the transition is set to be implemented within a depth range of 600 m, and each neighboring point has a depth increase of 0.2 km. A confining pressure gradient at is adopted to account for the depth difference. The two black arrows in Figure 16 indicate how the vertical P-wave velocity () and bulk density () move toward the illite loading curve if pure fluid expansion (downward) or fluid loss (rightward) happens in shale. The other three scattered points that do not fall along the pure fluid expansion/loss trend denote combinations of the two cases when some of the released interlayer water is preserved while the other escapes.
In the given example, the modeling results suggest that the P-wave velocity can decrease by 11% and density reduces little () due to the pure fluid expansion, and the increase of the P-wave velocity (2%) is less than the increase of density (6%) due to the pure fluid loss. In the case of the pure fluid expansion, the velocity decrease is significant because it is dominated by the excess free water transferred from the interlayer water; the density reduces a little because the porosity of shale is increased slightly by the increased pore pressure and constant confining pressure. By contrast, in the case of the pure fluid loss, the density increase is significant, which mainly results from the escaped water transferred from the interlayer water. Meanwhile, the fluid loss also leads to an increase of shale stiffness because of the loss of the soft component and pore water. Nevertheless, the P-wave velocity is slightly increased because the increased bulk density is counteracted by the increased stiffness. The evolution of the P-wave velocity and the density of shale in the transition zone is attributed to the hydration status of clay and the distribution of released interlayer water.
The pore pressure increase in megapascals is color coded. Two primary factors result in the increased pore pressure: the fluid expansion during the smectite-to-illite transition and the increased confining pressure. The modeling results show that compared with the initial point, the pore pressure is increased in all possible states but most significantly in the pure fluid expansion (up to 20 MPa). During the transition, the interlayer water is released as free water and may drain or be preserved in the rock. The more fluid that is trapped, the higher pore pressure that is generated according to equation 13. For the case of pure fluid loss, no excess pressure is generated by the fluid phase change. However, the increased confining pressure will contribute to the pore pressure increase (up to 14 MPa), and its pore pressure increase approximately equals the increased confining pressure because the Skempton’s coefficient is close to one for the rock saturated with water.
Figure 17 presents examples of calibrating our modeling results with well-logging data. The average quartz content (5% for well A and 14% for well F) of the shale within the smectite-to-illite transition zone from logging curves in Figure 9 is chosen as the modeling parameters. In Figure 17a, the upscaled P-wave velocity () and bulk density () are from well A in Figure 7a, and the pore pressures (MPa) are color-coded and approximated with the mud weights from its logging run. Our modeling results in the circles are calibrated with the data point at a depth of 2652 m where the velocity reversal starts. Pore pressures are estimated by adding the pore pressure at 2652 m to our modeled pore pressure increase. For well A, its velocity and density data trend can be interpreted using our modeling results. The data trend shifts downward from the dark-blue curve at shallow depths to the light-blue curve at a deeper section of the well, caused by the fluid expansion in the smectite-to-illite transition zone. With an increasing amount of interlayer water released and trapped, the P-wave velocity decreases, the density decreases a little, and the pore pressure increases. The transition zone spans approximately 500 m, corresponding to an increase of confining pressure. Then, the excess pore fluid tends to be squeezed out of pore space with the continuing compaction, leaving a trend of moving along the illite loading curve. The velocity and density data during reloading (deeper than 3150 m) match the modeling trend moving from the pure fluid expansion to the pure fluid loss case. The modeling points near the illite loading curve share the least amount of interlayer water and are modeled to lose the released water and close the soft pores gradually. The pore pressure of the reloading zone is higher than the one modeled losing the released water because the modeled results do not have a large depth increase as the logging data do.
Factors affecting the geophysical signatures of the smectite-to-illite transition
For the elastic characteristics of the trends I and II smectite-to-illite transition, though they have distinct signatures on density-slowness crossplot, they both allow “velocity trends deviating from the density trends” and will follow the illite compaction trend eventually. Their geophysical signatures depend on their geologic settings, and several factors might have impacts on the magnitude of unloading overpressure. Increased amounts of the initial smectite content, burial rate, or geothermal gradient can increase pore pressure. Accordingly, an increasing clay content leads to a decrease of permeability, which prevents the pore pressure from decreasing. Geologic time will determine the status of the clay diagenesis and pore pressure.
Figure 18 schematically illustrates that for two shales that have the same age, their pressure gradient development with respect to time can be different. In the past, there was no fluid volume increases in the two smectite-bearing shales. Because both of them continue to compact with depth, unloading caused by the smectite-to-illite transformation occurs. Their peak overpressure is determined by the initial smectite content and factors relating to leakage and can be reached at a similar time. Their pore-pressure relaxation time, during which overpressure reduces to normal hydrostatic pressure, depends on the permeabilities of overlying shaly formations if no lateral fluid migration occurs. We make an observation of the pore-pressure gradient at the present time marked by the dashed line. For the blue curve, trend I shale’s pore pressure is close to its peak overpressure and still needs more time to get equilibrium. However, for the red curve, a trend II shale has experienced its peak overpressure time and the pore pressure is close to restoring to hydrostatic pressure. Eventually, both shales’ pore-pressure gradients reach equilibrium when the pressure dissipates. We conceptually categorize the transition into one of fluid expansion, fluid loss, or a combination of these two to perform a first-order investigation of the impacts of smectite-to-illite transformation on the elastic responses of shale and pore pressure increase. Substantially, these cases involve the same chemical reaction but present different water distribution statuses due to various geologic settings.
Cautions for modeling development
Regarding the strategy in modeling development, several points are worth our attention. These are listed as following:
According to the reaction of the smectite-to-illite transformation (equation 1), comparing the constituents’ elastic moduli before and after the transformation, we approximate the evolution of elastic properties and density during the smectite-to-illite transition to a process of smectite clay dehydration (Burst, 1969). Because, first, the difference in elasticity between quartz and K-feldspar is much smaller than the elasticity difference between smectite and illite in addition to free water. Referring to Mavko et al. (2009), the solid grain densities for quartz and K-feldspar are approximately 2.65 and on average, and the bulk moduli for quartz and K-feldspar are both near 37 GPa, whereas the shear modulus of quartz is higher than that of K-feldspar. The reported grain density of smectite ranges from 2.0 to due to the uncertain hydration status, whereas illite’s grain density is larger than smectite’s and is in the range of 2.6– (Castagna et al., 1993; Katahara, 1996; Totten et al., 2002; Chitale, 2010; Hornbach and Manga, 2014). Moreover, for the clay structure, platelets of smectite and illite have a T-O-T structure and similar thickness. Last but not least, the shear rigidity of the interlayer water indicated by the dehydration/hydration process (Ebrahimi et al., 2012) is close to the interlayer shear modulus when the clay surfaces are charged (Kolstø and Holt, 2013), which is closer to the condition of the smectite-to-illite reaction.
From the perspective of basin modeling, the Arrhenius equation is used to characterize the reaction rate of the smectite-to-illite transition to be dependent on the formation temperature and burial time (Pytte and Reynolds, 1989). It is pertinent to emphasize that this diagenesis in shale, involving physical and chemical changes during millions of years of geologic time, is an irreversible and highly nonlinear process. The elastic properties of minerals at crystalline scale should vary continuously and complexly upon the smectite-to-illite transformation, instead of an abrupt way of dehydration as we simplify in the modeling steps. In our modeling, we use an effective property of clay mineral at the crystalline scale, basal layer spacing, to account for the averaged statistical results in the progress of transition. Although it is difficult to have direct measurements that indicate the basal layer spacing, the calibration of the bulk density and porosity from the logging data as well as the empirical loading curves of smectite-rich and illite-rich shales improves our confidence in the modeling results.
Cautions for modeling applicability
Regarding the applicability and predictive power of the established modeling framework, several points are worth attention as listed following:
Our model does not account for the porosity and pore-geometry change developed during compaction. A constant mechanical compacted porosity is set as an input parameter. Combining the mechanically compacted pores with its aspect ratio, a starting point without the smectite-to-illite transition is modeled and calibrated by the empirical compaction trend of smectite-rich shale. The initial clay basal layer spacing is constrained by the free water porosity, and we compute a series of points with decreasing basal spacing to speculate on how rock properties change during scenarios of fluid expansion, fluid loss, and their combination.
Crack-like soft pores are opened little by little with an increasing amount of released interlayer water during the smectite-to-illite transition. The effective pore compressibility upon pore-pressure change () will increase during the smectite-to-illite transition. However, solving the pore pressure increase in equation 13 uses a final regardless of the intermediate evolution of . Thus, the pore pressure increase is underestimated from this perspective.
Our model considers the pore-pressure generation and shale elastic properties at a passive basin condition, which does not exhibit much tectonic activity. If the largest principal stress is not in the vertical direction, the compaction behavior is different from that in a passive basin (Hauser et al., 2014). At a local scale, a salt dome may also change the stress distribution and therefore the pore-pressure accumulation, so the predictive power of the model under the effects of salt may be affected. Besides, the overburden stress gradient is assumed to be a constant in the model. Additional procedures may be required before applying the model to deep-water basins, where the overburden stress gradient increases faster with depth than some other geographical settings.
The analysis of the density-sonic velocity crossplot can be used to elucidate the pore-pressure evolution mechanism during the smectite-to-illite transition. We categorize two trends of smectite-to-illite transition from the wells with unloading overpressure. The density changes little while the velocity decreases for trend I, and the density increases while the velocity changes little for trend II. Trend I corresponds to the fluid expansion, in which the released interlayer water cannot escape out of pore space; trend II corresponds to the fluid loss, in which the released interlayer water drains out of the system. The occurrence of the two trends corresponds to different geologic settings: The formations overlying above trend I tend to have lower permeability or higher clay content than that above trend II. The local variation of clay content may cause heterogeneity of the pore-pressure distribution.
We also present a framework to model the elastic property evolution characteristics of shale within the smectite-to-illite transition zone. In particular, fluid expansion, fluid loss, and their combination are depicted through different bulk volume models. By using the proposed modeling strategy, our modeling result can quantitatively interpret the change of elastic properties and pore pressure. Modeling results are in agreement with the observed trends from the velocity and density data. The calibration of the model with logging data indicates that the modeling framework allows us to better understand the velocity and density variation as well as the pore pressure increase in the transition zone. The modeling framework presented provides a theoretical method to supplement conventional empirical models and quantitatively interpret how the shale velocity, density, and pore pressure evolve if smectite-to-illite transition occurs.
This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDA14010203) and the National Natural Science Foundation of China (grant no. 41874124). The authors also would like to thank the Fluids/DHI consortium for their financial support. Authors express their gratitude to F. Hilterman and Geokinetics for the use of their well-log library and D. Ebrahimi for the use of his calculated stiffness of Na-montmorillonite with various numbers of interlayer water. Authors would also like to thank K. Katahara, S. Tai, H. Yan, and H. Yu for their discussions.
DATA AND MATERIALS AVAILABILITY
All of well logging data are courtesy of Fred Hilterman. People can contact him to access the data.
INVERSION OF INTERLAYER WATER ELASTICITY
From the modeling results, we observe that the bulk and shear moduli of interlayer water as a function of basal spacing are not monotonic but fluctuate around 10 and 2 GPa. According to Figures 13a and 14, and are underestimated (errors of less than 15%), which reflects the deficiency when we use a method based on the Backus average, to model the horizontal stiffness of Na-montmorillonite from molecular simulation. This is caused by the Coulomb force between solid clay and interlayer water. Modeling results of and have small errors less than 1% and 5%, respectively, which gains our confidence when characterizing the vertical stiffness of smectite with this method. Hereby, we do not use and model because it usually has larger uncertainty than other elastic constants (Yan et al., 2015).
AVERAGING WITH PARTICLE ORIENTATION DISTRIBUTION
In a TI medium, it is sufficient to rotate the elastic tensor of crystalline aggregates by a series of discrete successive rotations and average the results. The polar angle varies from to with an interval of , and the azimuthal angle varies from 0 to with an increment of .