Deformation characteristics of a clayey interbed during fluid injection

Engineering Geology 183 (2014) 185–192

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

Technical Note

Deformation characteristics of a clayey interbed during fluid injection Xuejun Zhou ⁎, Thomas J. Burbey National Energy Technology Laboratory—Regional University Alliance (NETL-RUA), Department of Geosciences, Virginia Tech, United States

a r t i c l e

i n f o

Article history: Received 1 October 2013 Received in revised form 17 September 2014 Accepted 5 October 2014 Available online 14 October 2014 Keywords: Ground deformation Fluid injection Cam clay model Poroelasticity FEM Abaqus

a b s t r a c t Surface deformation due to fluid withdrawal has long been observed at the surface above aquifers and oil reservoirs. Uplift associated with fluid injection has also been observed. Although the processes of subsidence and uplift are reversible in a poroelastic setting, the presence of clayey interbeds can result in ground deformation behavior non-reversible because of their low permeability and nonlinear behavior. In this investigation a Cam clay model is used in conjunction with a poroelastic model to simulate the presence of a laterally extensive lens-shaped clayey interbed in a sandstone aquifer during fluid injection. Spatial and temporal changes associated with this interbed during aquifer pressurization are captured and can be clearly differentiated from the deformation of the surrounding poroelastic aquifer formation. Results for pore pressures and ground surface deformation patterns can be categorized into three distinct time intervals: an early time interval where the aquifer is pressurized but not the interbed leading to a lower region above the interbed; an intermediate time interval in which the interbed pressures slowly increase, approaching the pressure of the adjacent aquifer, leading to a potentially large surface uplift above the interbed; and a late time interval in which pressure equilibration is achieved between the interbed and aquifer and a highly non-symmetrical irregular surface deformation pattern results, which may be higher or lower than the background depending on the interbed characteristics. Reservoir configuration is found to be an important factor influencing ground deformation behavior, with more obvious deformation always occurring in a laterally confined aquifer as opposed to a laterally infinite aquifer. Published by Elsevier B.V.

1. Introduction Land deformation due to fluid withdrawal or injection has long been observed over many oil reservoirs and aquifers systems (Tolman and Poland, 1940; Yerkes and Castle, 1969; Vasco et al., 2001; Bell et al., 2002; Teatini et al., 2011a,b). As the prospect of CO2 sequestration progresses as an important environmental issue, gaining knowledge about aquifer-system response to waste storage becomes vitally important (Vasco et al., 2010). With the development of sophisticated surface monitoring techniques such as InSAR and GPS, careful monitoring of surface deformation signals during fluid injection can yield critical information about the geometric configuration and characterization of the host reservoir (Fielding et al., 1998; Vasco, 2004). For historical reasons, the literature associated with land subsidence due to fluid withdrawal is far more abundant than those associated with land uplift due to fluid injection. In a poroelastic hypothetical framework, these two processes are opposite to each other and reversible. Exceptions include formations with high clay content because clays often result in non-reversible and nonlinear deformation behavior.

⁎ Corresponding author at: Mewbourne School of Petroleum & Geological Engineering, The University of Oklahoma, 100 E. Boyd Street, SEC 1210, Norman, OK 73019, United States. E-mail addresses: [email protected], [email protected] (X. Zhou).

http://dx.doi.org/10.1016/j.enggeo.2014.10.001 0013-7952/Published by Elsevier B.V.

These systems are more difficult to evaluate than purely elastic systems because a hydrodynamic (time) lag response typically occurs in systems containing clayey interbeds due to their low permeability. Consequently, their response to an imposed confining stress is dependent on the preconsolidation stress history because their virgin compaction curve and rebounding curve often do not coincide. In this investigation an interbed refers to a relatively thin lensshaped unit composed of fine-grained clay material with a much lower permeability and relatively higher porosity than the surrounding host aquifer formation, which is composed of sandstone (Figure 1). The hydrodynamic lag has been observed during fluid withdrawal as the clayey interbeds drain more slowly than the surrounding coarser aquifer material. This leads to a time-delay between the extraction of the groundwater and the occurrence of land subsidence (Shearer, 1998). The time lag for consolidation can range widely from several days to many decades or even longer (Bell et al., 1986; Hoffmann et al., 2001; Li et al., 2006; Shi et al, 2007). This hydrodynamic lag effect may also be expected during fluid injection, but it is unclear whether the mechanics of the deformation response will be similar to that of groundwater withdrawal. The goal of this technical note is to investigate the hydrodynamic lag and the surface deformation response caused by clayey interbeds during fluid injection and to correlate the deformation response with the pore pressure evolution in both the clayey interbed and the aquifer.

186

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

and Bentheim sandstone, respectively (Benson et al., 2005; Louis et al., 2005). The mechanical behavior of stiff soil and clayey rocks during loading and unloading are similar in many aspects. This includes an initial stiff response until the normal compression line is reached, followed by an increase of the state boundary surface and followed then by a stiff response upon unloading. Hence, instances occur in which the behavior of clayey rocks can be better described within the theoretical framework of a critical state model of soil mechanics (Vukadin, 2007). A Cam clay model is used in this investigation to simulate the hydromechanical response of a clayey interbed during fluid injection. A short introduction on the mathematical framework for the Cam clay model follows. 3. Theoretical background Several different approaches can be used to mathematically characterize aquifer heterogeneity. One approach is to use Eshelby's inclusion theory (Eshelby, 1957; Mura, 1982; Segall and Fitzgerald, 1998), which is applicable for an aquifer embedded with clayey interbed as investigated in this paper. In such a system, pore pressures increase differently in the different hydrologic units during fluid injection due to the different permeability of each unit, which leads to differing stain histories of each unit during the pressure change. Because clayey formations behave differently than the coarser-grained aquifer formation, we simulate an interbed as an inclusion in an otherwise homogenous poroelastic aquifer formation using a Cam clay model. 3.1. Cam clay model

