Boletín de la Sociedad Geológica Mexicana Volumen 65, núm. 1, 2013, p. 123-136 http://dx.doi.org/10.18268/BSGM2013v65n1a10

Analysis and simulation of regional subsidence accompanying groundwater abstraction and compaction of susceptible aquifer systems in the USA

Devin L. Galloway^1,*, Michelle Sneed¹

¹ U.S. Geological Survey, Sacramento, CA, USA.

^* This email address is being protected from spambots. You need JavaScript enabled to view it.

Abstract

Regional aquifer-system compaction and land subsidence accompanying groundwater abstraction in susceptible aquifer systems in the USA is a challenge for managing groundwater resources and mitigating associated hazards. Developments in the assessment of regional subsidence provide more information to constrain analyses and simulation of aquifer-system compaction. Current popular approaches to simulating vertical aquifer-system deformation (compaction), such as those embodied in the aquitard drainage model and the MODFLOW subsidence packages, have proven useful from the perspective of regional groundwater resources assessment. However, these approaches inadequately address related local-scale hazards—ground ruptures and damages to engineered structures on the land surface arising from tensional stresses and strains accompanying groundwater abstraction. This paper presents a brief overview of the general approaches taken by the U.S. Geological Survey toward understanding aquifer-system compaction and subsidence with regard to a) identifying the affected aquifer systems; b) making regional assessments; c) analyzing the governing processes; and d) simulating historical and future groundwater flow and subsidence conditions. Limitations and shortcomings of these approaches, as well as future challenges also are discussed.

Keywords: Subsidence, aquifer-system compaction, aquitard drainage, MODFLOW, poroelastic deformation.

Resumen

Es un gran reto en los Estados Unidos administrar los recursos de agua subterránea evitando los daños que pueden ocasionarse por subsidencia regional provocada por la explotación de los acuíferos. Los avances en la evaluación de la subsidencia regional han dado mayor información para llevar a cabo detalladamente los análisis y simulaciones de la compactación o enjutamiento de acuíferos. Los modelos que actualmente se emplean para simular la deformación vertical del acuífero (compactación) provocada por la extracción de agua subterránea, como los basados en el modelo de acuitardo drenado y los paquetes computacionales de subsidencia, como el programa MODFLOW, han demostrado ser útiles para la evaluación de los recursos de agua subterránea. Sin embargo, estos enfoques tratan de manera inadecuada los riesgos relacionados con fracturas locales en el subsuelo, lo que genera daños a las estructuras ubicadas en superficie por la generación de esfuerzos de tensión y deformaciones asociadas a la extracción de agua. Este artículo presenta una reseña breve de los enfoques generales adoptados por el Servicio Geológico de los EE.UU. hacia el entendimiento de la compactación de los acuíferos y el hundimiento del subsuelo asociado con respecto a: a) la identificación de los sistemas de acuíferos afectados; b) la realización de evaluaciones regionales; c) el análisis de los procesos que gobiernan el comportamiento; y d ) la simulación del flujo de agua subterránea en el pasado y futuro, así como las condiciones de subsidencia. Asimismo, se discuten las limitaciones y deficiencias de estos enfoques y los retos en el futuro.

Palabras Clave: Subsidencia, compactación del acuifero, acuitardo drenado, MODFLOW, deformación poro-elástica.

1. Introduction

The abstraction of groundwater for agricultural and municipal-industrial water supplies is a leading cause of subsidence in the United States, and accounts for the majority of anthropogenic subsidence in the nation (Galloway et al., 1999). Much of the groundwater related subsidence is manifest in the semi-arid West where groundwater is in great demand and withdrawals since the turn of the 20^th century have critically lowered groundwater levels, raised effective stresses and compacted susceptible aquifer systems.

Though the phenomenon of compacting subsurface reservoirs owing to lowered fluid pressures has been recognized since the early 20^th century in oil fields in Texas and California, the study of cause-and-effect relations in regional aquifer-system compaction and subsidence were not begun until later. The approaches developed in pioneering studies principally in California during 1960s–80s invoked the Terzaghi (1925) principle of effective stress and focused on characterizing in situhydro-mechanical behaviors (vertical deformation) of aquitards. The approaches were decidedly hydrogeological and regional in scale with respect to groundwater flow and subsidence (Poland and Davis, 1969; Poland, 1984).

The early California studies were done by the United States Geological Survey (USGS) in the context of hazard assessment with respect to the design, development and protection of water conveyance infrastructures (such as the California Aqueduct) in the San Joaquin Valley (Figure 1; e.g., see Ireland et al., 1984), and coastal flood protection measures in the Santa Clara Valley (Poland and Ireland, 1988; Ingebritsen and Jones, 1999). Though mitigation of subsidence hazards remains a principal motivation for subsidence research, more often, the roles of aquifer-system compaction and the accompanying subsidence are being considered in regional groundwater resource assessments. The problems of groundwater depletion in response to groundwater abstractions include consideration of aquifer-system compaction in many developed groundwater basins. Managing groundwater resources in sustainable ways means balancing the consequences of subsidence hazards against the benefits of producing water supplies.

In this context, USGS subsidence studies have continued to focus on regional scale groundwater flow and storage depletion, arguably with more emphasis on hydrogeological than geomechanical engineering approaches. This paper presents a brief overview of the general approaches taken by the USGS toward understanding aquifer-system compaction and subsidence in the USA with regard to a) identifying the affected aquifer systems, b) making regional assessments, c) analyzing the governing processes, and d) simulating the historical and future groundwater flow and subsidence conditions. Limitations and shortcomings of these approaches, as well as future challenges also are discussed.

Figure 1. Approximate location of maximum measured subsidence (9 m) in the San Joaquin Valley, California (1925–77) attributed to aquifer-system compaction caused by groundwater abstraction. Signs on pole (1925, 1955) are positioned at approximate former elevations of land surface. Pictured is Dr. J.F. Poland; photograph by R.L. Ireland, USGS, ca. 1977.

2. Affected Aquifer Systems

Withdrawal of subsurface fluids from clastic sediments has permanently lowered the elevation of about 26000 km²of land in the 48 conterminous United States —an area of similar extent to the State of Massachusetts (Holzer and Galloway, 2005). Most of the subsidence attributed to subsurface fluid withdrawal is caused by groundwater abstraction. Permanent subsidence can occur when groundwater is removed by pumpage or drainage. The reduction of fluid pressure in the pores and cracks of aquifer systems, especially in unconsolidated clastic rocks, is inevitably accompanied by some deformation of the aquifer system. Because the granular structure —the “skeleton”— of the fluid-bearing rocks is more or less compliant, a shift in the balance of support for the overlying material causes the skeleton to deform slightly. Both the aquifers and aquitards that constitute the aquifer system undergo deformation, but to different degrees. Almost all the permanent subsidence in aquifer systems is attributable to the compaction of aquitards during the typically slow process of aquitard drainage (Tolman and Poland, 1940).

Figure 2 shows 51 selected areas of known subsidence in the USA and the associated principal aquifer systems where subsidence is attributed primarily or secondarily to groundwater or oil and gas abstractions. Subsidence in 47 of the areas is attributed primarily to groundwater abstraction. All but 4 of the 47 subsidence areas produce groundwater from one of six principal aquifer systems —the Basin and Range basin-fill aquifers, the California Coastal Basin aquifers, the Coastal lowlands aquifer system, the Central Valley (California) aquifer system, the Northern Atlantic Coastal Plain aquifer system, and the Rio Grande aquifer system.

Each of the affected principal aquifer systems comprises a large (> 100 m) thickness of Quaternary age unconsolidated deposits with a significant fraction of clay minerals distributed in variable-thickness (< 5 m), discontinuous interbedded layers. Most have a laterally extensive thick (> 10 m) confining unit generally comprising lacustrine deposits. In the western aquifer systems groundwater abstraction for irrigation has been the principal historical water use. Where supplies were available, surface water has been imported to mitigate effects of groundwater depletion. Where urbanization has succeeded agricultural land use, groundwater demand commonly has decreased initially and subsequently has increased steadily with population growth.

