Land Subsidence Due to Creep of the Gulf Coast Aquifer System in the Houston-Galveston Region

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.

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 m²/d and that the transmissivity of the Evangeline aquifer ranges from 915 to 4,575 m²/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]).

View largeDownload slide

View Large

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

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

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 m³/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.

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 Δt_d 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.

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 H_C, Evangeline aquifer thickness H_E, and Burkeville confining unit (clay) thickness H_B 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⁻⁵ to 3.94 × 10⁻³ (dimensionless), with an average of 1.38 × 10⁻³, at the 13 borehole extensometer locations. The Chicot aquifer C_αvalue is in the range of 8.74 × 10⁻⁵ to 2.55 × 10⁻³ 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⁻³ to 1.60 × 10⁻³, with an average of 1.48 × 10⁻³ 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⁻⁵ to 3.94 × 10⁻³.

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.

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.

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⁻⁴ to 3.94 × 10⁻³ at 12 out of 13 borehole extensometer sites, which matches a secondary consolidation coefficient C_α range of 1 × 10⁻⁴ to 5 × 10⁻³ found by Mesri (1973) for most soils from laboratory tests. Only one C_α value for clayey silt dominant aquitards is 8.72 × 10⁻⁵ 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).

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 _c(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⁻⁵ to 3.94 × 10⁻³, with an average of 1.38 × 10⁻³. 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⁻⁴ to 3.94 × 10⁻³ 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⁻⁴ to 5 × 10⁻³ 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.