Modeling of fast ground subsidence observed in southern Saskatchewan (Canada) during 2008–2011

Fast ground subsidence in southern Saskatchewan (Canada) between the city of Saskatoon and Rice Lake was observed with the RADARSAT-2 interferometric synthetic aperture radar (InSAR) during 2008–2011. We collected 23 ascending Multi-Look Fine 3 Far (MF3F) and 15 descending Standard 3 (S3) RADARSAT-2 images and performed time-series analysis utilizing Small Baseline Subset (SBAS) and Multidimensional SBAS (MSBAS) methodologies. We observed two well-defined circular regions located a few kilometers apart and subsiding with the nearly constant rate of about 10 cmyr −1. MSBAS decomposition revealed the presence of both vertical and horizontal ground displacements. For further analysis we selected two highly coherent interferograms spanning from November to December 2009 until April 2010 thanks to particularly favorable ground conditions that displayed superior coherence. We performed modeling and inversion assuming spherical and sill source models in order to determine the source location, depth and strength. The sill source model produced the smallest residual of 0.7 cmyr −1 applied to ascending interferograms and 0.9 cmyr−1 applied to descending interferograms. A residual of 1.0 cmyr−1 was achieved with the sill model when both ascending and descending interferograms were used. This model suggested sources located at 1.3 and 1.2 km depth with radius of 1.0 and 1.3 km for eastern and western areas, respectively. The spherical model suggested slightly shallower sources located at 0.9 and 0.8 km. We could not precisely identify the cause of this deformation, but the observed subsidence rate and source depth suggest mining-related origin. Topographic changes produced by this subsidence rate over a long time may produce shallow groundwater redistribution and flooding of agricultural lands.