Subsidence due to subsurface fluid extraction in coastal regions including those near Houston, Texas, and New Orleans, Louisiana (Coastal lowlands aquifer system), Long Beach and Santa Clara Valley (California Coastal Basin aquifers), and other areas in the Northern Atlantic Coastal Plain aquifer system plays a role in the relative rise of local mean sea level (LMSL). LMSL is affected by land movements attributed to natural subsidence, and sea-level changes attributed to eustasy. The observed rate of eustatic rise is 1 – 2 mm/yr (IPCC, 2001), of which global climate change likely is a significant contributor. Estimated average subsidence rates vary widely and are confounded by estimates of the relative (subsidence and uplift) contributions from isostasy (González and Törnqvist, 2006). Locally, during the latter part of the 20^thcentury much of the isostatic decline can be attributed to subsidence caused by the compaction of sediments owing to the abstraction of subsurface fluids; sediment compaction rates greater than 100 mm/yr have been measured in coastal basins while compaction rates attributed to other processes contributing to relative sea-level rise were nearly 1–2 orders of magnitude less.

Subsidence attributed to aquifer-system compaction in the USA generally is largest in magnitude and distribution in the western aquifer systems and in the Houston area of the Coastal lowlands aquifer system. The susceptibility of these aquifer systems to aquifer-system compaction is enhanced because of the large groundwater abstractions and the prevalence of fine-grained (clays and silts) deposits. By virtue of their larger compressibilities under virgin (inelastic) loading conditions, compaction of the fine-grained sediments is responsible for most of the subsidence.

Within each principal aquifer system, several to numerous named aquifer systems have been identified and each may comprise multiple individual groundwater subbasins for which various degrees of groundwater-resources assessments have been made. Most of these assessments have focused primarily on the quantity and quality of available groundwater with regard to the demand for groundwater, and secondarily on the expected consequences of developing the resource—subsidence. As such, the basin and (or) aquifer-system hydrogeologic boundaries have been used to define regional accounting units for computing groundwater budgets and formulating other groundwater flow and subsidence assessments.

Because subsidence accompanying aquifer-system compaction is the result of a transient coupled hydro-mechanical process, the assessment of regional subsidence in these affected basins has followed the assessment of regional groundwater flow. Thus, the USGS approach generally has been first to characterize the steady-state, predevelopment regional groundwater flow system and second to evaluate the transient, post-development hydro-mechanical stresses (abstractions and resulting fluid-pressure and effective stress changes) and the system responses.

3. Regional assessments

Regional assessments include geologic, geophysical, hydrogeologic, and geodetic reconnaissance and measurement components, and typically some hydrogeologic and geodetic monitoring components. The available data and information are assimilated, and new data needs are identified. Ideally, as feasible, new data are collected. The broad goals are to provide a sound basis to formulate a conceptual model of groundwater flow and aquifer-system compaction, and to provide sufficient spatial and temporal detail from which to base estimates of key hydrogeologic properties governing flow and subsidence. An emphasis is placed on the integration of data from different sources to develop the concepts, property characteristics, and observational data to constrain analyses and simulations.

3.1. Geology and Geophysics

Generally the goal is to produce a geologic framework model from which a hydrogeologic framework model is developed. The geologic framework model is an interpretive 3D representation of the subsurface geology of the study area. Structural, stratigraphic and lithologic information assembled from available geologic maps and cross-sections, borehole lithology logs, and surface (including electrical, seismic, magnetic and gravity) and borehole (including electrical, gamma and acoustic) geophysical surveys. An emphasis is placed on defining the distribution of coarse- and fine-grained younger Cenozoic deposits (Holocene and Pleistocene), the contacts between the younger and older alluvium (Pliocene) and the consolidated rocks that may constitute the basement of the groundwater flow system.

3.2. Hydrogeology

The hydrogeologic characterization generally begins with a definition of the surface-water drainage basin boundaries and acquisition of digital elevation model data. The source areas of natural recharge and discharge are identified, as well as any surface-water bodies that may interact with the groundwater flow system. The available groundwater-level information is plotted and a preliminary potentiometric surface map is used to infer general groundwater flow paths. The boundaries of the groundwater flow system are defined using the available geologic framework model. The hydrostratigraphy is defined, and typically involves combining individual geologic units into hydrogeologic units —constituting water-bearing aquifers and any intervening confining units. The available groundwater quality data are used to constrain the definition of individual aquifers in the system as well as the generalized flow paths. Thus, a preliminary conceptual hydrogeologic framework model is composed.

Available aquifer-system hydraulic property data are used to refine and further constrain the model. Generally, hydraulic property data is scarce for confining units and the distributed thin aquitards (interbeds) adjacent to and within the aquifers. Obtaining this information is one goal of new data collection —drilling, testing, monitoring. New data collection typically focuses on obtaining information in areas where data gaps exist. Sometimes, a preliminary numerical model of the flow system is used to identify sensitive, critical areas where additional data would improve constraints on the hydrogeologic framework model.

Hydrologic monitoring sites are identified taking advantage of existing monitoring installations (stream gages on critical streams; precipitation gages; weather stations; routinely measured groundwater monitoring networks in the basin) where possible. For groundwater-level monitoring the goals are a) to establish the generalized static relations between heads throughout the basin and between individual aquifers within the aquifer system; and b) to define the seasonal and interannual variability in heads in response to climatic and anthropogenic stresses.

Generally, based on the available information, two hydrologic budgets are formulated: one for the pre-development, steady-state flow system; and one for the post-development transient flow system. The steady-state budget attempts to balance natural recharge and discharge fluxes given the distribution of historical, pre-development heads. The transient budget accounts for alterations in the balance owing to groundwater abstractions, changes in land use and any alterations to natural recharge and discharge sources. These preliminary hydrologic budgets typically form the basis for distributing specified fluxes in groundwater flow models (discussed below in section 5).

Figure 2. Selected, known areas of land subsidence, owing primarily or secondarily to groundwater or oil and gas abstractions in the 48 conterminous USA, and associated aquifer systems (modified from Galloway et al., 2008).

3.3. Geodesy

A regional geodetic assessment is needed to define the temporal and spatial extent of land subsidence in the study area and, where possible, to define relations between aquifer-system compaction, subsidence and any governing hydrogeologic factors. Where historical subsidence maps are unavailable, ancillary or anecdotal information that suggests subsidence may be occurring often is useful and pertinent to regional-scale subsidence processes especially where the subsidence is subtle. Some types of ancillary information that have proven useful in identifying susceptible groundwater basins include increased incidence of damaged wells—protruding and (or) collapsed well casings; repeat adjustments to local geodetic controls; riverine or coastal flooding; conveyance and drainage problems; and ground failures (Galloway et al., 2008, p. 60–67).

Subsidence reconnaissance typically includes an examination of the available historic geodetic surveys in the study area. Special attention is paid to repeat adjustments of vertical control that may signal changes in the geodetic control networks that mask evidence of subsidence. Interferometric synthetic aperture radar (InSAR) has proven useful as a reconnaissance tool in subsidence related to groundwater abstractions (Galloway and Hoffmann, 2007). Ground-surface displacements since 1992 can be detected in radiometrically favorable landscapes at high spatial resolution (30 m) and good temporal resolution (about every 35 days) with a potentially high-resolution of displacement (± 5 mm). Figure 3 shows a reconnaissance InSAR image showing subsidence attributed to aquifer-system compaction in Antelope Valley (Mojave Desert), California (Galloway et al., 1998).

Figure 3. InSAR-detected subsidence (October 1993 to December 1995) and historical (1930–92) subsidence (Ikehara and Phillips 1994), Antelope Valley, California. (Beige-colored areas signify regions of decorrelation of the radar; black-colored areas signify regions of small-magnitude uplift).

