Ground Surface Response to Geothermal Drilling and the Following Counteractions in Staufen im Breisgau (Germany) Investigated by TerraSAR-X Time Series Analysis and Geophysical Modeling

Abstract

The city of Staufen im Breisgau in southwest Germany suffers from a localized land uplift, which has occurred in the past six years in relation to geothermal drilling activities in 2007. So far, severe damages at 269 buildings have been recorded. The chemical transformation of anhydrite and water to gypsum, resulting in a volume increase, has been attributed as the cause of the uplift. Previous studies provide knowledge on the spatio-temporal displacement evolution from 2008 through 2011 using leveling and spaceborne Synthetic Aperture Radar Interferometry (InSAR) measurements, but lack a detailed representation of vertical and horizontal displacement contributions as well as geophysical modeling. This study focuses not only on continued observation analysis from June 2011 through July 2013, but also on obtaining and evaluating horizontal displacements in Staufen based on combined analysis of TerraSAR-X satellite imagery from both ascending and descending orbits. Applying the Small BAseline Subset (SBAS) approach a deceleration of annual cumulative line of sight (LOS) uplift is observable from 13.8 cm ± 0.3 cm (July 2008–July 2009) to 3 cm ± 0.3 cm (July 2012–July 2013) within area of maximum deformation NNE of the drilling zone. Conducting displacement decomposition on ascending and descending data of a common period (October 2012 through July 2013) yields in an approximately symmetric east- and westward motion with maximum values approximately 1 cm and 1.4 cm, respectively. The joint inversion of ascending and descending InSAR data for the common period from October 2012 through July 2013 shows that a horizontal rectangular source with length, width and depth of 177 m ± 19 m, 69 m ± 15 m and 89 m ± 9 m, respectively, can satisfactorily model the observation. The amount of opening at depth shows a decrease in time by about 71% for the period 2011–2012 as compared to period 2008–2009.

1. Introduction

The uplift process occurring in Staufen im Breisgau (hereafter called Staufen), southwest Germany (Figure 1a), is a remarkable localized event with first signs of damages being recognized at the end of 2007. Comprehensive investigations [1,2,3,4,5] have already been conducted to investigate the spatial and temporal pattern of deformation area in the city and to evaluate the causes to which the ongoing subsurface processes can be attributed to.

Seven boreholes for geothermal probes were drilled in September 2007 for the purpose of thermal management of the city hall. This is supposed to have facilitated water contact with anhydrite (calcium sulfate) lenses that are deposited within the Keuper-Gypsum strata. Anhydrite was found at depth between 61.5 m and 126.1 m and between 73.5 m and 105.75 m below ground level by analyzing core samples of the well drillings EKB2 and BB3, respectively [1,2], both located close to the city hall (Figure 1b). The chemical transformation of anhydrite and water to gypsum results in a volume increase [1]. Investigations by Sass and Burbaum [4] show that the swelling pressure cannot be compensated by the overburden, leading to surface deformation and severe damages. Until now, 262 private and seven municipal buildings sustained substantial crack damage; one building was demolished completely in August 2013 as it was damaged beyond repair [6]. Continuous monitoring is essential as it provides invaluable data for local authorities and stakeholders to reliably track ground deformation and assess its associated hazards.

Terrestrial leveling measurements have been performed since January 2008 [1] and regular data acquisition for applying Synthetic Aperture Radar Interferometry (InSAR) techniques started in July 2008, both aiming to assess the deformation area in space and time. The first InSAR study using Synthetic Aperture Radar (SAR) imagery of the TerraSAR-X (TSX) satellite from July 2008 through January 2009 revealed a localized displacement field oriented NE-SW. Displacement rates of up to 12 mm/month (14 cm/y in the case of linear motion) in the satellite’s line of sight (LOS) were derived for this period [4]. A more intensive investigation by Lubitz et al. [3] using time series analysis of TSX data from July 2008 to May 2011 found maximum LOS velocities of up to ~12 cm/y and showed evidences for distinct horizontal motions of several centimeters by comparing leveling and InSAR results. Moreover, Lubitz et al. [3] found a slowdown of uplift in time, which they relate to the countermeasures, namely groundwater pumping and borehole sealing, that have been conducted since September 2009 in the region [1].

The displacement analysis presented by Lubitz et al. [3] lacks a detailed representation of horizontal and vertical contributions. Whereas leveling measurements provide only information on the vertical component of displacement field, InSAR analysis combines both the vertical and horizontal motion contributions into one dimension, which is the line of sight from the ground to the satellite. Annually performed terrestrial position measurements by using tachymeter (distance and angular measurements) show strong evidences for horizontal dislocations, e.g., a survey point approximately 70 m NNE of the drilling zone moved about 6.4 cm in a northwestern direction from October 2009 through September 2010 and further 7.8 cm until November 2012. However, those measurements yield in a relatively poor spatial resolution and do not capture the potentially high spatial variability of displacement field. Therefore, from October 2012 through July 2013 we started TSX data acquisition in Staufen from a descending track in addition to the regular data acquisitions from the ascending track with the aim of differentiation between horizontal and vertical motion for the common time interval. We analyzed the new dataset in this study and also used the observed signals to constrain parameters of a geophysical pressure source model at depth, which in turn can be used to reconstruct surface deformation field [7]. The added value of standard deformation source modeling here is to obtain a first-order model for physical understanding of the source of the observed deformation in Staufen. The estimation of source depth is of main interest and will be considered with respect to anhydrite occurrences as detected by core sample analysis of the EKB2 and BB3 boreholes.

