Novel Quantification of Shallow Sediment Compaction by GPS Interferometric Reflectometry and Implications for Flood Susceptibility

Estimates of flood susceptibility and land loss in the world's coastal regions depend on our knowledge of sea level rise (SLR) from increases in ocean mass and volume, as well as knowledge of vertical land motion. Conventional approaches to the latter include tide‐gauge and Global Positioning System (GPS) measurements relative to well‐anchored monuments few meters below the surface. However, in regions of rapid Holocene sedimentation, compaction of this material can add a significant component to the surface lowering. Unfortunately, this process has been difficult to quantify, especially for the shallowest material above the monument. Here we use a new technique, GPS interferometric reflectometry, to estimate the rate of this process in the Mississippi Delta and the eastern margin of the North Sea. We show that the rate of shallow compaction is comparable to or larger than the rate of global SLR, adding 35% and 65%, respectively, to the rate of relative SLR by 2100.


Introduction
The Global Positioning System (GPS) allows measurement of Earth surface displacements with submillimeter accuracy. This has improved our understanding of sea level rise (SLR) and coastal flood susceptibility by correcting relative sea level change obtained from tide gauges (measured relative to the nearby land surface) for the effects of nearby vertical land motion (VLM) (Becker et al., 2020;Chen et al., 2017;Karegar et al., 2017;Wöppelmann et al., 2009;Wöppelmann & Marcos, 2016). Geodetic GPS receivers as well as most tide gauge instruments are installed on the top of buildings or mounted to concrete pillars and masts driven to refusal, typically more than 1-2 m depth. Thus, GPS receivers and tide gauge instruments measure VLM that occurs beneath the base of building foundation or monument platform, ignoring deformation above that depth. This is particularity important for tide gauges since they measure sea level change with respect to the base of their building or monument, rather than the sea bottom. Long-term VLM and assessment of future flood risk and land loss may therefore be underestimated (e.g., Keogh & Törnqvist, 2019). The current paradigm is that GPS and tide gauges sense the same VLM (Bevis et al., 2002;Bouin & Wöppelmann, 2010;Hamlington et al., 2016;Santamaría-Gómez et al., 2012;Schöne et al., 2009;Wöppelmann et al., 2007), in part because shallow VLM is difficult to measure and has been thought to be negligible. Some recent studies have suggested that shallow compaction could be important in certain coastal areas, but quantitative approaches have been challenging (Jankowski et al., 2017;Keogh & Törnqvist, 2019). Here we suggest a new approach to this problem based on GPS interferometric reflectometry (GPS-IR) (Larson et al., 2009). We demonstrate the technique using data from the Mississippi Delta and the east coast of North Sea, some of the world's thickest Holocene age sediment deposits and hence experiencing high rates of relative SLR.

Geologic Background
Low-lying coastal areas are susceptible to accelerating rates of SLR, especially if rates of land subsidence are high. Flooding, wetland loss, saltwater intrusion, shoreline erosion and related economic losses are the main negative consequences (Alam, 1996;Ingebritsen & Galloway, 2014;Milliman & Haq, 1996;Shirzaei & Bürgmann, 2018). Natural compaction of Holocene age deposits, often amplified by peat soil oxidation due to surface water drainage, is one of the main processes that causes coastal subsidence. Subsidence rates can exceed 10-20 mm/yr depending on thickness, age and characteristics of the sedimentary sequence (Brain, 2016;Dixon et al., 2006;Törnqvist et al., 2008;Van Asselen et al., 2011). Sediment compaction occurs over a range of depths. Compaction of the upper 5-10 m of the Holocene substrate, here called shallow subsidence, contributes significantly to the total subsidence of coastal marsh soils (Brain et al., 2012;Long et al., 2006;Törnqvist et al., 2008). Models and data suggest that the rate of shallow subsidence can be greater than current rate of regional and global SLR. This is especially true in coastal plains and river deltas, where close to a billion people live (Gebremichael et al., 2018;Teatini et al., 2011). For example, the average shallow subsidence in the uppermost 5 m of Holocene strata in the Mississippi Delta has been estimated to be 6.4 ± 5.4 mm/yr (Jankowski et al., 2017). Despite its importance, measurements of present-day shallow subsidence are difficult to make, relying mostly on rod surface elevation table marker horizon stations (Jankowski et al., 2017;Webb et al., 2013).