Fig. 1. Regional and pore-scale illustrations of clayey interbeds within an aquifer system. (Modified after Leake, 1990).

2. Properties of clayey rock Clayey rocks (clays, claystone, shale, mudstone) represent an important constituent of sedimentary basin fill (Hildenbrand and Krooss, 2003). Many of these rocks are considered to be soft rocks that lie between hard soils and (moderately) strong rocks. The distinction between soil and rock is to a certain extent arbitrary because the processes of weathering and diagenesis are gradual (Nova, 2010). The boundaries between soils and rocks as set by the uniaxial unconfined compressive strengths vary by different organizations and researchers. These strength values can range from less than 1 MPa to that over 25 MPa (Bieniawski, 1973; IAEG, 1979; BSI, 1981; ISRM, 1982; Hawkins and Pinches, 1992). Laboratory testing on some clay or clayey rocks show that these rocks' permeability can range greatly from 1 nD to 1 mD (values on the lower end are more common), their porosity can range from 0.03 to 0.55, Young's modulus can range from 0.1 to 30 GPa and Poisson's ratio can range from 0.13 to 0.38 (Neuzil, 1994; Giraud and Rousset, 1996; Dewhurst et al., 1998; Horsrud et al., 1998; Gale et al., 2007; Gasparre et al., 2007; Hight et al., 2007; Loucks and Ruppel, 2007; Engelder et al., 2009; Sarker and Batzle, 2010; Yang and Aplin, 2010; Zhou et al., 2010; Fidler, 2011; Vermylen, 2011; Eshkalak et al., 2013).One may note that their Poisson's ratios that were acquired in the lab may be undrained ratios in some cases. Their permeability also exhibits a higher degree of anisotropy (Zhang, 2005). For example, the ratios of the permeability measured parallel to bedding over those measured perpendicular to bedding for the Wilcox shale are greater than 10 (Kwon et al., 2004), whereas this ratio is only 4 for Berea sandstone (Zoback and Byerlee, 1976) and 2.2 and 1.2 for Crab Orchard sandstone

The ability of the Cam clay model to accurately simulate the behavior of a clayey formation has been demonstrated in the past (e.g. Helwany, 2007). A Cam clay model is chosen because it can simulate the entire range of expected behaviors of plastic materials (clays) that include (1) a yield criterion that predicts whether the material will respond elastically or plastically in response to a loading increment, (2) a strain hardening rule that controls the shape of the stress–strain response during plastic straining, and (3) a plastic flow rule that determines the direction of the plastic strain increment caused by a stress increment (Roscoe et al., 1958; Alonso et al., 1990). The total volumetric strain can be decomposed as: rate in the clayey inclusion dεinclusion vol inclusion

dε vol

el

pl

¼ dεvol þ dεvol

ð1Þ

pl where dεel vol and dεvol represent elastic (recoverable) strain rate and plastic (non-recoverable) strain rate, respectively. The loading function retained to account for the yield surface of clays is expressed in the form (Coussy, 2004):

2f ðσ ; q; σ co Þ ¼


2 1 q2 1 2 σ− σ co þ 2 − σ co 2 4 M

ð2Þ

where σ is equivalent pressure stress (or mean effective stress), q is the deviatoric stress, σco is the maximum effective stress (or the preconsolidation stress) as used in Eq. (2), and M is a material constant. A typical loading–unloading curve for a Cam clay model is shown in Fig. 2. Two key parameters of the Cam clay model are identified in Fig. 2, i.e., w (logarithm plastic bulk modulus) and m (logarithm elastic bulk modulus). w is the slope of the normal consolidation line, defined by w = − de/d(ln σ), and m is the slope of the unloading line, defined by m = − deel/d(ln σ), where e is void ratio and eel is its elastic component. The value of m depends on the type of clay and usually ranges between 5% and 25% of the value of w (Nova, 2010). In this investigation, we focused on the coupled hydrological– mechanical processes of a clayey interbed during fluid injection. The

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

187

a radially infinite reservoir of an elastic diffusive solid, the pore pressure increment can be expressed as (Rudnicki, 1986; Wang, 2000): Q r2 pðr; t Þ ¼ E1 4ct 4πðkb=μ Þ

Fig. 2. Loading–unloading curve of a Cam clay model (after Jorgensen, 1980).

fluid injection will increase the pore pressure in the aquifer as well as within the interbed (which is dependent on both space and time), thus reducing the effective stress. Subsequently the clayey interbed will undergo an unloading process that differs from that of a reduction of compressive pressure as conducted in laboratory tests or observed from overburden erosion in geological conditions because the reduction in the effective stress upon fluid injection is dependent on an increase of the fluid pressure instead of a reduction of the total stress. 3.2. Fluid injection in a poroelasticity domain The simulation of the aquifer is based on poroelastic theory. The general constitutive equations for a continuous porous field govern both the deformation of the elastic porous medium and the diffusion of fluid through the pore space of the material (Rice and Cleary, 1976; Berryman and Milton, 1991): εi j ¼


 1 ν 3ðν u −ν Þ pδ σ i j− σ kk δi j þ 2G 1þν 2GBð1 þ νÞð1 þ ν u Þ i j