This paper is organized as follows: Section 2 is devoted to the data which we use in this research. In Section 3.1 and Section 3.2 we introduce briefly the InSAR methodology and the geophysical modeling, respectively. Results are presented in Section 4, separated into InSAR derived results (Section 4.1) and outcome of the source modeling (Section 4.2). Section 5 and Section 6 are devoted to the discussion and conclusion, respectively.

Figure 1. (a) Location of the city Staufen im Breisgau; (b) Sketch of well locations EKB2 and BB3 that are currently used for ground water pumping. An orthophoto was used as background (acquired by State Survey Office of Baden-Württemberg).

Figure 1. (a) Location of the city Staufen im Breisgau; (b) Sketch of well locations EKB2 and BB3 that are currently used for ground water pumping. An orthophoto was used as background (acquired by State Survey Office of Baden-Württemberg).

2. Data

All SAR imagery used in this study were acquired by the German TerraSAR-X satellite. The ascending dataset from the relative orbit 116 (Stripmap mode, beam strip_011, HH polarized) covers the period from 2 June 2011 through 11 July 2013 and hence follows the time period investigated by the previous study [3]. Descending Stripmap data from the relative orbit 154 (beam strip_009R, HH polarized) were acquired from 1 October 2012 through 14 July 2013 and have an overlap period of 276 days (approximately nine months) with the ascending data. Larger data gaps in the winter periods 2011/2012 and 2012/2013 occurred due to mission conflicts with the TerraSAR-X add-on for Digital Elevation Measurement (TanDEM-X) data acquisition for the WorldDEM^TM.

As mentioned earlier, leveling measurements have been performed since 2008. However, for comparison with InSAR results in this study, we use only those leveling measurements that were conducted from 23 May 2011 through 15 July 2013. Until November 2012, the surveys were repeated at a two-month interval, followed by quarterly measurements in 2013. The reduction of repeated measurements was justified with the deceleration of the uplift and in consequence with the reduced capability of leveling instruments to capture displacements at small temporal intervals [8]. Additionally, the horizontal dislocation measurements were performed at several survey points once a year on 2 October 2009, 28 September 2010, 26 September 2011, 11 September 2012 and 7 October 2013.

3. Methodology

3.1. InSAR

Interferometric time series analysis of the SAR images was done using the Small BAseline Subset (SBAS) algorithm as implemented in the StaMPS/MTI (Stanford Method for Persistent Scatterer/Multi-Temporal InSAR [9]) software. Differential interferograms were generated with DORIS (delft object oriented radar interferometric software [10]) by using a Light Detection and Ranging (LiDAR) Digital Elevation Model (DEM) with 1 m spatial resolution (from the State Survey Office of Baden-Württemberg) for topographic correction. The SBAS method in StaMPS identifies single-look slowly-decorrelating filtered phase (SDFP) pixels directly from full resolution SAR images [9] and therein differs from SBAS methods of Berardino et al. [11] and Lanari et al. [12]. For interferogram generation of the ascending data, the perpendicular and temporal baselines were constrained to be below 450 m and 250 days, respectively. In the case of the descending dataset, we could decrease the temporal baseline below 100 days as the image acquisitions show a denser temporal distribution, which reduces the problem of network connectivity. The networks of the interferograms with small spatial and temporal baselines of the ascending and descending images of this study are shown in Figure 2a,b, respectively. The network of the ascending dataset consists of 113 interferograms based on 30 SAR images that cover about two years, and 21 descending SAR images of a nine-month period are the basis for 82 interferograms of the descending network.

3.2. Source Modeling

InSAR derived surface displacement inversion can be used for estimating geophysical source models at depth that cause surface deformation [13,14,15]. In this paper, we use a finite rectangular opening-mode dislocation source in homogeneous elastic half-space [7] to simulate the observed uplift in Staufen. This model yields a convenient first-order approximation to evaluate complex source characteristics of deformation. However, it should be noted that this model is not aiming at a detailed and realistic representation of the Earth’s situation, since it does not include effects due to local structures, surface topography, crustal layering, lateral inhomogeneity or effects due to an obliquely layered medium.

Figure 2. The network of Small BAseline Subset (SBAS) interferograms used for the time series analysis of ascending (a) and descending (b) TerraSAR-X images. Red dots show the individual Synthetic Aperture Radar (SAR) images and green lines represent the image pairs for interferogram generation.

Figure 2. The network of Small BAseline Subset (SBAS) interferograms used for the time series analysis of ascending (a) and descending (b) TerraSAR-X images. Red dots show the individual Synthetic Aperture Radar (SAR) images and green lines represent the image pairs for interferogram generation.

Several geometrical parameters define the source. These are: length, width, depth, strike (orientation in clockwise direction), dip, midpoint coordinates of the upper edge as well as the amount of uniform opening. Relating these unknown source parameters $(x)$ and known InSAR displacement results in [16]:

$$d=G(x)+ e ,$$

(1)

where d is LOS displacement, G the Green’s function and $e$ the observation error, yields in a non-linear inverse problem, which we solved for the parameter vector $(x)$ by using a Genetic Algorithm [17] with Bayesian probabilistic approach as described in [18].

To improve the efficiency of the inversion we selected a priori bounds for the parameters that are considered to be optimized. For example, exploration activities in Staufen, summarized in [1], locate the swelling area between 61.5 m and 99.5 m depth, which could be used to constrain the model’s boundaries of the source depth. We have decided to select 0–1 km as boundary conditions for depths to avoid an overly restrictive setting. Further a priori bounds for source parameters (length, width and strike) were selected based on displacement maps (spatial dimension and orientation) from the previous study by Lubitz et al. [3]. The chosen boundaries are listed in Table 1, forming the parameter framework of the applied Genetic Algorithm that starts the first generation by using a population of 500 random source models. In order to find robust source parameters, 500 generations were then created.

