From northwest to southeast, the HGR includes Grimes County (with a high elevation of close to 122 m a.m.s.l.), Montgomery County, Waller County, Harris County, and Galveston County (with a low elevation of 0 m along the coast of the Gulf of Mexico). The Gulf Coast aquifer system consists of three Quaternary and Tertiary aquifers and two Tertiary confining units composed of laterally discontinuous deposits of gravel, sand, silt, and clay. The aquifer units in the HGR are the Chicot, Evangeline, and Jasper aquifers (Kasmarek, 2013). The youngest and uppermost Quaternary aquifer, the Chicot aquifer, consists of Holocene- and Pleistocene-age sediments; the underlying Tertiary Evangeline aquifer consists of Pliocene- and Miocene-age sediments; and the oldest and most deeply buried Tertiary aquifer, the Jasper aquifer, consists of Miocene-age sediments. Stratigraphically positioned between the Evangeline and Japer aquifers is the Tertiary Burkeville confining unit of Miocene age sediments, which restricts groundwater flow between the Evangeline and Jasper aquifers. The lowermost unit of the Gulf Coast Tertiary aquifer system is the Miocene-age Catahoula confining system, which includes the Catahoula Sandstone. The Catahoula confining system consists of sands in the upper section and clay and tuff interbedded with sand in the lower section.

The Gulf Coast aquifer system consists of hydrogeologic units that dip and thicken from northwest to southeast; the geologic units 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). There is no confining unit between the Chicot and Evangeline aquifers; therefore, the aquifers are hydraulically connected, which allows groundwater flow between the aquifers. Because of this hydraulic connection, water-level changes that occur in one aquifer can affect water levels in the adjoining aquifer (Kasmarek et al., 2015). Supporting evidence of 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, which indicate that the areas where water levels have risen or declined are approximately spatially coincident (Carr et al., 1985). Hydraulic properties of the Chicot aquifer do not differ appreciably from the hydrogeologically similar Evangeline aquifer, but can be differentiated on the basis of hydraulic conductivity. The transmissivity of the Chicot aquifer ranges from 915 to 7625 m² d⁻¹ and the transmissivity of the Evangeline aquifer ranges from 915 to 4575 m² d⁻¹ (Meyer and Carr, 1979). The geologic units that make up the Chicot aquifer outcrop and extend inland from the Gulf of Mexico and terminates at the most northern updip limit of these units. The recharge rate across the outcrop area ranged from 6.35 to 177.8 mm yr⁻¹ (Noble et al., 1996). 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. In the outcrop and updip areas of the Jasper aquifer, the aquifer can be differentiated from the Evangeline aquifer on the basis of the depths to water below land-surface datum, which are shallower (closer to 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 Burkeville confining unit.

Almost all the permanent subsidence of a compressible aquifer system occurs due to the irreversible (or inelastic/nonrecoverable) compression or consolidation of aquitards through a slow (delayed) process of aquitard drainage. This concept, which was labeled “the aquitard drainage model” by Helm (Helm, 1984), has formed the theoretical basis of many successful subsidence investigations (Riley, 1969). The relation between changes in groundwater levels and compression of the aquifer system is based on the principle of effective stress first proposed by Terzaghi (1925). By this principle, when groundwater levels decrease due to discharges from an aquifer system, under a constant total load the support previously provided by the pore-fluid pressure will be transferred to the skeleton of the aquitard. Namely, the change in pore fluid pressure will be converted to effective stress on the skeleton of aquitards, which in turn causes the aquifer system's compaction. If the current effective stress is larger than the preconsolidation stress, the compaction is nonrecoverable (inelastic). In contrast, if effective stress is less than the preconsolidation stress, the compaction is recoverable (elastic). Conversely, when groundwater levels increase due to recharge to the aquifer system, the support previously provided by the skeleton is reduced and the change in effective stress is shifted onto the pore fluid, which results in the aquifer system's elastic expansion (Galloway et al., 1999). Therefore, in a general case the primary consolidation or compaction (s_p) can consist of two components: inelastic compaction (s_p−v) related to nonrecoverable specific skeletal storage (S_skv) of the aquitard(s) or confining unit(s) and elastic compaction (s_p−e) associated with recoverable specific skeletal storage (S_ske) of aquitard(s) or confining unit(s) and sand layers in one aquifer system. Therefore, we have the following equation $s_{\mathrm{p}}=s_{\mathrm{p}-\mathrm{v}}+s_{\mathrm{p}-\mathrm{e}}$, which yields the primary compaction rate as defined by Eq. (1)

where ${\dot{s}}_{\mathrm{s}\left(t\right)}$ denotes rate of secondary consolidation.

Figure 1Groundwater levels in Chicot and Evangeline aquifers and land subsidence at extensometer site Southwest in the Houston-Galveston region. (Modified from Liu et al., 2019).

Download

The S_skv values are 2 to 3 orders of magnitude larger than the values of (S_ske). This leads to inelastic compaction dominating land subsidence when it happens. Based on Terzaghi's consolidation theory, both ${\dot{s}}_{\mathrm{p}-\mathrm{v}}$ and ${\dot{s}}_{\mathrm{p}-\mathrm{e}}$ can be considered to be approximately zero when their consolidation degrees reach 99.4 % where the time factor, T_v ($=\mathrm{\Delta }t/{{\mathit{\tau }}^{\prime }}_{\mathrm{0}}$, where Δt is real time (T) and ${{\mathit{\tau }}^{\prime }}_{\mathrm{0}}$ is Terzaghi's time constant), equals 2.

