Tunneling- and dewatering-induced rapid differential ground rebound and delayed subsidence measured by InSAR in an urban environment

Abstract. During the excavation of the Alaskan Way Viaduct replacement tunnel in Seattle, Washington, a 17.5 m diameter tunnel boring machine (TBM) nicknamed “Big Bertha” was damaged after encountering unexpected subsurface conditions. Significant dewatering of multiple aquifers was required to reach the TBM for repairs. Groundwater drawdown and soil consolidation associated with dewatering created a 0.4  km2 region of initial subsidence with maximum vertical settlements exceeding 2.5 cm between August and December 2014. Dewatering wells remained operational until January 2016 and likely contributed to observed groundwater drawdown in areas outside the region of initial subsidence. To determine how an urban landscape with complex and poorly constrained geologic and hydrologic conditions responds to an extended period of dewatering within multiple aquifers, the rate, duration, spatial extent, and magnitude of dewatering-related displacements were analyzed by combining three paths of Sentinel-1 interferometric synthetic aperture radar data spanning November 2014 to October 2019 into a time series of vertical surface deformation using the minimum acceleration algorithm. Our results show that post-dewatering ground rebound within this complex hydrogeologic system occurred at faster rates and with more significant spatial deformation variability than initial subsidence, reaching rates of up to 17  cm/year coupled with potentially hazardous differential rebounds across short distances. In addition, prolonged groundwater pumping at depths greater than 60 m appears to have induced delayed subsidence over a larger area of ∼20  km2, reaching magnitudes of up to 3 cm and lasting for over 3 years after the cessation of pumping.


Introduction
Subsurface excavations in an urban environment often encounter saturated soils, especially in cities along coastal areas with shallow groundwater tables.When groundwater cannot feasibly be diverted by subsurface physical barriers (e.g., steel sheet piling or artificial ground freezing), underground construction projects commonly use dewatering wells to lower the groundwater table temporarily.Groundwater pumping reduces pore pressure and increases effective stress within subsurface strata, subsequently inducing subsidence in compressible geology. 1 Dewatering creates artificial and often highly dynamic groundwater systems, 2 especially in naturally complex geologic and hydrologic conditions.In urban settings, tunneling-induced subsidence also poses a significant hazard to the overlying infrastructure, especially when it occurs non-uniformly over short distances, resulting in differential settlement.Tunnelinginduced ground subsidence has been studied by researchers using various methods, such as in situ monitoring, analytical or numerical modeling, interferometric synthetic radar (InSAR) monitoring, and subsidence prediction by machine learning, e.g., Refs.3-8.Potential hazards induced by post-dewatering rebound are less often considered and rarely have they been studied with detailed spatial geodetic data in an urban setting.
Groundwater monitoring techniques generally rely on piezometer networks to measure spatial and temporal changes in groundwater levels.However, piezometer-based monitoring networks only provide sparsely distributed point information, which can be insufficient for characterizing the spatiotemporal changes in the groundwater regime accurately in complicated geological and hydrological subsurface settings.9][20] InSAR techniques have successfully been applied to measure deformation induced by tunneling activities in both hard rock [21][22][23] and soft rock or sediment rich environments. 6,24,25Studies have shown that groundwater recharge can cause surface uplift, the duration of which depends on factors such as geology, degree of aquifer depletion, and groundwater pumping rates, 9,15,26,27 though these studies often rely on measurements from individual interferograms and have not investigated tunneling related rebound in urban settings.
To understand how drained soils spatially respond to recharging in an urban environment with complex geology, surface deformation in downtown Seattle (Fig. 1) was estimated over  20 and locations of access shafts in which dewatering wells were installed, nearby GPS station, borehole used to characterize subsurface geology, and piezometer network used to monitor groundwater levels.
a 5-year period from 2014 to 2019, during which a network of dewatering wells at multiple depths was sequentially initiated and decommissioned.This quantitative assessment of postdewatering rebound identifies potential hazards associated with a short-term rapid ground rebound and shows the complex, long-term surficial response of multiple deforming aquifers within an artificially perturbed groundwater regime.

