The compaction measurements of Quaternary and Tertiary Gulf Coast aquifer system sediments in the Houston-Galveston region (TX) show spatially variable compression of 0.08 to 8.49 mm/yr because of geohistorical overburden pressure when groundwater levels in the aquifer system were stable after about the year 2000. An aquifer-system creep equation is developed for evaluating this variable compression, with a thickness-weighted average creep coefficient based on Taylor's (1942) secondary consolidation theory. The temporal variation of aquifer system creep can be neglected in a short-term observation period (such as a decade) after a long-term creep period (such as over 1,000 years) in geohistory. The creep coefficient of the Gulf Coast aquifer system is found to be in a range of 8.74 × 10−5 to 3.94 × 10−3 (dimensionless), with an average of 1.38 × 10−3. Moreover, for silty clay or clay-dominant aquitards in the Gulf Coast aquifer system the creep coefficient value varies in the range of 2.21 × 10−4 to 3.94 × 10−3, which is consistent with values found by Mesri (1973) for most soils, which vary in the range of creep coefficient, 1 × 10−4 to 5 × 10−3. Land subsidence due to secondary consolidation of the Gulf Coast aquifer system is estimated to be 0.04 to 4.33 m in the 20th century and is projected to be 0.01 to 0.64 m in the 21st century at the 13 borehole extensometer locations in the Houston-Galveston region. The significant creep should be considered in the relative sea level rise, in addition to tectonic subsidence and primary consolidation.

The Houston-Galveston region (HGR) in Texas, comprising Harris, Galveston, Fort Bend, Montgomery, Brazoria, Liberty, San Jacinto, Walker, Grimes, Waller, and Chambers Counties (Figure 1), has been one of the largest areas of land subsidence (LS) in the United States. It was in the early 1900s that the Houston area began to see the first true signs of human-induced LS―initially attributed to the extraction of oil gas from beneath the surface (Pratt and Johnson, 1926). The Houston area has been subsiding because of the combined effects of groundwater withdrawal, hydrocarbon extraction, salt dome movement, and faulting (Qu et al., 2015). By 1979, 3 m of subsidence had occurred in the HGR, and approximately 8,288 km2 of land had subsided more than 0.3 m (Coplin and Galloway, 1999; Kasmarek et al., 2015). Over southeastern Harris County, the maximum subsidence reached 4 m during the 1915–1917 and 2001 periods (Kasmarek et al., 2009). LS caused by fluid withdrawals was first documented in the HGR in the Goose Creek oil field in southeastern Harris County (Pratt and Johnson, 1926). Subsidence continued throughout the 20th century as a result of groundwater withdrawal that depressurized the major aquifers in this area, thereby resulting in the compaction of the aquifer sediment (Kasmarek et al., 2010).
Figure 1.

Hydrogeologic section A–A′ of the Gulf Coast aquifer system in Grimes, Montgomery, Harris, and Galveston Counties, TX (modified from Braun et al. [2019] and Baker [1979, 1986]).

Figure 1.

Hydrogeologic section A–A′ of the Gulf Coast aquifer system in Grimes, Montgomery, Harris, and Galveston Counties, TX (modified from Braun et al. [2019] and Baker [1979, 1986]).

Historically, groundwater has been the primary water source for industrial use, municipal supply, and irrigation, and groundwater use in the HGR had sharply increased for a few decades to meet the needs of the rapidly growing population (Seifert and Drabek, 2006). In addition, the complex geologic setting, laterally diverse subsurface hydrological units, regional faults, hydrocarbon extraction, and salt dome movement in the HGR (Coplin and Galloway, 1999) have made it difficult to characterize the source(s) of the observed LS (Qu et al., 2015).

A network of 13 discrete borehole extensometers was installed within the HGR to monitor groundwater-level changes and measure accumulated clay compaction to better understand the extent and magnitude of the regional subsidence in the 1970s (Kasmarek et al., 2010). It was identified that most of the subsidence in the HGR has occurred as a direct result of groundwater withdrawals that depressurized and dewatered the Chicot and Evangeline aquifers, thereby causing compaction of the aquifer sediments (Galloway et al., 1999; Kasmarek, 2013; and Kasmarek et al., 2015). Historically, groundwater withdrawn from the Chicot, Evangeline, and Jasper aquifers had been the primary source of water for municipal supply, commercial and industrial use, and irrigation in the HGR since the early 1970s (Kasmarek, 2013; HGSD, 2017). This is the “primary consolidation” of the unconsolidated Quaternary and semi-consolidated Tertiary aquifer systems (Liu et al., 2020) due to groundwater withdrawal. In the study area, sand layers are more transmissive and less compressible than are fine-grained clay and silt layers, and these sand layers depressurized more rapidly than did the clay and silt layers. When groundwater withdrawing rates change, pressure equilibrium is to reestablish more rapidly in the sand layers than that in the clay and silt layers.

The amount of compaction of the sand layers is usually minor when compared to that of the clay and silt layers (Trahan, 1982; Galloway et al., 1999). Because most compaction of subsurface sediments is inelastic and largely permanent only a small amount of rebound of the land-surface elevation can occur because of unloading or increase of pore water pressure. While the compaction of one thin clay and silt layer typically will not cause a measurable decrease in the land-surface altitude, a measurable amount of subsidence can occur when an aquifer system comprises numerous stratigraphic sequences of sand layers and clay and silt layers (e.g., characteristic of the Gulf Coast aquifer system) that are subjected to depressurizing and compaction (Gabrysch and Bonnett, 1975). Groundwater-level fluctuations are measured by the U.S. Geological Survey (USGS) in over 700 wells in an 11-county area annually in the HGR to develop a regional depiction of groundwater levels. The cumulative compaction in the Chicot and Evangeline aquifers is measured at the 13 borehole extensometer stations in this region, with data collection extending from 1973 to 1980. Compaction measurements of the 13 borehole extensometers recorded through about the year 2000 corroborate primary consolidation analysis with monitored groundwater-level changes (Liu et al., 2019). However, a nearly constant rate of subsidence with spatial variation of 0.08 to 8.49 mm/yr emerged while the groundwater levels in the two aquifers were being managed to be stable in trend since about 2000 after groundwater-level elevations had recovered to −42.4 m in 2010 (from −97.1 m in 1990) in the Chicot aquifer and to −56.0 m in 2005 (from −125.0 m in 1984) in the Evangeline aquifer, respectively (Liu et al., 2019). This gives rise to the following question: Is the nearly constant rate of consolidation of 0.08 to 8.49 mm/yr still attributed to the primary consolidation due to groundwater withdrawal?