Because of the weight of the overburden and the inelastic compaction characteristics of the clay layers, about 90 % of the compaction is permanent (Gabrysch and Bonnett, 1975). Three main sedimentation stages are defined with respect to the degree of self-weight consolidation as: the clarification regime, zone-settling regime, and compression regime (Fitch, 1983). The above Quaternary and Tertiary aquifer systems are still in the third compression stage. This compression was referred to as secondary consolidation (creep) by Taylor (1942) or as “self-weight consolidation” by Been and Sills (1981). Therefore, it is assumed that secondary consolidation (creep) exists in the upper unconsolidated aquifer systems due to geohistorical overburden pressure. For an unconsolidated sediment layer with an initial thickness of H (L), the secondary consolidation s_s(t) can be approximated by Eq. (2)

where C_α is the dimensionless coefficient of secondary compression of the sediment layer, t₁ is an initial reference time for secondary compression, t is time larger than or equal to t₁. Equation (2) yields ${\dot{s}}_{\mathrm{s}\left(t\right)}=(C_{\mathit{\alpha }}H/\ln \mathrm{10})\frac{\mathrm{1}}{t}$ by taking the derivative with respect to time t for subsidence rate. The percentage decrease (D_S) of ${\dot{s}}_{\mathrm{s}\left(t\right)}$ from t to t+Δt can be derived as

with $D_{\mathrm{S}}=\mathrm{100}\left[{\dot{s}}_{\mathrm{s}\left(t\right)}-{\dot{s}}_{\mathrm{s}(t+\mathrm{\Delta }t)}\right]/{\dot{s}}_{\mathrm{s}\left(t\right)}$ as D_s approaches zero when t≫Δt, which implies that ${\dot{s}}_{\mathrm{s}}\approx$ a constant. In other words, the changing value of ${\dot{s}}_{\mathrm{s}}$ over the Δt period is negligible and can be ignored. This negligibly variable rate is called a pseudo-constant rate of secondary consolidation (Liu et al., 2019). For example, if a period (Δt) is considered to be 10 years, 990, 1990, and 9990 years are needed for specified subsidence rate decrease percentages of 1.0 %, 0.5 %, and 0.1 %, respectively. The secondary consolidation rate ${\dot{s}}_{\mathrm{s}}$ is a pseudo-constant if 1.0 %, 0.5 %, and 0.1 % subsidence rate changes are considered as negligible. The secondary consolidation for the Quaternary and Tertiary sediments can be considered to have been more than 1000 years since the youngest and uppermost sediments of the Holocene Chicot aquifer were formed in the Greenlandian Age (4200 to 8200 years ago) and the Northgrippian Age (8200 to 11 700 years ago).

Therefore a bulk subsidence rate ${\dot{s}}_{\left(t\right)}$ from aquifer-system compaction can be the sum of primary inelastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{v}\left(t\right)}$, primary elastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{e}\left(t\right)}$ and secondary compaction rate ${\dot{s}}_{\mathrm{s}\left(t\right)}$, i.e.,

Equation (4) (Liu et al., 2019) is employed to analyze extensometer-measured compaction rate for the three components in response to groundwater level changes in aquifers. Three distinct compaction characteristics for the three components must be correctly applied in this analysis: inelastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{v}\left(t\right)}$ is 10 to over 100 times larger than elastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{e}\left(t\right)}$ when groundwater levels are lower than preconsolidation pressure head (Epstein, 1987; Hanson, 1989; Helm, 1978; Liu and Helm, 2008; Sneed and Galloway, 2000); elastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{e}\left(t\right)}$ can be negative (land rebounding) while inelastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{v}\left(t\right)}$ decreases rapidly but is never negative when groundwater is recovering; and secondary compaction rate ${\dot{s}}_{\mathrm{s}\left(t\right)}$ does not change in response to changes in groundwater levels.