Geology and Hydrogeology
The greater Seattle area is built on successive glacial, fluvial, and lacustrine deposits truncated by unconformities and deformed through faulting and folding due to a complicated geomorphic and tectonic history. 28The geologic units include Holocene nonglacial deposits, Pleistocene younger glacial deposits, and older and nonglacial deposits.Postglacial bog and recessional lake deposits frequently contain highly compressible peat and organic-rich salt deposits in locations throughout Seattle. 20,28,29Many peat deposits in the area are now covered by fill and occur as discontinuous lenses, in which case mapping is not plausible.Subsurface materials show significant variability in spatial extent and thickness, leading to a complicated groundwater regime.Figure 2 illustrates the schematic geologic cross-section along the tunnel alignment from Blanchard Street (north) to South Massachusetts Street (south) in Downtown Seattle.
At least two aquifers are known to exist beneath Pioneer Square and the surrounding areas in downtown Seattle, including an unconfined aquifer within surficial deposits and a confined aquifer composed primarily of sand and gravel. 20,30,31These aquifers are separated by an aquitard layer composed of cohesive clay and silt that can be found at multiple depths from 3 to 65 m below the surface that range from 1 to 30 m in thickness. 30Data from the deepest borehole available within the study area (Table 1) provide a generalization of the aquifer system.The maximum depths of the confined aquifer are poorly constrained.Regional groundwater studies indicate a confining depth of at least 245 m, whereas the thickness of unconsolidated deposits composing the Puget Sound aquifer system could reach up to 730 m or more in the downtown Seattle area. 29,312 Dewatering Timeline and Initial Subsidence A 37 m deep access shaft was constructed from October 2014 to January 2015, to reach the damaged tunnel boring machine (TBM) (Fig. 1).A series of dewatering wells were installed to lower the groundwater table to excavate the access shaft and gain access to the TBM.20,30,32,33 A total of 15 dewatering wells were installed in and immediately around the access shaft and screened at multiple depths to 61 m, dewatering both the unconfined and confined aquifers.Four internal access shaft wells screened to 46 m began collectively pumping 55 to 110 m 3 ∕day on October 2, 2014.By October 29, seven other shallow wells installed around the outside the shaft within a series of Settlement Mitigation Piles (SESMP) were pumping a total of 320 to 650 m 3 ∕day.Four deep wells (61 m depth) began pumping on November 8 at a combined rate of 3500 to 4100 m 3 ∕day.30,34 Within 2 months of the initiation of the deep well pumping, piezometric elevation recorded at a 115 m deep steam plant well located ∼850 m northwest of the access shaft (Fig. 1) fell by over 7.5 m. 30 The wells within the access shaft were shut off by February 19, 2015, whereas the remaining shallow wells and all deep wells continued pumping until they were sequentially shut off between January 22 and January 24, 2016.35 Dewatering and subsequent consolidation were shown to directly cause the formation of an elliptical subsidence feature in downtown Seattle, reaching ∼0.4 km 2 in size.Previously published RADARSAT-2 spotlight data processed via the multidimensional small baseline subset (MSBAS) technique 36 measured vertical deformation reaching a maximum of 3.5 cm between August 2014 and August 2015, with maximum subsidence rates of 10 cm∕year occurring between August 31 and December 29, 2014.20,37 Damaged utilities as a result of this event include a 400 m stretch of 20-in.diameter water line piping beginning 200 m north of the access shaft location.34