The time dependency on the volume change of clay observed in one-dimensional compression has been an active research area for decades. For practical applications, the volume changes over time are arbitrarily divided into “primary consolidation” and “secondary consolidation,” though the latter is widely known to occur in the entire compression process. In the case of primary consolidation, the volume changes are mainly governed by the changes in effective stress (hydrodynamic effect). The volume changes occurring after the primary consolidation or during the secondary consolidation are dominated by the viscous behavior of the clay and silt system. Within a compressible unit during the secondary consolidation the effective stress remains relatively constant after the pore water pressure reaches equilibrium for a given boundary condition. Most research activities have focused on modeling the time-dependent or viscous response during secondary consolidation observed in the laboratory or in the field. Upon loading, the hydrodynamic and viscous effects occur simultaneously. It is very challenging to determine those two separate effects experimentally or analytically. In addition, study on basic mechanisms contributing to the intrinsic viscous behavior that occurs in creep compression is still limited. It was suggested that the recent nearly constant rate of consolidation of 0.08 to 8.49 mm/yr within the Gulf Coast aquifer system could be attributed to secondary consolidation in sedimentation due to geohistorical overburden pressure or self-weight (Liu et al., 2019; Liu, Li, Fasullo, et al., 2020; and Liu, Li, Fang, et al., 2020). The goal of this article is to present an equation with an equivalent creep coefficient and to apply it to quantify creep deformation of the Gulf Coast aquifer system in the 20th and 21st centuries, respectively, based on 13 borehole extensometers’ compaction measurements in the HGR.

From northwest to southeast, the HGR investigated in this article includes Grimes County (with an elevation close to 122 m), Montgomery County, Waller County, Harris County, and Galveston County (with an elevation from 0 to 15 m) along the coast of Gulf of Mexico (Figure 1). The three primary Quaternary and Tertiary aquifers in the Gulf Coast aquifer system in the HGR study area are the Chicot, Evangeline, and Jasper aquifers (Figures 1 and 2), which comprise laterally discontinuous deposits of gravel, sand, silt, and clay. The youngest and uppermost Quaternary Chicot aquifer consists of Holocene- and Pleistocene-age sediments; the underlying Tertiary Evangeline aquifer comprises Pliocene- and Miocene-age sediments; and the oldest and most deeply buried Tertiary Jasper aquifer consists of Miocene-age sediments (Figures 1 and 2) (Baker, 1979, 1986). The Burkeville confining unit between the Evangeline and Jasper aquifers consists of Miocene Fleming Formation Lagarto clay. The Miocene-age Catahoula confining system, which includes the Catahoula Sandstone, is the lowermost unit of the Gulf Coast Tertiary aquifer system. The Catahoula confining system consists of sands in the upper section and clay and tuff interbedded with sand in the lower section.
Figure 2.

Geologic and hydrogeologic units of the Gulf Coast aquifer system in the Houston-Galveston region study area, TX (modified by Kasmarek et al. [2015] from Sellards et al. [1932]; Baker [1979]; Meyer and Carr [1979]).

Figure 2.

Geologic and hydrogeologic units of the Gulf Coast aquifer system in the Houston-Galveston region study area, TX (modified by Kasmarek et al. [2015] from Sellards et al. [1932]; Baker [1979]; Meyer and Carr [1979]).

Numerous authors have contributed to the body of knowledge and understanding of the complex stratigraphic and hydrogeologic relations of the Gulf Coast aquifer system in the HGR study area (Figure 2). Using this information, a series of groundwater flow models were created, the most recent being that of Kasmarek (2013) on hydrogeology and simulation of groundwater flow and land-surface subsidence in the Northern part of the Gulf Coast aquifer system, Texas, 1891–2009 (HAGM, 2013). These models provide an evaluative tool that can be used by water-resource managers to help regulate and conserve the important natural water resource of the aquifer system.

The percentage of clay and other fine-grained clastic material generally increases with depth downdip (Baker, 1979). Over time, geologic and hydrologic processes created accretionary sediment wedges (stacked sequences of sediments) that are more than 2,318 m thick at the coast (Figure 1) (Chowdhury and Turco, 2006). The sediments composing the Gulf Coast aquifer system were deposited by fluvial-deltaic processes and subsequently were eroded and redeposited (re-worked) by worldwide episodic changes in sea level (eustacy) that occurred because of oscillations between glacial and interglacial climate conditions (Lambeck et al., 2002). The Gulf Coast aquifer system consists of hydrogeologic units that dip and thicken from northwest to southeast (Figure 1); the aquifers thus crop out in bands inland from and approximately parallel to the coast and become progressively more deeply buried and confined toward the coast (Kasmarek, 2013). The Burkeville confining unit restricts groundwater flow between the Evangeline and Jasper aquifers by being stratigraphically positioned between the Evangeline and Jasper aquifers (Figure 1).