Table 1. Lower and upper boundaries used in the inversion.

Boundaries	Length (km)	Width (km)	Depth (km)	Dip (°)	Strike (°)
Lower Boundary	0	0	0	−20	0
Upper Boundary	0.5	0.5	1	0	90

Table 1. Lower and upper boundaries used in the inversion.

Boundaries	Length (km)	Width (km)	Depth (km)	Dip (°)	Strike (°)
Lower Boundary	0	0	0	−20	0
Upper Boundary	0.5	0.5	1	0	90

4. Results

The multi-temporal analysis of TSX SAR images from both ascending and descending orbits by applying the SBAS method is the main basis of the present study. The results provide not only information on the spatial and temporal evolution of the surface displacement, but also serve as important input for the source modeling.

4.1. InSAR Derived Motion

Figure 3a,b present the cumulative LOS displacement obtained from ascending data acquired between 2 June 2011 and 11 July 2013 and from the descending dataset for the period 1 October 2012 through 14 July 2013, respectively. It shows an elliptical-shaped uplift bowl oriented NE-SW. The LOS displacement values are with respect to a reference area marked by a red star. For consistency, the reference area was chosen at the same position as for the data analysis in [3] and is located within a non-deforming area. The drilling area has been highlighted with a black circle. For comparison with previous study, a red dashed contour has been added representing the 1 cm uplift contour of the cumulative LOS displacement of the 2008–2011 period from [3]. Symbols for reference point, drilling zone and 1 cm boundary of the 2008–2011 SBAS result apply for all figures in this paper when needed.

By looking at the red contour in Figure 3a, we observe that compared to 2008–2011 the displacement pattern of 2011–2013 period has not changed significantly in its spatial extent. The local cumulative LOS uplift maximum for 2011–2013 is 8.5 cm, which is located approximately 50 m northward the drilling zone. Figure 3b shows a similar motion pattern derived from the descending dataset as compared to the ascending dataset. For this nine-month observation period the maximum cumulative uplift is 2.4 cm, which is located directly east of the drilling site. In comparison with the ascending data of both the 2008–2011 (red dashed contour) and the 2011–2013 (Figure 3a) periods, the descending result shows further expansion of displacement field to the east.

4.2. Modeled Source Parameters

For the geophysical modeling we considered not only the SAR data processed in this study (period June 2011 through July 2013), but also the SBAS results of the period July 2008 through May 2011 as presented in [3]. Investigating different periods of surface ground motion helps us to better assess the spatial and temporal variability of the deformation source at depth. Figure 4 provides an overview about the best-fit modeling results (left column) in comparison to the observed LOS displacements (middle column) for different periods. The 2008–2011 period is presented at four different intervals: from 22 July 2008 through 12 December 2008 (five months, Figure 4a), from 22 July 2008 through 11 August 2009 (eleven months, Figure 4b), from 22 July 2008 through 24 May 2010 (22 months, Figure 4c) and from 22 July 2008 through 22 May 2011 (34 months, Figure 4d). The 2011–2013 period is shown at two intervals: from 2 June 2011 through 21 June 2012 (12 months, Figure 4e) and from 2 June 2011 through 11 July 2013 (25 months, Figure 4f). Finally, modeling of the descending result is shown in Figure 4g covering the period from 1 October 2012 through 14 July 2013 (nine months).

Figure 3. (a) Cumulative line of sight (LOS) displacement in cm derived by using the Small BAseline Subsets (SBAS) approach based on the ascending data from 2 June 2011 through 11 July 2013. The red rectangle highlights the area of maximum deformation that is statistically analyzed in Section 5.1. (b) Cumulative LOS pattern based on the descending data from 1 October 2012 through 14 July 2013. The red star shows the reference area to which the displacement refers, the black circle indicates the drilling site and the red contour represents the 1 cm boundary of the ascending LOS displacement from 2008–2011 [3] for comparison.

Figure 3. (a) Cumulative line of sight (LOS) displacement in cm derived by using the Small BAseline Subsets (SBAS) approach based on the ascending data from 2 June 2011 through 11 July 2013. The red rectangle highlights the area of maximum deformation that is statistically analyzed in Section 5.1. (b) Cumulative LOS pattern based on the descending data from 1 October 2012 through 14 July 2013. The red star shows the reference area to which the displacement refers, the black circle indicates the drilling site and the red contour represents the 1 cm boundary of the ascending LOS displacement from 2008–2011 [3] for comparison.

Figure 4. Modeled surface displacement (first column) due to a rectangular finite dislocation source for different periods of the ascending data of [3] in (a)–(d) and of those processed in this study in (e)–(f). In 4(g) the descending outcome is shown. Corresponding source characteristics are listed in Table 2. Observed line of sight (LOS) cumulative displacements and residuals (difference between observation and model prediction) are shown in column two and three, respectively. The outline of the estimated rectangular source is shown by a dashed contour exemplarily for the corresponding period in 4(a), left column, and the black circle indicates the drilling field location.

Figure 4. Modeled surface displacement (first column) due to a rectangular finite dislocation source for different periods of the ascending data of [3] in (a)–(d) and of those processed in this study in (e)–(f). In 4(g) the descending outcome is shown. Corresponding source characteristics are listed in Table 2. Observed line of sight (LOS) cumulative displacements and residuals (difference between observation and model prediction) are shown in column two and three, respectively. The outline of the estimated rectangular source is shown by a dashed contour exemplarily for the corresponding period in 4(a), left column, and the black circle indicates the drilling field location.