Methods
Measurements of ground deformation associated with dewatering near the access shaft were acquired by processing one ascending path (137) and two descending paths (115 and 13) of Sentinel-1 data.All available acquisitions before October 19, 2019, were used (Table 2).9][40] GAMMA software 41 was used to create a series of single master interferograms from dual vertically (VV) polarized single look complexes (SLCs).A single master date was chosen for each data path.To reduce temporal decorrelation, the master scene dates near the middle of the InSAR time series were selected.Wrapped Wnuk, Zhou, and Gutierrez: Tunneling-and dewatering-induced rapid. . .interferograms were created using the scenes directly before and after each master scene date to visually inspect the pair-wise interferograms for turbulent atmospheric effects.After master scenes were selected and verified, slave images were co-registered to the geometry of the chosen master scenes, with topographic phase contribution removed via the 3D elevation program (3DEP) 1/9 arc-second resolution digital elevation model (DEM). 42AMMA's interferometric point target analysis was used to create a time series of persistent scatterer (PS) pixels.PS pixels (Fig. 3) were identified using the mean to standard-deviation ratio (MSR) of the radar intensity backscatter.A threshold MSR value of 1.1 was used to select pixels with low temporal variability and high-intensity backscatters that provide stable phase returns over the period of interest.SLC values at PS pixel locations were extracted and written to a point data stack.A phase reference point, to which all measured phase changes are relative, was then selected for all three paths of data (Fig. 3).The reference area was visually selected to be on top of a large, flat-roofed building, contain PS pixels from all three data paths, and located outside the initial subsidence patch.Initial differential interferograms were formed and then analyzed in the temporal domain to measure the linear dependence between the perpendicular baseline and phase for each PS pixel, which was used to calculate DEM corrections.These corrections were subsequently used to update perpendicular baseline values.Linear deformation rates were then calculated for each pixel.The residual phase, considered to be a deviation from the linear rate, remains present in the data and is composed of atmospheric contamination from secondary scenes, baseline measurement errors, non-linear deformation, and random noise.Because the atmospheric phase contribution from the master scene correlates with space and time, it was extracted using a low-pass temporal filter of the residual phase and then removed from all interferograms.The process was then iterated using the corrected DEM, updated perpendicular baselines, and atmosphere-corrected interferograms.Final deformation values were produced via a phase unwrapping algorithm based on minimum cost flow techniques applied to a Delaunay triangular network.
After LOS deformation time series were constructed for all input datasets, values of PS pixels (Fig. 3) were interpolated to a common 10 m gridded surface to reduce noise from erroneously selected PS pixels and make the data more accessible for post-processing.PS time-series were combined in a post-processing stage using the minimum acceleration (MinA) algorithm 43 to estimate vertical and east-west horizontal deformation from multiple LOS time series. 44,45Input LOS displacement time series were connected to the unknown velocities among consecutive time acquisitions, creating a system of underdetermined independent linear equations: where d is the deformation velocity, ϑ is the radar incidence angle, and φ is the LOS azimuth of the n LOS displacement time series.Many InSAR sensors, including Sentinel-1, are side-looking systems onboard satellites traveling along polar orbits, preventing LOS measurements outside a primarily E-W orientation that obscures the horizontal deformation in the N-S direction. 46herefore, N-S deformation velocities were assumed to be zero when applying Eq. ( 1) for the purposes of this study.Equations were solved in a least squares sense by applying first-order Tikhonov regularization in combination with the L-curve method. 47,48Regularization was enforced by applying a range of regularization coefficients to the input series of linear equations, creating a series of potential solutions ranging from solutions that fit the data without minimizing acceleration to solutions that minimize acceleration but not the residual error (Fig. 4).All potential solutions plot along an L shaped curve, with the optimum solution located at the bottom left corner.The solution found at the corner of the L-curve is optimized because it provides the most precisely balanced solution between minimizing the solution error and minimizing the displacement acceleration.The regularization portion of the analysis was performed using the regularization tools MATLAB package. 49Regularized displacement velocities were then independently time-integrated to extract a single deformation time series.
To assess the algorithm accuracy, results from the Sentinel-1 MinA algorithm output were compared with previously published, temporally overlapping RADARSAT-2 vertical displacement data calculated using the MSBAS algorithm. 36,37As the MSBAS approach consists of generating and simultaneously processing a full sequence of multiple-platform differential SAR interferograms, the primary difference between the algorithms lies in the fact that the MinA approach is a post-processing step in which completed LOS time series are used as inputs.MinA output displacements were additionally compared to available GPS data fixed to the North American plate, processed by and downloaded from the Nevada Geodetic Laboratory 50 for further validation.

Continued Subsidence and Rebound
The first available Sentinel-1 image for the area of interest was acquired on November 6, 2014, 35 days after the initiation of dewatering in the shallower wells but before the deep wells began pumping on November 8.During the period in which all dewatering wells were active, widespread subsidence was observed throughout the pioneer square area [Fig.5(a)].In February 2015, subsidence rates decreased coincidentally with the shutoff of all wells within the access shaft on February 19 (Fig. 6).Vertical deformation during this time period shows continued subsidence with scattered patches of uplift [Fig.5(b)].After the cessation of groundwater pumping by January 24, 2016, an immediate widespread ground rebound was observed in the area immediately around the access shaft [Figs.5(c)-7] in temporal coincidence with observed groundwater recharge (Fig. 8).A maximum total rebound of nearly 3 cm was observed at point PD between January 12 and August 5, 2016, with a maximum uplift rate of 17 cm/year sustained from January 29 to March 14.Rebound varied greatly over short distances during this period, with differential rebound magnitudes exceeding 2 cm between points PA and PB in the months after the deep wells stopped operating (Fig. 6).
The mapped surface rebound shows that areas of greatest rebound are not found in areas of greatest initial subsidence and that most subsided areas experienced some degree of permanent deformation (Fig. 7).This is especially apparent at point PA, which lies in a region that initially subsided more than 3 cm but experienced little to no rebound (Fig. 6).This non-uniform surface rebound created differential deformation of magnitudes greater than 2 cm across distances shorter than 100 m in some locations for a period of time.