There is no significant confining unit between the Chicot and Evangeline aquifers; therefore, the aquifers are hydraulically connected, which allows groundwater flow between the aquifers (Figure 1). Because of this hydraulic connection, water-level changes occurring in one aquifer can affect water levels in the adjoining aquifer (Kasmarek and Robinson, 2004).Supporting evidence for the interaction of groundwater flow between the Chicot and Evangeline aquifers is demonstrated by comparing the two long-term (1977–2015) water-level–change maps (Kasmarek et al., 2015), indicating that the areas in which water levels have declined or risen are approximately spatially coincident. Hydraulic properties of the Chicot aquifer do not differ appreciably from those of the hydrogeologically similar Evangeline aquifer but can be differentiated based on hydraulic conductivity (Carr et al., 1985). From aquifer test data, Meyer and Carr (1979) estimated that the transmissivity of the Chicot aquifer ranges from 915 to 7,625 m2/d and that the transmissivity of the Evangeline aquifer ranges from 915 to 4,575 m2/d. The Chicot aquifer outcrops and extends inland from the Gulf of Mexico coast and terminates at the most northern updip limit of the aquifer. Proceeding updip and inland of the Chicot aquifer, the older hydrogeologic units of the Evangeline aquifer, the Burkeville confining unit, and the Jasper aquifer sequentially outcrop (Figure 1). The aquifer in the outcrop and updip areas of the Jasper aquifer can be differentiated from the Evangeline aquifer based on the depths to water below land-surface datum, which are shallower (closer to the land surface) in the Jasper aquifer compared to those in the Evangeline aquifer. Additionally, in the downdip parts of the aquifer system, the Jasper aquifer can be differentiated from the Evangeline aquifer on the basis of stratigraphic position relative to the elevation of the Burkeville confining unit (Figure 1). Table 1 illustrates clay vertical hydraulic conductivity and inelastic and elastic-specific storage values with burial depth estimated porosity by Kelley and Deeds (2019).

Table 1.

Estimated porosity, specific compressibility, and specific storage of clay beds as a function of depth of burial (from Kelley and Deeds [2019]).

Figure 3 shows the 11 borehole extensometer station locations in the HGR (Kasmarek and Lanning-Rush, 2003). There are two borehole extensometers (shallow and deep) at the Baytown and Clear Lake stations, respectively. Each of the other nine stations has only one borehole extensometer. In total, there are 13 borehole extensometers. Detailed information on the scientific theory, construction, and operation of borehole extensometers is presented in Gabrysch (1984). Five borehole extensometers were installed in Harris (four) and Galveston (one) counties and began recording compaction data in July 1973: East End, Baytown Shallow, Baytown Deep, and Seabrook in Harris County and Texas City in Galveston County (see Figure 3). The borehole extensometer Johnson Space Center was installed in 1962 in Harris County to record compaction. Additional borehole extensometers were added to the network at four locations in Harris County during the 1974–1976 period: Addicks in 1974, Pasadena in 1975, and Clear Lake Shallow and Clear Lake Deep in 1976 (see Figure 3). In 1980, the final three borehole extensometers were installed in Harris County: Lake Houston, Northeast, and Southwest (see Figure 3). Since activation or installation that has taken place between 1973 and 1980, compaction data have been constantly recorded and periodically collected about every 28 days at the 13 borehole extensometers, thereby providing site-specific rates of compaction that are accurate to within 0.3 mm (Kasmarek et al., 2015).
Figure 3.

Location of borehole extensometer sites, Houston Galveston region, TX (modified from Kasmarek and Lanning-Rush [2003]). (Note: The location of wells LJ-65-21-229 and LJ-65-21-227 is at the same location of borehole extensometer Southwest).

Figure 3.

Location of borehole extensometer sites, Houston Galveston region, TX (modified from Kasmarek and Lanning-Rush [2003]). (Note: The location of wells LJ-65-21-229 and LJ-65-21-227 is at the same location of borehole extensometer Southwest).

Subsidence and Its Components

In this article, subsidence denotes an observed field cumulative land subsidence, S(t), of a compressible aquifer system with a total thickness of H from one borehole extensometer. S(t) consists of two components: primary and secondary/creep consolidations, Sp(t) and Sc(t), respectively as the below Eq. 1:
formula
(1)
Subsidence rate can be written in Eq. 2 from Eq. 1:
formula
(2)
where (t) = dS/dt, forumlap(t) = dSp/dt, and forumlac(t) = dSc/dt. The secondary consolidation rate forumlac can be identified from a subsidence curve with time when the primary consolidation rate forumlap(t) becomes zero [i.e., (t) = forumlac(t)]. If forumlap−v and forumlap−e denote subsidence rates of primary consolidation due to changes in inelastic (virgin) storativity of aquitards and elastic storativity of sand layers and aquitards, respectively, in the aquifer system, Eq. 2 can be further written as the following Eq. 3:
formula
(3)

Each term in Eq. 3 will be applied to the qualitative analysis or interpretation using the subsidence data measured from fields.

Time for Completion of Primary Consolidation