In general, the model solutions for all individual periods fit observations well, confirmed by residual images (observation minus model) in the right column of Figure 4. However, noticeable residuals can be observed for datasets covering long periods (Figure 4c,d). Interestingly, later observations of similar time intervals (Figure 4f) show a significantly reduced magnitude of displacement, but a similar residual pattern as in Figure 4b,c,d, although of lower values. Several possibilities may be attributed to the residuals, which are discussed in Section 5.

The estimated source parameters corresponding to best-fit models in Figure 4 are summarized in Table 2. Similarities in the results of different periods need to be considered with care, since they are based on different observation geometries and separated processing.

Table 2. Source parameters from inversion of different Synthetic Aperture Radar Interferometry (InSAR) datasets for different time periods as presented in Figure 4. The standard deviation has been computed from conditional posterior probability density functions (PDF) as described in [18]. Values shown without standard deviation resulted from unimodal posterior PDFs.

Source Parameters	Ascending						Descending
	Jul 2008–Dec 2008	Jul 2008–Aug 2009	Jul 2008–May 2010	Jul 2008–May 2011	Jun 2011–Jun 2012	Jun 2011–Jul 2013	Oct 2012–Jul 2013
Length (m)	175 ± 27	173 ± 15	171 ± 14	164 ± 12	171 ± 32	174 ± 24	167 ± 26
Width (m)	51 ± 8	65 ± 9	69 ± 2	71	64 ± 14	64 ± 12	101 ± 12
Depth (m)	63 ± 12	49	46	45	51 ± 4	52 ± 2	48 ± 23
Dip (°)	−1 ± 4	0 ± 3	0 ± 2	0 ± 2	0 ± 5	0 ± 4	0 ± 6
Strike (°)	39 ± 9	38 ± 5	39 ± 5	38 ± 4	34 ± 10	33 ± 8	41 ± 12
Opening (cm)	15	28	42	56	8	15	2
Opening rate (cm/month)	3	2.6	1.9	1.7	0.7	0.6	0.2

Table 2. Source parameters from inversion of different Synthetic Aperture Radar Interferometry (InSAR) datasets for different time periods as presented in Figure 4. The standard deviation has been computed from conditional posterior probability density functions (PDF) as described in [18]. Values shown without standard deviation resulted from unimodal posterior PDFs.

Source Parameters	Ascending						Descending
	Jul 2008–Dec 2008	Jul 2008–Aug 2009	Jul 2008–May 2010	Jul 2008–May 2011	Jun 2011–Jun 2012	Jun 2011–Jul 2013	Oct 2012–Jul 2013
Length (m)	175 ± 27	173 ± 15	171 ± 14	164 ± 12	171 ± 32	174 ± 24	167 ± 26
Width (m)	51 ± 8	65 ± 9	69 ± 2	71	64 ± 14	64 ± 12	101 ± 12
Depth (m)	63 ± 12	49	46	45	51 ± 4	52 ± 2	48 ± 23
Dip (°)	−1 ± 4	0 ± 3	0 ± 2	0 ± 2	0 ± 5	0 ± 4	0 ± 6
Strike (°)	39 ± 9	38 ± 5	39 ± 5	38 ± 4	34 ± 10	33 ± 8	41 ± 12
Opening (cm)	15	28	42	56	8	15	2
Opening rate (cm/month)	3	2.6	1.9	1.7	0.7	0.6	0.2

Based on the ascending data, the dislocation source is found to be located between 45–63 m depth with a length ranging between 164 m and 175 m and a width ranging between 51 m and 71 m. The estimated parameters from the descending data lie within the range of the ascending results, except for the width, which is 101 m. As Table 2 shows, the dip angle is found to be zero or close to zero, representing a horizontal plane, which is oriented about 38° from the north (strike). Modeled strike values vary by 4°–5° for different datasets and periods. Except for the amount of opening, other parameters are found to be approximately similar in all inversions.

A joint inversion of the ascending and descending data covering a nine month period from October 2012 through July 2013 yields in similar parameter values for length, width, dip and strike (Table 3) as the individual inversions of ascending data only. Interestingly, from this inversion the source is estimated deeper (89 m ± 9 m) as compared to the single-data inversions, showing that combined data analysis enhances the geophysical modeling. The opening value from the combined inversion is three times larger than of descending data inversion only, resulting in a rate similar to that from ascending data inversion for June 2011–June 2012 and June 2011–July 2013 periods. The modeled and observed surface displacements as well as the residuals are presented in Figure 5.

As seen in Table 2 a relative opening of 56 cm was estimated within the first three years of SAR observation in Staufen, reaching to 15 cm for the following two years (2011–2013), and 2 cm for the nine month covered by descending data and 6 cm from joint inversion (Table 3), respectively. Considering the rate of opening (last row of Table 2 and Table 3), a decrease in time is clearly observable.

Table 3. Estimated source parameters from joint inversion of ascending and descending data for the period from October 2012 through July 2013.

Source Parameters	Oct 2012–Jul 2013
Length (m)	177 ± 19
Width (m)	69 ± 15
Depth (m)	89 ± 9
Dip (°)	0 ± 1
Strike (°)	37 ± 5
Opening (cm)	6
Opening Rate (cm/month)	0.7

Table 3. Estimated source parameters from joint inversion of ascending and descending data for the period from October 2012 through July 2013.

Source Parameters	Oct 2012–Jul 2013
Length (m)	177 ± 19
Width (m)	69 ± 15
Depth (m)	89 ± 9
Dip (°)	0 ± 1
Strike (°)	37 ± 5
Opening (cm)	6
Opening Rate (cm/month)	0.7

Figure 5. Modeled surface displacements (first column), observed line of sight (LOS) cumulative displacements (second column) and residuals (third column) from joint inversion of ascending and descending data for the period from 9 October 2012 through 11 July 2013. The black circle in the model-images indicates the location of the drilling site.