Technique
Geodetic GPS antennas receive two kinds of signals: strong direct signals from the GPS satellites and weak reflected signals from the surrounding environment. Conventional carrier phase positioning uses the direct signals to determine three-dimensional estimates of antenna position. Here, we use the precise point positioning algorithm (Zumberge et al., 1997) implemented in GipsyX software to produce daily time series of antenna positions relative to the International GNSS Service frame 2014 (IGS14)  (see supporting information Text S1). The vertical component of position change is attributed to displacement that occurs beneath the structural foundation of the GPS monument and reflects deep VLM ( Figure 1a). The reflected signals travel a longer path than the direct signals, reaching the antenna later and interfering with the direct signals. The interference signals are recorded as the receiver-generated signal-to-noise ratio (SNR). The power spectral densities of detrended SNR data include peaks with corresponding frequencies, which are linearly dependent on the reflector height (Larson et al., 2009) (see supporting information Text S2). The GPS-IR method has previously been used to derive snow depth (Larson et al., 2009;Siegfried et al., 2017) and sea level changes (Larson et al., 2013). Here, we use the GPS-IR method for the first time to measure height changes attributed to ground surface changes related to shallow displacements that occur within the shallow layer between the surface and the base of the GPS monument ( Figure 1a).
Resolving shallow displacements for sites with significant local vegetation is challenging because apparent reflector height changes may be affected by changes in vegetation height associated with seasonal fluctuations and/or variable rainfall (Small et al., 2010). It is important to apply masks on elevation angles and azimuths to avoid reflections from vegetation while maintaining an adequate number of observables. We estimated linear trends in shallow and deep VLM time series using the Median Interannual Difference Adjusted for Skewness algorithm (MIDAS), rate uncertainties following Allan Variance of the Rate method (AVR), and tested the statistical significance of rates with a nonparametric modified Mann-Kendall's test (see Text S3).
Interferometric Synthetic Aperture Radar (InSAR) is another potential technique for measuring shallow sediment compaction and can provide important spatial detail. It is most successful in urban and suburban areas where strong radar scatterers (e.g., buildings) exist, obviating the problem of vegetation, which can reduce phase coherence (Bekaert et al., 2017;Dixon et al., 2006;Higgins, 2016;Shirzaei & Bürgmann, 2018;Teatini et al., 2011). However, the best radar scatterers in these areas tend to be buildings anchored at depth; hence, these may underestimate shallow compaction. It is also important to recall that InSAR is a relative technique, requiring the use of one or more reference points, such as the GPS technique described here.

Thickness of Holocene Sediments
Many coastal plains host river deltas that formed in the late to mid Holocene (Stanley & Warne, 1994) above Pleistocene units that may themselves host poorly compressed sediment layers. Three-dimensional GIS data sets of regional Holocene geology have recently been compiled to create comprehensive map of the Holocene-Pleistocene surface of the Louisiana coastal plain (Heinrich et al., 2015)   Louisiana is composed of a very thick and high-porosity Holocene sediment wedge, tied to the balance between rates of sea level change, land subsidence, and sediment supply and deposition. Holocene deposition begins north of Lake Pontchartrain at latitude 30.5°N, increasing seaward from a thin veneer of~1 m to a layer exceeding 100 m thickness at the coastline, and >150 m at the shelf margin ( Figure 1b).
The Eastern margin of the North Sea has complex sedimentological and morphological processes, influenced by tidal mechanisms of the North Sea, rivers estuaries, and a sandy coastal barrier, bounded by the cliff coast of northern France in the south and the tip of Denmark in the north. Holocene sediments span a 50 km wide strip along the coast with westward and northward thickening. Holocene deposits can exceed 30 m in the Rhine Delta, central Netherlands coast and along the Wadden Sea barrier islands ( Figure 1c). The Holocene sequence in the Rhine Delta is composed of river fluvial and tidal basin deposits, thickening westward to the coastline. The sediment wedge on the central coast of Netherlands deviates from the Rhine Delta by sea ingression and peat brook (Beets & van der Spek, 2000). The Wadden Sea is one of largest uninterrupted natural coastal barrier systems in the world, consisting of barrier islands, tidal plains, and salt marshes, with most sediments imported from the North Sea. The thickness of these Holocene wedges reaches to 50 m along the barrier islands in the Dutch Wadden Sea (Figure 1c), diminishing northward along the German and Danish coasts.