A survey of the existing historical benchmarks combined with any new benchmarks is made using GPS or spirit leveling. Advantage is taken of the ancillary and other reconnaissance information such as InSAR to design the network to optimize coverage of affected areas and to ensure stable geodetic control. The survey results are used to make current subsidence maps, and form the basis for comparison with future repeat surveys of the network (e.g., Sneed et al., 2001, 2002; Sneed and Brandt, 2007). Figure 4 shows an example from the Coachella Valley, California.

Figure 4. Areas of subsidence, 26 October 2003 through 12 June 2005, as shown by subsidence contours interpreted from InSAR, Coachella Valley, California (modified from Sneed and Brandt, 2007, fig. 8b). The three new GPS stations were sited in subsidence affected areas based on previous reconnaissance InSAR.

The new geodetic and subsidence information can be used with a suite of InSAR or persistent scatterer interferograms (PSI) to develop ground-truthed time-series displacement maps for the period 1992–present for study areas where suitable images and favorable radiometric conditions exist (Galloway and Hoffmann, 2007).

For more precise, accurate and direct measurements of aquifer-system compaction at local scales, borehole extensometers are used. Borehole extensometers measure the change in vertical distance in the interval between the land surface and a reference point or “subsurface benchmark” at the bottom of a deep borehole (Riley, 1986). If the subsurface benchmark is established below the base of the compacting aquifer system the extensometer can be used as the stable reference or starting point for local geodetic surveys. The deformation history generated by a borehole extensometer provides the basis for stress-strain analysis and modeling that is used to estimate and constrain the material properties governing compaction of the aquifer system (see the Analyses and Simulation sections, 4 and 5).

Several types of early borehole extensometer designs are presented in Poland (1984). The anchored cable and pipe (free standing and counterweighted) extensometers have been used widely in a number of successful subsidence investigations. More recently, dual-stage counterweighted pipe extensometers (Figure 5) have been used to measure compaction simultaneously in two depth intervals (Heywood, 1993; Metzger et al., 2002). Counterweighted pipe extensometers are capable of measurement resolutions of 0.01 - 0.1 mm (Riley, 1969).

Figure 5. Schematic of dual-stage counterweighted borehole extensometer (modified from Galloway et al., 2008, fig. 2.16A).

Multiple position borehole extensometers (MPBXs) that incorporate markers anchored to the formation borehole have been used effectively to monitor subsidence caused by groundwater abstraction. Magnetic markers have been used in the Republic of China (ROC) (Liu et al., 2004; Hwang et al., 2008) to compute vertical displacements in boreholes using repeat borehole logging with magnetic sensors on calibrated lines or tapes to measure temporal changes in marker positions. This method is capable of monitoring ten to several tens of marker positions in a single borehole at measurement resolutions of about 1 - 2 mm over depths of several hundred meters.

Horizontal displacement occurs in aquifer systems in response to pumping and seasonal recharge/discharge stresses (Wolff, 1970; Carpenter, 1993; Helm, 1994; Hsieh, 1996; Bawden et al., 2001; Burbey, 2001a, 2001b; Li, 2007a, 2007b, 2007c). However, historically, in areas of subsidence attributed to fluid abstraction, horizontal displacement of the land surface has been measured at only a few places.

Several kinds of horizontal extensometers are used to measure horizontal ground motion at earth fissures caused by changes in groundwater levels (Carpenter, 1993). Buried horizontal extensometers constructed of quartz tubes or invar wires are useful when precise, continuous measurements are required on a scale of 3 - 30 m. Tape extensometers measure changes across inter-monument distances up to 30 m with repeatability of approximately 0.3 mm.

GPS has the ability to measure both horizontal and vertical movements. Uplift and subsidence associated with managed subsurface fluid production (injection and extraction) is accompanied by measurable horizontal movements in the Earth’s crust. If the points are directly on the margin of the subsidence/uplift feature, then the ratio of vertical to horizontal motion may be nearly 1:1 (Bawden et al., 2001). For example, in the San Gabriel Valley in southern California, groundwater pumping pulled nearby continuous GPS (CGPS) stations inward toward the region of maximum drawdown (Figure 6), and record rainfall in the same region in the winter-spring 2005 produced more than 4 cm of uplift with greater than 1 cm of radial outward motion of the nearby CGPS station (King et al., 2007).

Figure 6. Horizontal GPS displacement vectors for the Upper San Gabriel Valley, California. Groundwater levels in the valley declined about 3.5 m between May and October 1999 corresponding with approximately 12 mm of land subsidence (Bawden, 2002). The neighboring continuous GPS stations “vyas” and “lphs” are pulled inward towards the zone of maximum subsidence which generally corresponds with the drawdown cone.

Other promising future applications of InSAR measurements may prove useful in evaluating horizontal deformation in aquifer systems (e.g., Burbey, 2001a, 2001b, 2002, 2005; Hoffmann and Zebker, 2003).

Ground-based tripod light detection and ranging (T-LiDAR) is a portable remote sensing instrument that uses an infrared laser to scan the landscape and generate very detailed (centimeter to sub-centimeter) and accurate 3D digital models of the scanned target at distances from 2 to 2000 m. The 3D positional accuracies and laser spot spacing are a function of laser beam divergence and angular step-width, where individual point position accuracies of ± 4 mm at 100 m are common. A full 3D image is obtained by scanning a target from multiple directions to characterize all sides of the target area and to minimize shadowing. T-LiDAR scans imaged from different vantage points are aligned and combined through an algorithm that computes a best-fit surface through the individual points in each scan and then minimizes the misfit between common surfaces in each scan.

Changes in the position of the land surface or a structural feature can be obtained through differencing of precisely aligned T-LiDAR images collected at different times, known as differential LiDAR (Dif-LiDAR). There are two approaches for Dif-LiDAR analysis: absolute and relative. Absolute Dif-LiDAR measures land-surface change by differencing two LiDAR datasets that are georeferenced using GPS. This approach produces the best information of land-surface change, but requires additional field equipment and time for GPS data collection and processing and may be unnecessary for all scientific applications. Alternatively, relative Dif-LiDAR applies the same best-fit surface-matching algorithm used in the alignment of an individual scan to common “stable” regions outside of the area of interest that is changing. This approach is ideal for resolving very detailed spatial changes within a well-defined deformation zone or imaging change in a subset of a larger dataset where absolute positioning is not required. T-LiDAR is an ideal technique for measuring spatial and temporal changes in regions that are actively deforming, but the technique may be too labor intensive for characterizing subsidence features at scales greater than a few square kilometers.

4. Analyses

Analysis of the geologic, hydrogeologic and geodetic information is best achieved using an integrated approach. The spatial distributions of susceptible (compressible) aquifer-system deposits and the spatial and temporal distributions of groundwater discharge and recharge stresses, groundwater-level variations, aquifer-system compaction and land-surface displacements are used to determine whether regional groundwater abstraction causes aquifer-system compaction and land subsidence. If so, the assessment information is used to formulate a conceptual model of groundwater flow and land subsidence. A simple approach that has been widely applied to regional groundwater flow and subsidence problems is based on the aquitard drainage model, which uses conventional groundwater flow theory and two principles of consolidation to describe the relations between fluid pressure, intergranular stress and fluid flow: 1) the principle of effective stress (Terzaghi, 1925),

(1a)

where

p is pore-fluid pressure,

and are components of the effective stress and total stress tensors of order two, respectively,

i and j for i = 1 to 3 and j= 1 to 3, represent the Cartesian coordinates x, y and z, respectively, and

is the Kronecker delta, where

,

which, for 1D vertical stress problems reduces to

where σ'_zz and σ_zz are the vertical effective and total stresses, respectively (Figure 7); and 2) the theory of hydrodynamic lag, which describes the delay in draining aquitards and is explained later in this section. This concept has formed the theoretical basis of many successful subsidence investigations associated with depressuring porous media (Helm, 1984a; Holzer, 1998).

