Land subsidence contributions to relative sea level rise at tide gauge Galveston Pier 21, Texas

Relative sea level rise at tide gauge Galveston Pier 21, Texas, is the combination of absolute sea level rise and land subsidence. We estimate subsidence rates of 3.53 mm/a during 1909–1937, 6.08 mm/a during 1937–1983, and 3.51 mm/a since 1983. Subsidence attributed to aquifer-system compaction accompanying groundwater extraction contributed as much as 85% of the 0.7 m relative sea level rise since 1909, and an additional 1.9 m is projected by 2100, with contributions from land subsidence declining from 30 to 10% over the projection interval. We estimate a uniform absolute sea level rise rate of 1.10 mm ± 0.19/a in the Gulf of Mexico during 1909–1992 and its acceleration of 0.270 mm/a2 at Galveston Pier 21 since 1992. This acceleration is 87% of the value for the highest scenario of global mean sea level rise. Results indicate that evaluating this extreme scenario would be valid for resource-management and flood-hazard-mitigation strategies for coastal communities in the Gulf of Mexico, especially those affected by subsidence.

thickness of 229 to 305 m 12 . The sandstone near SG32 is shown in Supplementary Fig. S8. In East Texas and along the Gulf Coast to the Rio Grande River, the Yegua Formation overlies the Cook Mt. Formation, and is overlain by the Caddell Formation. SCnBR is negligible in the over semiconsolidated strata (T OC + C OC ). Minimal groundwater has been developed from this minor aquifer in Texas and no significant groundwater-level decline has been measured by the Texas Water Development Board (TWDB) 13 . Thus, SPC in the strata is negligible. Therefore, land subsidence (LS) measured from GPS station SG32 is attributed solely to LSBR.
Outcropped over semi-consolidated Tertiary strata (Calvert Bluff Formation) at GPS station LDBT: The Calvert Bluff Formation underlying GPS station LDBT 107 km southwest of GPS station SG32 is mostly 14 . The sandstone is medium to fine grained, well sorted, crossbedded, and lenticular. This formation is about 305 m thick. Owing to the over semi-consolidated strata (T OC + C OC ) SCnBR is negligible. Because no groundwater has been developed from the strata in the region, SPC is negligible. Thus, LS measured from the paired reference GPS station LDBT is attributed to LSBR.