Δm f ¼

B¼


 3ρ f ðνu −νÞ 3 σ kk þ p 2GBð1 þ νÞð1 þ ν u Þ B

1=K−1=K s   1=K−1=K s þ n 1=K f −1=K n

ð3Þ

ð4Þ

ð5Þ

where δij is the Kronecker delta, B is Skempton's B coefficient as defined by Eq. (5), νu and ν are undrained and drained Poisson's ratio, respectively, G is shear modulus, p is pore pressure, K is drained bulk modulus of rock, Ks is solid grain's bulk modulus, 1/Kn is the unjacketed pore compressibility, Kf is the pore fluid bulk modulus, and n is porosity. Δmf is the change of fluid mass. The hydro-mechanical coupling is accomplished through the governing equation of mass conservation for the infiltrating pore fluid:   ∂qx =∂x þ ∂qy =∂y þ ∂qz =∂z þ ∂m f =∂t ¼ 0

ð6Þ

∂qx/∂x represents the divergence of the water flux along a coordinate direction x. When an injection source occurs as a continuous line source in

! ð7Þ

where k is the intrinsic permeability and relates with hydraulic conductivity κ by k = κμ/ρfg, μ is the viscosity of fluid, the exponential integral E1 is “the Well Function” in hydrogeology, r is the distance from the injection well, b is the confined aquifer thickness, c is the hydraulic diffusivity, and Q is the injection flow rate. This equation only applies to a horizontal well at very early times, before the pore pressure front reaches the top and bottom of the reservoir layer. The pore pressure behavior in the horizontally layered system as considered in this note is more complicated and a numerical approach is needed. Analytical and numerical methods as applied to the effects of fluid withdrawal on land subsidence have been studied in the past by considering different boundary conditions (Geertsma, 1966, 1973; Gambolati and Freeze, 1973; Gambolati, 1974; Safai and Pinder, 1980; Helm, 1994; Burbey, 2006). In this paper, we assume a single phase flow condition. A single-phase flow is also a good approximation for the pore pressure evolution for a multiphase flow scenario over a long period of time (Zhou et al., 2009; Chang et al., 2013). However, a multiphase flow simulation is always desired for a detailed research associated with CO2 sequestration. 4. Conceptual model Abaqus (Hibbitt et al., 1998) was selected as the numerical tool for simulating our interbed-aquifer system (Figure 1) as it is a sophisticated finite-element software package capable of handling the hydromechanical coupling problem of ground deformation due to fluid injection. Abaqus also contains a large material library for poroelastic and Cam clay models. A plane strain model with dimensions of 13,000 m × 300 m (Figure 3) is developed to simulate our interbed-aquifer system. The real problem is more likely to be a three dimensional problem. However, a 2D model has the advantage to cover a large extent with high resolutions for the purpose of demonstration. A 2D plane strain model may also be applicable to some real cases if the injection occurs through a horizontal well. The overburden has a thickness of 200 m, and the aquifer has a thickness of 100 m. The interbed is located at the center of the aquifer and has a flat hexagon (lens) shape with a thickness of 20 m at its thickest extent and extending for a length for 800 m. The interbed thins at both ends (from 20 m to 0 m over a distance of 100 m along each end) and extends a total of 1000 m, which is long enough to avoid the stress arching effect. We monitor the ground deformation and aquifer/interbed behavior through the survey lines AA′, BB′ and CC′ (Figure 3). The injection wells is located at B, which is at the half depth of the aquifer. AA′ is on the ground surface; BB′ is a horizontal line through the interbed in the aquifer and CC′ is a vertical line through the middle of the interbed and the aquifer (1500 m from the injection well). AA′ and BB′ extend to a horizontal distance of 3000 m from the injection well. The boundary conditions are represented by rollers on both the left and right sides of the model and the bottom of the aquifer is fixed in the vertical direction. The input parameters for a baseline analysis are summarized in Table 1 for the overburden, aquifer and interbed. Several simplifications were made to both restrict the complexity of our model conceptualization and for the purpose of better identifying individual unit behavior associated with fluid injection. The first assumption is that a single interbed is used. In reality an aquifer system may contain numerous interbeds and their distribution can also be highly complex. However, a single interbed allows us to understand such formation's behavior more clearly and in a straightforward manner. Additionally a single elastic property is assumed for the entire

188

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

Fig. 3. Schematic model configuration.

overburden. In reality an overburden may consist of many different layers and some layers may have properties similar to that of interbed and can also absorb the injected fluid (Zhou et al., 2009). In this study, the overburden serves a single purpose of adding a gravity weight to the aquifer and can reflect the ground deformation after fluid injection. Thirdly, the hydraulic conductivity of the interbed is likely to be different in its vertical and horizontal directions; however, as these parameter values are all supposed to be much lower than that of the aquifer, only one hydraulic conductivity value is assigned to the interbed based on the understanding of a clay rich rock (Neuzil and Pollock, 1983; Neuzil, 1994). We note that it is only a trivial extension to include anisotropic permeability. In a poroelastic model as used in Abaqus, Young's modulus E relates with the logarithmic bulk modulus m, Poisson's ratio ν and elastic tensile strength pel t by (Hibbitt et al., 1998):    3ð1−2ν Þ 1 þ eel  el E¼ σ þ pt m