Figure 7. Hydrogeologic units of an idealized aquifer system and interface stresses between a confined aquifer and an aquitard (modified from Galloway et al., 1999, p. 10).

The standard diffusion equation for 3D transient saturated groundwater flow in a confined aquifer can be expressed by

,

(2)

where

his hydraulic head,

is specific storage,

K is hydraulic conductivity,

W is volumetric flux of sources and (or) sinks of water, and

t is time.

Equation 2 can be generalized for aquifer systems comprising aquifers and aquitards (e.g. Figure 7) where the parameters and  and K represent properties of the aquifers, aquitards or the aquifer system.

The aquifer specific storage defined by Jacob (1940) assumed 1-D vertical deformation only (see Cooper, 1966; Gambolati, 1973) with the storage changes proportional to head through the specific storage S_s (Cooper, 1966):

(3)

or expressed another way (Riley, 1969),

(4)

where

ᾱ is the vertical skeletal compressibility,

β_w is the compressibility of water,

n is porosity,

is the unit weight of water,

is the skeletal specific storage, where

= ρ_wgᾱ and

  is the water specific storage, where

  = ρ_wgnβ_w.

The vertical skeletal compressibility, ᾱ, can be defined based on the vertical effective stress and vertical strain by

,

(5)

where is the change in thickness of a control volume with initial thickness of a deformable geologic unit-compaction is defined as a positive . This is the matrix compressibility used in standard formulations for transient saturated groundwater flow (equation 2). Two skeletal compressibilities can be further defined: 1) ᾱ_e for the elastic range of stress where , (where is the previous maximum effective stress); and 2) ᾱ_v for the virgin or inelastic range of stress where . Thus, two skeletal specific storages can be defined based on equation 4, and , for the elastic and inelastic ranges of stress, respectively.

For effective stress changes below the previous maximum effective stress (preconsolidation stress threshold), the compaction or expansion of both aquitards and aquifers is elastic —that is, approximately proportional to the change in effective stress over a moderate range in stress, and fully recoverable if the stress reverts to the initial condition. For effective stress changes above the preconsolidation stress threshold, the “virgin” compaction of aquitards is chiefly inelastic —that is, not fully recoverable upon decrease in effective stress. This virgin compaction includes a recoverable elastic component that is small when compared to the inelastic component. The inelastic and elastic compaction is roughly proportional to the change in logarithm of effective stress. In contrast to aquitards, the compaction of aquifers is chiefly elastic but it may include an inelastic component.

The aquitard drainage modelis based on the generality that in developed aquifer systems, aquitards provide water and aquifers transmit it to wells —that flow in aquitards principally is vertical owing to the large contrast in hydraulic conductivity between aquitards and aquifers (generally greater than two orders of magnitude). As such, a 1D form of Equation 2 can be used to describe groundwater flow in aquitards:

where

′ denotes material properties of the aquitard, and

K'_z is the aquitard vertical hydraulic conductivity.

For virgin compaction, recognizing that .

Terzaghi (1925) developed an analytical solution for an equivalent of Equation 6 to simulate the equilibration of head in a saturated clay sample with uniform initial head where only vertical flow is permitted, in response to a specified instantaneous step change in head at the top and bottom of the sample. This process describes the theory of hydrodynamic lag(consolidation) which was extended to the analysis of aquitard drainage (Riley, 1969), and subsequently to the simulation of aquitard drainage (Helm, 1975, 1976).

The theory describes the delay in draining aquitards when heads are lowered in adjacent aquifers, as well as any residual compaction of the aquitards that may continue long after the heads are initially lowered. The application of the hydrodynamic consolidation theory of soil mechanics to aquifer-system compaction has been summarized lucidly by Riley (1969, p. 425–426).

For a doubly-draining aquitard (or discontinuous interbed), the same aquifer head changes are assumed to occur at the upper and lower surfaces. The time constant for this type of aquitard is (Riley, 1969; Helm, 1975, p. 470):

,

(7)

where

b' is the thickness of the aquitard,

for the elastic range of stress, where is the elastic specific storage, and for the inelastic range of stress, where is the inelastic specific storage.

A dimensionless time factor, T_D, is defined as

.

(8)