Results
We derived time series of shallow VLM at six GPS sites in the Mississippi Delta and six stations along the eastern margin of the North Sea. Azimuth and elevation angle masks were imposed to isolate ground reflections with minimal vegetation cover from additional sources of reflection (see supporting information Figure S1 and Table S1). The GPS sites in Grand Isle (GRIS) and Boothville-Venice (BVHS) along the coast of Louisiana are anchored within the thickest unconsolidated Holocene sediments in the Mississippi Delta (~60 m thickness) at a depth of about 20 m from the surface (Figure 1b). The GPS site in Terschelling Island (TERS) along the coast of Northern Netherlands is anchored within one of the thickest Holocene sequences in Europe (31 m thickness) at a depth of about 15 m from the surface (Figure 1c). Rate analysis of shallow VLM indicates that highest subsidence rates occur at sites where the Holocene sediments are substantial. We observe a shallow rate of −5.6 ± 0.6 mm/yr at GRIS and −4.6 ± 2.7 mm/yr at BVHS, reflecting shallow compaction at depths above 20 m. We observe a corresponding rate of −4.1 ± 1.8 mm/yr at TERS, reflecting shallow compaction at depths above 15 m. Stations with very shallow foundation depths (smaller than 1-2 m) show no significant shallow compaction (FSHS, MSIN, and TGBF) ( Figure 2).
One way to validate the shallow rates derived from GPS-IR is to consider stations underlain directly by stable Pleistocene sediment; another is to see if the shallow subsidence rates are consistent with independently derived values from the Rod Surface-Elevation Table-Marker Horizon method (RSET-MH) (Cahoon et al., 2002;Webb et al., 2013) (see supporting information Text S4). Station ALNB in Alabama and station ESBC in Denmark are anchored in consolidated Pleistocene age sediment at depth of 2 and 10 m from the surface, respectively, and are not affected by shallow compaction (Figure 2a). We compared our GPS-derived shallow rates with the rates measured from RSET-MH technique in the Mississippi Delta (Figures 3a, 3b, and S3). A RSET-MH site consists of vertical stainless steel rod attached to a benchmark hammered typically 10-25 m into the substrate, which quantifies surface elevation change with respect to the bottom of benchmark, and an artificial layer of white sand on the ground surface (a marker horizon) which measures vertical sediment accretion. The RSET-MH thus provides net surface elevation change above the base of benchmark, which is attributed to shallow subsurface processes (Jankowski et al., 2017). Shallow rates from the RSET-MH show qualitatively good agreement within uncertainty to our estimates of shallow rates from GPS-IR.
For GPS stations anchored in Holocene age sediments, the deep VLM rates reflect compaction resulting from both lower Holocene and pre-Holocene sediments ( Figure S10), as well as deeper processes including glacial isostatic forebulge collapse (Wickert et al., 2019), Holocene sedimentary isostatic adjustment (Kuchar et al., 2018), growth fault movement (Shen et al., 2017), and anthropogenic movements due to fluid withdrawal (Dokka, 2011). The deep subsidence rate is highest on the coast of the Mississippi Delta in south Louisiana: 6.3 ± 0.4 mm/yr (GRIS), 6.8 ± 0.4 mm/yr (LMCN), and 4.6 ± 0.6 mm/yr (BVHS). The highest deep subsidence rate in the European study area is 1.0 ± 0.5 mm/yr (TERS). The deep subsidence rates are significantly higher in the Mississippi Delta compared to the Rhine Delta and the east coast of the 10.1029/2020GL087807