ð8Þ

Table 1 Material properties for overburden, aquifer, and interbed. Overburden Density (kg/m3) Young's modulus (E) Aquifer Density (kg/m3) Void ratio Fluid specific weight (kN/ m3) Hydraulic conductivity (m/ s) Interbed Density (kg/m3)

2000 95 MPa

Poisson's ratio

0.25

2000 0.23 9.81

Porous elasticity Log bulk modulus (m) Poisson's ratio

0.040 0.20

1 × 10−5

Elastic tensile strength(pel t )

7.0 MPa

1900

Fluid specific weight(kN/ m3) Initial void ratio (e)

9.81

Hydraulic conductivity (m/ s) Porous elasticity Log bulk modulus (m)

1 × 10

Poisson's ratio Elastic tensile strength (pel t )

0.25 0

0.0102

−12

Plasticity Log plastic bulk modulus (w) Stress ratio, M Initial yield surface size (kPa) Wet yield surface size a Flow stress rate b

where eel is the elastic change of the void ratio. A careful selection of the input parameters for poroelasticity allows for significant flexibility in simulating the sedimentary rock's poromechanical behavior so that it can be represented, for example, as Indiana limestone and Berea sandstone (Zhou and Burbey, 2012). Young's modulus of a clayey interbed is sensitive to the confining stress and the interbed formation tends to be stiffer at higher confining stresses, which agrees with field observations. Young's modulus of the clayey interbed is estimated as 0.3 GPa in this investigation. The stress within the aquifer must be equal to the lithostatic stress prior to fluid injection. The initial stresses are geostatic: the vertical stress in the aquifer is equal to the weight of the overburden, and the horizontal stress is equal to 0.63 of the value of the vertical stress. The ratio of 0.60 is a typical lower bound for horizontal stress to ensure a stability state (Zoback, 2007). We increased the overburden weight in a stepwise manner until the total overburden weight was achieved. The void ratio of the clayey interbed was decreased from 0.60 to 0.537 during this initial in-situ stress installation procedure. The model is allowed to be drained over sufficient time and consolidation is allowed to occur. Next, we inject fluid at a constant flow rate of 1.0 × 10−3 m3/s per meter normal to the model plane over 1.21 × 108 s (1406 days or 3.85 years). Here we consider two different scenarios. The first (Case A) is a large (“infinite”) reservoir (essentially no pore pressure changes along the distant lateral boundary) and the second (Case B) is a semiconfined reservoir with a permeability boundary located 3000 m from the injection well. In the first scenario we assign a constant pore pressure at the distant lateral boundary (13 km away from the injection

0.60

0.110 1.0 5.0 1.0 1.0

a The wet yield surface size is a parameter that defines the size of the yield surface on the “wet” side of the critical state. It is used to modify the shape of the yield surface on the “wet” side of the critical state, so the elliptic arc on the “wet” side of the critical state has a different shape from the one on the “dry” side. b Flow stress rate is the ratio of the flow stress in triaxial tension to the flow stress in triaxial compression and determines the shape of the yield surface in the plane of principal deviatoric stresses (Hibbitt et al., 1998).

Fig. 4. Pore pressure boundary conditions as indicated by point B′ (3000 m away from injection point) during fluid injection for both cases (inset is the figure in logarithmic time scale).

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

189

constant pore pressure boundary will occur once leaking commences. The threshold of pore pressure increase is determined by many factors including the minimum confining stress, the rock properties, fracture distributions, etc. For the purpose of comparing results between the large reservoir (Case A) and a small semi-confined reservoir (Case B), we apply a constant pore pressure of 4.13 × 105 Pa on day 4 along the interface of the permeable and less permeable regions (Line DD′ in Figure 3) until the end of the fluid injection. Fig. 4 shows the pore pressure change measured at B′ for Case B. Clearly, pore pressure levels for Cases A and B are significantly different as summarized in Fig. 4. Consequently, simulated ground deformations are also very different even if the injection schemes are the same for the two cases. 5. Numerical testing results Fig. 5. Pore pressure as measured at 1500 m away from the injection point in the aquifer and the midpoint of the interbed (time in log scale) during fluid injection.

well). The pore pressure change as measured at B′ (3 km away from the injection well) will exhibit logarithmic behavior (Case A in Figure 4), i.e., the pore pressure at early time will increase sharply and then level off after a long period of time. This agrees with the behavior as determined by Eq. (7) with r = 3000 m. In the second scenario (referred as Case B hereafter), we reduce the permeability in the aquifer region east of line DD′ in Fig. 3 from 10−5 m/s to 10−9 m/s to simulate a semi-confined reservoir that extends outward beginning at a distance 3000 m from the injection well. The pore pressure will increase linearly upon a constant injection rate in such a configuration because of the impedance effect of the low permeability region. However, the pore pressure cannot be increased arbitrarily high either in reality or in numerical modeling. In reality, this is due to the leaking effect that can occur at high pore pressures, which can be caused by the increased permeability of adjacent formations or reopening of some preexisting fractures, etc.; in numerical modeling high pore pressures can result in ubiquitous tension failures at early time steps. A reasonable assumption is that the pore pressure will increase linearly at early time periods and then level off because a