Terzaghi (1925) developed an analytical solution to simulate the equilibration of pore water pressure in an individual saturated clay layer (or aquitard), with a uniform initial pore water pressure where only vertical flow is permitted, in response to a specified instantaneous step change (note: practically gradual or step-by-step change) in the hydraulic head at the bottom and top of its surrounding upper and lower aquifers, respectively. This process elucidates the theory of pore water dissipation process (consolidation), which was extended to the analysis of aquitard/confining unit drainage (Riley, 1969) in an aquifer system and subsequently to the simulation of aquitard/confining unit drainage. This concept, which was referred to as “the aquitard drainage model” by Helm (1984), has formed the theoretical basis of many successful subsidence investigations (Holzer, 1998). For the doubly-draining clay layer, the same aquifer hydraulic head changes are assumed to occur at the upper and lower surfaces. A time constant (not infinity) for all aquitards in the aquifers system is τ′0 (=Ssk (b0 equiv/2)2/K′, where Ssk, b0 equiv, and K′ denote vertical skeletal-specific storage (dimensionless), equivalent thickness (L), and vertical hydraulic conductivity (L/T) of the aquitards, respectively) (Riley, 1969). The primary consolidation degree in the aquifer system reaches 93.1 percent, 99.4 percent, and 100 percent when time factor Tv (=Δt/τ′0) (where Δt is real time or Terzaghi time (T), since pore pressure is immediately reduced to be lower than pre-consolidation pressure for activation of inelastic consolidation) equals 1, 2, and ∞, respectively. The term “no-delay interbeds” in Hoffmann et al. (2003) is used to denote the interbeds for which τ′0 is short compared to the time steps used in the MODFLOW simulation, while the primary consolidation degree is 93.1 percent. The primary consolidation is considered fully completed when the time factor Tv reaches 2 with consolidation of 99.4 percent (Liu and Li, 2016). (t) ≈ forumlac(t) when time is larger than Δt and forumlap(t) ≈ 0 based on Eq. 2. In practice, the inelastic consolidation rate forumlap−v in Eq. 3 is dominant, with delay when pore pressure in aquitard(s) within the aquifer system and in confining bed(s) is kept lower than its/their preconsolidation pressure, since the magnitude of inelastic-specific skeletal storage is about 2 or 3 magnitude orders larger than that of elastic-specific storage. The deformation of the aquifer system performs elastically because the aquifer system's sand storativity is usually much larger than the aquitard/confining layer's storativity after the inelastic deformation ends. Therefore, the trend elastic deformation disappears when groundwater level approaches stability for a long-term period in trend after inelastic consolidation ends, while the secondary consolidation (creep) of the aquifer system appears. The time between the historical lowest groundwater level and the end of inelastic compaction of the Gulf Coast aquifer system is denoted by the inelastic delay time Δtd in this article. Δtd can be determined by analyzing the relationship between measured compaction and water level of the Gulf Cost aquifer system. The inelastic delay time factor TV−v can be estimated from Δtd using the following Eq. 4:
formula
(4)
The consolidation degree is larger than 99.4 percent when TV−v is equal to 2 if there are lowering and/or recovery periods of groundwater level before and after the lowest groundwater level. The equivalent aquitard thickness of the Gulf Coast aquifer system in the HGR is estimated to be 3.4 m to 3.7 m in Figure 4 based on the fact that the Δtd values are 27 years and 32 years, respectively. For instance, Figure 4B indicates that if (a) Δtd = 32 years, (b) an aquitard equivalent thickness = 3.7 m, and (c) the values of K′ and S′skv in Table 1 are applied to the burial depth of 0–900 m, TV−v will be larger than 2, which means that the consolidation degree of the aquitard will be larger than 99.4 percent. Furthermore, based on Figure 4B, if the equivalent thickness increases to 6 m, only those aquitards located within 0–300-m burial depth can reach more than 99.4 percent consolidation degree (TV−v > 2) in 32 years. Thus, Eq. 4 can be a tool for the first-hand estimate to evaluate if the process of primary consolidation is completed.
Figure 4.

Time factors (TV−v) for inelastic compaction of aquitards in the Gulf Coast aquifer system estimated by Eq. 4 with K′ and Sskv values in Table 1 and assumed equivalent aquitard thickness b0 equiv. (A) b0 equiv of the Gulf Coast aquifer system is probably equal to 3.4 m for Δtd = 27 years, and (B) b0 equiv is probably equal to 3.7 m for Δtd = 32 years.

Figure 4.

Time factors (TV−v) for inelastic compaction of aquitards in the Gulf Coast aquifer system estimated by Eq. 4 with K′ and Sskv values in Table 1 and assumed equivalent aquitard thickness b0 equiv. (A) b0 equiv of the Gulf Coast aquifer system is probably equal to 3.4 m for Δtd = 27 years, and (B) b0 equiv is probably equal to 3.7 m for Δtd = 32 years.

Secondary Consolidation Equation of Gulf Coast Aquifer System

As a result of the weight of the overburden and the inelastic compaction characteristics of the clay layers, about 90 percent of the compaction is permanent (Gabrysch and Bonnett, 1975) during loading. Three main sedimentation stages are defined with respect to the concentration degree in self-weight consolidation (Tong et al., 2011): the clarification regime (suspension), zone-settling regime (settling without consolidation), and compression regime (consolidation) (Fitch, 1983). The above Quaternary and Tertiary aquifer systems are still in the third compression stage. This compression was referred to as “secondary consolidation” by Taylor and Merchant (1940) or as “self-weight consolidation” by Been and Sills (1981). Let Hi and Cαi denote the thickness (L) and the secondary consolidation coefficient (dimensionless) of soil layer i, which includes any individual sand, silt, and clay layers within a compressible Gulf Coast aquifer system with a total soil layer number N. The compaction Sci(t) of each individual aquitard layer, i, is assumed to be in a secondary consolidation due to geohistorical overburden pressure and follows Taylor's (1942) creep equation Sci(t) = CαiHi log t/t0 (Taylor, 1942), where t and t0 represent creep time and tt0. If a bulk creep deformation Sc(t) in Eq. 1 is measured from a borehole extensometer for a compressible aquifer system, Sc can be written as follows:
formula
(5)
where forumla and forumla, which is the aquifer system thickness. So, Cα denotes a secondary consolidation coefficient of an aquifer system. The Cα value of the aquifer system can be dominated by that of aquitards or confining layers when compared to the negligible Cα of sand layers.

Pseudo-Constant Rate of Secondary Consolidation

Creep rate forumlac(t) = (CαH/ln10) 1/t follows Eq. 5 by taking the derivative with respect to time t. The decrease percentage (DS) of forumlac(t) from t to t + Δt was derived with [forumlac(t) − forumlac(t + Δt)]/forumlac(t), as follows:
formula
(6)
DS approaches zero when t ≫ Δt, which implies that forumlac(t) ≈ a constant when t ≫ Δt. Put differently, the changing value of forumlac(t) over the Δt period (such as 10- or 20-year observation period) is difficult to be identified and can be ignored. This negligibly variable rate is referred to as the “pseudo-constant rate of secondary consolidation” (Liu et al., 2019). Figure 5 shows that creep subsidence rate decrease percentage (DS) changes with the secondary consolidation time (t) for each given time period (Δt) of 5, 10, 20, 30, 40, and 50 years. Table 2 displays how many years are needed for specified subsidence rate percentage changes of 1.0 percent, 0.5 percent, and 0.1 percent in one given period. For example, if an observation period (Δt) is 10 years, 990; 1,990; and 9,990 years are needed for specified subsidence rate decrease percentages of 1.0 percent, 0.5 percent, and 0.1 percent, respectively. The secondary consolidation rate forumlac(t) during an observation period (Δt) is considered a pseudo-constant if t ≫ Δt. For example, over a 14-year observation during the 2003–2017 period (Δt = 14 year) and 1,000-year creep (t = 1,000), based on Eq. 6, DS = 1.38 percent, which means a value of forumlac(t) = 3.87 mm/yr will drop by 1.38 percent and reduce to 3.82 mm/yr over a 1,000-year creep (see Figure5). If this drop is ignored, forumlac(t) is approximately considered as a constant. This approximation is also applied to other forumlac(t) values, with a range changing from 0.08 to 8.49 mm/yr in the HGR (Liu et al., 2019).
Figure 5.