Delayed Subsidence
In the months after the deepest wells began pumping over 3500 m 3 ∕day, areas outside the initial subsidence event also began to subside (Fig. 9).This post-rebound, delayed subsidence was more widespread than the initial subsidence event.The subsided region includes the steam plant,  located at point PV, where piezometric elevation in a 115 m deep well decreased by 7.5 m in the weeks following the onset of deep well pumping (Fig. 8).Despite the rapid piezometric response, subsidence at the steam plant was not observed until April 2015 (Fig. 10), 5 months after deep well pumping began.From April 10, 2015, to April 9, 2019, point PV subsided a total of 2.4 cm at an average rate of 0.6 cm∕year.Points PX-PZ were all observed to reach a maximum  settlement before rebounding by ∼1 cm.The timing of the rebound is not uniform, as point PZ reached maximum settlement on October 23, 2018, whereas points PX and PY continued to subside until July 29, 2019.Points PQ-PT begin subsiding at similar times and rates to points PV-PZ but reach maximum subsidence sooner, as early as October 2016.Points PV-PZ continued subsiding, whereas points PQ-PT began rebounding in late 2018 and 2019, making a full recovery to their original elevations.
5 Discussion and Conclusions

Data Comparison
Vertical displacements calculated from Sentinel-1 data using the MinA algorithm are generally very similar to those calculated from RADARSAT-2 spotlight data (Fig. 11) using the MSBAS algorithm. 37Datasets from both RADARSAT-2 and Sentinel-1 possess C-band sensor wavelengths (∼5.5 cm); however, the RADARSAT-2 spotlight data possess a greater spatial resolution (1.6 × 0.8 m) than that of the Sentinel-1 data (2.3 × 14.1 m).The spatial resolution of InSAR data determines the area over which reflecting objects, or scatterers, can be identified.Therefore, a finer spatial resolution can potentially identify scatterers that coarser resolutions may not, creating the potential for two datasets of the same wavelength to measure reflections from different scatterers that do not necessarily exhibit the same deformation pattern.Despite the use of separate input datasets, the small differences in deformation output indicate that the MinA and MSBAS algorithms are comparable in terms of the ability to extract vertical deformations from combined ascending and descending InSAR data.
To further assess the accuracy of vertical and horizontal displacements derived from Sentinel-1 data, they were compared to available GPS data 51 within the study area (Fig. 12) and assume a common starting elevation.Differences in vertical elevation change measured by the two geodetic methods vary slightly over time.These differences can be attributed to several potential influencing factors.Small offsets between the location of the InSAR PS pixels and the GPS antenna can lead to measurement discrepancies due to the spatial differences in deformation magnitude or different targets of measurement (i.e., measuring the ground surface elevation changes versus that of the top of a building).Larger differences observed during the summer months can be attributed to seasonal elevation changes that are observed in the vertical GPS data but are not as strong in the InSAR time series.Error in the InSAR measurements is present and likely due to improper weighting during the optimization process, as well as the presence of residual phase that includes tropospheric effects from slave scenes, baseline measurement errors, non-linear deformation, and random noise.Deformation estimates shown here could likely be improved through the application of additional post-processing tools such as the Toolbox for Reducing Atmospheric InSAR noise 52 to reduce atmospheric contamination.Overall, both horizontal and vertical displacements derived from the InSAR time series are comparable in terms of the timing and magnitude of observed deformation.