As injection begins the pore pressures in the aquifer and the interbed will be increase; however, a hydrodynamic lag effect in the interbed occurs because of its much lower permeability than the aquifer (Figure 5). Because of the interbed's low permeability, the pore pressure remains almost unchanged inside the interbed during the early stages of injection. The pore pressure in the interbed lags behind that of the aquifer until enough time has allowed for equilibrium of the interbed with the aquifer formation. The pore pressures follow different pathways in the aquifer and interbed for both cases. The pore pressure distribution can be characterized by three distinct time intervals (Figure 5); the early time period (b 2e5 s or 55 h), intermediate time period (2e5–1e7 s or 55 h–115 days) and the late time period (N1e7 s or 115 days–3.85 years); at early time, pore pressures in the aquifer increase rapidly, while little change occurs in the interbed; during intermediate time pore pressure in the aquifer has neared equilibrium in response to the lateral boundary condition, while pressures just begin to increase in the middle of the interbed; at the late time, pore pressures in both domains converge and tend to be stabilized even as injection continues. In fact, we even observe minor pore pressure increase (b 200 Pa) during the early time interval for both cases. The interbed pore pressure is increased at such an early time not because of incoming fluid but

Fig. 6. Pore pressure profiles along the vertical survey line CC′ (Figure 3) at different time steps for Case A (left) and Case B (right). A curve of an earlier time always appears at a lower level than a later time. The last two curves at time steps of 1.15e7 and 2.15e7 s cluster together for both cases.

190

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

Fig. 7. Comparison of void ratio change over effective stress for both cases. Note that the inelastic compaction portion of the curve corresponds to the geologic time-scale initialization of the model, while the elastic unloading portion corresponds to the much short injection phase.

because the interbed is squeezed by the expanding aquifer. This conclusion is verified by porosity changes occurring in both units. Fig. 6 shows the pore-pressure distribution in the interbed over time along the vertical survey line for both cases. Pore pressure gradients are much higher near the interface of the aquifer and interbed and the larger gradients occur in Case B rather than in Case A.

Increased pore pressures result in decreased effective stresses, which is an unloading process. Fig. 7 shows the simulated stress–strain curves prior to fluid injection and the swelling curves during fluid injection as measured in the interbed for both cases. The relation between the change in void ratio (strain) and change in effective stress resembles the curves shown in Fig. 2. The compaction curves for both cases coincide because the loading processes are the same; however, the rebounding curve of Case B is much longer than that of Case A. Perhaps most importantly, ground deformation behavior is significantly impacted by the underground reservoir permeability boundary conditions (Zhou and Burbey, 2014) and the distribution of embedded inhomogeneous inclusions. This can be clearly visualized by Fig. 8. The region above the interbed has deviated significantly from the curvature without the interbed. Overall the ground deformation in Case B is more significant than that of Case A for all the time intervals. The ground deformation behavior can be characterized by three time intervals. During early injection time (b2e5 s or 55 h), smaller uplift occurs above the region over the interbed and differential uplift between the interbed and aquifer reaches their maximum by the end of this time interval; during intermediate time (from 55 h to day 115) this same region begins to uplift; and during the late time interval (day 115–3.85 years) the rate of ground deformation slows and equilibrium is reached. In the end, the region above the interbed has experienced considerably more uplift than the adjacent areas in this baseline analysis. The regions above the edges of the interbed experience the greatest differential uplift, which may pose potential problems to infrastructures on the surface. The final ground deformation is determined by both the aquifer and interbed properties. Consequently, we constrain our interest to the interbed only. We note that the behavior of a clayey formation is very sensitive to its clay content and the type of clay. Some types of clay such as montmorillonite are much more expansive than others, such as illite.

Fig. 8. Ground uplift during fluid injection for Case A (left) and Case B (right). A curve of an earlier time always appears at a lower level than a later time. The last two curves at time 1.15e7 s and 2.15e7 s for Case A cluster together.

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

The immediate economic loss due to the expansive clay problem exceeds 1 billion US dollars annually in China (Shi et al., 2002), and this number could also be very high for the USA. There is evidence that shallow buried clayey formations are generally more expansive than those that are deeply buried (Weaver, 1989), but these scenarios may be very complex and site-specific. The ground deformation over long time periods is mainly determined by the logarithm elastic bulk modulus of the clayey interbed (“m” in Figure 2). The slope of the unloading line as measured on shale from the Gulf of Mexico is one or two orders of magnitude less than that of the consolidation line; however, the uncertainty of the slope is obviously quite high (Chang and Zoback, 2010). Fig. 9 shows a sensitivity analysis performed on the logarithm of the elastic bulk modulus. One can see that over long injection times the magnitude of ground deformation is mainly governed by this parameter and the region above the clayey interbed may either be much higher or even still much lower than its neighboring regions depending on whether the formation is (very) expansive or not. A more detailed study on the formation behavior correlated with clay type and content is beyond the scope of this paper, and further research is needed. 6. Conclusions As coupled processes, fluid flow and deformation of different formations form a set of separate domains, which cannot be analyzed separately from each other (Gutierrez and Lewis, 2002). For a clayey layer several meters to several tens of meters thick the time delay with respect to the surrounding aquifer can be years, decades or even centuries (Ireland et al., 1984). We simulate a laterally extensive lens-shaped clayey interbed embedded in a sandstone aquifer by using a Cam clay model in conjunction with a poroelastic model. We consider two pore pressure boundary conditions, (1) an “infinitely” extended large aquifer and (2) a small semi-confined reservoir. Spatial and temporal changes associated with this interbed during fluid injection are captured. The evolution of pore pressure distribution in both the aquifer and interbed can be classified into three time intervals, as follows: 1. During the early time period (up to several days in this analysis), the pore pressure increases quickly in the aquifer but only minimally changes near the margins of the interbed. The surface deformation signal results in a lower region above the interbed while its neighboring regions with no interbed all begin to uplift after injection commences. The pore pressure changes in the interbed at this time are