such that for T_D equals unity, and τ' is the time required to attain about 93 percent of the ultimate consolidation. Time constants can be computed for both the elastic (τ'_e ) and inelastic (τ'_v ) stress ranges. Detailed development of the hydrodynamic consolidation theory for soils summarized above for aquitards is given by Scott (1963, p. 162–197).

A key component of the analyses is to determine whether any observable subsidence can be attributed to virgin or inelastic compaction of the aquifer system. Aquifer-system deformation in the elastic range of stress typically is small and reversible. Seasonal, reversible land surface displacements (subsidence and uplift) of a few cm are typical for many alluvial aquifer systems with significant fractions of fine-grained deposits. Generally, for larger magnitude subsidence that accumulates over annual periods in the presence of declining groundwater levels, inelastic compaction (largely irreversible) is suspected. This analysis involves an evaluation of the preconsolidation stresses in the aquifer system. Paired historical groundwater-level and subsidence time series have been used to estimate the groundwater level threshold corresponding to the onset of inelastic compaction in the aquifer system (Holzer, 1981; Sneed and Galloway, 2000). The threshold groundwater level is a measure of the initial, predevelopment preconsolidation stress. Subsequent groundwater level lows accompanying groundwater abstraction represent new values of the preconsolidation stress. Water-level fluctuations above previous lows may cause measureable elastic compaction, though depending on the stress history, this deformation may be superimposed on inelastic compaction occurring in relatively thick fine-grained deposits with sufficiently large time constants.

An analysis of regional aquifer-system compaction may include estimates of elastic and inelastic skeletal storage properties made from observable, coincident vertical land-surface or aquifer-system displacements, and changes in aquifer-system hydraulic heads. The approach follows from the relations for and ᾱ given in Equations 4 and 5 for σ'_zzexpressed in terms of equivalent hydraulic head, and can be written as

where ^*denotes material properties of the aquifer system (aquitards and aquifers).

For the case where the total stress is constant Equation 9 reduces to

.

(10)

For examples of estimates based on aquifer-system compaction time series measured in borehole extensometers see Riley (1969), Hanson (1989), and Sneed and Galloway (2000). For an example based on land-surface displacements measured using InSAR see Hoffmann et al. (2001).

5. Simulation

Methods to simulate compaction in aquifer systems based on concepts of aquitard drainage were developed by Helm (1972, 1975, 1976, 1978), Witherspoon and Freeze (1972), Gambolati and Freeze (1973), Narasimhan and Witherspoon (1977), and Neuman et al. (1982). The numerical formulation of the aquitard drainage model in the 1D (vertical) model COMPAC (Helm, 1984b, 1986) led directly to important simulation tools and powerful predictive techniques for land subsidence caused by water-level fluctuations within a confined aquifer system. COMPAC has been used widely to analyze compaction, estimate critical aquifer-system parameters, and predict future subsidence at borehole extensometer sites for which there are detailed records of water-level changes, compaction and (or) subsidence (e.g. Epstein, 1987; Hanson, 1989; Pope and Burbey, 2003, 2004; Liu and Helm, 2008a, 2008b).

Regional aquifer-system compaction and land subsidence simulation modules were incorporated into earlier versions of the MODFLOW groundwater flow model (McDonald and Harbaugh, 1988; Harbaugh and McDonald, 1996) as Interbed Storage packages 1, 2 and 3—IBS1, IBS2, IBS3 (Leake, 1990, 1991; Leake and Prudic, 1991). The approach implements the concept of interbed storage whereby discontinuous aquitards (termed interbeds) interspersed within aquifers can deform vertically in response to head changes in the adjacent aquifer. The deformation of the interbeds is coupled to the groundwater flow equation through interbed storage. The formulations are based on the aquitard drainage model and individually address different processes (no delayed drainage, delayed drainage, and variable geostatic load and stress dependent hydraulic properties with no delayed drainage). These packages were subsequently upgraded for newer versions of MODFLOW (Harbaugh et al., 2000; Harbaugh, 2005) as the SUB package (Hoffmann et al., 2003b) which combines the functionalities of IBS1 and 2, and the SUB-WT package (Leake and Galloway, 2007) which updates the IBS3 package.

Each of the MODFLOW subsidence packages partially couples subsidence (the sum of aquifer-system compaction) to groundwater flow (Equation 2) through the source term, W, and the skeletal specific storage term (Equations 9 and 10). Each of the formulations simulate constant hydraulic conductivity, i.e. stress-dependent hydraulic conductivity is not simulated in the present packages. For each of the packages 3D flow is simulated in aquifers and can be simulated in confining units. Flow in interbeds is not simulated in SUB-WT and in SUB where no-delayed drainage is specified, 1D vertical flow can be simulated in interbeds using the SUB package for systems of interbeds where delayed drainage is specified. The SUB package simulates constant total stress and constant skeletal specific storage and Equation 10 describes the relation between changes in head and compaction, and skeletal specific storage. The SUB-WT package simulates variable total stress and simulates aquitard skeletal specific storage () as a function of effective stress by solving a nonlinear form of Equation 6 where and thus is a function of head (Leake and Galloway, 2007). The relation between effective stress and skeletal specific storage is described using

where e₀ is the initial void ratio and C is the compression index that takes on values of C_c and C_r for σ'_zz greater than and less than or equal to, respectively, the preconsolidation stress σ'_max. Similarly, for in Equations 10 and 11,

The formulation in Equation 11 is based on the relation where C, and thus , varies proportionally with (Leake and Galloway, 2007).

Delayed drainage is simulated explicitly in the SUB package (delay option) by numerical solution of Equation 6 for head in a system of delay interbeds, and using Equation 10 to compute compaction at discrete intervals in the system of interbeds. If delayed drainage is not simulated (SUB package—no delay option) compaction occurs instantaneously with changes in head in the aquifer and its interbeds. This approach is valid for thin interbeds in which heads can equilibrate readily (or within a model time step) with head changes in the surrounding aquifer. It is worth noting that delayed drainage can be simulated with the SUB package using the no-delay option by specifying multiple vertically adjacent model layers comprising interbeds. Delayed drainage is not simulated in the SUB-WT package

MODFLOW has been used widely by the USGS and others mostly in Arizona, California and Texas to simulate regional groundwater flow, aquifer-system compaction, and land subsidence (e.g. Hanson et al., 1990, 2003, 2004; Hanson and Benedict, 1994; Wilson and Gorelick, 1996; Galloway et al., 1998; Kasmarek and Strom, 2002; Hoffmann et al., 2003a; Leighton and Phillips, 2003; Faunt, 2009; Leake and Galloway, 2010). Figure 8 shows simulated results of an inverse model of subsidence in Antelope Valley, California, using the SUB package. MODFLOW has been used in 1D vertical mode, similar to COMPAC, to simulate compaction measured at borehole extensometer sites for which there are detailed records of water-level changes (e.g. Sneed and Galloway, 2000; Pavelko, 2003; Sneed, 2008).

Figure 8. Simulated subsidence and drawdown, Antelope Valley, California, 1996-99, constrained by historical terrestrial geodetic measurements and InSAR (Hoffmann et al., 2003a). a) Model-derived compaction time constants and aquifer-system inelastic skeletal storage coefficients, where . b) Simulated subsidence and data residuals (simulated minus InSAR-derived). c) Simulated drawdowns and kriged measured drawdowns (from Galloway and Hoffmann, 2007).

6. Limitations and challenges

The chief limitation of the approach described above for analyzing and simulating aquifer-system compaction using the standard groundwater flow equation and the classic storage coefficient is the assumption of 1D vertical stress and strain inherent in the formulation of the skeletal specific storage coefficient. From the perspective of groundwater flow and accounting for the volume of groundwater released from or taken into storage this simplifying assumption is not very limiting. Gambolati et al. (2000) demonstrate that in sedimentary basins the results of 3D fully coupled hydro-mechanical models and 3D flow coupled with 1D subsidence models are practically indistinguishable at any time of practical interest, with the computational cost of the former, however, much higher than the latter.

From a hazards perspective regarding the 3D stresses and strains affecting formation of earth fissures, motions on surface faults and engineered structures on the land surface, the 1D deformation approach is limiting and is inappropriate for evaluating these local effects. Because groundwater sustainability issues increasingly are focusing on consequences of groundwater depletion, these hazard issues are becoming more important in groundwater resources management strategies (Alley et al., 1999). To address these concerns more rigorous and robust approaches are needed to evaluate 3D stresses and strains accompanying groundwater abstraction. Such modeling approaches based on finite-element formulations of the Biot (1941) poroelastic coupled fluid flow and mechanical deformation are being used to address these local effects (e.g. Hernandez-Marin and Burbey, 2009). However, the application of these methods in hydrogeology is currently the subject of ongoing research and not widely accessible to hydrogeologists, costly in terms of commercial software procurement and computational burden, and generally inapplicable in current forms to complex heterogeneous regional groundwater flow problems. Though, in sedimentary basins, a decoupled approach using a 3-D flow model followed by a 3-D mechanical model can be a very good trade-off between a 1D consolidation model and a fully coupled 3D Biot model. In numerical terms there are definite advantages in using two symmetric positive definite models in place of one large symmetric but indefinite model such as the Biot model that may also display non trivial problems of numerical ill-conditioning as shown by Ferronato et al. (2001).

The principal challenge seems to be one of scale as the significant horizontal components of deformation in these systems generally are focused in key local areas where heterogeneities contribute to large lateral variations in vertical stresses, e.g. near pumping wells, the edges of groundwater basins, or near faults or significant facies changes internal to the aquifer system. One interesting new approach might include embedded modeling approaches where local-scale, Biot-type 3D models could be embedded in a simpler regional-scale model. Recent advances in local grid refinement for MODFLOW (Mehl et al., 2006) and extension to embedding finite-element models in regional finite-difference models (Dickinson et al., 2007) suggest that an embedded approach is promising.

Another limitation of the regional MODFLOW approach is the inability presently to model aquifer systems that simultaneously exhibit features simulated by the SUB package (time-dependent drainage and compaction of thick aquitards) and the SUB-WT package (changing geostatic stress and stress-dependent aquitard storage properties). Incorporation of these capabilities into a single subsidence simulator would enhance flexibility and facilitate applications to commonly encountered developed aquifer systems with a water table aquifer overlying a confined aquifer.

Acknowledgments

This work was supported by the U.S. Geological Survey’s Federal Cooperative and Groundwater Resources Programs. Thoughtful colleague reviews by Stanley A. Leake and Jason P. Pope (both USGS) improved the manuscript. Gerald W. Bawden (USGS) provided helpful suggestions on T-LiDAR.

References

Alley, W.M., Reilly, T.E., Franke, O.L., 1999, Sustainability of ground-water resources: Denver, Colorado, United StatesU.S. Geological Survey, Circular 1186, 79 p. http://pubs.usgs.gov/circ/circ1186/index.html, cited 16 Sept 2010.

Bawden, G.W., 2002, Optimizing GPS arrays to image both tectonic and anthropogenic deformation: Eos Transactions, American Geophysical Union (Fall Meeting Supplement), 83, abstract G22A-11.

Bawden G.W., Thatcher W., Stein R.S., Hudnut K.W., Peltzer G., 2001, Tectonic contraction across Los Angeles after removal of groundwater pumping effects: Nature, 412, 812-815.

Biot, M.A., 1941, General theory of three-dimensional consolidation: Journal of Applied Physics, 12, 155-164.

