Remarkable Urban Uplift in Staufen im Breisgau, Germany: Observations from TerraSAR-X InSAR and Leveling from 2008 to 2011

Abstract

As geothermal energy is of increasing importance as a renewable energy source, there is a high demand for comprehensive studies to prevent failure during implementation, as is the case in Staufen im Breisgau, Germany. The drilling of seven wells for the geothermal heating of the city hall in 2007 is thought to have disturbed the existing hydro-geological system in the complex structured transition zone of the Upper Rhine Graben and the Schwarzwald massif. This event has led to uplift, related to the transformation of anhydrite to gypsum, which affects the infrastructure of the city centre via the generation of large cracks. This study focuses on the application of the InSAR Small Baseline Subset (SBAS) approach using 50 X-band radar images from the German TerraSAR-X satellite (TSX) to map the spatial and temporal patterns of the deformation field in detail. X-band InSAR time series analysis for the three-year time period from July 2008 through May 2011 indicates maximum velocities of ∼12 cm/yr in the line of sight (LOS) direction, from the ground to the satellite, approximately 50 m northeast of the drilling field. In comparison with leveling data for the same time period, TSX data analysis better delineates the border of the deformation area, and it is able to map the amount of deformation associated with different parts of the city. Moreover, this comparison indicates contributions of horizontal motion, as is expected for uplift patterns.

1. Introduction

Staufen im Breisgau (hereafter referred to as Staufen) is a small city in southwestern Germany, approximately 50 km north of Basel and 120 km south of Karlsruhe (Figure 1). It is located at the eastern transition zone of the Upper Rhine Graben to the Schwarzwald massif, and their interaction controls the regional hydrology and geology [1]. A system of numerous faults has created a complex mosaic of blocks with varying strata sequences. Figure 1 shows the regional lineaments [2,3].

As part of the renovation of the historical city hall of Staufen and its neighboring building in 2006–2007, renewable shallow geothermal energy was selected as the new heating source. The region of the Upper Rhine Graben has a high potential for this technology due to its high subsurface temperature gradients [4]. Borehole drilling to a depth of 140 m was permitted in Staufen after an assessment of the geological and hydrological conditions. In September 2007, seven boreholes for geothermal double-U-tube probes were drilled, adjacent to the city hall, up to the permitted depth of 140 m [5]. Their locations are presented in Figure 2(a). At the end of the year, hairline cracks appeared in the city hall and in other neighboring buildings. The cracks (Figure 2(b)) increased rapidly in size and number. As of November 2011, 262 private buildings and seven public buildings have shown severe damage [6].

The deformation has been attributed to uplift that is caused by the disturbance of the existing subsurface hydro-geological system below the historical city centre. This disturbance is thought to have been triggered by the drilling, which most likely created connections between different groundwater horizons, above and below the strata of Keuper, and deposited anhydrite lenses. The mineral anhydrite (calcium sulphate) transforms into gypsum through contact with water and leads to a subsurface volume increase. As the drilling project has led to severe damage of buildings and some urban infrastructures, the project has since been stopped.

Deformation rates of up to 1 cm/month [5,7] were measured by conventional leveling surveys from February 2008 through February 2009. They have been carried out at regular intervals since January 2008, shortly after the first hairline cracks appeared, to map the spatial extent and magnitude of the surface deformation. Gradually, the number of survey points was increased to better capture the displacement field. The locations of the points are presented in Figure 3, color-coded by the date of their first measurement to show the expansion of the survey network. Between March 2009 and May 2010, the survey interval was decreased to two weeks for a subset of points located in the uplift centre and marked with black dots in Figure 3.

Although highly accurate, conventional leveling surveys are a time-consuming, labor-intensive, and expensive technique of providing adequately dense mapping of the spatio-temporal pattern of surface deformation. Spaceborne Synthetic Aperture Radar (SAR) sensors allow high-resolution monitoring of the earth’s surface at regular intervals and over large areas. The revisit time of the German TerraSAR-X (TSX) satellite for a single location, and with the same viewing geometry, is 11 days, and the spatial coverage of images recorded in Stripmap mode (resolution of approximately 3 m), which are used in this study, is approximately 30 km in range and 50 km in azimuth.