Fig. 9. Ground deformation curves of different values of swelling slope m (or logarithm elastic bulk modulus) for Case A during fluid injection after a period of 3.85 years.

191

due to the mechanical impacts of the adjacent aquifer layers. Their expansion compacts the interbed as evidenced by the porosity reduction of the clays. 2. During intermediate time (from several days to several months in this analysis) the pore pressure in the aquifer has reached equilibrium with the boundary conditions and the rate of pore pressure change in the interbed exceeds that of the aquifer. Uplift of the region above the interbed begins to occur and the greatest differential changes in deformation occur near the edges of the interbed. The result may be larger or smaller uplift depending on the poromechanical properties of the aquifer and interbed. 3. During late time (from several months to several years in this analysis) the pore pressure in the aquifer and interbed equilibrates. The resulting ground deformation signature will be determined by the swelling index of the interbed and poroelastic property of the aquifer. The differential uplift is governed by the poroelastic property of the aquifer and the swelling potential of the interbed. Regions above the clayey interbed will deform in different time scales and an interbed with a higher content of expansive clays could cause significant differential uplift. The region above the flank of a clayey interbed shows varied tilting behavior during fluid injection, which may be the most problematic region for any infrastructures on the surface. On the other hand, delayed uplift may also indicate the existence of a clayey interbed, and the magnitude of uplift can yield critical information about the rebounding (stress–strain) slope (or the logarithm elastic bulk modulus) of the formation from this study. References Alonso, E.E., Gens, S.A., Josa, A., 1990. A constitutive model for partially saturated soils. Geotechnique 40, 405–430. Bell, F.G., Cripps, J.C., Culshaw, M.G., 1986. A review of the engineering behaviour of soils and rocks with respect to groundwater. Geochem. Soc. Lond. Eng. Geol. Spec. Publ. 3, 1–23. Bell, J.W., Amelung, F., Ramelli, A.R., Blewitt, G., 2002. Land subsidence in Las Vegas, Nevada, 1935–2000: new geodetic data show evolution, revised spatial patterns, and reduced rates. Environ. Eng. Geosci. VIII (3), 155–174. Benson, P.M., Meredith, P.G., Platzman, E.S., White, R.E., 2005. Pore fabric shape anisotropy in porous sandstone and its relation to elastic and permeability anisotropy under isostatic pressure. Int. J. Rock Mech. Min. Sci. 42, 890–899. http://dx.doi.org/10.1016/ j.ijrmms.2005.05.003. Berryman, J.G., Milton, G.W., 1991. Exact results for generalized Gassmann's equations in composite porous media with two constituents. Geophysics 56, 1950–1960. Bieniawski, Z.T., 1973. Engineering classification of jointed rock masses. Trans. S. Afr. Inst. Civ. Eng. 15, 335–344. British Standard Institution, B.S.I., 1981. Code of Practice for Site Investigation, BS 5930. HMSO, London. Burbey, T.J., 2006. Three-dimensional deformation and strain induced by municipal pumping, Part 2: numerical analysis. J. Hydrol. 330 (3–4), 422–434. Chang, C., Zoback, M.D., 2010. Viscous creep in room-dried unconsolidated Gulf of Mexico shale (II): development of a viscoplasticity model. J. Pet. Sci. Eng. 72, 50–55. Chang, K.W., Hesse, M.A., Nicot, J.-P., 2013. Reduction of lateral pressure propagation due to dissipation into ambient mudrocks during geological carbon dioxide storage. Water Resour. Res. 49, 2573–2588. Coussy, O., 2004. Poromechanics. John Wiley & Sons Ltd., West Sussex, England 0-47084920-7. Dewhurst, D.N., Aplin, A.C., Sarda, J.-P., Yang, Y., 1998. Compaction-driven evolution of porosity and permeability in natural mudstones: an experimental study. J. Geophys. Res. 103 (B1), 651–661. Engelder, T., Lash, G.G., Uzcategui, R., 2009. Joint sets that enhance production from Middle and Upper Devonian gas shales of the Appalachian Basin. Am. Assoc. Pet. Geol. Bull. 93, 857–889. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. A Math. Phys. Sci. 241 (1226), 376–396. Eshkalak, M.O., Mohaghegh, S.D., Esmaili, S., 2013. Synthetic, Geomechanical Logs for Marcellus Shale. SPE-163690-MS. SPE Digital Energy Conference and Exhibition held in The Woodlands, Texas, USA, 5–7 March 2013. Fidler, L.J., 2011. Natural fracture characterization of the New Albany Shale, Illinois Basin, United States(MS thesis) University of Texas at Austin. Fielding, E.J., Blom, R.G., Goldstein, R.M., 1998. Rapid subsidence over oil fields measured by SAR interferometry. Geophys. Res. Lett. 25 (17), 3215–3218. Gale, J., Reed, R.M., Holder, J., 2007. Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments. Am. Assoc. Pet. Geol. Bull. 91 (4), 603–622. Gambolati, G., 1974. Second-order theory of flow in three-dimensional deforming media. Water Resour. Res. 10, 1217–1228.