Deformation
An analysis of post-subsidence surface deformation in downtown Seattle reveals complex hydromechanical coupling of groundwater with multiple aquifers.Subsidence is observed until the shutoff of access pit wells on February 19, 2015 [Fig.5(a)].Shutoff of the shallow wells initiated a period of opposing deformational processes [Fig.5(b)], where coarse-grained deposits within the upper, unconfined aquifer began to elastically rebound and fine-grained deposits continued to compress because of time-dependent consolidation behavior, e.g., Refs.53-55, in combination with compaction of the confined aquifer, as deep well pumping rates remained constant during this time.After the remaining wells were shut off, instantaneous and spatially heterogeneous uplift was observed at rates as high as 17 cm/year [Fig.5(c)].Rebound rates were observed to exceed those of initial subsidence, likely because the remaining wells were shut off over 2 days in January 2016 after being sequentially turned on over several weeks between October and November 2014, allowing for the rapid recovery of groundwater table elevation (Fig. 8).
The mapped surface rebound shows that areas of greatest rebound are not found in areas of greatest initial subsidence and that most subsided areas experienced some degree of permanent deformation.This is especially apparent at point PA, which lies in a region that initially subsided more than 3 cm but experienced little rebound.As it is well known that clay-rich deposits can be orders of magnitude more compressible than sand when drained, 1,56,57 the majority of permanent deformation can be attributed to the compaction of low-permeability aquitard layers.The magnitude of permanent deformation therefore depends on the depth and thickness of these clay-rich deposits, which are known to vary greatly within the subsurface, [28][29][30] resulting in spatially variable ground rebound (Figs. 6 and 7).][60] A 5-month lag is observed between the start of deep well pumping and the onset of subsidence at point PV, located at the steam plant where groundwater table elevations were lowered by 7.5 m in the months after deep well pumping began (Figs. 8 and 9).Delayed subsidence rates in areas outside the initial subsidence event, including the steam plant, were less than 1 cm/year but continued in some locations until July 2019, over 3 years after the final wells were shut off.This widespread lagging subsidence is likely caused by time-dependent consolidation of finegrained deposits within the confined aquifer, the full extent, depth, and thickness of which are poorly constrained.Point PZ subsided over 2 cm by October, 2018, and is located 2.7 km away from the access shaft and dewatering wells, suggesting that delayed subsidence extends beyond the area of initial rapid subsidence.A small degree of rebound is observed at points PV-PZ, likely due to an elastic response within coarse-grained layers.Some areas, including points PQ-PT, are observed to fully rebound.This indicates that sediments in these locations are primarily composed of coarse-grained material, allowing for a full elastic recovery in response to the restored pore water pressure.The spatial extent of delayed subsidence spans an area of land ∼20 km 2 in size (Fig. 9).
In conclusion, data show that hydro-mechanically coupled rebound can be highly unpredictable and reach rates in excess of 17 cm∕year.Rapid deformation events pose a hazard to the overlying infrastructure in urban environments underlain by highly variable geology, suggesting that more cautious approaches could be used when decommissioning dewatering wells in such areas.Delayed subsidence magnitudes of up to 3 cm were observed over 3 km away from dewatering wells and lasted for over 3 years after they were shut off, indicating that the deepest wells had accessed groundwater from the upper reaches of the confined Puget Sound aquifer.Areas impacted by delayed subsidence display a full continuum of elastic behavior, where areas composed of coarser-grained sediments make a full elastic rebound and areas composed of finergrained deposits undergo extended subsidence and less rebound.The spatial variability of rebound and duration of post-dewatering surface deformation are well captured by highresolution SAR data that can monitor spatially and temporally dynamic groundwater interactions with subsurface strata in greater detail than other point-based monitoring systems.

Fig. 1
Fig. 1 (a) Location and elevation of the study area.(b) Contours of subsidence previously measured by Samsonov et al.20 and locations of access shafts in which dewatering wells were installed, nearby GPS station, borehole used to characterize subsurface geology, and piezometer network used to monitor groundwater levels.

Fig. 3
Fig. 3 Map showing the location of PS pixels and the reference area used for analysis.

Fig. 4
Fig. 4 Illustration of possible solutions along an L-curve.Modified from Ref. 6.

Fig. 8
Fig. 8 Change in piezometric groundwater elevation over time.Piezometer locations as shown in Fig. 1.

Fig. 10
Fig. 10 Vertical deformation time series of points PV-PZ as shown in Fig. 9.

Fig. 12
Fig. 12 Vertical and horizontal deformation spanning the period of interest from GPS station V102 and MinA displacements derived from Sentinel-1 data at the same location.

Table 2
Input Sentinel-1 data used for the InSAR time series analysis.

Table 1
Lithological descriptions with depths from borehole 1.