Geophysical Research Letters
KAREGAR ET AL.  Uncertainties in the RSET-MH rates are estimated through the linear regression and are scaled by factor 3 to approximate time-correlated noise (e.g., Mao et al., 1999). The GPS rates without error bar are predicated shallow rate based on linear regression model of the correlation between total VLM rates and thickness of Holocene age sediments at GPS sites (Figure 3d). (c) Relationship between deep VLM rates and Holocene sediment thickness. (d) Total VLM rates (sum of shallow and deep VLM rates) against the Holocene deposits thickness. (e) Shallow VLM rates versus the depth of the GPS antenna. Error bars are 1 σ. The gray shaded area is 95% confidence interval of the regression line estimated based on 100,000 realizations of the bootstrap method accounting for uncertainty in VLM rates (see supporting information Text S5). The red line is a fitting line based on mean of bootstrapped regression parameters. Inset are probability density functions for Pearson moment-product correlation coefficients estimated from bootstrap results. The red vertical dashed lines indicate 95% bootstrapped confidence interval for correlation coefficient.
North Sea (Figure 2), mainly reflecting the greater thickness of Holocene sediment in the former, and different GPS foundation depths and proximity to glacial isostatic forebulge collapse. The ALNB site in the Mississippi Delta and the LETO, VLIS, ESBC, and most likely TGDA sites on the eastern margin of the North Sea are anchored on the stable Pleistocene surface and experience relatively slow rates of deep subsidence (less than 1.0 mm/yr) (Figure 2). The deep subsidence rates on these Pleistocene-anchored sites agree with three recent glacial isostatic adjustment (GIA) models (see supporting information Table S2), which would imply that GIA forebulge collapse is the dominant process contributing to deep subsidence at these sites. However, station MARY in East New Orleans records uplift presumably resulting from deep groundwater recharge that began in 2016 ( Figure S4 and below).
VLM in the Mississippi Delta near metropolitan New Orleans, Louisiana, is affected by groundwater withdrawal and human intervention (Dokka, 2011). Station MARY is located at the NASA Michoud facility in Eastern New Orleans where the highest subsidence rates were reported in previous studies . This site is grounded at a depth of 2000 m below the surface and is known as one of the deepest anchored GPS sites in the world. Using conventional GPS positioning, we estimate displacements that occur deeper than 2,000 m. The long-term deep vertical rate is +2.2 ± 0.7 mm/yr and represents the positive poroelastic response of ground surface to the ongoing aquifer recharge that began in 2016 ( Figure S4). The approximate 2.5 m antenna height used at station MARY allows sensing of the first Fresnel zone (ellipse located along the satellite ground track) with maximum dimensions of 4 m by 60 m. However, the vegetation at this site is intense and is classified as short-grass steppe. Thus, the reflector height changes include the compound effects of ground surface motion and vegetation growth ( Figure S5). We applied a strict mask on elevation angles and azimuth to limit the reflection data to regions with minimum vegetation cover ( Figures S5 and S6). We observe a shallow subsidence rate of 2.5 ± 1.3 mm/yr, presumably mainly driven by compaction of 16 m Holocene sediment at this site.