Figure 5. Modeled surface displacements (first column), observed line of sight (LOS) cumulative displacements (second column) and residuals (third column) from joint inversion of ascending and descending data for the period from 9 October 2012 through 11 July 2013. The black circle in the model-images indicates the location of the drilling site.

Our inversion ran with a Poisson ratio of 0.25. The Okada dislocation model that we have applied here does not include the effects of layering. The influence of material properties can be relevant, in particular in shallow, inhomogeneous and geologically complex regions as they can have an effect on the spatial extent of the surface deformation [19,20,21,22]. A study by Pritchard and Simons [23] analyzing dislocation models in elastic half-space suggests that the material properties seem to be of secondary relevance for source depth estimation as compared to other geometry parameters [23]. We tested our inversion with varying Poisson ratios for the dislocation source (0.09, 0.1, 0.13, 0.15 and 0.3) and found deviations of up to 5% for the source parameters listed in Table 2 and Table 3.

5. Discussion

In the following, we first discuss the temporal evolution of uplift, and then assess horizontal motions by evidences from comparing time series of SBAS and leveling, from combination of ascending and descending data as well as from source modeling. Finally, we discuss our modeling results.

5.1. Uplift Deceleration in Response to Counteractions

Monitoring of temporal variations in the displacement patterns allows us to evaluate the response of ground surface to countermeasures (borehole sealing and groundwater pumping) implemented in autumn 2009 and spring 2011 [2]. For this, we statistically analyze an area of 30 m by 30 m within the region of maximum LOS displacement about 50 m NNE of the drilling zone (red square in Figure 3a). At this location, we observe a cumulative LOS uplift of 13.7 cm ± 0.3 cm (mean value within the defined area) from July 2008–July 2009, which then was reduced to cumulative LOS uplift of 3 cm ± 0.3 cm four years later from July 2012–July 2013 (Table 4). The deceleration in uplift is significantly observable when comparing the periods 2009/2010 with 2010/2011 as well as the periods 2010/2011 with 2011/2012. A correlation with the permanent groundwater pumping at the new well BB3 since spring 2011 in addition to well EKB2 seems reasonable. When considering the decrease in the amount of opening (Table 2), we find compatible indications of the effectiveness of countermeasures.

Table 4. Line of sight (LOS) uplift statistics within area of maximum displacement (red square in Figure 3a) for different periods. The first three columns refer to the Small BAseline Subsets (SBAS) analysis of Lubitz et al. [3] (24 slowly-decorrelating filtered phase (SDFP) pixels within the defined area) and the last two represent the period under investigation of this study (27 SDFP pixels within the defined area).

Statistics	22 July 2008–31 July 2009	31 July 2009–18 July 2010	18 July 2010–22 May 2011	5 July 2011–13 July 2012	13 July 2012–11 July 2013
Minimum (cm)	12.8	10.3	6.4	3.7	2.3
Maximum (cm)	14.2	11.7	7.8	4.8	3.5
Mean (cm)	13.7	11.1	6.9	4.01	3
Median (cm)	13.8	11	6.8	4.1	3
Standard deviation (cm)	0.3	0.3	0.3	0.2	0.3

Table 4. Line of sight (LOS) uplift statistics within area of maximum displacement (red square in Figure 3a) for different periods. The first three columns refer to the Small BAseline Subsets (SBAS) analysis of Lubitz et al. [3] (24 slowly-decorrelating filtered phase (SDFP) pixels within the defined area) and the last two represent the period under investigation of this study (27 SDFP pixels within the defined area).

Statistics	22 July 2008–31 July 2009	31 July 2009–18 July 2010	18 July 2010–22 May 2011	5 July 2011–13 July 2012	13 July 2012–11 July 2013
Minimum (cm)	12.8	10.3	6.4	3.7	2.3
Maximum (cm)	14.2	11.7	7.8	4.8	3.5
Mean (cm)	13.7	11.1	6.9	4.01	3
Median (cm)	13.8	11	6.8	4.1	3
Standard deviation (cm)	0.3	0.3	0.3	0.2	0.3

5.2. Horizontal Displacement Estimation

5.2.1. Evidence from Comparing Time Series of SBAS and Leveling Measurements

Comparisons between ascending SBAS results and time series of leveling measurements for the 2011–2013 period at 13 selected locations within or in close distance to the uplift area are plotted in Figure 6. Nine locations (A, B, C, G, H, I, J, K, M) have been selected from previous study conducted for the 2008–2011 period [3] in order to better assess the temporal evolution of displacement. Three points have been selected to investigate the northeastern area (D, E, F) and one point has been selected south of the drilling zone (L). Similar to our earlier study, we compared directly vertical leveling measurements with LOS SBAS displacements to assess the contribution of horizontal motion in the deformation field. For each presented location the vertical displacements as measured by leveling are shown as red triangles and the mean LOS displacements as derived from SBAS are presented as blue dots. Mean values and error bars of SBAS results were computed from a window size of 7.5 m by 7.5 m with leveling point as center location.

The SBAS results show good consistency with leveling measurements obtained for points outside the displacement field (location A) and some locations in its interior (e.g., location M with approximately linear uplift). However, some clear discrepancies are observable at some locations. For example, at locations I, K and L (southeastern part of uplift area) significant deviations of the SBAS result from the leveling motion evolution are observed from the beginning of the investigated period (e.g., up to 2.4 cm for point I). Generally, contributions of eastward, northward and southward motion in addition to uplift yield in smaller LOS values than leveling measurements for observations from an ascending orbit. Such cases can be seen for locations I, K and L, which are located on the southeastern part of the graben block (Figure 7). The situation at point K is striking as the declining trend is not due to downward movement, but rather due to contribution of eastward and/or southward motion to overall LOS.