By applying time series methods of SAR Interferometry (InSAR), in particular, the Persistent Scatterer Interferometry (PSI) or the Small Baseline Subset (SBAS) approach, surface displacements at centimeter- to millimeter-level accuracies can be measured over time [8–11]. Both methods have shown their applicability in different fields of natural hazards and anthropogenically induced events [12–15]. The first InSAR-based observations of ground uplift in Staufen [16] are based on the analysis of nine TSX Stripmap images from July 2008 through January 2009 using differential repeat-pass interferometry. Further monitoring of the deformation field is essential for a better understanding of the on-going subsurface processes, in particular to determine if the rates of motion vary, and is helpful for the evaluation of counter measures. Hence, in this study, we present the result of the SBAS time series analysis of 50 TSX images acquired over Staufen from July 2008 to May 2011. Comparison with data from leveling surveys for the same period allows the evaluation of this multi-temporal InSAR method. We are also interested in what additional information the remotely sensed deformation data contain with respect to the leveling surveys.

This paper is organized as follows. Section two provides information about the available SAR data. The details of the SBAS methodology as implemented in the StaMPS software [17] are introduced in section three. Section four presents the results of the time series processing. In section five, we compare the SBAS product with the leveling data, consider the aspects of horizontal motion and discuss the result with respect to the geological system and against the background of the subsurface processes. The conclusion and summary are given in section 6.

2. Data

For our study, data from the German TSX satellite over a time period of almost three years—from 22 July 2008 through 22 May 2011—are available. The SAR images were acquired in Stripmap mode in an ascending orbit path direction (relative orbit 116, beam strip_011) with HH polarisation. The orbit path and image frame are schematically shown in Figure 1. The incidence angle ranges from 37.9° in the near range to 40.6° in the far range within the SAR image. The ground resolution is approximately 2.8 m in the range direction and approximately 3.3 m in the azimuth direction. In total, 53 TSX-images were acquired. Three images of the data stack, namely those from 12 January 2010, 14 February 2010 and 27 November 2010, display snow coverage leading to strong decorrelation for interferometric processing with respect to the other images. Hence, only the remaining 50 images were used for data processing.

A timeline of the SAR image acquisitions and leveling surveys is presented in Figure 4. It shows a high temporal sampling, in particular between March 2009 and June 2010, for both datasets with approximately similar periodicity: 11 days for the spaceborne data versus 14 days for the terrestrial recordings. This sampling allows the development of deformation over time to be resolved in detail and our remote sensing result to be compared with the terrestrial measurements.

3. Methodology

Classical repeat-pass SAR interferometry is a suitable tool for detecting surface displacement and monitoring coherent areas [18–21]. Interferometric coherence is a pixel-based measure of the signal similarity between two SAR images and ranges from 0 (total decorrelation) to 1 (total correlation). Temporal signal decorrelation is the limiting factor for InSAR, particularly in rural and vegetated areas due to the progressive change of the scattering properties [22] caused by vegetation growth, erosion, or harvesting. The InSAR time series methods of PSI and SBAS overcome this limitation. PSI focuses on identifying only coherent pixels dominated by a single scatterer mechanism that are imaged in every recorded SAR scene, while SBAS identifies those pixels showing only slow decorrelation over short time intervals [17]. Both techniques work better in urban areas because man-made objects, which often cause high-energy signal echoes, are highly coherent in time.

Due to strong uplift at a rate of several cm/yr in Staufen (inferred preliminary from leveling measurements), we applied the SBAS method as implemented in the StaMPS/MTI (Stanford Method for Persistent Scatterer/Multi-Temporal InSAR) software for the time series analysis. The SBAS approach refers to a network of individual differential interferograms (multi-master) that take into account minimized temporal, perpendicular and Doppler baselines to minimize decorrelation effects. In terms of high displacement rates, using multi-master interferograms with short temporal baselines reduces the number of phase difference ambiguities to be resolved, potentially allowing a more reliable phase unwrapping than the single-master PSI technique without a priori temporal model phase unwrapping.