Percentage change of secondary consolidation subsidence rate with time from Eq. 6 for Δt = 10, 20, 30, 40, and 50 years, respectively.

Figure 5.

Percentage change of secondary consolidation subsidence rate with time from Eq. 6 for Δt = 10, 20, 30, 40, and 50 years, respectively.

Table 2.

Time of the secondary consolidation needed for a specified subsidence rate decrease in given periods (modified from Liu et al. [2019]).

Estimate of Secondary Consolidation Subsidence

From Eq. 5, the secondary consolidation subsidence ΔSc from time t1 to time t2 can be estimated by the following Eq. 7:
formula
(7)

With the data from fields, the coefficient Cα can be calibrated using Eq. 7.

Roughly Steady Groundwater Withdrawal since 2001

Artificial primary consolidation first occurred in the early 1900s in areas where relatively large volumes of groundwater, oil, and gas were extracted. Primary consolidation continued throughout the 20th century, due primarily to groundwater pumpage (Galloway et al., 1999; Kasmarek, 2013). Groundwater withdrawal in the HGR experienced the following four periods (see Figure 6C): (1) a slight withdrawal rate of about 0.19 million m3/d from 1891 to 1930 for 40 years; (2) increasing withdrawal rates for 46 years from 0.57 million m3/d in 1931 to 4.28 million m3/d in 1976, with an average growth of 0.05 million m3/d per year; Near Texas City, the withdrawal of groundwater for public supply and industry caused more than 0.5 m of subsidence between 1906 and 1943 (Galloway et al., 1999); (3) decreasing withdrawal rates for 25 years from 4.28 million m3/d in 1976 to 3.03 million m3/d in 2001; and (4) roughly steady withdrawal rates of around 3 million m3/d from 2001 to date for more than 17 years.
Figure 6.

Groundwater withdrawal and groundwater-level fluctuations in the Houston-Galveston region: (A) Measured and simulated groundwater level in the Chicot aquifer (modified from Kasmarek [2013]); (B) Measured and simulated groundwater level in Evangeline aquifer (modified from Kasmarek [2013]); and (C) Groundwater withdrawal (data from Kasmarek [2013] after Liu et al. [2019]).

Figure 6.

Groundwater withdrawal and groundwater-level fluctuations in the Houston-Galveston region: (A) Measured and simulated groundwater level in the Chicot aquifer (modified from Kasmarek [2013]); (B) Measured and simulated groundwater level in Evangeline aquifer (modified from Kasmarek [2013]); and (C) Groundwater withdrawal (data from Kasmarek [2013] after Liu et al. [2019]).

Groundwater-Level Change with Groundwater Withdrawal

The lowering of groundwater level in the HGR due to groundwater withdrawal from 1891 to 1976 (Figure 6C) was observed and simulated by the USGS using steel tape measurements and MODFLOW software. Based on their simulation results in Figure 6A and B, from 1891 to 1900, the groundwater levels were about 21.35 m and 9.15 m in Chicot and Evangeline aquifers, respectively. This would be close to the status observed during the pre-development of groundwater before 1891. From 1901 to 1930, groundwater levels were lowered to 8.2 m and 5.2 m in Chicot and Evangeline aquifers, respectively. Based on measured results in Figure 6A and B, in 1976, the lowest groundwater levels corresponding to the highest regional groundwater pumpage of 4.28 million m3/d were −82 m and −86 m in Chicot and Evangeline aquifers, respectively; then, groundwater levels were raised about 41.5 m to −40 m in 2008 for Chicot aquifer and about 43 m to −43 m in 2008 for Evangeline aquifer. After 2008, groundwater levels have been roughly stable.

The lowest groundwater-level depression cones for Chicot and Evangeline aquifers in 1976 likely are very close to the observed lowest depression cones in Figure 7 A and B, respectively, which were found by the USGS Subsidence Reviewer. Through checking groundwater-level data monitored and published by USGS, the lowest groundwater levels were −97.1 m on 12 January 1990 (compared to 16 January 1984) at Well LJ-65-21-150 for the Chicot aquifer and −125.0 m on 9 January 1984 (compared to 9 January 1978) at Well LJ-65-13-904 for the Evangeline Aquifer, respectively. The maximum drawdown caused by groundwater withdrawal in 1976 was estimated to be about 112 m for the Chicot aquifer and 134 m for the Evangeline aquifer, respectively, in the HGR. Based on the two wells, groundwater levels were raised by about 54.6 m to −42.4 m in 2010 for the Chicot aquifer and about 64.4 m to −56.0 m in 2005 for the Evangeline aquifer. At the two wells, these aquifers’ groundwater levels were roughly stable in trend after 2010 and 2005, respectively.
Figure 7.