Comparison between InSAR and leveling measurement at point E indicates contributions of westward and/or northward motion in addition to the uplift, resulting in slightly larger LOS values than leveling measurements.

The other time series plots at locations B, C, D, F, G, H, and J show an interesting feature. First, their corresponding SBAS and leveling observations agree well until the winter period of 2011/2012, then a shift in spring 2012 is observed in which the SBAS derived motions show smaller magnitudes than leveling. No unwrapping problems have been recognized at those locations during processing, which may account for the observable shift. Among these points, B, C, G, H and J were analyzed before by [3], for which larger SBAS magnitudes than leveling for the period 2008–2011 were observed, an indication for the westward horizontal motion contributions. An explanation for the shift in the current uplift trend observable at these locations may be related to a transient deformation caused e.g., by reduction of westward motion or/and increase in northward motion since spring 2012. This observation cannot be verified against current ground reference measurements because of their poorer temporal resolution; the dislocation measurements were performed only annually in autumn season. Interestingly, all affected locations are at the horst block (Figure 7), which may indicate an influence of the fault system on the change of the motion rates. Complex situations of surface deformation bulges in interaction with faults can lead to complex effects on magnitude and direction of displacements [24]. Whether the observed shift reflects such complexity or even a reactivation of the fault system requires further investigations.

Figure 6. Time series of displacements at selected locations. Red triangles correspond to leveling measurement and blue dots with error bars to Small BAseline Subsets (SBAS) results using the ascending data from 2011–2013 time period. The locations are highlighted as black dots in the centered image (enlarged subset of the displacement field of Figure 3a).

Figure 6. Time series of displacements at selected locations. Red triangles correspond to leveling measurement and blue dots with error bars to Small BAseline Subsets (SBAS) results using the ascending data from 2011–2013 time period. The locations are highlighted as black dots in the centered image (enlarged subset of the displacement field of Figure 3a).

Figure 7. Indication of horizontal motions in approximately east- (red) and westward (blue) direction that is revealed by subtracting interpolated cumulative Small BAselines Subset (SBAS) line of sight (LOS) displacement and interpolated leveling results of the 2008–2011 period from [3]. Ascending orbit path and eastward look-direction (range) are displayed. Green filled pentagons show locations of survey points, for which annual dislocation measurements are performed. Arrows display horizontal motions for the period from 26 September 2011 through 7 October 2013. Local faults at 50 m depth (brown line features) as detected by the State Office for Geology, Resources and Mining (LGRB) were added. The locations, whose time series are presented in Figure 6, are marked with white filled cycles.

Figure 7. Indication of horizontal motions in approximately east- (red) and westward (blue) direction that is revealed by subtracting interpolated cumulative Small BAselines Subset (SBAS) line of sight (LOS) displacement and interpolated leveling results of the 2008–2011 period from [3]. Ascending orbit path and eastward look-direction (range) are displayed. Green filled pentagons show locations of survey points, for which annual dislocation measurements are performed. Arrows display horizontal motions for the period from 26 September 2011 through 7 October 2013. Local faults at 50 m depth (brown line features) as detected by the State Office for Geology, Resources and Mining (LGRB) were added. The locations, whose time series are presented in Figure 6, are marked with white filled cycles.

5.2.2. Evidence from Ascending and Descending Data Processing

Comparing Figure 3a,b, we observe a further extension of the displacement area to the east in the descending result. This observation can be explained in part by the different acquisition geometries of the two orbit paths (eastward and westward look direction, look angles of 39.7° and 36.1° for ascending and descending data, respectively) and by the orientation of the surface objects that scatter the radar signal back to the antenna. The influence of the object orientation may be evidenced by an increased number of SDFP pixels in the descending result, although the shorter time span (9 months) and in consequence a smaller descending dataset is also a reasonable explanation for a higher amount of SDFP pixels. Approximately 1.66 times more SDFP pixels were found by the time series analysis in the descending data than in the ascending data (21,474 versus 12,943).

Although different acquisition geometries cause small variations in the pattern of displacement obtained from using ascending and descending data, the overall pattern in both results is similar to each other. This demonstrates that the vertical process is dominant in Staufen. However, the contribution from horizontal displacement as discussed in the previous section make slight differences between ascending and descending results.

Subtraction of interpolated SBAS LOS results from interpolated leveling data for the 2008–2011 period (Figure 7) shows a pattern that is divided into an approximately westward (blue) and an eastward motion field (red) separated along a roughly NE-SW axis. This orientation agrees with the overall displacement field orientation and represents the consequent horizontal motion of a surface bulge. Interestingly, the westward motion occurs at the horst and the eastward motion at the graben block. Dislocation measurements at several survey points confirm this pattern for the later period from 26 September 2011 through 7 October 2013 as highlighted with arrows in Figure 7, showing the difference map between leveling and InSAR data from the ascending orbit.

In order to assess the pattern observed in Figure 7, we also compute horizontal motions by combining LOS displacements from ascending and descending SAR data. The common period covers approximately nine months, from 9 October 2012 through 11 July 2013. The magnitudes of SBAS derived displacements are similar for both datasets, with maximum of 2.4 cm for descending and 2.8 cm for ascending measurements. Data gaps in both ascending and descending results occur due to temporal decorrelation of the signal in vegetated areas (not shown here). To avoid noisiness or misinterpretation due to interpolation, we decided to consider only the measured point information. Hence, we have chosen the ascending SDFP locations to be the locations of interest and therefore computed a descending value at those positions by taking the mean of all descending SDFP values within a search window of 30 m by 30 m with the ascending SDFP location being the center point.