Figure 5 shows the network of the small baseline interferograms used for the SBAS analysis of surface deformation in Staufen. The interferograms were computed from TSX single-look complex (SLC) images by using DORIS (Delft object oriented radar interferometric software [23]). Prior focusing was not necessary, as the TSX data are delivered in SLC format. Initially, 289 interferograms were processed, but only 249 interferograms were used in the SBAS-processing due to the removal of the three SAR images with strong snow cover.

As differential interferograms are required for displacement analysis, we removed the topographic information from the interferometric phase by transforming a Digital Elevation Model (DEM) into phase values and subtracting it from the interferograms. We used a high-resolution LiDAR DEM with 1 m resolution (provided by the State Survey Office of Baden-Württemberg) for the topographic phase correction and geocoding of interferograms.

In the SBAS implementation of StaMPS/MTI, time series analysis is only conducted for those pixels for which the phase, after filtering of the azimuth and range spectra, shows slow decorrelation over short time intervals, the so-called slowly decorrelating filtered phase (SDFP) pixels [17]. The previous filtering aims to reduce decorrelation effects due to geometry and non-overlapping Doppler frequencies. A computational cost-effective initial selection of the SDFP candidates was performed at the highest resolution (single-look image) by setting a threshold (0.6) for the amplitude difference dispersion, which is the standard deviation of the amplitude difference between the master and slave divided by the mean amplitude. Afterwards, an iteratively conducted phase analysis yielded the phase stability estimation of each candidate. In the StaMPS algorithm, after the removal of topography and earth flatting, spatial correlation was assumed for the differential interferometric phase residuals of SDFP pixels related to deformation of the earth’s surface, the spatially correlated height error, the atmospheric delay, and orbit inaccuracies. These contributions were estimated by the bandpass filtering of surrounding pixels. The statistical analysis of the spatially correlated terms, the spatially uncorrelated (look-angle error) terms, and the amplitude difference dispersion were used for the final SDFP pixel selection. This method estimates displacement time series without prior consideration of a temporal model for the deformation and thus allows the derivation of temporally varying deformation processes.

Three-dimensional phase unwrapping was applied on the final sets of SDFP pixels by using a statistical cost flow algorithm that works for both single-master (PSI) and multi-master (SBAS) methods [24]. The time series of displacement for each SDFP pixel was derived by least-squares inversion. The reader is referred to [17] and [24–26] for further details on the SBAS algorithm, the unwrapping method, and StaMPS in general.

4. Results

The multi-temporal SBAS processing of the available TSX dataset reveals a clear NE-SW elliptical-shaped deformation pattern, as presented in Figure 6. As the motion is towards the satellite, it is specified as uplift. Figure 6 shows the cumulative uplift on 22 May 2011, after three years of SAR observation. The displacement field extends approximately 290 m NE-SW and approximately 190 m NW-SE. A large part of the deformation field has risen more than 25 cm (dark blue to violet) in the line of sight (LOS), and the maximum three-year cumulative uplift of 33.5 cm is located approximately 50 m northeast of the drilling field.

The temporal evolution of the uplift at three-month intervals is illustrated in Figure 7, which apparently shows a smooth increase in uplift over time. However, a closer look at the differences in deformation at the selected dates (Figure 8) reveals that with time, the uplift rate decreased since autumn/winter 2009. For example, comparing the uplift rates of a SDFP pixel located within the main deformation area (47.882313°N, 7.732191°E) for the same season in different years, we observe a decrease by 43 %; 4.12 cm between 29 October 2008 and 25 January 25 2009, compared to 2.35 cm between 27 October 2009 and 23 January 2010.

5. Discussion

5.1. Leveling and InSAR Comparison