Discussion and Conclusions
The shallow subsidence rate derived from the GPS-IR technique depends on two factors: (1) the thickness and compressibility of Holocene sediments above the base of benchmark and (2) the depth of GPS monument. The subsidence rates in coastal plains generally correlate with the thickness of Holocene sediments (Edwards, 2006;Horton et al., 2013;Jankowski et al., 2017;Karegar et al., 2015;Mazzotti et al., 2009) . To discern the impacts of these factors, we used a nonparametric bootstrap statistical correlation analysis, taking into account uncertainty of GPS rates (see supporting information Text S5). Station MARY was excluded from this analysis due to effects from groundwater recharge. The deep subsidence rates derived from GPS show moderate correlation with sediment thickness with a broad (95%) confidence interval (Figure 3c). GPS stations anchored in Holocene sediment record the contribution of compaction occurring in sediment below their anchoring depths, as well as contributions from deeper processes. Accounting for both shallow and deep VLMs (total VLM), we find a strong linear relationship between the rate of total VLM and sediment thickness with a narrow confidence interval (Figure 3d). The higher dependence of total VLM rates on Holocene thickness suggests that (i) shallow VLM as derived by GPS-IR is mainly caused by sediment compaction; (ii) conventional GPS positioning underestimates the VLM since it misses shallow sediment compaction; and (iii) regions with thicker and younger sediments experience faster subsidence. The trend of higher total subsidence rates with thicker Holocene sediment sections appears to be true from both alluvial and fluvial fans, although our data are limited.
The shallow subsidence rates from GPS-IR are inversely correlated with foundation depth of the GPS stations (where known) (Figure 3d). Stations anchored at greater depths show faster shallow subsidence rates compared to stations anchored at shallower depth, due to the greater volume of compressible sediments above the anchor ( Figure S10). However, the confidence bound is rather large because some stations (e.g., ALNB, MSIN, and ESBC) are grounded on consolidated Pleistocene age sediments with minimal compaction. These findings also support the idea that present-day subsidence rates in the Mississippi Delta and along the Eastern margin of the North Sea are mainly caused by compaction of Holocene strata. It should be mentioned that the GPS anchoring depth is not always recorded in the metadata since many of these sites are primarily installed for surveying applications. We recommend that foundation depth be recorded for future GPS installations.
Given the empirical relation between the total rate of VLM and Holocene sediment thickness (Figure 3d), we can predict the value of total VLM at three sites in the Mississippi Delta where GPS-IR was not applicable due to environmental considerations. Conventional positioning allows a precise estimate of deep VLM. Thus, we subtracted the measured deep rates of VLM from predicted total VLM to estimate shallow rates at these sites. The predicted shallow rates can then be compared with independent measurements from the closest RSET-MH sites (grayed stations in Figure 3b). For stations near Cocodrie (LMCN, 0369) and in the Shell Beach area (SBCH, 4548), the GPS and RSET-MH instruments are anchored at approximately the same depth from the ground surface (Tables S3 and S4). The predicted shallow rates from our empirical model agree well with measured rates from RSET-MH within their uncertainty. However, RSET-MH data show slightly greater shallow rates for sites close to the ENG5 GPS station south of New Orleans. The GPS foundation at ENG5 station is 3 m deep whereas RSET-MH sites (3641 and 3664) are anchored at a depth of about 25 m, thereby recoding a higher compaction rate. Unfortunately, it is not possible to evaluate this empirical model for stations along the North Sea, because independent measurements of shallow rates are not available at reasonably close distances to GPS sites. Nevertheless, our empirical analysis suggests that knowledge of Holocene sediment thickness and the deep VLM rate from conventional GPS positioning is sufficient to determine the rate of shallow subsidence in areas where GPS-IR is not applicable due to ground cover near the antenna.
The current deep-seated subsidence rate of Mississippi Delta varies from 4.6 ± 0.6 to 6.8 ± 0.4 mm/yr based on three southern delta GPS stations. The rates of VLM calculated by GPS-IR here suggest signification shallow subsidence rates, up to 5.6 mm/yr in this region, nearly 3 times faster than local absolute SLR (with respect to the center of mass of the Earth) recorded by tide gauges since 1924 (Karegar et al., 2015) and 3-4 times greater than twentieth century rate of global mean SLR (1.1 to 2 mm/yr; Dangendorf et al., 2017). Total land subsidence rates (sum of shallow and deep subsidence rates; Table S4) here thus range from 9.2 ± 2.8 to 11.9 ± 0.7 mm/yr, implying rates of relative SLR of up to 15 mm/yr. For the eastern coast of the North Sea and the Rhine Delta, the deep subsidence rate is smaller than 1 mm/yr but the shallow subsidence rate is as high as 4.1 mm/yr, varying locally with sediment thickness, again exacerbating rates of relative SLR.
The contribution of VLM to sea level projection is generally limited to GIA models (Grinsted et al., 2015;IPCC, 2014;Jackson & Jevrejeva, 2016;Slangen et al., 2014) or is locally accounted for by using tide gauges records (Kopp et al., 2014). Shallow VLM has not typically been incorporated into sea level projection due to lack of direct observations. Our results have important implications for sea level projection and estimation of future flood susceptibility and land loss in deltas and other coastal regions underlain by thick sequences of Holocene sediment. Probabilistic projections of relative SLR between the years 2000 and 2100 predict a very likely (5th-95th percentiles) rise of 1.1-2.0 m in GRIS (Southern Louisiana), under representative concentration pathway (RCP) 8.5 global warming scenario (Kopp et al., 2014). The shallow subsidence in GRIS is comparable to expected rate of SLR resulting from individual contributing processes (Kopp et al., 2014;Jackson & Jevrejeva, 2016), adding 0.56 m (median, 50% probability) rise by 2100. We estimate a very likely range of 1.6 and 2.7 m of relative SLR by 2100 under RCP 8.5 scenario ( Figure S9). Station TERS (Northern Netherlands) is expected to experience a 0.66 m (medium) SLR by 2100, with a very likely range of 0.29 to 1.12 m. The shallow subsidence at TERS results in the highest SLR projection along the margin of the North Sea, exceeding 1.0 m (medium) by 2100. The lower and upper bounds of very likely ranges of RCP 8.5 projected SLR is expected to be 0.4 and 1.92 m, respectively. In areas of increased sediment deposition, such as Grand Isle and Terschelling Island, rapid natural shallow subsidence adds 35% and 65% (95th percentile) respectively to the rate of relative SLR by 2100 ( Figure S9). While these findings imply that shallow subsidence is an important contributor to sea level projection, we note that these estimates are conservative in terms of flood hazard in the sense that we do not consider short time scale water-level variability due to storm events, tidal flooding, and precipitation.