192

X. Zhou, T.J. Burbey / Engineering Geology 183 (2014) 185–192

Gambolati, G., Freeze, R.A., 1973. Mathematical simulation of the subsidence of Venice. Water Resour. Res. 9, 721–733. Gasparre, A., Nishimura, S., Minh, N.A., Coop, M.R., Jardine, R.J., 2007. The stiffness of natural London Clay. Geotechnique 57 (1), 33–47. Geertsma, J., 1966. Problems of rock mechanics in petroleum production engineering. Proc. Congr. Int. Soc. Rock Mech. 1st, 585–594. Geertsma, J., 1973. Land subsidence above compacting oil and gas reservoirs. J. Pet. Technol. 25, 734–744. Giraud, A., Rousset, G., 1996. Time-dependent behavior of deep clays. Eng. Geol. 41, 181–195. Gutierrez, M.S., Lewis, R.W., 2002. Coupling of fluid flow and deformation in underground formations. J. Eng. Mech. 128, 779–787. Hawkins, A.B., Pinches, G.M., 1992. Engineering description of mudrocks. Q. J. Eng. Geol. 25, 17–30. Helm, D.C., 1994. Horizontal aquifer movement in a Theis–Thiem confined system. Water Resour. Res. 30, 953–964. Helwany, S., 2007. Applied Soil Mechanics. John Wiley & Sons Inc., New Jersey, USA. Hibbitt, Karlsson, Sorensen, Inc, 1998. ABAQUS, Standard Users Manual. Hibbitt, Karlsson & Sorensen Inc., Pawtucket, Rhode Island. Hight, D.W., Gasparre, A., Nishimura, S., Minh, N.A., Jardine, R.J., Coop, M.R., 2007. Characteristics of the London Clay from the Terminal 5 site at Heathrow Airport. Geotechnique 57 (1), 3–18. Hildenbrand, A., Krooss, B.M., 2003. CO2 migration processes in argillaceous rocks: pressure-driven volume flow and diffusion. J. Geochem. Explor. 78–79, 169–172. Hoffmann, J., Galloway, D.L., Zebker, H.A., Amelung, F., 2001. Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry. Water Resour. Res. 37, 1551–1566. Horsrud, P., Sonstebo, E.F., Boe, R., 1998. Mechanical and petrophysical properties of North Sea Shales. Int. J. Rock Mech. Min. Sci. 35 (8), 1009–1020. International Association of Engineering Geology IAEG, 1979. Report of the commission on engineering geological mapping. Bull. IAEG 19, 364–371. International Society for Rock Mechanics ISRM, 1982. Suggested Methods: Rock Characterization, Testing and Monitoring. In: Brown, E.T. (Ed.), (Oxford). Ireland, R.L., Poland, J.F., Riley, F.S., 1984. Land subsidence in the San Joaquin Valley, California as of 1980. US Geol. Surv. Prof. Pap. 437-I, 93. Jorgensen, D.G., 1980. Relationships between basic solid-engineering equations and basic groundwater flow equations. USGS Water Supply paper 2064. USGS, Washington, USA. Kwon, O., Kronenberg, A.K., Gamgi, A.F., Johnson, B., Herbert, B.E., 2004. Permeability of illite-bearing shale: 1. Anisotropy and effects of clay content and loading. J. Geophys. Res. 109, B10206. http://dx.doi.org/10.1029/2004JB003055. Leake, S.A., 1990. Interbed storage changes and compaction in models of regional groundwater flow. Water Resour. Res. 26 (9), 1939–1950. Li, C., Tang, X., Ma, T., 2006. Land subsidence caused by groundwater exploitation in the Hangzhou-Jiaxing-Huzhou Plain, China. J. Hydrogeol. 14, 1652–1665. Loucks, R.G., Ruppel, S.C., 2007. Mississippian Barnett Shale: lithofacies and depositional setting of a deep-water shale-gas succession in the Forth Worth Basin, Texas. AAPG Bull. 91 (4), 579–601. Louis, L., David, C., Metz, V., Robion, P., Menendz, B., Kissel, C., 2005. Microstructural control on the anisotropy of elastic and transport properties in undeformed sandstones. Int. J. Rock Mech. Min. Sci. 42 (7–8), 911–923. Mura, T., 1982. Micromechanics of Defects in Solids. Martinus Nijhoff, The Hague (494 pp.). Neuzil, C.E., 1994. How permeable are clays and shales? Water Resour. Res. 30 (2), 145–150. Neuzil, C.E., Pollock, D.W., 1983. Erosional unloading and fluid pressure in hydraulically tight rocks. J. Geol. 91 (2), 179–193. Nova, R., 2010. Soil Mechanics. John Wiley & Sons, Inc., Hoboken, New Jersey, USA. Rice, J.R., Cleary, M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous-media with compressible constituents. Rev. Geophys. 14 (2), 227–241. Roscoe, K.H., Schofield, A.N., Wroth, C.P., 1958. On the yielding of soils. Geotechnique 8, 22–53.