To evaluate the result of the TSX SBAS processing of Staufen, we compare it with leveling measurements that cover approximately the same period as the SAR imagery (Figure 4). The cumulative vertical motion, measured from 19 August 2008 through 23 May 2011, of each survey point (triangles) is overlaid on the SBAS result (circles) in Figure 9 and shown in the same color scale as the LOS displacement of the SDFP pixels. In comparison, the SBAS method represents the surface displacement at a higher level of detail due to the higher density of pixels in the InSAR analysis. Dense point clouds are an advantage of InSAR time series methods, especially if applied to urban areas. Man-made objects usually guarantee a high coherence over time. For vegetated areas, as in the upper right corner of Figure 9, no SDFP pixels were found, demonstrating the decorrelation-related limitation of the method. Although not yet performed, leveling surveys could fill such gaps and provide a complementary analysis to the InSAR and existing leveling data.

The boundary of the SBAS-derived uplift ellipse agrees well with the leveling measurements, as well as the amount of displacement of points at the border. Differences on the order of several centimeters occur in the inner part of the displacement area. To determine the cause of deformation differences, the displacement time series of several survey points (red square in Figure 9) was compared with that of nearby SDFP pixels, which were intentionally selected from different parts of the deformation field. The selected points are labeled with letters (A–I) from west to east for simplicity, but for reconstruction purposes, the reader is referred to the official survey point IDs provided on the left side of Figure 9.

The time series plotted in Figure 10 show approximately uniform (linear) motion for both leveling and InSAR measurements, but the displacement rate decreases slightly over time, which generally coincides with the areal results presented in Figure 8. The leveling points and SDFP pixels in Figure 10 are marked with triangles and circles, respectively. Corresponding points are presented in the same color. As leveling and InSAR data acquisitions began on different dates, the time series have been adjusted to the same start date via linear interpolation. The time series start on 22 July 2008 (first TSX image), but the data for some points are shown only from 20 January 2009 due to the gradually increasing number of leveling marks.

As shown in Figure 10, the selected points show different amounts of uplift depending on their location. Leveling point A and the corresponding selected SDFP pixel (SBAS A) are located at the border of the deformation field (Figure 9), at which both types of measurements coincide (Figure 10(b)). Among the selected point pairs, pair E presents the largest cumulative uplift of 28.8 cm (leveling) and 31.2 cm (SBAS). This point pair is located north of, but close to, the drilling field. For all other points, the general trend is similar for both leveling and SBAS-derived time series, but the time series show discrepancies in magnitude between leveling and InSAR. According to those differences, Figure 10 is divided into three sub-charts: (a) shows the time series of points with larger SBAS values than leveling values, (b) gives those with small discrepancies between InSAR values and leveling values and (c) presents those with larger leveling values than SBAS values.

5.2. Implications of Horizontal Motion

The discrepancy between InSAR and leveling points, as depicted in Figure 10, can be attributed to the contribution of horizontal displacement. Figure 11 provides a simplified scheme to visualize the “positive” or “negative” contribution of horizontal motion to the overall LOS displacement.

In the case of horizontal motion toward the satellite, its contribution is added to the vertical motion (uplift), resulting in a LOS displacement of higher magnitude than measured by leveling (Figure 11(a) and time series in Figure 10(a)). All point pairs showing this type of difference in their time series are suggested to also experience horizontal motion towards the west. Point pairs showing no discrepancy in magnitude likely represent only vertical motion. The opposite case of Figure 10(a) is shown in Figure 10(c), where leveling measurements of selected points show larger magnitudes of motion than those based on SBAS. In the case of horizontal motion away from the ascending satellite (here, toward the northeast), the uplift and horizontal components have opposing impacts on the LOS motion (Figure 11(b)). This results in a reduced LOS displacement compared to the leveling value. Depending on the incidence angle, the magnitude of uplift, and the magnitude and direction of horizontal motion, this effect can theoretically lead to an increase in slant range (see Figure 11(b), t1 versus t1’) and can be misleadingly interpreted as subsidence. Hence, a displacement field interpretation of a one-dimensional LOS result is limited, and ascending and descending data for the same period are required. Work is currently underway to precisely assess the amount of horizontal deformation in Staufen using ground- and satellite-based geodetic measurements [28,29].