From 1931 to 1976, groundwater levels decreased as groundwater withdrawal rates increased from 0.57 to 4.3 million m³ d⁻¹, causing pressure heads in aquitards of the Chicot and Evangeline aquifers to continually decline (Kasmarek, 2013). In response to reductions in groundwater withdrawal rates from 4.3 to 3.0 million m³ d⁻¹ between 1976 and 2001 groundwater-levels rebounded and have stabilized since then. Only the Chicot aquifer and Evangeline aquifers are involved with the 11 extensometers except North East and Southwest. Of the 13 extensometers in operation in the HGR, six measure subsidence in the Chicot aquifer and the other seven measure subsidence in both the Chicot and Evangeline aquifers at their respective locations. A total of 19 groundwater level wells in the two aquifers at or near the 11 extensometer station locations were employed to provide supporting evidence to be used in conjunction with the extensometer data in order to analyze changes in inelastic and elastic compaction and secondary consolidation in response to trends in groundwater-level fluctuations based on the methodology (Eq. 4). As one example, Fig. 1 shows the results based on analysis of observed land subsidence at borehole extensometer Southwest, which measures compaction in both the Chicot and Evangeline aquifers from 17 June 1980 to 28 December 2017 and groundwater levels at monitoring piezometers LJ-65-21-229 (completed in the Chicot aquifer) and LJ-65-21-227 (completed in the Evangeline aquifer) from 4 April 1980 to 1 October 2018. The preconsolidation hydraulic head was set to be −21.35 m in HAGM model (Kasmarek, 2013) for the aquitards within the two aquifers. Six periods were divided based on variable rate characteristics of elastic, inelastic, and creep compaction corresponding to groundwater level change. From Fig. 1 during Period I (4 April 1980–16 August 1990), groundwater levels in the Chicot and Evangeline aquifers ranged from −65 to −85 m and from −92 to −121 m, respectively. During period I, both wells exhibited water-levels much lower than the required −21.35 m of the initial uniform preconsolidation pressure head, thus initiating inelastic compaction which continued due to delay even after historical lowest water levels were recorded in this period. Thus, the inelastic compaction from the aquitards dominated the subsidence at this location: ${\dot{s}}_{\left(t\right)}$ was 46.92 mm yr⁻¹ from 1980 to 1987 then decreased to 31.46 mm yr⁻¹ from 1988 to 1990. The subsidence characteristics during Period I would be ${\dot{s}}_{\mathrm{p}-\mathrm{v}}\gg {\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}>\mathrm{0}$. During Period II (16 August 1990 to 25 March 1993), groundwater levels in the Chicot and Evangeline aquifers increased from −83 to −60 m and from −110 to −82 m, respectively. The 23 and 28 m groundwater level rise caused a land rebounding rate of 14.9 mm yr⁻¹. Thus, elastic rebounding of the two aquifers dominated the deformation at this location. The subsidence characteristics during Period II would be ${\dot{s}}_{\mathrm{p}-\mathrm{e}}<\mathrm{0},{\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{s}}>\mathrm{0}$ and $\dot{s}={\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}<\mathrm{0}$. During Period III (25 March 1993 to 22 January 1998), groundwater levels in the Chicot and Evangeline aquifers were further raised about 11 m to reach −49 m and about 9 to −72 m, respectively. The additional 9 to 11 m in groundwater level recovery did not cause land surface to rebound further, although the elastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{e}}$ was less than zero. The trend in the subsidence rate approached approximately zero, which implies ${\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{s}}>\mathrm{0}$ and ${\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{s}}\approx -{\dot{s}}_{\mathrm{p}-\mathrm{e}}$ from Eq. (4). Thus, the elastic rebounding of the two aquifers approximately offset the combination of inelastic compaction and secondary consolidation at this location. The subsidence characteristics during Period III would be ${\dot{s}}_{\mathrm{p}-\mathrm{e}}<\mathrm{0},{\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{s}}>\mathrm{0}$ and $\dot{s}={\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}\approx \mathrm{0}$. During Period IV (22 January 1998 to 20 September 2000), groundwater levels in the Chicot and Evangeline aquifers were lowered about 17 to −65 m and about 13 to −82 m, respectively. The 13 to 17 m maximum groundwater level lowering caused land subsidence in trend and the elastic compaction rate ${\dot{s}}_{\mathrm{p}-\mathrm{e}}$ is larger than zero. The inelastic consolidation from aquitards within the two aquifers continued for more than about 21 years with a decreasing rate which approached zero during this period, since the lowest historical regional groundwater levels occurred in response to the period of maximum groundwater withdrawals from 1977 to 1984. The subsidence characteristics in Period IV would be ${\dot{s}}_{\mathrm{p}-\mathrm{e}}>\mathrm{0},{\dot{s}}_{\mathrm{p}-\mathrm{v}}\to \left(\text{approaches}\right)\mathrm{0}$ and $\dot{s}={\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}>\mathrm{0}$. During Period V (20 September 2000 to 18 September 2003), groundwater levels in the Chicot and Evangeline aquifers rose an additional 5 to −50 m and about 10 to −70 m, respectively. The 5 to 10 m groundwater level rise caused neither further land rebounding nor significant subsidence. This happened only when inelastic compaction ceased $({\dot{s}}_{\mathrm{p}-\mathrm{v}}\approx \mathrm{0}$) and when elastic rebounding offset the secondary consolidation $({\dot{s}}_{\mathrm{s}}\approx -{\dot{s}}_{\mathrm{p}-\mathrm{e}}$). Thus, it would appear that the delay in compaction from inelastic specific skeletal storage of aquitards (Helm, 1984) within the Chicot and Evangeline aquifers at extensometer site Southwest ceased during or before 2000. The subsidence characteristics in Period V would be ${\dot{s}}_{\mathrm{p}-\mathrm{e}}<\mathrm{0},{\dot{s}}_{\mathrm{p}-\mathrm{v}}\approx \mathrm{0}$ and $\dot{s}={\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}\approx \mathrm{0}$. During the last Period VI (18 September 2003 to 28 December 2017), groundwater levels in the Chicot and Evangeline aquifers exhibited an almost stable trend of $\mathrm{1.13}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹ (Fig. 1) (0.03 m yr⁻¹) and $\mathrm{4.59}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹ (Fig. 1) (0.14 m yr⁻¹), respectively. This leads to the conclusion that the trend in elastic compaction can be considered negligible $({\dot{s}}_{\mathrm{p}-\mathrm{e}}\approx \mathrm{0}$). Only secondary consolidation emerged $({\dot{s}}_{\mathrm{s}}>\mathrm{0}$) since both ${\dot{s}}_{\mathrm{p}-\mathrm{e}}\approx \mathrm{0}\text{and}{\dot{s}}_{\mathrm{p}-\mathrm{v}}\approx \mathrm{0}$. Thus the subsidence characteristics in Period VI would be ${\dot{s}}_{\mathrm{p}-\mathrm{e}}\approx \mathrm{0},{\dot{s}}_{\mathrm{p}-\mathrm{v}}\approx \mathrm{0}$, therefore $\dot{s}={\dot{s}}_{\mathrm{p}-\mathrm{v}}+{\dot{s}}_{\mathrm{p}-\mathrm{e}}+{\dot{s}}_{\mathrm{s}}\approx {\dot{s}}_{\mathrm{s}}$, which equals 0.0106 mm d⁻¹ (3.87 mm yr⁻¹) (Fig. 1).