Measured groundwater-level depression cones (A) in Chicot aquifer (lowest: −97.14 m on 12 January 1990 compared to 16 January 1984) and (B) in Evangeline aquifer (lowest: −125.04 m on 9 January 1984 compared to 9 January 1978) (from USGS shape files on https://pubs.usgs.gov/ds/793/).

Figure 7.

Measured groundwater-level depression cones (A) in Chicot aquifer (lowest: −97.14 m on 12 January 1990 compared to 16 January 1984) and (B) in Evangeline aquifer (lowest: −125.04 m on 9 January 1984 compared to 9 January 1978) (from USGS shape files on https://pubs.usgs.gov/ds/793/).

Identification of Secondary Consolidation of the Gulf Coast Aquifer System

Theoretically, secondary consolidation fully appears in trend when groundwater levels are stable in trend, while inelastic primary consolidation disappears in all compressible aquifer systems. Appendix A (https://www.aegweb.org/e-eg-supplements) shows how an appearance period of the pseudo-constant secondary compaction forumlac was identified in the Chicot and Evangeline aquifers from 1980 to 2017 from the total compaction measured in the two Gulf Coast aquifer systems at borehole extensometer Southwest (see Figure 3 for location) in the HGR through empirical analysis of forumlap−v and forumlap−e variations and disappearances in six different periods of groundwater-level change (Figure 8). Thereafter, the period of the pseudo-constant secondary compaction forumlac in the Gulf Coast aquifer system for each of other 12 borehole extensometers was identified by application of this same methodology (shown in Figure 9) (Liu et al., 2019). Only the Chicot aquifer and Evangeline aquifers are involved with those 11 borehole extensometers, except Texas City and Baytown Shallow in Figure 3. A total of 19 groundwater-level wells (Table 3) in the two aquifers at or near the 11 locations of the 13 borehole extensometers were employed in the analysis of inelastic and elastic compaction variations forumlap−v and forumlap−e corresponding to trends of groundwater-level fluctuations for identifying the pseudo-constant periods of the secondary consolidation based on the methodology in section METHODOLOGY.
Figure 8.

Primary consolidation due to groundwater withdrawal ended in about 2000, and secondary consolidation became apparent at borehole extensometer Southwest in Houston when groundwater-level trends were stable. (I) Inelastic subsidence dominated by primary consolidation; (II) Subsidence dominated by elastic rebound due to unloading; (III) Subsidence due to primary consolidation (delay compaction) and creep offset by rebound during unloading; (IV) Subsidence caused by both primary consolidation and secondary consolidation, but it is primarily elastic compaction, because reloading stress is less than the historical one; (V) Subsidence due to secondary consolidation offset by rebound during unloading; and (VI) Subsidence dominated by secondary consolidation represented with a red trend line. Slope of the trend line (forumlac) is 0.0106 mm/d or 3.87 mm/a in a pseudo-constant rate; 427.76 mm in the logarithmically linear equation shows a value for 0.4343 CbαH in Eq. 5 (modified from Liu et al. [2019]).

Figure 8.

Primary consolidation due to groundwater withdrawal ended in about 2000, and secondary consolidation became apparent at borehole extensometer Southwest in Houston when groundwater-level trends were stable. (I) Inelastic subsidence dominated by primary consolidation; (II) Subsidence dominated by elastic rebound due to unloading; (III) Subsidence due to primary consolidation (delay compaction) and creep offset by rebound during unloading; (IV) Subsidence caused by both primary consolidation and secondary consolidation, but it is primarily elastic compaction, because reloading stress is less than the historical one; (V) Subsidence due to secondary consolidation offset by rebound during unloading; and (VI) Subsidence dominated by secondary consolidation represented with a red trend line. Slope of the trend line (forumlac) is 0.0106 mm/d or 3.87 mm/a in a pseudo-constant rate; 427.76 mm in the logarithmically linear equation shows a value for 0.4343 CbαH in Eq. 5 (modified from Liu et al. [2019]).

Figure 9.

Secondary consolidation (creep) (Sc) simulated with linear and logarithmically linear trends from observed cumulative compactions at 13 borehole extensometer locations in the Houston-Galveston region (data source: USGS). The secondary consolidation period for each site is given in Table 3. The slope values in the linear and logarithmically linear trend equations are in mm/d. The slope values in the logarithmically linear trend equations are dimensionless, which are employed to compute secondary consolidation (creep) coefficient values (see columns 7 and 8 in Table 4). (A) Extensometers Addicks, Texas City, Seabrook, ClearlakeShallow and Johnson Space Center; (B) Extensometers Northeast, Lake Houston, EastEnd, ClearLakeDeep and Southwest; and (C) Pasadena, BaytownShallow and BaytownDeep.

Figure 9.

Secondary consolidation (creep) (Sc) simulated with linear and logarithmically linear trends from observed cumulative compactions at 13 borehole extensometer locations in the Houston-Galveston region (data source: USGS). The secondary consolidation period for each site is given in Table 3. The slope values in the linear and logarithmically linear trend equations are in mm/d. The slope values in the logarithmically linear trend equations are dimensionless, which are employed to compute secondary consolidation (creep) coefficient values (see columns 7 and 8 in Table 4). (A) Extensometers Addicks, Texas City, Seabrook, ClearlakeShallow and Johnson Space Center; (B) Extensometers Northeast, Lake Houston, EastEnd, ClearLakeDeep and Southwest; and (C) Pasadena, BaytownShallow and BaytownDeep.

Table 3.

Secondary consolidation appearance periods based on groundwater-level trend at or near borehole extensometer sites in the Houston-Galveston region (modified from Liu et al. [2019]).

Table 4.

Identification of secondary consolidation (creep) coefficient values for the Gulf Coast aquifer system in the Houston-Galveston region.

Columns 5 and 6 in Table 3 illustrate the starting date and the ending date, respectively, of the full appearance of secondary consolidation at each borehole extensometer site at each of the 13 borehole extensometer sites (Southwest; Texas City; Seabrook; Johnson Space Center; Clear Lake Shallow and Deep; Baytown Shallow and Deep; Addicks; East End; Northeast; Pasadena; and Lake Houston) based on monitored groundwater-level data from wells near each site. The slope of the groundwater-level trend, such as from Figure 8, was given in column 7 in m/d and in column 8 in m/yr during the appearance period at each site. All other 10 starting dates (or inelastic delay compaction end dates) from 7 January 2004 through 25 January 2008 commenced after the starting date (18 September 2003) of the appearance of secondary consolidation at borehole extensometer site Southwest (Table 3). The inelastic compaction delay time Δtd relative to 1976 (the lowest groundwater time due to groundwater withdrawal) can be determined to be 27–32 years. All the periods are after the inelastic compaction ceased within period (V) from 20 September 2000 to 18 September 2003 identified in Figure 8. The groundwater-level trend values from 0.21 to 0.39 m/yr in column 8 in Table 3 are a bit larger for borehole extensometer sites East End, Northeast, and Pasadena, but the groundwater-level difference between the starting date and the ending date of the corresponding secondary consolidation appearance period is very small (0.02 to .014 m). This small groundwater-level difference signifies that the cumulative elastic deformation approaches zero during that period.

Previously, the reason for the huge subsidence fluctuations at Baytown Shallow and Deep and Pasadena from 2010 to 2014 (seeFigure 9C) was not well understood. Water usage in the area has not increased. The only other physical factor is reaction along faults in the area, which was detected by mapping ground deformation with multi-temporal InSAR (Qu et al., 2015). However, if this was indeed true, all borehole extensometer locations would be affected similarly (email communication with USGS hydrologist, Jason Ramage, 2018). Secondary compaction at all other 10 borehole extensometer sites seems to have dominated land subsidence without primary inelastic consolidation, occurring since around 2003. The secondary consolidation will continue for a very long period into the future, since inelastic compaction was fully completed in 2003, if the current groundwater management is kept with the stability of groundwater level.

Secondary Consolidation Coefficient of the Gulf Coast Aquifer System

Logarithmically linear simulations of secondary consolidation at the 13 borehole extensometers in the HGR since 1973 are depicted in Figure 9 with black solid trend lines. The equations for each logarithmically linear trend line simulated using Microsoft Excel show the slope values of the logarithmically linear trend line in millimeters. The slope values of 9.256 to 939.8 mm approximately represent 0.4343 CαH in Eq. 5 and are shown in column 7 in Table 4 . Column 6 in Table 4 shows all pseudo-constant secondary consolidation values from Table 3. The Chicot aquifer thickness HC, Evangeline aquifer thickness HE, and Burkeville confining unit (clay) thickness HB involved by borehole extensometer are given in columns 2, 3, and 4, respectively, in Table 4. The total thickness H of a Gulf Coast Aquifer system involved in one borehole extensometer is shown in column 5. Meanwhile, the values of the secondary consolidation coefficient Cα are displayed in column 8, as derived from values in column 7, with total thickness H values in column 5. The creep coefficient Cα value of the Gulf Coast aquifer system is in a range of 8.74 × 10−5 to 3.94 × 10−3 (dimensionless), with an average of 1.38 × 10−3, at the 13 borehole extensometer locations. The Chicot aquifer Cαvalue is in the range of 8.74 × 10−5 to 2.55 × 10−3 at borehole extensometers Texas City and Baytown Shallow. The Cα value of the Chicot and Evangeline aquifers in addition to the Burkeville confining unit is in the range of 1.38 × 10−3 to 1.60 × 10−3, with an average of 1.48 × 10−3 at borehole extensometers Southwest and Northeast. At the other nine borehole extensometer locations, the Cα value of the Chicot and Evangeline aquifers is in the range of 2.21 × 10−5 to 3.94 × 10−3.

Secondary Consolidation Subsidence in the 20Th and 21St Centuries, Respectively

The secondary consolidation subsidence of the involved Gulf Coast aquifer system in the 20th century is estimated to be 0.04 to 4.33 m (see column 9 in Table 4) at the 13 borehole extensometers. The highest creep of 4.33 m in the 20th century is at Addicts because the dominant aquitards are clay, with the highest secondary consolidation coefficient of 0.00394, which likely led to the low-lying areas near Addicts with similar geology. The secondary consolidation subsidence of the involved Gulf Coast aquifer system in the 21st century (2001–2100) is estimated to be 0.01 to 0.64 m (see column 10 in Table 4) at the 13 borehole extensometers. The creep will continue with a much smaller rate after 2100.

The Period of Primary Inelastic Consolidation

From Figure 6A and B, it can be said that the inelastic consolidation started in about 1942 during pumping period 2 (1931–1976), when the groundwater level was lowered below the pre-consolidation pressure level of −21.5 m almost at the same time in both Chicot (Well LJ-65-14-912) and Evangeline (Well LJ-65-22-618) aquifers (Kasmarek, 2013). This inelastic compaction was completed before a time from 18 September 2003 to 25 January 2008 (Table 3). Based on groundwater level in the two wells in 2000, the updated lowest pre-consolidation pressure level would be −50.3 m in the Chicot aquifer and −51.2 m in the Evangeline aquifer, respectively. Therefore, the primary consolidation is elastic or recoverable, since inelastic compaction was completed if the groundwater level in the Chicot aquifer and in the Evangeline aquifer here is higher than −50.3 m and −51.2 m, respectively. The period of primary inelastic consolidation period existed for about 61–66 years, approximately from 1942.

Dominant Aquitard(s) Inferred from Secondary Consolidation Coefficient Values

The values of secondary consolidation coefficient (Cα) found in Table 4 may be understood to represent characteristics of dominant aquitard(s) rather than sands in the Gulf Coast aquifer system. Magnitude orders of −3, −4, and −5 would apparently represent aquitard(s) dominated by clay, silty clay, and clayey silt, respectively, in the Gulf Coast aquifer system. The inferred dominant aquitard(s) in the Gulf Coast aquifer system at the 13 borehole extensometer locations is/are given in column 11 in Table 4. Most aquitards are clay or silty clay. The Cα value for silty clay or clay-dominant aquitards from column 8 in Table 4 is in the range of 2.21 × 10−4 to 3.94 × 10−3 at 12 out of 13 borehole extensometer sites, which matches a secondary consolidation coefficient Cα range of 1 × 10−4 to 5 × 10−3 found by Mesri (1973) for most soils from laboratory tests. Only one Cα value for clayey silt dominant aquitards is 8.72 × 10−5 at borehole extensometer site Texas City (Table 4), with aquitard(s) thickness percentage of 20.8 percent in the Chicot aquifer (Kasmarek and Robinson, 2004).

Significance of Secondary Consolidation in Gauged Sea Level Rise

The secondary consolidation existing in the Gulf Coast aquifer system in the HGR is found from the 13 borehole extensometer compaction measurements. The pseudo-constant creep rate during a very short period (such as a decade) after a long-term geological history of sedimentation (such as 1,000 years or longer) played an important role in this finding (Liu et al., 2019). The secondary consolidation coefficient values found in Table 4 help describe slowly variable creep rate with equation forumlac(t) = (CαH/ln10) 1/t. However, the secondary consolidation can be concealed by a significant primary inelastic consolidation for about 61 years from 1942 to 2003 due to huge groundwater withdrawals from the aquifer systems in the HGR. The secondary consolidation can fully appear while the inelastic consolidation is completed for the historical maximum pre-consolidation pressure in the Gulf Coast aquifer system, and elastic deformation approximately approaches zero in trend. Therefore, because of groundwater withdrawal, the secondary consolidation of the Gulf aquifer system is one land subsidence component, in addition to tectonic subsidence and primary consolidation subsidence. Tide gauge Galveston Pier 21 (see location in Figure 3) has the longest tide records (for 112 years since 1909) among the 25 gauges along the Gulf Coast of Mexico. Sea level in Galveston Bay along the Gulf Coast of Mexico has risen about 71 cm, with a mean rate of 6.51 mm/yr at tide gauge Galveston Pier 21 since 1909, based on National Oceanic and Atmospheric Administration (NOAA) estimation. This mean sea level rise (SLR) rate is 3.8 times larger than the global mean SLR rate of 1.7 mm/yr (Parris et al., 2012; Walsh et al., 2014). Tebaldi et al. (2012) and Zervas et al. (2013) estimated land subsidence rate values at some NOAA tide gauge stations by using the difference between relative sea level rise rate to global mean sea level rise of 1.70 mm/yr. For example, the land subsidence rate at the tide gauge Galveston Pier 21 (Figure 3) in the HGR was estimated to be 4.72 mm/yr (Zervas et al., 2013), which should include secondary consolidation subsidence of the Gulf Coast aquifer system at Galveston Pier 21 in addition to regional basement rock subsidence due to tectonics and primary consolidation subsidence due to groundwater withdrawal (Liu, Li, Fasullo, et al., 2020).

In this article, the land subsidence caused by compaction of the Gulf Coast aquifer system and measured by borehole extensometers in the HGR is assumed to be the sum of primary consolidation due to groundwater withdrawal and secondary consolidation due to geo-historical overburden pressure. The primary consolidation comprises inelastic and elastic components caused by groundwater-level or pore water pressure-lowering due to groundwater withdrawal from the compressible aquifer systems. The inelastic or non-recoverable consolidation is induced by the lowed pore water pressure head in aquitard(s) or confining units when it is lower than its pre-consolidation pressure. Inelastic consolidation dominates land subsidence when it happens because inelastic-specific skeletal storage value is about 2-3 magnitude orders larger than elastic-specific skeletal storage value. The secondary consolidation may exist, at least, for more than 1,000 years in the Gulf Coast aquifer system. In this article, the secondary consolidation coefficient of the Gulf Coast aquifer system is defined with a thickness-weighted average of each individual creep coefficient of aquitards in the aquifer system, based on Taylor's (1942) secondary consolidation theory. The rate of secondary consolidation behaves as a pseudo-constant characteristic, especially if it has elapsed over 1,000 years since the youngest and uppermost sediments of the Holocene Chicot aquifer were formed in the Greenlandian Age (4,200–8,200 years ago) and the Northgrippian Age (8,200–11,700 years ago), respectively.

Significant inelastic primary consolidation likely started around 1942, when the groundwater levels in the Chicot and Evangeline aquifers were below the pre-consolidation pressure of −21.35 m. According to Terzaghi's theory, any primary consolidation can be more than 99 percent completed when its consolidation time factor reaches 2 if there is no updated pre-consolidation pressure within aquitards or confining units. The inelastic primary consolidation of the Gulf Coast aquifer system in the HGR was completed from 2003 to 2008. Then, the inelastic primary consolidation disappeared or tended to disappear when the aquifer groundwater level exhibited stability, while the pseudo-constant secondary consolidation fully appeared. The equivalent thickness of aquitards in the Gulf Coast aquifer system is estimated to be 3.4 to 3.7 m in this article.

Thirteen borehole extensometers in the HGR were built to monitor the compaction of the Quaternary and Tertiary Gulf Coast aquifer since 1973. The pseudo-constant secondary consolidation rate was identified to be 0.08 to 8.49 mm/yr from the borehole extensometer compaction data. The secondary consolidation coefficient of the Gulf Coast aquifer system is found to be in a range of 8.74 × 10−5 to 3.94 × 10−3, with an average of 1.38 × 10−3. It is inferred that magnitude orders of −3, −4, and −5 in the secondary consolidation coefficient values would represent aquitard(s) dominated by clay, silty clay, and clayey silt, respectively, in the Gulf Coast aquifer system. The secondary consolidation coefficient value range of 2.21 × 10−4 to 3.94 × 10−3 for silt clay or clay-dominant aquitards in the Gulf Coast aquifer system matches well the individual secondary consolidation coefficient value range of 1 × 10−4 to 5 × 10−3 found by Mesri (1973) for most soils from laboratory tests. The secondary consolidation subsidence of the involved Gulf Coast aquifer system is estimated to be 0.04 to 4.33 m in the 20th century and is projected to be 0.01 to 0.64 m in the 21st century at the 13 borehole extensometers in the HGR. In addition to regional tectonic subsidence and land subsidence due to groundwater withdrawal, the significant creep subsidence of the Gulf Coast aquifer system in the 20th and 21st centuries suggests a need to consider the secondary consolidation subsidence due to geohistorical overburden pressure as a new factor of relative sea-level rise along the Gulf Coast of Mexico.

The compaction data and groundwater data were provided by USGS (https://txpub.usgs.gov/houston_subsidence/home/). Direct requests for these materials may be made to the provider, as indicated in the “Acknowledgments” section.

This research is supported by the National Science Foundation (NSF) grant 1832065 entitled “Identification of urban flood impacts caused by land subsidence and sea-level rise in the Houston-Galveston region” and Maryland Port Administration (MPA) grant 51831 “Use of dredged material to protect low-lying areas in the Chesapeake Bay.” The authors express their gratitude to Jason Ramage (USGS) and John Ellis (USGS) for accessing and explaining data and reviewing an early version of the manuscript and to William Mike Chrismer for his help in data collection, as well as assistance from Tranell Griffin with GIS mapping and Ermei Liu with mapping support. The authors appreciate two journal peer reviewers, J. M. Sharp and an anonymous reviewer, for their insightful and constructive comments that led to an improved manuscript.

Open Access Article