5.3. Geological Considerations

The region surrounding the Upper Rhine Graben is a complex structured system of blocks and faults dominated by the extensional tectonics of the Upper Rhine Graben and its interaction with the Schwarzwald massif. The main fracture zones attributed to the Upper Rhine Graben strike NE-SW, whereas those of the Schwarzwald crystalline primarily strike E-W or NW-SE (Figure 1). The subsurface beneath the city centre of Staufen reflects this typical fault formation of the transition zone. Several horst and graben blocks are oriented NE-SW and are cross-bordered by faults striking NW-SE [5].

Potential natural causes for the uplift in Staufen are unlikely, not only due to the strong temporal and spatial correlation of the deformation field with the drilling activities but also due to a lack of evidence pointing to other causes. No earthquakes have been recorded in this particular period [5].

Different exploration activities have been performed by the National Agency for Geology, Resources and Mining (LGRB) to obtain information about surface and subsurface geological structures and to achieve a detailed understanding of the processes that have led to the observed surface deformation. These activities include analysis of samples from boreholes EWS 1, EWS 2, and EWS 3 (Figure 2(a)), geological exploration by the drilling of wells EKB 1 and EKB 2 (location shown in Figure 2(a) in blue), two-dimensional seismic surveys along five profiles in 2009, and temperature profiling in all seven wells. Based on the geophysical measurements (seismic surveys) and drilling samples, the locations of faults were identified and several blocks were characterized [5].

According to the LGRB investigations, the drilling area is located on a graben block, and the core sample analysis revealed the typical stratigraphic sequence of the region: upper Muschelkalk, lower Keuper, middle Keuper (Schilfsandstein and Keuper-Gypsum), and quaternary sediments. Anhydrite (calcium sulphate) and swellable clay minerals (mainly corrensite) were found in the non-lixiviated Keuper-Gypsum layer of the middle Keuper. Four hydraulically disconnected groundwater levels were also discovered in the graben block in which the drilling was performed. All boreholes dip to the southeast, and several of them cross faults that separate the graben block into smaller pieces [5]. Borehole EWS 7 reaches the hydraulically separated groundwater aquifer of the lower Keuper. During the drilling in September 2007, artesian groundwater was reported. It is suggested that water found a path through the insufficiently sealed annular space of one or more boreholes to anhydrite lenses. Temperature anomalies were detected in all boreholes except EWS 1 (maximum of 0.7 °C in EWS 7 between 12 February 2009 and 17 October 2009), which may indicate anhydrite swelling due to the exothermic character of the transformation into gypsum. This chemical transformation is accompanied by a volume increase that could theoretically reach 60%, if there is sufficient water [30]. Due to the pressure change, new cracks may occur into all directions within the Keuper-Gypsum, potentially forming additional paths for water to travel to other anhydrite lenses. The swelling is located between 61 m and 99 m depth [5], where the overlying rock mass is too small to hamper the swelling [16]. As the quantities and distribution of anhydrite and groundwater influx are unknown, the complex spatio-temporal evolution of subsurface processes cannot be described in detail. Nevertheless, the observed uplift at the surface can be considered an effect of the described chemical reactions.

The existing geological and hydro-geological system may have an influence on the spatial pattern of the uplift in Staufen [31,32]. The localized elliptical shape of the deformation field is oriented NE-SW (Figure 6), and thus a significant structural control by the Upper Rhine Graben related horst- and graben-blocks can be assumed. In Figure 12(a), the SBAS-derived cumulative uplift pattern from 22 May 2011 is overlaid with the geological structures derived from the geophysical explorations of the LGRB and from the 1:25,000 geological map to investigate potential spatial correlations. To better assess the deformation field with respect to the geology, we extracted two profiles across the area of uplift: profile A (NW-SE) and profile B (NE-SW). As observed in Figure 12(b) and 12(c), both profiles show convex-downward curvature as they reach their area of maximum uplift. For Profile B, no measurement points are available in the vegetated area in the northeast, which reduces the mapping of deformation there.