Introduction
Ground subsidence, a downward motion of the earth surface in relation to its position at an earlier time, can be produced by natural and anthropogenic phenomena.Among the natural causes are the earthquakes and more general fault motion (Beavan et al., 2010(Beavan et al., , 2011;;Wen et al., 2013), volcanic thermal and pressurized sources (Samsonov and d'Oreye, 2012), karst processes (Klimchouk, 2005), consolidation of sediments (Mazzotti et al., 2009) and permafrost thawing (Short et al., 2011;Chen et al., 2012).The anthropogenic ground subsidence is usually caused by mining (Gourmelen et al., 2007;Samsonov et al., 2013a), groundwater (Bawden et al., 2001;Samsonov et al., 2010Samsonov et al., , 2011a;;González andFernández, 2011), oil andgas extraction (Mayuga andAllen, 1969;Geertsma, 1973), and heavy weight loading (Samsonov et al., 2014).The temporal pattern of subsidence can vary from very slow, occurring over a long period of time (Samsonov and d'Oreye, 2012;Samsonov et al., 2013a), to nearly abrupt with no pre-or post-event deformation (Beavan et al., 2010(Beavan et al., , 2011)).Subsidence produced by earthquakes, mining collapse and sinkholes generated by karst processes is usually spatially and temporally localized, with no detectable preevent motion; other processes usually occur over an extended period of time and can be detected by ground-based or spaceborne measurements.Ground uplift or heave, the process opposite to ground subsidence, is also observed but with a rare frequency and usually associated with natural volcanic activities (Biggs et al., 2010;Amelung et al., 2000;Fialko and Pearse, 2012) or anthropogenic liquid, gas or steam injection (Teatini et al., 2011;Samsonov et al., 2013a).
Three-dimensional high-frequency measurements of ground motion by a global position system (GPS) receiver installed at a subsiding site can map the progress of subsidence as a function of time (Hofmann-Wellenhof et al., 2001).However, due to high equipment and installation costs GPS measurements often cannot achieve sufficient spatial resolution necessary for initial detection of subsidence, as well as for determining its present spatial extent (Lanari et al., 2004;Gourmelen et al., 2007).In contrast space-borne interferometric synthetic aperture radar (In-SAR) measurements are particularly capable of mapping ground deformation with a very high spatial resolution over a large area, with a high precision and a moderate to low cost (Sansosti et al., 2010;Cigna et al., 2011;Hooper et al., 2012;Lubitz et al., 2013).
An InSAR product or interferogram is computed from two SAR images acquired by a space-borne sensor at two different times and captures the occurred line-of-sight (LOS) displacements (Massonnet and Feigl, 1998;Rosen et al., 2000).InSAR processing while being complicated is well understood and can be easily performed using one of many available commercial or open source software packages.Repeatedly acquired SAR images can be used to produce stacks of interferograms that can be further presented as a time series, similar to those produced by GPS but with significantly lower temporal resolution equal to the sensor revisit time.Typical spatial resolution of SAR sensors is 1-20 m over 10 km × 10 km to 300 km × 300 km area, and temporal resolution is 1-41 days, with 24-35 days being most common.Various modern space-borne sensors, satellites and satellite constellations, are capable of acquiring SAR data with various characteristics that can be adjusted to a particular application.For example, the Canadian RADARSAT-2 satellite can acquire spotlight beam data with range-azimuth resolution of 1.6 m × 0.8 m and coverage of 18 km × 8 km and wide beam data with range-azimuth resolution of 13.5 m × 7.7 m and coverage of 150 km × 150 km, but many other beam configurations are available depending on the user requirement.SAR data polarization, acquisition incidence angle, and left or right looking geometry can also be easily manipulated.The air-borne SAR sensors are also available, but their operational cost is significantly higher and processing methodology is significantly more complex, which somewhat limits their accessibility for regular users.
In this paper we present analysis of relatively fast ground subsidence observed in southern Saskatchewan (Canada) in a sparsely populated agricultural area located between the city of Saskatoon and Rice Lake.The subsidence consists of two nearly circular regions 1-2 km in diameter and separated by about 1 km in the northeast to southwest direction.Initial study using the Multidimensional Small Baseline Subset (MSBAS) method (Samsonov and d'Oreye, 2012) revealed the presence of both fast vertical (∼10 cm yr −1 ) and moderate horizontal (∼4 cm yr −1 ) motion occurring with a nearly constant rate during 2008-2011 (Samsonov et al., 2013b).In this work, we perform modeling and inversion of the In-SAR results in order to determine the source parameters, such as location, depth, strength and shape, and discuss potential causes of subsidence.Additionally we investigate the impact of subsidence on the environment at the local scale in relation to potential increase of flooding hazard of the agricultural lands overlaying this subsidence.Geologically, this region is characterized by an approximately 80-100 m succession of Tertiary to Quaternary age sediments of glacial origin (till, clay, silt, sand, and gravel) deposited over Cretaceous age basement rocks (Saskatchewan Geological Survey, 2013a).The terrain surface area is mainly covered by agricultural fields, with extensive water bodies recharged by precipitation, snowmelt and shallow groundwater.

Data processing
For mapping ground deformation we collected RADARSAT-2 SAR data from two different beams: ascending Multi-Look Fine 3 Far (MF3F) and descending Standard 3 (S3).The SAR data coverage and satellite acquisition parameters are presented in Fig. 1  Each SAR data set was processed independently with GAMMA software (Wegmuller and Werner, 1997) in the following way.A single master for each set was selected, and remaining images were re-sampled into the master geometry.All possible differential interferograms with perpendicular baseline less than 400 m were generated while applying multilooking 7 × 7 and 3 × 6 for MF3F and S3 interferograms respectively.The topographic phase was removed using the 20 m resolution digital elevation model (DEM) provided by GeoBase (2013) with an approximate reported elevation accuracy of 5 m.Interferograms were filtered using the adaptive filtering with filtering function based on local fringe spectrum (Goldstein and Werner, 1998) and unwrapped using the minimum cost flow algorithm (Costantini, 1998).A residual orbital ramp was observed in many interferograms, and baseline refinement was performed to correct it.The minor interpolation of each interferogram was performed in order to increase the coverage reduced by decorrelation.The  corrected interferograms were geocoded to the DEM grid, and the sub-region of interest outlined in brown in Fig. 1 was re-sampled to a common latitude/longitude grid with GMT scripts (http://gmt.soest.hawaii.edu/).On average the coherence of interferograms was very low, which can be explained by seasonal changes in the land cover conditions, including snow coverage and flooding.Nevertheless, for the Small Baseline Subset (SBAS) processing we could select 29 MF3F and 23 S3 interferograms with the mean coherence after filtering above 0.5.We applied SBAS algorithm described in Samsonov et al. (2011b) that solves for deformation rates and the residual topographic noise simultaneously.No additional filtering was applied, and time series were reconstructed by integration of deformation rates.
The time series for two regions of fastest subsidence (from here on called western and eastern) are presented in Fig. 2a  and b.These results demonstrate that during the observation period the rate of deformation is nearly constant and not affected by seasonal groundwater fluctuation.
In order to analyze the spatial extent of this subsidence, we applied the Multidimensional Small Baseline Subset (MS-BAS) method (Samsonov and d'Oreye, 2012).This method integrates ascending and descending InSAR data and produces two-dimensional time series and rates of ground de- formation.During MSBAS processing we applied 2-D highpass spatial filtering to remove a plane from the interferograms caused by minor errors in baseline estimation and by long wave-length atmospheric noise.Calibration was performed by subtracting a constant offset from each interferogram calculated against a single reference region that was assumed to be stable.The topographic correction was applied, and the regularization parameter was selected equal to 0.25.The 1-D temporal and 2-D spatial low-pass filtering were not performed.The inversion was performed, and vertical and horizontal east-west time series and linear deformation rates were computed.We present here only the vertical time series in Fig. 3 as horizontal deformation rate and time series for selected points were already discussed in Samsonov et al. (2013b) and shown to reach during 2008-2011 nearly 10 cm.In Fig. 3 it can be seen that the spatial extent of deformation increases nearly linearly with time.
Most analyzed interferograms are at least partially incoherent, which contributes to the inaccurate phase unwrapping and reduction in precision during SBAS and MS-BAS processing.Therefore, for modeling and inversion purposes instead of using linear deformation rates that can be of moderate precision due to severe decorrelation, we selected only two highly coherent interferograms.The ascending MF3F 20091115-20100408 (Fig. 2c) and descending S3 20091205-20100428 (Fig. 2d) interferograms show very high coherence probably due to particular ground conditions, partially frozen soil with no snow cover or foliage.

Modeling approach
The shape of the deformation pattern (ascending and descending interferograms and MSBAS) indicates that the causative source of the deformation could involve a localized volume at some depth.In addition, the approximately circular deformation pattern also could be consistent with volume/pressure change inside a source with axisymmetric geometry.Now, to interpret quantitatively the deformation process responsible for the observed ground displacements, we consider a simplified representation of the subsurface beneath the Rice Lake area.We select an idealized earth structure.First, we assume that the earth is flat, which is a good approximation for the low relief southern Saskatchewan region.Second, the earth is treated as an elastic body, for which its elastic properties are homogenous (uniform across the media) and isotropic (no dependence on orientation).Therefore, for the modeling interpretation we consider only uniform elastic isotropic half-space deformation source solutions.
Two different analytical solutions are tested, a spherical and a sill (lens-like) source.A spherical point source at depth can simulate ground displacements caused by volume/pressure changes of magma or fluids (Mogi, 1958;Samsonov et al., 2010).A second solution calculates surface deformation due to hydrostatic pressure change inside a horizontal circular fracture ("penny-shaped") in an elastic, homogeneous half-space (Sun, 1969;Fialko et al., 2001).Both models could be appropriate representation of the unknown deformation process at depth.Using one of the models, we can express the relationship between surface deformation and source parameters as where d indicates the surface displacement data, the observation error, m the source parameters, and G is Green's function that relates d to m.We adopt a non-linear global inversion technique, bounded simulated annealing, to calculate the best-fitting model m, which finds the minima of the misfit function in a least squares sense (González et al., 2010).All the interferograms were downsampled using a quadtree scheme (Jonsson et al., 2002) to reduce the computational burden while preserving the observed deformation patterns.To avoid solutions trapped in local minima, we re-started the non-linear inversion 100 times with a homogeneously distributed random set of initial model parameters, with the same bounding limits as the model parameters.We used an improved version of the method presented in González et al. (2010), which introduces realistic noise at each iteration (González et al., 2013).We simulated the random noise using an exponential and cosine theoretical covariance function, where r is the distance between data points, σ 2 the phase noise variance, and a and b coefficients that control the correlation lag.The correlation parameters were estimated for each selected interferogram away from the deformation zone.Finally, we determined the best-fitting model parameters as the mean of the 100 best-fitting set of modeled parameters.

Results
For the modeling, we select the most coherent interferograms from the ascending and descending passes.In ascending pass, we select an interferogram spanning the period 20091115-20100408 from the Multi-Look Fine beam mode (MF3F).The unwrapped phase is converted to ground displacement and scaled to corresponding annual rates (cm yr −1 ).Then, the data are spatially reduced but keep the statistically significant part of the deformation signal using a quadtree algorithm.The quadtree algorithm reduces the ascending interferogram to 361 data points.This reduced set is used as input for the inversion.During the inversion, noise is simulated based on the estimated parameters on the full-resolution interferogram far from the deformation area, with σ = 0.50871cm yr −1 , a = 0.41988 km −1 , and b = 4.8e −8 km −1 .Before inversion and quadtree spatial reduction, a residual orbital trend is estimated using a robust regression algorithm to reduce the effect of the presence of strong localized deformation.We model the two subsidence signals east of Rice Lake by estimating the source parameters of two spherical sources.Namely, we inverted to obtain four parameters for each source model (latitude, longitude, depth and volume change).All parameters could vary within reasonably large bounds, for horizontal coordinates ±1.5 km around the maximum deformation point, depth could vary from 0 to 10 km deep, and for volume change ( V ) a range of values from 0 to −2 km 3 are permitted.The results of the modeled ascending interferogram are presented in Fig. 4.
A second model is based on two sill sources, for which we inverted to obtain five parameters for each source model (latitude, longitude, depth, radius and volume change).Again, all parameters could vary within reasonably large bounds.For horizontal coordinates ±1.5 km around the maximum deformation point, depth could vary from 0 to 5 km deep, radius between 0.1 and 2 km, and for volume change ( V ) values from 0 to −2 km 3 are permitted.The results of the modeled ascending interferogram are presented in Fig. 5.
For the descending pass, we select an interferogram spanning a similar period as the ascending pair.The selected period is 20091205-20100428 from the standard beam mode  (S3).The unwrapped phase is also converted to ground displacement and scaled to annual rates (cm yr −1 ).In addition, the data are spatially reduced using the quadtree algorithm leaving 385 data points as input for the inversion.During the inversion, noise is simulated based on the estimated parameters on the full-resolution interferogram far from the deformation area, with σ = 0.60259cm yr −1 , a = 0.024478km −1 , and b = 0.06812km −1 .One can note that the noise correlation is stronger in the descending interferogram than in the ascending one.
A set of two models composed by two spherical and two sill sources are estimated using the descending data.Bounds on the estimated model parameters are the same as for the case of the ascending data.Results of the best-fitting models are presented in Figs. 6 and 7 for spherical and sill source models, respectively.
Finally, we jointly invert for the very same set of models (two spherical and two sill sources) using the ascending and descending data simultaneously.Spherical source model results are presented in Fig. 8, and the results for the corresponding sill source model are shown in Fig. 9.In addition, numerical values for the best-fitting source models are presented in Table 2 In order to compare model performances objectively, we applied F test statistics (Stein and Gordon, 1984;González and Fernández, 2011) to the three inversion scenarios (ascending, descending and joint).To pass the test, the empirically computed F ratio (F e ) should be larger than the theoretical F ratio (F t ) computed based on F probability density distribution for the appropriate number of degrees of freedom.For all scenarios we observed that improvement is statistically significant at a 99 % confidence level (F e > F t ): ascending F e = 48.9> F t = 4.7, descending F e = 81.9> F t = 4.7; joint F e = 402.7 > F t = 4.6.Therefore, the sill model is statistically preferable to the spherical model.Although the actual shape of underground mine is unknown, it seems logical to assume that its vertical dimension (height) is smaller than the horizontal dimension (radius or width and length), resembling a sill-shaped model.

Discussion
Modeling and inversion of InSAR measurements allowed us to determine parameters of the sources responsible for the observed subsidence.Both spherical and disk (sill) models suggested deep sources, but the spherical model produced consistently shallower sources located at about 0.7-1.0km, while sill model suggested sources located at 1.1-1.4km depth.The variability of depth estimation during the Monte Carlo simulations was insignificant, excluding a few outliers, and observed discrepancy is believed to be purely model related.The smallest root mean squared error (RMSE) of 0.7 cm yr −1 was achieved in case of a sill model and ascending data.From Fig. 2c-d it is apparent that precision Table 2. Best-fit InSAR inversion results for two shallow spherical models, and two shallow sill models (source type).Models are constrained using as input data: subsets of ascending, descending and combined.E stands for east source, and W for west source.Longitude and latitude are given in decimal degrees, depth and radius in km, and volume change ( V ) in km 3 .Model misfit as root mean squared error (RMSE) in cm yr −1 .and resolution of the ascending MF3F interferogram used for modeling is significantly higher than precision and resolution of the descending data.Therefore, dependence of RMSE on data type is expected.When both ascending and descending interferograms were used, the rms values were slightly higher (1.0-2.1 cm yr −1 ) than for each data set independently, but no outliers were produced (Fig. 10), which suggests the robustness and stability of the obtained parameters.Due to the high quality of InSAR data used for modeling and inversion, precision of determined parameters estimated as a standard  deviation was high, and it was degraded to a reasonable from practical point of view values.For example, source depth was determined with precision of a few meters, and it was degraded to 0.1 km (see Table 2).Time series of subsidence presented in Figs.2a, b and 3, which span nearly 3 yr of InSAR observation during 2008-2011, do not show any seasonal pattern expected in the case of groundwater-related deformation.We do not observe increase of subsidence rate during relatively dry season in the late summer and early fall and groundwater recharge due to snowmelt in spring.The steady motion rate revealed by InSAR and substantial source depth produced by modeling suggest that the cause of subsidence is likely not related to groundwater withdrawal as was suggested in our preliminary study.Similarly gas-and oil-extraction-related causes can also be excluded since these natural resources are not available in this region (Information Services Corporation of Saskatchewan, 2013).An alternative hypothesis might be due to active mining operations.

Data set
Potash deposits are spread throughout this region and extensively mined.Saskatchewan's rich potash deposits are found in the prairie formation, which lies at depths greater than 1000 m beneath much of southern Saskatchewan (http: //www.ir.gov.sk.ca/Potash).In particular, around the Rice Lake area the two most important extraction levels are the Patience Lake and Belle Plain members with up to 30 and 20 m in thickness, respectively (Saskatchewan Geological Survey, 2013b).As reported by the Saskatchewan Geological Survey (2013a), two active potash mining sites are located in the broad area of Rice Lake.Moreover, it is reported that mining depth at the potash mining site located a few kilometers east of the observed subsidence is 1021 m (Potash-Corp, 2013), and the mining method is conventional underground.Since the mining depth closely matches the source depth determined by modeling and inversion, we speculate that observed subsidence is likely mining related.Removal of potash by mining produces redistribution of stress in the overlying media, which in turn manifest as subsidence at the surface.However, none of the mines are directly located at the deformation areas but at 8 km south and 12 km east, respectively.Therefore, we cannot be certain that these mines are the actual cause of the deformation.Potash (potash salt) reserves in Saskatchewan are some of the largest in the world.Some estimates suggest that existing reserves would be able to supply the world demand for 100 yr at current levels.Therefore if the origin is the potash mining, the monitoring of the ground deformation could be a potentially useful tool for its safe operation.
Figure 11 shows the topographic relief and the land cover of the area undergoing subsidence.It can be seen that this agricultural region is flat with only a few meters of relief.The lowest elevation feature is co-located with the center of the western subsidence, and a moderate elevation gradient is observed in the east-west direction.This region is also characterized by the presence of multiple water bodies recharged by precipitation, snowmelt and shallow groundwater.Subsidence ongoing for a sufficiently long period of time will produce redistribution of shallow groundwater and surface water towards the centers of subsidence, and the relative groundwater level will become shallower.These effects in turn could increase susceptibility of the region to seasonal and weatherrelated flooding.

Conclusions
Ground subsidence observed with RADARSAT-2 InSAR during 2008-2011 is studied in detail.Time-series analysis suggests that the subsidence rate and the spatial extent expansion are approximately constant in time.The modeling and inversion revealed two sources located at depths 0.7-1.0km assuming a spherical model and 1.1-1.4km assuming a sill model.These depths match the reported depth of potash mining operations active in the region and located a few kilometers east of the observed subsidence.Therefore, the mining-related cause of subsidence is considered most likely.If the origin is the potash mining, monitoring of the ground deformation could be a potentially useful tool for its safe operation.Assessment of surface topography and distribution of surface water bodies suggest that subsidence occurring for a long period of time may increase susceptibility of this region to seasonal and weather-related flooding.

Fig. 1 .
Fig. 1.RADARSAT-2 Multi-Look Fine 3 Far (MF3F) and Standard 3 (S3) data sets used in this study are outlined in black.Region of interest (with coordinates longitude 252.88-253.025and latitude 52.04-52.12degrees) is outlined in brown.Selected roads are shown as black thin lines and lakes as blue polygons.In background GeoBase (2013) 20 m resolution digital elevation model of southern Saskatchewan, Canada.

Table 1 .
RADARSAT-2 Multi-Look Fine 3 Far (MF3F) and Standard 3 (S3) data used in this study: time span (in YYYYMMDD format), azimuth α and incidence θ angles, number of available SAR images N, and number of used interferograms M.

Fig. 2 .
Fig. 2. Cumulative line-of-sight deformation time series calculated from RADARSAT-2 MF3F (a) and S3 (b) beams with SBAS technique (Samsonov et al., 2011b) for points of maximum subsidence.Vertical dotted lines in (a)-(b) mark SAR images, for which wrapped differential interferograms are shown in (c)-(d).Due to particular ground conditions during acquisition of these images, corresponding interferograms display superior coherence.Region of interest with longitude 252.88-253.025and latitude 52.04-52.12degrees is shown.

Fig. 3 .
Fig. 3. Cumulative vertical ground displacements calculated from RADARSAT-2 MF3F and S3 data with MSBAS technique (Samsonov and d'Oreye, 2012) for period 14 December 2008-21 August 2011 for region of interest with coordinates longitude 252.88-253.025and latitude 52.04-52.12degrees.Displacements were clipped for clarity to maximum values [−20, 20] cm; color scale is shown in first image.Selected roads are shown as black thin lines and lakes as blue polygons.

Fig. 4 .
Fig. 4. (a) Spatial subset of observed ground deformation rate in ascending interferogram (20091115-20100408); (b) simulated deformation rate based on best-fitting parameters for two spherical sources; (c) residual between (a) and (b); (d) profile between points A and B (roughly west-east) across two subsiding areas, in black observed, in red simulated model and in green residual.

Fig. 5 .
Fig. 5. (a) Spatial subset of observed ground deformation rate in ascending interferogram (20091115-20100408); (b) simulated deformation rate based on best-fitting parameters for two sill sources; (c) residual between (a) and (b); (d) profile between points A and B (roughly west-east) across two subsiding areas, in black observed, in red simulated model and in green residual.

Fig. 6 .
Fig. 6.(a) Spatial subset of observed ground deformation rate in descending interferogram (20091205-20100428); (b) simulated deformation rate based on best-fitting parameters for two spherical sources; (c) residual between (a) and (b); (d) profile between points A and B (roughly west-east) across two subsiding areas, in black observed, in red simulated model and in green residual.

Fig. 7 .
Fig. 7. (a) Spatial subset of observed ground deformation rate in descending interferogram (20091205-20100428); (b) simulated deformation rate based on best-fitting parameters for two sill sources; (c) residual between (a) and (b); (d) profile between points A and B (roughly west-east) across two subsiding areas, in black observed, in red simulated model and in green residual.

Fig. 8 .
Fig. 8. (a) Ascending deformation rate (20091115-20100408); (b) simulated deformation rate based on joint inverted best-fitting parameters for two spherical sources; (c) residual between (a) and (b); (d) descending deformation rate (20091205-20100428); (e) simulated deformation rate based on joint inverted best-fitting parameters for two spherical sources; (f) residual between (d) and (e); (g) and (h) profiles A-B across two deformation areas for the ascending and descending, respectively (black -observed, red -model and green -residual).

Fig. 9 .
Fig. 9. (a) Ascending deformation rate (20091115-20100408); (b) simulated deformation rate based on joint inverted best-fitting parameters for two sill sources; (c) residual between (a) and (b); (d) descending deformation rate (20091205-20100428); (e) simulated deformation rate based on joint inverted best-fitting parameters for two sill sources; (f) residual between (d) and (e); (g) and (h) profiles A-B across two deformation areas for the ascending and descending, respectively (black -observed, red -model and green -residual).

Fig. 10 .
Fig. 10.Profiles across line A-B showing best-fitting source parameters of the 100 inverted sources for each data set.Spherical sources are indicated with orange circles, while sill sources are red oblate ellipses.(a) Best-fitting sources constrained with ascending data; (b) with descending data, and (c) joint ascending and descending inversion.Note that all best-fitting sources are distributed over a very small area (indicating that interferogram noise was indeed relatively small), except for some outliers in panels (a) and (b).

Fig. 11 .
Fig. 11.Relief map (a) and 7 May 2009 MF3F intensity (b) with 2 cm contours of subsidence measured during 14 December 2008-21 August 2011 plotted over.Subsidence coincides with low relief agricultural land susceptible to flooding.