The secondary consolidation rate value range of 0.08 to 8.49 mm yr⁻¹ was found by applying the above analysis to 10 other extensometer sites: Texas City, Seabrook, Johnson Space Center and Clear Lake (which shares two groundwater wells), Baytown, Addicks, East End, Northeast, Pasadena, and Lake Houston, each of which has monitored groundwater level data from wells at or near the sites. The temporal variation of the pseudo-constant secondary consolidation could not be easily identified from the current compaction observations in the HGR because its emerging period is still relatively short at just 2 to 15 years. However, the spatial variation of secondary consolidation rate depends on the extensometer location. The reason for this spatial variation would include total thickness of compressible aquifer systems, clay cumulative thickness percentage, individual aquitard/confining unit thickness, and overburden pressure history.

The compaction data and groundwater data are accessible on https://txpub.usgs.gov/houston_subsidence/home/ (USGS, 2020, last access: 22 March 2020).

YL conceived this study, developed the methodology, conducted data curation, analysis and simulation, as well as wrote the manuscript. JL contributed to the methodology development. MR and TG contributed to data curation. JL and ZNF contributed to review and revision of the manuscript.

The authors declare that they have no conflict of interest.

This article is part of the special issue “TISOLS: the Tenth International Symposium On Land Subsidence – living with subsidence”. It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17–21 May 2021.

This research has been supported by the The U.S. National Science Foundation (grant no. 1832065).