Horizontal motions were obtained by conducting the LOS decomposition [25,26] that takes into account the satellite heading (α) and look angles (θ) of the ascending (α_asc = 350.05°, θ_asc = 39.7°) and descending (α_desc = 190.45°, θ_desc = 36.1°) acquisitions. Applying [25]:

$$d_{LOS} = \text{cos}(\theta )d_V - \text{sin}(\theta )\text{cos}(\alpha )d_E + \text{sin}(\theta ) \text{sin}(\alpha )d_N$$

(2)

the sensitivity decomposition of LOS displacement for ascending data is [0.77 −0.63 −0.11] [d_V, d_E, d_N]^T and for descending data [0.81 0.58 −0.11] [d_V, d_E, d_N]^T. It shows that the east-west (d_E) motion sensitivity for both orbit paths is considerably large compared to the north-south (d_N) sensitivity and of similar relevance as the vertical motion (d_V). Due to near-polar orbit and right-looking viewing geometry of the sensor, north-south sensitivity is lowest in InSAR measurements. Therefore, we only focus on decomposing east-west and vertical motion by solving the linear system of equations. The results d_V and d_E are presented in Figure 8a,b, respectively.

Having a closer look at the east-west component (Figure 8b), the motion pattern is divided into two opposed parts following the expectations of the surface bulge’s expansion. Oriented in NE-SW direction, the deformation area is expected to experience horizontal motion westwards for the northwestern part and an eastward horizontal motion for the southeastern part. The magnitudes reach up to 1 cm for the eastward motion (positive values) and 1.4 cm for the westward motion (negative values). Annually performed dislocation measurements at survey points revealed maximal values for eastward and westward motion of 1 cm and 3 cm, respectively, from 9 November 2012 through 7 October 2013. Comparing these measurements with decomposed maxima (slightly different periods), we find comparable values for eastward motion and underestimation of westward motion (Figure 8b) by applying LOS decomposition. Although not comparable according to the magnitudes, the spatial distribution of horizontal motions of the nine-month-period is similar to the pattern derived by subtracting leveling information from 2008–2011 ascending SBAS result as shown in Figure 7.

It is worth noting that applying a forward model by using the best-fit Okada solution for the period 2 June 2011 through 11 July 2013 (source parameters from Table 2) provides a modeled surface displacement field that also can be decomposed into all three directions and hence help for better understanding of horizontal motions. For comparison with the LOS decomposed results, the vertical and east-west components are shown in Figure 8c,d, respectively. Although the modeled motion patterns represent a longer period, a pronounced similarity can be observed between Figure 8a,c,b,d, respectively, which provides another evidence for the contribution of horizontal displacement to the overall deformation phenomenon operative in Staufen.

Figure 8. (a) Vertical and (b) horizontal east-west motion contributions derived from line of sight (LOS) decomposition of ascending and descending SBAS results for the period 9 October 2012 through 11 July 2013. Estimated vertical (c) and east-west (d) motion from Okada modeling for the period 2 June 2011 through 11 July 2013 for comparison purposes. The black circle indicates the drilling site.

Figure 8. (a) Vertical and (b) horizontal east-west motion contributions derived from line of sight (LOS) decomposition of ascending and descending SBAS results for the period 9 October 2012 through 11 July 2013. Estimated vertical (c) and east-west (d) motion from Okada modeling for the period 2 June 2011 through 11 July 2013 for comparison purposes. The black circle indicates the drilling site.