Burbey, T.J., 2001a, Stress-strain analyses for aquifer-system characterization: Ground Water, 39, 128-136.

Burbey, T.J., 2001b, Storage coefficient revisited—Is purely vertical strain a good approximation?: Ground Water, 39, 458-464.

Burbey, T.J., 2002, The influence of faults in basin-fill deposits on land subsidence, Las Vegas, Nevada, USA: Hydrogeology Journal, 10, 525-538.

Burbey, T.J., 2005, Use of vertical and horizontal deformation data with inverse models to quantify parameters during aquifer testing, in Zhang A., Gong S., Carbognin L., Johnson A.I. (eds.), Land Subsidence: Proceedings of the 7^thInternational Symposium on Land Subsidence: Shanghai, Shanghai Scientific & Technical Publishers 560-569.

Carpenter, M.C., 1993, Earth-fissure movements associated with fluctuations in ground-water levels near the Picacho Mountains, south-central Arizona, 1980–84: Washington, D.C., United StatesU.S. Geological Survey, Professional Paper 497-H, 49 p. http://pubs.er.usgs.gov/usgspubs/pp/pp497H, cited 6 Sept 2009.

Cooper, H.H., Jr., 1966, The equation of groundwater flow in fixed and deforming coordinates: Journal of Geophysical Research, 71, 4785-4790.

Dickinson, J.E., James, S.C., Mehl, S.W., Hill, M.C., Leake, S.A., Zyvoloski, G.A., Faunt, C.C., Eddebbarh, A.A., 2007, A new ghost-node method for linking different models and initial investigations of heterogeneity and nonmatching grids: Advances in Water Resources, 30, 1722-1736, doi:10.1016/j.advwatres.2007.01.004.

Epstein, V.J., 1987, Hydrologic and geologic factors affecting land subsidence near Eloy, Arizona: Tucson, Arizona, U.S. Geological Survey, Water-Resources Investigations Report 87-4143, 28 p. http://pubs.er.usgs.gov/usgspubs/wri/wri874143, cited 14 Sept 2010.

Faunt, C.C. (ed.), 2009, Ground-water availability of California’s the Central Valley Aquifer, California: Sacramento, California, U.S.Geological Survey, Professional Paper 1766, 225 p. http://pubs.usgs.gov/pp/1766/, cited 15 Sept 2010.

Ferronato M., Gambolati G., Teatini P., 2001, III-conditioning of finite element poroelasticity equations: International Journal of Solids & Structures, 38, 5995-6014.

Galloway, D.L., Bawden, G.W., Leake, S.A., Honegger D.G., 2008, Land subsidence hazards, inBaum, R. L., Galloway, D.L., Harp, E.L., (eds), Landslide and land subsidence hazards to pipelines: U.S. Geological Survey Open-File Report 2008-1164, chapter 2, http://pubs.usgs.gov/of/2008/1164/, cited 4 Sept 2010.

Galloway, D.L., Hoffmann, J., 2007, The application of satellite differential SAR interferometry-derived ground displacements in hydrogeology: Hydrogeology Journal, 15, 133-154, doi: 10.1007/s10040-006-0121-5.

Galloway, D.L., Hudnut, K.W., Ingebritsen, S.E., Phillips, S.P., Peltzer, G., Rogez, F., Rosen, P.A., 1998, Detection of aquifer system compaction and land subsidence using interferometric synthetic aperture radar, Antelope Valley, Mojave Desert, California: Water Resources Research, 34, 2573-2585.

Galloway, D., Jones, D.R., Ingebritsen, S.E. (eds), 1999, Land subsidence in the United States: Reston Virginia, United States Geological Survey, 177 p. U.S. Geological Survey Circular 1182, http://pubs.usgs.gov/circ/circ1182/, cited 4 Sept 2010.

Gambolati G., 1973, Equation for one-dimensional vertical flow of groundwater: 1. The rigorous theory: Water Resources Research, 9, 1022-1028.

Gambolati, G., Freeze, R.A., 1973, Mathematical simulation of the subsidence of Venice: 1. Theory: Water Resources Research, 9, 721-733.

Gambolati G., Baú, D., Teatini P., Ferronato M., 2000, Importance of poroelastic coupling in dynamically active aquifers of the Po river basin, Italy: Water Resources Research, 36, 2443-2459.

González, J.L., Törnqvist, T.E., 2006, Coastal Louisiana in crisis—Subsidence or sea level rise?: Eos, Transactions, American Geophysical Union, 87, 493-498.

Hanson, R.T., 1989, Aquifer-system compaction, Tucson Basin and Avra Valley, Arizona: Tucson, Arizona, U.S. Geological Survey, Water-Resources Investigations Report 88-4172, 69 p. http://pubs.er.usgs.gov/usgspubs/wri/wri884172, cited 12 Sept 2010.

Hanson, R.T., Anderson, S.R., Pool, D.R., 1990, Simulation of ground-water flow and potential land subsidence, Avra Valley, Arizona: U.S. Geological Survey Water-Resources Investigations Report 90-4178, http://pubs.er.usgs.gov/usgspubs/wri/wri904178, cited 15 Sept 2010.

Hanson, R.T., Benedict, J.F., 1994, Simulation of ground-water flow and potential land subsidence, upper Santa Cruz Basin, Arizona: U.S.United States Geological Survey, Water-Resources Investigations Report 93-4196, 47 p. http://pubs.er.usgs.gov/usgspubs/wri/wri934196, cited 15 Sept 2010.

Hanson, R.T., Martin, P., Koczot, K.M., 2003, Simulation of ground-water/surface-water flow in the Santa Clara—Calleguas ground-water basin, Ventura County, California: Sacramento, California, United States Geological Survey, Water–Resources Investigations Report 02-4136, 214 p.

Hanson, R.T., Li, Z., Faunt, C.C., 2004, Documentation of the Santa Clara Valley regional ground-water/surface-water flow model, Santa Clara County, California: Sacramento, California, United States Geological Survey, Scientific Investigations Report 2004-5231, 75 p.

Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model—The ground-water flow process, in U.S. Geological Survey Techniques and Methods: Reston, Virginia, United States Geological Survey.

Harbaugh, A.W., McDonald, M.G., 1996, User’s documentation for MODFLOW-96, an update to the U.S. Geological Survey modular finite-difference ground-water flow model: Reston, Virginia, United States Geological Survey, Open-File Report 96-485, 56 p.

Harbaugh, A.W., Banta, E.R., Hill, M.C., McDonald, M.G., 2000, MODFLOW-2000, the U.S. Geological Survey modular ground-water model—User guide to modularization concepts and the ground-water flow process: Reston, Virginia, United States Geological Survey, Open-File Report 00-92, 121 p.

Helm, D.C., 1972, Simulation of aquitard compaction due to changes in stress (abstract): Eos, Transactions, American Geophysical Union, 53, 979.

Helm, D.C., 1975, One-dimensional simulation of aquifer system compaction near Pixley, Calif. 1. Constant parameters: Water Resources Research, 11, 465-478.

Helm, D.C., 1976, One-dimensional simulation of aquifer system compaction near Pixley, Calif. 2. Stress-dependent parameters: Water Resources Research, 12, 375-391.

Helm, D.C., 1978, Field verification of a one-dimensional mathematical model for transient compaction and expansion of a confined aquifer system, Verification of mathematical and physical models in hydraulic engineering, in Proceedings of the 26^thHydraulic Division Specialty Conference: College Park, Maryland, American Society of Civil Engineers, 189-196.

Helm, D.C., 1984a, Field-based computational techniques for predicting subsidence due to fluid withdrawal, inHolzer, T.L. (ed.), Man-induced Land Subsidence: Reviews in Engineering Geology, 6, 1-22.