Rudnicki, J.W., 1986. Fluid mass sources and point forces in linear elastic diffusive solids. Mech. Mater. 5, 383–393. Safai, N.M., Pinder, G., 1980. Vertical and horizontal land deformation due to fluid withdrawal. Int. J. Numer. Anal. Methods Geomech. 4, 131–142. Sarker, R., Batzle, M., 2010. Anisotropic elastic moduli of the Mancos B shale—an experimental study. SEG Annual Meeting, October 17–22, 2010, Denver, Colorado. Segall, P., Fitzgerald, S.D., 1998. A note on induced stress changes in hydrocarbon and geothermal reservoirs. Tectonophysics 289 (1–3), 117–128. Shearer, T.R., 1998. A numerical model to calculate land subsidence, applied at Hangu in China. Eng. Geol. 49, 85–93. Shi, B., Jiang, H., Liu, Z., Fang, H.-Y., 2002. Engineering geological characteristics of expansive soils in China. Eng. Geol. 67, 63–71. Shi, X.Q., Xue, Y.Q., Ye, S.J., Wu, J.C., Zhang, Y., Yu, J., 2007. Characterization of land subsidence induced by groundwater withdrawals in Su-Xi-Chang area, China. Environ. Geol. 52, 27–40. Teatini, P., Gambolati, G., Ferronato, M., Settari, A., Walters, D., 2011a. Land uplift due to subsurface fluid injection. J. Geodyn. 51, 1–16. Teatini, P., Castelletto, N., Ferronato, M., Gambolati, G., Janna, C., Cairo, E., Marzorati, D., Colombo, D., Ferretti, A., Bagliani, A., Bottazzi, F., 2011b. Geomechanical response to seasonal gas storage in depleted reservoirs: a case study in the Po River basin, Italy. J. Geophys. Res. 116, F02002. http://dx.doi.org/10.1029/2010JF001793. Tolman, C.E., Poland, J.F., 1940. Groundwater infiltration, and ground surface recession in Santa Clara Valley, Santa Clara County, California. EOS Trans. AGU 21, 23–34. Vasco, D.W., 2004. Estimation of flow properties using surface deformation and head data: a trajectory-based approach. Water Resour. Res. 40, W10104. http://dx.doi. org/10.1029/2004WR003272. Vasco, D.W., Karasaki, K., Kishida, K., 2001. A coupled inversion of pressure and surface displacement. Water Resour. Res. 37 (12). http://dx.doi.org/10.1029/2001WR000391 [issn: 0043-1397]. Vasco, D.W., Rucci, A., Ferretti, A., Novali, F., Bissell, R.C., Ringrose, P.S., Mathieson, A.S., Wright, I.W., 2010. Satellite-based measurements of surface deformation reveal fluid flow associated with the geological storage of carbon dioxide. Geophys. Res. Lett. 37, L03303. http://dx.doi.org/10.1029/2009GL041544. Vermylen, J.P., 2011. Geomechanical studies of the Barnett shale, Texas, USA(Ph.D dissertation) Stanford University. Vukadin, V., 2007. Modeling of the stress–strain behavior in hard soils and soft rocks. Acta Geotech. Slov. 2, 5–15. Wang, H.F., 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press. Weaver, C.E., 1989. Clays, Muds, and Shales. Developments in Sedimentology. 44. Elsevier Sciences Publishers, B.V. New York, USA. Yang, Y., Aplin, A.C., 2010. A permeability–porosity relationship for mudstones. Mar. Pet. Geol. 27 (8), 1692–1697. Yerkes, R.F., Castle, R.O., 1969. Surface deformation associated with oil and gas field operation in the United States. In: Tison, L.J. (Ed.), Land SubsidenceIASH Publ. 88 vol. 1. Int. Assoc. of Sci. Hydrol., Louvain, Belgium, pp. 55–64. Zhang, L., 2005. Engineering Properties of Rocks. Elsevier Ltd., Oxford, UK. Zhou, X.J., Burbey, T.J., 2012. A FEM approach to measure Skempton's B coefficient for supercritical CO2-saturated rock. Environ. Eng. Geosci. 18 (4), 343–355. Zhou, X.J., Burbey, T.J., 2014. How horizontal surface deformation during fluid injection correlates to reservoir permeability setting. Environ. Eng. Geosci. 20 (3), 305–320. Zhou, Q., Birkholzer, J.T., Mehnert, E., Lin, Y.F., Zhang, K., 2009. Modeling basin‐and plume‐ scale processes of CO2 storage for full‐scale deployment. Ground Water 48 (4), 494–514. Zhou, X.J., Zeng, Z., Liu, H., 2010. Laboratory testing on Pierre shale for CO2 sequestration under clayey cap rock, paper ARMA10-107. Proc. 44th U.S. and 5th U.S.–Canada Rock Mech Symposium, Salt Lake City, UT, USA. June 27–30 (12 pp.). Zoback, M.D., 2007. Reservoir Geomechanics. Cambridge University Press, NY. Zoback, M.D., Byerlee, J.D., 1976. Effect of high-pressure deformation on the permeability of Ottawa Sand. AAPG Bull. 60, 1531.