Although there is the potential issue of small-scale sinkholes due to water-driven dissolution of the recently formed gypsum (karstification) [3,4], we currently do not see any evidence in the new data that was analyzed in the present work. The InSAR results from both ascending and descending orbits do not show any subsidence, when considered individually and in combined analysis (vertical motion component d_V based on LOS decomposition (Figure 8a).

5.3. Model Evaluation

Although the source model used in this study is simple, the estimated source characteristics (Table 2) appear to infer well some parameters obtained independently with geophysical investigations. The estimated source depth of 51 m ± 6 m (mean of inversion results of ascending data) from single-data inversion represents a shallow source located several meters above the upper anhydrite layer at 61.5 m depth, as obtained from core sample analyses of the well drilling EKB2, which is directly located within the drilling field of the geothermal probes (Figure 1b). The analysis of the excavation well BB3, which is located about 30 m east of the city hall revealed the upper anhydrite layer to be at 73.5 m depth [2]. By considering spatial variations of the depth of the upper anhydrite layer and the simplicity of the model, the estimated depth of the dislocation source represents the opening process to be operative close to the upper anhydrite layer based on single-data inversion (Figure 9d). The combined inversion of ascending and descending data leads to a source depth estimation of 89 m ± 9 m, which represents a depth within the interval of swellable region (61.5 m and 99.5 m below ground level) as assumed by LGRB [1]. Discussion about the estimated depth of a modeled source requires considering the trade-off between source strength and its depth. A deep and strong source can lead to the same surface deformation as a weak and shallow source [23]. Fixing the depth of the source to the upper anhydrite layer of the EKB2 well sample (61 m below ground level) for single-data inversion leads to an increase of source opening of 102% for the period 2008–2011 (56 cm vs. 113 cm). However, this will not affect the decrease in time for the amount of relative opening as shown in Table 2.

The situation, however, is quite complex. As anhydrite occurs as lenses, one could conclude multiple sources in combined action leading to the observed signal at the surface. This is hard to model as the deformation is localized and in consequence a single source model, despite being simple, was considered as an adequate first-order representation of the subsurface behavior.

The occurrence of residuals between modeled surface displacements due to a single pressure source and SBAS observations (Figure 4, right column) may partly indicate the described source complexity. Moreover, as mentioned in Section 3.2, the applied model does not include factors such as local geological structures, surface topography or crustal layering, which all can account for residual occurrence. Surface topography in terms of steep slopes can significantly influence modeled source geometry and depth leading to misinterpretation of deformation as Cayol and Cornet [27] have described for volcano deformations. As our study area is generally flat, this impact has not been considered in our source estimation. Influencing topographic effects can be related to urban features, which we considered by using a high-resolution LiDAR-DEM to remove the topographic phase contribution from the signal.

Figure 9a shows an orthophoto of the historical city center of Staufen overprinted with local faults located 50 m below ground level that were derived from LGRB investigations. To better investigate the residual pattern derived from modeling, we plot the residuals of the period from 22 July 2008 through 22 May 2011 in Figure 9b and show a profile in NW-SE direction across the residuals with cross sections of the faults (Figure 9c). The residual plot and profile show spatial correlation with geological blocks that are separated by faults, in particular of the graben block S and SE of the drilling zone (see profile between 225 m and 250 m distance).

Figure 9. (a) Orthophoto of the city center of Staufen (acquired by the State Survey Office Baden-Württemberg) overprinted with faults at 50 m below ground level derived from investigations of the State Office for Geology, Resources and Mining (LGRB). (b) Residuals between observation and model for the period from 22 July 2008 through 22 May 2011 (Figure 4d) were added to the orthophoto that is shown in 9(a). The black circle shows the location of the drilling area. (c) Model, observation and residual values (in cm) for the period from 22 July 2008 through 22 May 2011 along the profile AB that is presented in 9(b) with cross sections of the faults. (d) Sketch of estimated source depth from Synthetic Aperture Radar Interferometry (InSAR) inversion (single and joint inversion) compared to upper and lower anhydrite layer boundaries known from core sample analysis of the boreholes EKB2 and BB3. The location is not with respect to the profile, but related to the corresponding geological block.

Figure 9. (a) Orthophoto of the city center of Staufen (acquired by the State Survey Office Baden-Württemberg) overprinted with faults at 50 m below ground level derived from investigations of the State Office for Geology, Resources and Mining (LGRB). (b) Residuals between observation and model for the period from 22 July 2008 through 22 May 2011 (Figure 4d) were added to the orthophoto that is shown in 9(a). The black circle shows the location of the drilling area. (c) Model, observation and residual values (in cm) for the period from 22 July 2008 through 22 May 2011 along the profile AB that is presented in 9(b) with cross sections of the faults. (d) Sketch of estimated source depth from Synthetic Aperture Radar Interferometry (InSAR) inversion (single and joint inversion) compared to upper and lower anhydrite layer boundaries known from core sample analysis of the boreholes EKB2 and BB3. The location is not with respect to the profile, but related to the corresponding geological block.

Heterogeneous geological conditions as for example different anhydrite concentrations and varying rock masses can cause different swelling pressures at the surface [5], which cannot be predicted by our single source model. In addition, correlation with building locations and orientations, in particular NW of the drilling zone, are observable when comparing Figure 9a,b. Different building foundations and loads react differently to subsurface opening process, and hence different building stresses can influence the detectable surface deformation. For further information the reader is referred to [5]. Both heterogeneous geology (complex source situation) and building stress may explain the residuals between model and observation. Moreover, the assumed elastic character of the model is restricted here, in that the damages at the buildings are permanent and consequently non-elastic.

As there exists a potential relation between residuals and geological features, we recommend a detailed seismic monitoring or a transect of permanent GPS stations across the faults (from NW to SE) to assess their potential kinematics. The development of a complex hydro-mechanical model including not only different strata (material properties), faults and ground water levels, but also the possibility of several interacting swelling sources (anhydrite lenses) is suggested as potential future work.

6. Conclusions

In this study, we have focused on extending the knowledge about the remarkable urban uplift in Staufen by a combined analysis of TerraSAR-X images from both ascending and descending orbits. Interferometric time series analysis of spaceborne Synthetic Aperture Radar (InSAR) imagery confirms deceleration of uplift as measured with terrestrial surveying techniques [1,2]. Within the area of maximum deformation NNE of the drilling zone, we found a reduction of the annual cumulative line of sight uplift of up to 78% when comparing the periods July 2008–July 2009 and July 2012–July 2013.

Although uplift is the dominant process, the decomposed horizontal motions for the common nine-month period of ascending and descending data (October 2012 through July 2013) demonstrate significant east- and westward motions (up to 1.4 cm in westward and up to 1 cm in eastward direction) in an approximately symmetric pattern along the horst-graben boundary. Inversion of ascending and descending InSAR data individually suggest a shallow horizontal rectangular source with length ranging between 164 m and 175 m, width ranging between 51 m and 101 m and depth between 45 m and 63 m. Performing a joint inversion of both datasets of the period October 2012 through July 2013 results in length, width and depth of 177 m ± 19 m, 69 m ± 15 m and 89 m ± 9 m, respectively. Although the single-data inversion provides a good fit to the observation, the source depth resulted by joint inversion agrees more closely with the depth of swellable region inferred from independent core sample analysis. The strength of the source (the amount of relative opening) varies in time and shows a decrease by a factor of 3.5 when comparing the periods before counteraction in autumn 2009 and after permanent implementation in spring 2011. The amount of opening at depth has reduced by 71% for period 2011/2012 as compared to period 2008/2009.

Although the model adopted in this paper is a well-fitted first-order approximation, residuals between model and observation indicate potential influences of varying building stresses and fault kinematics. More geophysical investigations are required for their assessment and to also better resolve for additional deformation sources that may account for the complex subsurface situation.