Helm, D.C., 1984b, Latrobe Valley subsidence predictions: The modeling of time-dependent ground movement due to groundwater withdrawal: Joint Report of Fuel Department and Design Engineering and Environment Department, State Electricity Commission of Victoria, Melbourne, 2 vols.

Helm, D.C., 1986, COMPAC: A field-tested model to simulate and predict subsidence due to fluid withdrawal: Australasian Geomechanics Computing Newsletter, 10, 18-20.

Helm, D.C., 1994, Horizontal aquifer movement in a Theis-Thiem confined aquifer system: Water Resources Research, 30, 953-964.

Hernandez-Marin, M., Burbey, T.J., 2009, The role of faulting on surface deformation patterns from pumping induced ground-water flow: Hydrogeology Journal, 17, 1859-1875, doi: 10.1007/s10040-009-0501-8.

Heywood, C.E., 1993, Monitoring aquifer compaction and land subsidence due to ground-water withdrawal in the El Paso, Texas–Juarez, Chihuahua area, inPrince, K.R., Galloway, D.L., Leake, S.A. (eds.), U.S. Geological Survey Subsidence Interest Group Conference, Edwards Air Force Base, Antelope Valley, California, 1992—Abstracts and summary: Reston, Virginia, United States Geological Survey, Open-File Report 94-532, 14-17.

Hoffmann, J, Zebker, H.A., 2003, Prospecting for horizontal surface displacements in Antelope Valley, California, using satellite radar interferometry: Journal of Geophysical Research - Earth Surface, 108, doi:10.1029/2003JF000055.

Hoffmann J., Zebker, H.A., Galloway, D.L., Amelung, F., 2001, Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry: Water Resources Research, 37, 1551-1566.

Hoffmann, J., Galloway, D.L., Zebker, H.A., 2003a, Inverse modeling of interbed storage parameters using land subsidence observations, Antelope Valley, California: Water Resources Research, 39, 1031-1042, doi: 10.1029/2001WR001252.

Hoffmann, J., Leake, S.A., Galloway, D.L., Wilson, A.M., 2003b, MODFLOW-2000 ground-water model—User guide to the subsidence and aquifer-system compaction (SUB) package: Tucson, Arizona, U.S. United States Geological Survey, Open-File Report 03-233, 46 p.http://pubs.usgs.gov/of/2003/ofr03-233/, cited 14 Sept 2010.

Holzer, T.L., 1981, Preconsolidation stress of aquifer systems in areas of induced land subsidence: Water Resources Research, 17, 693-704.

Holzer, T.L., 1998, History of the aquitard-drainage model in land subsidence case studies and current research, inBorchers JW (ed,), Land Subsidence Case Studies and Current Research. Proceedings of the Dr. Joseph F. Poland Symposium on Land Subsidence, 1995, Sacramento, California: Belmont, California, Star Publishing, Association of Engineering Geologists, 7-12.

Holzer, T.L., Galloway, D.L., 2005, Impacts of land subsidence caused by withdrawal of underground fluids in the United States, in Ehlen, J., Haneberg, W.C., Larson, R.A. (eds.), Humans as Geologic Agents: Reviews in Engineering Geology, 16, 87-99, doi: 10.1130/2005.4016(08).

Hsieh, P.A., 1996, Deformation-induced changes in hydraulic head during ground-water withdrawal: Ground Water, 36, 1082-1089.

Hwang, C., Hung, W.C., Liu, C.H., 2008, Results of geodetic and geotechnical monitoring of subsidence for Taiwan High Speed Rail operation: Natural Hazards, 47, 1-16, doi: 10.1007/s11069-007-9211-5.

Ikehara, M.E., Phillips, S.P., 1994, Determination of land subsidence related to ground-water level declines using global positioning system and leveling surveys in Antelope Valley, Los Angeles and Kern Counties, California, 1992: Reston, Virginia, United States Geological Survey, Water-Resources Investigations Report 94-4184, 101 p.

Ingebritsen, S.E., Jones, D.R., 1999, Santa Clara Valley, California—A case of arrested subsidence, inGalloway, D.L., Jones, D.R., Ingebritsen, S.E. (eds.), Land Subsidence in the United States: Reston, Virginia, United States Geological Survey, Circular 1182, 15-22.

International Panel on Climate Change (IPCC), 2001, Climate change 2001—the scientific basis: Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate: Cambridge and New York, Cambridge University Press, 881 p.

Ireland, R.L., Poland, J.F., Riley, F.S., 1984, Land subsidence in the San Joaquin Valley, California as of 1980: Washington, D.C., United States Geological Survey, Professional Paper 437-I, 93 p.

Jacob, C.E., 1940, On the flow of water in an elastic artesian aquifer: Transactions, American Geophysical Union, 21, 574-586.

Kasmarek, M.C., Strom, E.W., 2002, Hydrogeology and simulation of ground-water flow and land-surface subsidence in the Chicot and Evangeline aquifers, Houston-Galveston region, Texas: Austin, Texas, United States Geological Survey, Water-Resources Investigations Report 02-4022, 68 p. http://pubs.er.usgs.gov/usgspubs/wri/wri024022, cited 15 Sept 2010.

King, N.E., Argus, D., Langbein. J., Agnew, D.C., Bawden, G., Dollar, R.S., Liu, Z., Galloway, D., Reichard, E., Yong, A., Webb, F.H., Bock, Y., Stark, K., Barseghian, D., 2007, Space geodetic observation of expansion of the San Gabriel Valley, California, aquifer system, during heavy rainfall in winter 2004–2005: Journal of Geophysical Research, 112, doi:10.1029/2006JB004448.

Leake, S.A., 1990, Interbed storage changes and compaction in models of regional ground-water flow: Water Resources Research, 26, 1939-1950.

Leake, S.A., 1991, Simulation of vertical compaction in models of regional ground-water flow, inJohnson, A.I. (ed.), Land Subsidence: Proceedings of the Fourth International Symposium on Land Subsidence: Houston, Texas. International Association of Hydraulic Sciences, 565-574.

Leake, S.A., Galloway, D.L., 2007, MODFLOW ground-water model—User guide to the Subsidence and Aquifer-System Compaction Package (SUB-WT) for Water-Table Aquifers: U.S. Geological Survey Techniques and Methods: Reston, Virginia, United States Geological Survey, Techniques and Methods 6–A23, 42 p 6–A23, http://pubs.usgs.gov/tm/2007/06A23/, cited 14 Sept 2010.

Leake, S.A., Galloway, D.L., 2010, Use of the SUB-WT Package for MODFLOW to simulate aquifer-system compaction in Antelope Valley, California, USA, inCarreón-Freyre, D., Cerca, M., Galloway, D.L. (eds.), Land Subsidence, Associated Hazards and the Role of Natural Resources Development: Proceedings of the Eighth International Symposium on Land Subsidence: Queretaro, Mexico, International Association of Hydraulic Sciences, 61-67.

Leake, S. A., Prudic, D. E., 1991, Documentation of a computer program to simulate aquifer-system compaction using the modular finite-difference ground-water flow model, inTechniques of Water-Resources Investigations: Reston, Virginia, United States Geological Survey, 68 p.

Leighton, D.A., Phillips, S.P., 2003, Simulation of ground-water flow and land subsidence in the Antelope Valley ground-water basin, California: Sacramento, California, United States Geological Survey, Water-Resources Investigations Report 03-4016, 107 p.

Li, J., 2007a, Transient radial movement of a confined leaky aquifer due to variable well flow rates: Journal of Hydrology, 333, 542-553.

Li, J., 2007b, Analysis of radial movement of an unconfined leaky aquifer due to pumping and injection: Hydrogeology Journal, 15, 1063-1076.

Li, J., 2007c, Analysis of radial movement of a confined aquifer due to pumping and injection: Hydrogeology Journal, 15, 445-458.