Piecewise equation of sea level rise (SLR) at tide gauge Galveston Pier 21
Future GMSLR has been projected using the quadratic equation E global (t)=0.0017t+b g t 2 , where 0.0017 m/a refers to the trend for the GMSLR determined from observations from 1900 to 1992, with an acceleration-related constant b g [m/a 2 ] determined to be 0.0, 2.71×10 -5 , 8.71×10 -5 , and 1.56×10 -4 m/a 2 for the lowest, intermediate-low, intermediate-high, and highest GMSLR scenarios, respectively 15,16 . The highest scenario of GMSLR by 2100 is derived from a combination of estimated ocean warming from the IPCC Fourth Assessment Report (AR4), GMSLR projections and a calculation of the maximum possible glacier and ice sheet loss by the end of the century 15 . The intermediate-high scenario is based on an average of the high end of semi-empirical, GMSLR projections 15 . Semi-empirical projections utilize statistical relations between observed global sea level change, including recent ice sheet loss and air temperature. The intermediate-low scenario is based on the upper end of IPCC AR4 GMSLR projections resulting from climate models using the B1 emissions scenario 15 . The lowest scenario is based on a linear extrapolation of the historical SLR rate derived from tide gauge records beginning in 1900 (1.7 mm/a) 15 . Note, bg = 0.5ag where ag denotes GMSLR acceleration [m/a 2 ] with values of 0.0, 5.42×10 -5 , 1,74×10 -4 , and 3.12×10 -4 m/a 2 for the lowest, intermediate-low, intermediate-high, and highest scenarios of GMSLR, respectively. By comparison, SLR acceleration in the Chesapeake Bay, U.S. was estimated to be about 0.05−0.10 mm/a 2 (5.0×10 -5 − 1×10 -4 m/a 2 ) from the 1930s to 2011 17 . The historical and estimated future magnitudes of RSLR are two additional factors for scientists, engineers and planners to consider when devising sustainable solutions involving adaptive engineering, such as seawalls [18][19][20] , levee systems [21][22][23] , and pile-elevated building foundations 24 for future coastal settlements in the U.S. A set of piecewise equations representing a time-dependent model of RSLR at a tide gauge was developed as follows: RSLR(t)=(a r +s BR )(t-t 0 )+0.4343C H ln(t/t 0 ) +C, t∈(t 0 ,t 1 ] (S1-1) RSLR(t)=(a r +s BR )(t-t 0 )+0.4343C H ln(t t 0 ⁄ ) +p l (t-t 1 )+C, t∈(t 1 ,t 2 ] (S1-2) RSLR(t)=(a r +s BR )(t-t 0 )+0.4343C H ln(t t 0 ⁄ ) +p l (t 2 -t 1 )+C, t∈(t 2 ,t 3 ] (S1-3) RSLR(t)=(a r +s BR )(t-t 0 )+0.4343C H ln(t t 0 ⁄ ) +p l (t 2 -t 1 )+0.5a g (t-t 3 ) 2 +C, t∈( t 3 , 2100] (S1-4) where a r and s BR denote a regional uniform ASLR rate [ 26 to optimize the fit to observed RSLR. Values of t 1 and t 2 were determined by analyzing RSLR data with information of regional and local land subsidence, subsurface fluid withdrawal, groundwater level and subsurface fluid-flow simulation. The value of 1992 for t 3 is from Parris et al. 15 and represents the start of the period of ASLR acceleration which extends to 2100 in this analysis. Note, after 1983, the term p l (t 1 -t 1 ) in equations (S1-3-4) is constant and represents the SPC contribution for subperiods beyond 1983. The supplementary equation set (S1) for a particular tide gauge reflects the fact that RSLR varies along the U.S. coastline 27 . The U.S. Army Corps of Engineers (USACE) 16 has incorporated sea level change into civil works programs using a single equation that is similar to supplementary equation (S1-4) in Regulation No. 1100-2-8162.
Land subsidence from bedrock systems (LSBR) In this study, LSBR is defined as the portion of land subsidence (LS) attributed to TS and SCBR in bedrock systems. LSBR can be measured at a tide gauge's paired reference GPS station that is anchored on bedrock (pre-Cretaceous) or over semi-consolidated Tertiary and Cretaceous strata for which SPC and SCnBR are negligible. The TS is assumed to be caused by comprehensive intraplate or interplate tectonic activities [28][29][30] , which include regional and local faulting, and glacial isostatic adjustment (GIA) 30 . SCBR represents subsidence due to creep of bedrock systems 31,32 . For purposes of this study, the LSBR trend in terms of an annual rate is assumed to be constant for a specific station during the relatively short time period (human time scale) represented by these analyses. Global GPS height data from 2,567 GPS stations on land are available from the Jet Propulsion Laboratory (JPL) 33 Fig. S4C) representing a rate of LS equal to LSBR of 2.67 mm/a. This value was used to represent LSBR at reference GPS station TXGA for tide gauge Galveston Pier 21 as described in the main article.

Regional ASLR estimation
For purposes of this study, the tectonic conditions in a coastal region were assumed to be identical in a small zone with similar strata. For tide gauges and their paired reference GPS stations where SPC and SCnBR are negligible, LSBR was used to estimate ASLR in a coastal region because those tide gauges were assumed to measure height changes attributed solely to ASLR and LSBR. For example, tide gauge Cedar Key and its reference GPS station XCTY (Figs. 1 & 2) 55 km distant are established over the same outcropped over semi-consolidated (see Table 1) Tertiary limestone (T OC ) 7 , where superscript OC denotes an overconsolidated stress status (σ 0 ' <σ c ' ), where σ 0 ' and σ c ' and denote current and historically maximum stresses, for example in Fig. S3) 25 . The NOAA measured trend in relative sea level referred to as RSLR in this paper is 2.13 mm/a (Fig.  Supplementary S1B) 35 , while the NOAA measured LS in trend (LSBR) is 0.88 mm/a (± 0.43 mm/a standard deviation) ( Supplementary Fig. S4A) 34 . At this gauge, RSLR comprises only ASLR and LSBR, as LSnBR (SPC and SCnBR) in the T OC and the Cretaceous strata C OC ( Fig. 2  No significant trends in groundwater-level decline were evident in the available water-level records (1961-94). Therefore, SPC in the Tertiary over semi-consolidated Floridan limestone aquifer system at Cedar Key and Cross City is considered to be negligible, and as discussed in the previous section, SCnBR is also considered to be negligible in these strata.

Subsidence due to primary compaction (SPC)
SPC is the compaction of compressible aquifer systems caused by subsurface fluid withdrawal. SPC in this paper includes the elastic and inelastic (virgin) deformation of aquitards (fine-grained deposits [clays and silts] with low permeability) and the elastic deformation of aquifers (coarsegrained deposits with moderate to high permeability) in aquifer systems. An analytical solution was developed 38,39 to simulate SPC, a coupled compaction and fluid-flow process, in a water saturated, doubly draining clay layer. The clay layer has a uniform initial pore-fluid pressure wherein only vertical water flow is permitted, and identical, instantaneous step changes in hydraulic head (or equivalent fluid pressure) occur at the bottom and top of the clay layer 38,39 . This process, based on changes in effective stress resulting from changes in pore-fluid pressure, describes the equilibration of fluid-pressure and resulting compaction, which was subsequently extended to the analysis 40  ⁄ , where ∆t is the change in time since the initial step change in hydraulic head at the upper and lower boundary of the aquitard, equals 1, 2 and ∞, respectively. SPC is deemed fully completed when T v reaches 2 with 99.4% of the ultimate compaction realized 46 . Thus, a SPC period (∆t) can be theoretically estimated from the above expression for T v with specified values for τ 0 ' and T v . In practice, generally SPC is dominated by inelastic deformation which occurs when the pore-fluid pressure declines result in increased effective stresses in aquitard(s) and confining unit(s) that are greater than the historical maximum effective stress (i.e., σ o ' >σ c ' ), typically defined by the previous minimum pore-fluid pressure in those units. The deformation of the coarse-grained deposits constituting the aquifers generally proceeds only elastically for both decreases and increases in aquifer pore-fluid pressure. The elastic skeletal specific storage of the aquifers governs their deformation and is usually much smaller than S sk ' 47 . Elastic deformation of both aquifers and aquitard/confining units occurs during periods of pore-fluid pressure increase (or groundwater level recovery; for example, see period II in Fig. 4) 47 . For the aquitard drainage model, inelastic compaction is governed by S sk ' which is typically about 1-3 magnitude orders larger than the elastic skeletal specific storage of the aquifers and the aquitards/confining units [48][49][50][51] . Notably, SPC owing to inelastic compaction of aquitards and confining units can result in appreciable land subsidence and thus, increased RSLR. In the Houston-Galveston region it was demonstrated that SPC began in about 1937 and proceeded for at least the next 63 years (∆t) 47 . SPC stopped at many borehole extensometer sites between about 2000 and 2004 (Fig. 4).
Subsidence due to creep of non-bedrock aquifer systems (SCnBR) SCnBR represents the deformation caused by creep behavior of sedimentary materials under a constant load. Due to the weight of the overburden (geostatic stress) and the inelastic compaction characteristics of the aquitards/confining units, about 90 percent of the deformation is permanent 52 . With regard to the degree of self-weight compression, three main sedimentation stages are defined: the clarification regime, zone-settling regime, and compression regime 53 . Quaternary, Tertiary and Cretaceous aquifer systems with a stress condition of σ o ' >σ c ' , remain in the compression regime and experience compaction under self-weight, which has been referred to as creep 25,54 . The path A-B in Fig. S3 shows SCnBR at a constant historical maximum effective stress σ c ' due to the overburden (geostatic) stress. For an unconsolidated/semi-consolidated sediment layer with an initial thickness of H [L], SCnBR can be approximated using S c (t)=C α Hlog(t t 0 ⁄ ), which is employed in supplementary equations (S1-1-4) for the variable SCnBR rate, where t 0 denotes an initial reference time for compression of creep, and t signifies time larger than or equal to t 0 25 . Taking the derivative with respect to time t of the expression above for S c gives Ṡc= C α H (t ln 10) ⁄ , the subsidence rate. The decreased percentage (D S ) of Ṡc from t to t+∆t was derived using [Ṡc(t)-Ṡc(t+∆t)]/Ṡc(t) as follows 47 : For t≫∆t, where ∆t can be a short observing period such as 10-20 years, D S approaches zero which implies that Ṡc is approximately a constant. In other words, the changing value of Ṡc over the ∆t period can be ignored. For example, if a current observation period (∆t) is considered to be 10 years, 990, 1990, and 9990 years of SCnBR are needed to achieve decreased percentages of 1.0, 0.5 and 0.1% of the specified subsidence rate, respectively 47 . This negligibly-variable SCnBR rate 47 is used to estimate SPC.

Negligible SCnBR in over semi-consolidated Tertiary and Cretaceous strata
SCnBR is assumed to be negligible in Tertiary (T OC ) and Cretaceous (C OC ) strata shown in Figs. 1, 2 and Table 1. This is because the current effective stress σ ′ in T OC or C OC is less than its historical maximum effective stress σ c ' (Supplementary Fig. S3) owing to the outcropping or uplifting of the strata and the geological removal of the overlying Quaternary sediments Qh and Qp, which is referred as an overconsolidated condition 25 . Experimental investigations show that sediment creep rate with an overconsolidated stress status is significantly slower than that without 55 . In addition, the primary factor is that when compared to the age of geological strata our human observation is in such a short time period that the change in SCnBR is insignificant, which has also been discussed after supplementary equation S2. For instance, the length C-D in Supplementary Fig. S3