As shown in Figure 12(a,b), the southeastern boundary of the deformation field seems to be controlled by one of the NE-SW striking faults, which are derived from geophysical measurements of the LGRB. Due to the limitation of temporal decorrelation in vegetated areas, the InSAR data do not provide enough information on the northeastern part of the deformation area (Figure 12 (a,c)). Hence, the possible bounding of the uplift field on its eastern edge by the assumed “Roter-Berg” Fault, which strikes NWSE, cannot be confirmed nor excluded by the existing data. Leveling measurements from a dense network of survey points in this area could fill such a data gap and would show the complementary benefit of the two measurement techniques. The current geophysical knowledge is not sufficient to explain the spatial limit of the deformation field to the north and west as detected by applying InSAR. Evidence for bounding faults are not yet obtained, but cannot be excluded from the considerations.

Counter measures against the on-going subsurface swelling have been applied, namely the backfilling of the wells since November 2009 and the pumping of water since September 2009 to reduce the water inflow to the swellable strata [33]. Exploration well EKB 2 is used for this purpose, with a constant pumping rate of 1.9 L/s. Since March 2011, this has been supported by an additional well southeast of the drilling field that pumps at a rate of 2.6 L/s [28]. Figure 8 and both the SBAS and leveling time series provided in Figure 10 show a deceleration of uplift that can likely be attributed to these remediation activities. Another possible scenario, as described by [34], may be the natural evolution of the geochemical process of transforming anhydrite to gypsum, which can stop under certain external conditions (depending on the temperature and salinity of groundwater and the stiffness of the surrounding rock mass). Although theoretically possible, this scenario cannot be evaluated because the amount of anhydrite, gypsum, and water that has entered the system is unknown. An estimation of the time until complete transformation, the total volume change, and hence the final amount of uplift is not possible. Alternatively, a reduced rate of anhydrite transformation due to the isolation of anhydrite from water by the recently created gypsum crust, which acts as low-permeability barrier [37], could explain the deceleration of observed uplift in Figure 8 and the time series in Figure 10. However, the effectiveness of counter measures does seem reasonable.

A potential future scenario may arise: the dissolution of the recently formed gypsum in water (karstification), leading to a volume loss. The risk involved is in this the structural collapse of possible subsurface cavities [35], which may lead to small-scale sinkholes [16], thus causing further damage. SBAS time series analysis has the potential of capturing early signs of subsidence caused by gypsum dissolution due to more densely sampled points than leveling measurements and the small image acquisition interval of 11 days. Considering the current spatio-temporal results (Figure 6 and Figure 10), indications of such a process are not currently visible. Continued monitoring of the area is therefore strongly recommended.

6. Conclusion

In this study, we have demonstrated the capability of X-band InSAR time series analysis to record surface displacements occurring at rapid rates in Staufen, Germany. Because this method considers interferograms of small temporal baselines, it is suitable for detecting and correctly unwrapping the rapid surface uplift in Staufen, which occurs at rates greater than 10 cm/yr. The spatially limited deformation field is precisely reflected in the SBAS time series; its spatial extent and the trend coincide with the leveling measurements. Missing survey points and the incapability of the InSAR technique to overcome temporal decorrelations in vegetated areas result in a lack of information northeast of the uplift area. The availability and analysis of both remote sensing and terrestrial survey data indicate horizontal velocity contributions of several centimeters. As neither detailed terrestrial measurements of horizontal deformation, nor proper TSX data acquisitions from a descending orbit are available for the observation period, the full three-dimensional velocity field cannot be determined. Accordingly, the evaluation and validation of the InSAR-based deformation field at Staufen by using leveling data is only partially possible. Future data acquisitions in ascending and descending orbits, in combination with recent GNSS measurements, enable the discrimination between horizontal and vertical displacement. The combination of these various datasets together with either numerical or analytical geophysical modeling could provide additional details for geological interpretations and lead to a better understanding of the complex subsurface processes.