Liu, C.H., Pan, Y.W., Liao, J.J., Hung, W.C., 2004, Estimating coefficients of volume compressibility from compression of strata and piezometric changes in multiaquifer system in west Taiwan: Engineering Geology, 75, 33-47, doi: 10.1016/j.enggeo.2004.04.007.

Liu, Y., Helm, D.C., 2008a, Inverse procedure for calibrating parameters that control land subsidence caused by subsurface fluid withdrawal: 1. Methods. Water Resources Research, 44, doi:10.1029/2007WR006605.

Liu, Y., Helm, D.C., 2008b, Inverse procedure for calibrating parameters that control land subsidence caused by subsurface fluid withdrawal: 2. Field application. Water Resources Research, 44, doi:10.1029/2007WR006606.

McDonald, M.G., Harbaugh, A.W., 1988, A modular three-dimensional finite-difference ground-water flow model, inU.S. Geological Survey Techniques of Water-Resources Investigations: Reston, Virginia, U.S.United States Geological Survey,, book 6, chap. A1, http://pubs.usgs.gov/twri/twri6a1/, cited 14 Sept 2010 586 p.

Mehl, S., Hill, M.C., Leake, S.A., 2006, Comparison of local grid refinement methods for MODFLOW: Ground Water, 44, 792-796, doi:10.1111/j.1745-6584.2006.00192.x.

Metzger, L.F., Ikehara, M.E., Howle, J.F., 2002, Vertical-deformation, water-level, microgravity, geodetic, water-chemistry, and flow-rate data collected during injection, storage, and recovery tests at Lancaster, Antelope Valley, California, September 1995 through September 1998: Sacramento, California, United States Geological Survey, Open-File Report 01-414, 149 p.http://pubs.usgs.gov/of/2001/ofr01414/, cited 6 Sept 2010.

Narasimhan, T.N., Witherspoon, P.A., 1977, Numerical model for land subsidence in shallow groundwater systems, in: Proceedings of the Second International Symposium on Land Subsidence, 1976: Anaheim, California, International Association of Hydraulic SciencesDec 1976, IAHS Publication 121, , 133–-144143., http://iahs.info/redbooks/a121/iahs_121_0133.pdf, cited 14 Sept 2010.

Neuman, S.P., Preller, C., Narasimhan, T.N., 1982, Adaptive explicit-implicit quasi three-dimensional finite element model of flow and subsidence in multiaquifer systems: Water Resources Research, 18, 1551-1561.

Pavelko, M.T., 2003, Estimates of hydraulic properties from a one-dimensional numerical model of vertical aquifer-system deformation, Lorenzi Site, Las Vegas, Nevada: Carson City, Nevada, U.S. United States Geological Survey, Water-Resources Investigations Report 03-4083, 36 p.http://pubs.usgs.gov/wri/wri034083/index.html, cited 15 Sept 2010.

Poland, J.F. (ed), 1984, Guidebook To to Studies Of of Land Subsidence Due To to Ground-Water Withdrawal: Paris, UNESCO, Studies and reports in hydrology 40, http://wwwrcamnl.wr.usgs.gov/rgws/Unesco/, cited 4 Sept. 2010305 p.

Poland, J.F., Davis, G.H., 1969, Land subsidence due to withdrawal of fluids, inVarnes, D.J., Kiersch, G. (eds.): Reviews in Engineering Geology, 2, 187-269.

Poland, J.F., Ireland, R.L., 1988, Land subsidence in the Santa Clara Valley, California, as of 1982: Reston, Virginia, United States Geological Survey, Professional Paper 497-F, 61 p. http://pubs.er.usgs.gov/usgspubs/pp/pp497F, cited 4 Sept. 2010.

Pope, J.P., Burbey, T.J., 2003, Characterization and modeling of land subsidence due to ground-water withdrawals from the confined aquifers of the Virginia Coastal Plain, inPrince, K.R., Galloway, D.L. (eds.), U.S. Geological Survey Subsidence Interest Group Conference: Proceedings of the technical meeting, 2001: Galveston, Texas, United States 27–29 Nov 2001, U.S. Geological Survey, Open-File Report 03-308, 49-–56, http://pubs.usgs.gov/of/2003/ofr03-308/, cited 14 Sept 2010.

Pope, J.P., Burbey, T.J., 2004, Multiple-aquifer characterization from single borehole extensometer records: Ground Water, 42, 45-58.

Riley, F.S., 1969, Analysis of borehole extensometer data from central California, inTison, L.J. (ed.), Land subsidenceSubsidence: Proceedings of the Tokyo Symposium: Tokyo, International Association of Hydrological Sciences, Sept 1969, IAHS Publication 88, 423-431, http://iahs.info/redbooks/a088/088047.pdf, cited 6 Sept 2010.

Riley, F.S., 1986, Developments in borehole extensometry, inJohnson, A.I., Carbognin, L., Ubertini, L. (eds.), Land subsidenceSubsidence: Proceedings of the Third International Symposium on Land Subsidence, 1984:, Venice, Italy, International Association of Hydrological SciencesMar 1984, IAHS Publication 151, 169-–186, http://iahs.info/redbooks/a151/iahs_151_0169.pdf, cited 6 Sept 2010.

Scott, R.F., 1963, Principles of Soil Mechanics: Palo Alto, California, Addison-Wesley, 550 p.

Sneed, M., 2008, Aquifer-System Compaction and Land Subsidence Data and Simulations, the Holly Site, Edwards Air Force Base, California: Sacramento, California, California State University-Sacramento, M.S. thesis, 40 p.

Sneed, M., Brandt, J.T., 2007, Detection and measurement of land subsidence using global positioning system surveying and interferometric synthetic aperture radar, Coachella Valley, California, 1996-2005: Reston, Virginia, United States Geological Survey, Scientific Investigations Report 2007-5251, 31 p. http://pubs.usgs.gov/sir/2007/5251/, cited 17 Sept 2010.

Sneed, M., Galloway, D.L., 2000, Aquifer-system compaction and land subsidence: measurements, analyses, and simulations―the Holly site, Edwards Air Force Base, Antelope Valley, California: Reston, Virginia, United States Geological Survey, Water-Resources Investigation Report 00-4015, 65 p. http://pubs.usgs.gov/wri/2000/wri004015/, cited 12 Sept 2010.

Sneed, M., Ikehara, M.E., Galloway, D.L., Amelung, F., 2001, Detection and measurement of land subsidence using global positioning system and interferometric synthetic aperture radar, Coachella Valley, California, 1996-98: Sacramento, California, U.S. Geological Survey, Water-Resources Investigations Report 01-4193, 21 p. http://pubs.usgs.gov/wri/wri014193/, cited 17 Sept 2010.

Sneed, M., Stork, S.V., Ikehara, M.E., 2002, Detection and measurement of land subsidence using global positioning system and interferometric synthetic aperture radar, Coachella Valley, California, 1998-2000: Sacramento, Califoria, United States Geological Survey, Water-Resources Investigations Report 02-4239, 29 p. http://pubs.usgs.gov/wri/wri024239/, cited 17 Sept 2010.

Terzaghi, K., 1925, Settlement and consolidation of clay: Engineering News-Record, 95, 874-878.

Tolman, C.F., Poland, J.F., 1940, Ground-water, salt-water infiltration, and ground-surface recession in Santa Clara Valley, Santa Clara County, California: Transactions, American Geophysical Union, 1, 23-35.

Wilson, A.M., Gorelick, S.. 1996, The effects of pulsed pumping on land subsidence in the Santa Clara Valley, California: Journal of Hydrology, 174, 375-396.

Witherspoon, P.A., Freeze, R.A., 1972, The role of aquitards in a multiple-aquifer system: Geotimes, 17(4), 22-24.

Wolff, R.G., 1970, Relationship between horizontal strain near a well and reverse water-level fluctuations: Water Resources Research, 6, 1721-1728.

Manuscript received: December 18, 2010.
Corrected manuscript received: April 2, 2012.
Manuscript accepted: April 2, 2012.