Perspective and Outlook
Our findings represent the first attempt to examine the GPS-IR technique to quantify shallow subsidence rates in actively subsiding coastal plains, which have important implications for the interpretation of tide gauge records in terms of long-term relative SLR and flood hazard. We suggest that the GPS-IR method should be used to quantify shallow subsidence and correct the rate of relative SLR at colocated tide gauges. The number of GPS sites that can be used to retrieve the reflector height is currently limited by two main 10.1029/2020GL087807 Geophysical Research Letters factors. First, since the primary aim of existing geodetic GPS networks is accurate positioning for geophysical studies and/or survey engineering based on pseudorange and carrier phase measurements, the SNR observable are not always archived in metadata associated with sites. We recommend recording SNR data for GPS-IR applications. Second, the GPS antenna must be located close to a fairly planar natural surface (preferably bare ground), so that with minimum imposed azimuth and elevation angle masks, desirable reflections from the area of interest are obtained.
In addition to correcting the rate of relative SLR, direct measurement of shallow VLM can constrain sediment compaction models on the evolving state of accumulating sediments, potentially improving the reliability of compaction rate prediction Meckel et al., 2006). Furthermore, time series of shallow displacement can be used to assess whether natural subsidence of deltaic sediments and SLR is compensated by sediment accumulation. For example, construction of artificial levees, canals and dams reduces river flooding and sediment flow onto many deltaic plains (Syvitski et al., 2009). The GPS-IR technique used here for measuring displacement above the base of monument or building foundations has potential application for engineering purposes such as monitoring foundation stability.