A model-backfeed deformation estimation method for revealing 20-year surface dynamics of the Groningen gas field using multi-platform SAR imagery International Journal of Applied Earth Observations and Geoinformation

The Groningen gas field, the largest natural gas field in Europe, was discovered in 1959 and started its production in 1963. The Earth surface above it experienced subsidence over the past six decades because of gas extraction activities. To accurately reveal this surface movement with satellite SAR data, our study first proposes and demonstrates a model-backfeed (MBF) InSAR deformation estimation method to improve InSAR deformation time series modeling. This method allowed us to include a priori knowledge and to iteratively optimize functional and stochastic models. Using this method and employing a spatio-temporal SAR data integration method based upon Monte Carlo and Multiple Hypothesis Testing methods, we retrieved the 20-year subsidence history of the Groningen gas field by integrating 32 ERS-1/2, 68 Envisat and 82 Radarsat-2 SAR images. The results show that the maximum cumulative surface subsidence in this gas field has been as much as 15 . 5 cm between 1995 and 2015. In terms of precision and accuracy, our MBF method offered a better result than the standard Multi-Temporal InSAR analysis method: the Ensemble Coherence increased by 10% – 19% and Spatio-Temporal Con- sistency decreased by 2% – 20%. In terms of accuracy, our results better concur with the external GNSS reference observations. We further show that the spatio-temporal SAR data integration method has better links with multi- platform SAR data if the uncertainties of the InSAR geolocation and temporal deformation are included. The study demonstrates that the MBF method optimized the estimation of deformation parameters and mitigated unwrapping errors.


Introduction
The Groningen gas field was discovered in 1959. At that time it was the largest natural gas field in the world. Its discovery was followed by an increasing gas demand from the Netherlands, Germany, Belgium and Northern France since 1960 (Richter et al., 2020). During its six decades of exploration, the Groningen surface area has shown a gradual continuous subsidence of approximately, 7.5 mm yr − 1 between 1970 and 2003 (Schoustra, 2004), 7 mm yr − 1 between 2003 and 2007 (Ketelaar, 2008) and 8 mm yr − 1 between 2010 and 2015 (Van Eijs and Van der Wal, 2017). Mulder and Perey (2018) and Bourne and Oates (2020) revealed the driving mechanisms for these surface subsidence in this gas field. They hypothesize that surface deformation responds to a decrease in pore pressure when gas is pumped from the sandstone layer. As the decreased pressure can not support the weight of the layers on top, it results in surface subsidence that compresses the underlying layers. As subsidence pertaining to gas extraction activities could be long-lasting, it is important to reveal the gradual subsidence in the long term.
In practice, a priori knowledge on signal smoothness and temporal behavior of ground targets can facilitate Persistent Scatterer Interferometry (PSI) processing. For instance, signal smoothness in space is often postulated, if InSAR observations are spatially correlated (Ferretti et al., 2001;Kampes, 2005). A deformation model is constructed when temporal behavior of targets is (partially) known (Hooper et al., 2004;. Apparently, the targets are untraceable if they do not hold such a postulation and/or do not apply to a pre-defined model (Aobpaet et al., 2013). Sousa et al. (2010) showed that accurate deformation modelling secures the quality of PSI results. As a follow-up, realistic and often complicated deformation models have been proposed and tested. For example, Zhang et al. (2012) and Hu et al. (2019) have included the linear deformation model into their MT-InSAR methodologies, while Chaussard et al. (2014) retrieved land subsidence with significant linear trends in central Mexico. The deformation component related to thermal expansion, precipitation and seasonal variation has been parameterized and demonstrated by , Ma et al. (2019), Özer et al. (2019), Zhao et al. (2021), while Van Leijen (2014) introduced a single-breakpoint model for earthquake applications. Many of these studies are case driven and therefore selection of deformation model has not yet been automated. The aim of this study is to automate the deformation modelling for each PS. It does so by developing an algorithm that embraces a range of linear and nonlinear models for testing and systematically identifies the best deformation model.
In the past years, MT-InSAR has been used to monitor the subsidence over the Groningen gas field. Ketelaar (2008) used the multi-track ERS and single-track Envisat SAR data to retrieve the maximum vertical subsidence with 8 mm yr − 1 during 1992-2003 and 7 mm yr − 1 during 2003-2007. Beers (2018 identified 6 mm yr − 1 maximum deformation velocity in satellite Line-Of-Sight (LOS) direction during 2009-2016 using Radarsat-2 SAR data. Van Eijs and Van der Wal (2017) presented an InSAR subsidence map with an 8 mm yr − 1 maximum deformation velocity in the vertical direction during 2010-2012 and 2013-2015, also using Radarsat-2 SAR data. To unravel the long-term (e.g. up to two decades) surface dynamics in this field, multi-platform SAR data integration is recommended, but it was not yet exploited in those studies.
This study focused on developing and demonstrating a modelbackfeed InSAR deformation estimation method, the so-called MBF method, to generate more reliable deformation time series estimations. It also further expands and demonstrates a spatio-temporal SAR integration method for revealing 20-year surface dynamics in the Groningen gas field using the 32 ERS-1/2 (1995( ), 68 Envisat (2003( -2010 and 82 Radarsat-2 (2009-2015) SAR images. The paper is structured as follows. Section 2 describes the methods and procedures of modelbackfeed InSAR deformation estimation, spatio-temporal data integration and quality control. Section 3 introduces the geology of the Groningen gas field and the data used. Section 4 displays the procedure of data processing for this study. Section 5 is devoted to the experimental results, followed by the discussion and the conclusion in Section 6 and Section 7. The workflow of the standard SAR data processing, model-backfeed deformation estimation method, spatio-temporal data integration and quality control.

Methods
This section first reviews a standard SAR data processing in Section 2.1 and then introduces the model-backfeed InSAR deformation estimation method in Section 2.2. Next, a spatio-temporal SAR data integration method is introduced in Section 2.3 and Section 2.4. Section 2.5 presents the method to evaluate the precision and accuracy of the results.

Standard SAR data processing
We depict the processing workflow in Fig. 1 and start with reviewing the standard SAR data processing. For every SAR mission, one can separately conduct interferogram generation processing using a conventional Differential InSAR technique (DInSAR). These interferograms from different SAR missions are imported separately into an MT-InSAR processing tool, such as the Delft implementation of PSI (DePSI), for InSAR time series analysis (Van Leijen, 2014).
The Persistent Scatterers (PS), including Persistent Scatterer Candidates (PSC) and Persistent Scatterer Potential Points (PSP), can be selected by thresholding the Amplitude Dispersion Index (ADI) (Ferretti et al., 2001). As traditionally applied in geodetic networks, PSC as the first-order PS are assumed to be of the highest precision with a relatively low ADI (e.g., ADI ⩽0.2). These PSC form a sparse initial network, which models and estimates errors from, e.g., the Atmospheric Phase Screen (APS) and orbital uncertainty. The number of PSC usually accounts for less than 5% of the total amount of PS (PSC and PSP) for similar cases as ours. By definition, the PSP are less stable over time with a slightly larger ADI (e.g. ⩽0.4). The PSP are used to form the densification network during the PSP densification step.
In the densification step, every PSP is connected to a specific number of the nearby PSC for temporal and spatial unwrapping. Hence, the PSP densification depends upon the setup of the initial PSC network. And because of the high proportion PSP among all PS (e.g., about 95% for our case), its quality influences the quality and quantity of the final PS results in, for instance, PS time series. For temporal unwrapping, the temporal behavior can be pre-defined and assumed to be linear. It could be inappropriate, however, for relatively low-quality PSP, which are more likely to have a complicated nonlinear temporal behavior. For spatial unwrapping, a PSP is linked with multiple PSC (e.g., three) to build the double arc phase. If the phase closure of the double arc phase loop is not zero, the PSP in the loop is treated as an alleged low-quality PSP and it will be discarded. As such, the PSP quality is secured, whereas the density of PSP is decreased. This may limit the practicality of PSI, especially in rural areas. To make a trade-off between the quantity and quality of PSP, one can decide to reduce the number of arcs between PSP and PSC by defining a smaller arc total number threshold (e.g., one or two), even though the correctness of unwrapping results could be reduced. Details about DePSI are found in Van Leijen (2014).

Model-backfeed InSAR deformation estimation
The MBF method includes two steps. On a point-by-point basis, it first applies post-InSAR time series analysis results to detect the best deformation model and to obtain the corresponding posterior variance (see Section 2.2.1). Next, it reconstructs the functional and stochastic models and reintroduces these models back to phase unwrapping stage for re-estimation and adjustment of the deformation, topographic, atmospheric and noise phase parameters (see Section 2.2.2). As the best deformation models are inversely fed into the phase unwrapping, we name it the MBF method.

Best deformation model determination of PS
In ths past, a library of canonical functions has been developed to supply a series of physically realistic deformation models. These are helpful to detect the most probable deformation model among all options and well describes the movement pattern of every PS . These potential deformation models consider real world influential factors like, e.g., temperature, periodic signal and a sudden event. Here we illustrate five representative canonical deformation functions: linear function of time, temperature-related function, single breakpoint function, Heaviside function and periodic function. One can combine these five functions to form all potential deformation models, as shown in Table 1 and Appendix A. Assuming that the default deformation model M 0 is a linear function of time, other potential deformation models can be treated as the alternative models, denoted as M J (J ∈ [1, 2, 3, 4, 5, 6]) (cf. Table 1

Pointwise update for PS functional and stochastic model and deformation re-estimation
Using the detected best deformation model M B and the corresponding posterior variance σ B , the functional and stochastic model are updated on a point-by-point basis. For a SAR data stack with acquisition period is T, the k th interferogram formed by two SLCs at t 1 and t 2 (t 2 >t 1 and t 1 ,t 2 ∈ T), an arc phase (ϕ k pq ) formed by two pixels is named as double-difference phase. The p indicates PSC registry and q indicates PSP registry. The functional model of the double-difference phase is written as where unw{} is the unwrapping operator, a k pq is the integer ambiguity of arc phase (a k pq ∈ Z), ϕ k pq,H is the phase of relative height difference, ϕ k pq,defo is the phase of relative displacement, ϕ k pq,atmo is the phase of relative atmospheric delay, ϕ k pq,n is the relatively unmodeled noise phase (e.g., thermal noise, coregistration errors), λ is the wavelength, B k ⊥,q is Table 1 The composition of the potential models. The binary index 1/0 indicates either the canonical function included or not.

Canonical functions Models
the perpendicular baseline, θ q is the local incidence angle, R q is the slant range, ΔH k pq is the relatively residual height, s is the master atmospheric delay in meters. Further, M k pq,B is the selected best deformation model from all possible models listed in Eq. (5), in which A B and x B are the corresponding design matrix and unknown parameters, respectively. The corresponding stochastic model of the double-difference phase is mainly shaped by the unmodeled atmospheric noise, deformation noise and scattering/thermal noise and is expressed as where Q pq,atmo , Q pq,defo and Q pq,n are the variance-covariance matrices of unmodeled atmospheric phase noise, deformation noise and thermal and scattering noise, respectively. As PS, e.g., point p and q in this case, are usually the pixel with a high signal-to-noise ratio (SNR) and relatively stable over time. The thermal and scattering noise of PS is marginal especially when the computing double-difference phases. Q pq,atmo and Q pq,defo can be initially pre-defined empirically and are further adjusted according to least-squares estimation. The posterior variance σ B (see Section 2.2.1) can be used to fine-tune Q pq,defo after MHT. After establishing the functional and stochastic models, the ambiguities of the double-difference phase are calculated, leaving the spatial reference and densification networks unchanged. The Integer Least-Squares (ILS), Integer Bootstrapping (IB) and Ambiguity Function (AF) methods are among the most popular methods that are often applied for ambiguity calculation (Van Leijen, 2014). After unwrapping with e.g., the IB method, the absolute (unwrapped) interferometric phase of q in the PSP registry at the kth interferogram, can be obtained by (1) and Φ k p represents the absolute interferometric phase of the PSC point p, which was derived during the unwrapping in the reference network. With this step, the deformation values of all PS are updated and re-estimated.

Spatial integration
The proposed model-backfeed InSAR deformation estimation method in Section 2.2 can offer the improved deformation estimation results. We are therefore now able to obtain a more reliable deformation time series estimation of PS. Following the workflow of Fig. 1, we move to spatial integration, starting with identifying tie point pairs (TPPs) in space. Note that a TPP is composed of multi-PS (e.g., two) captured by different SAR missions that represents a common ground object. We propose a deformation time series post-processing method  to identify TPPs in space. This method considers the difference in spatial resolution, observation direction and polarization, for which the multi-platform SAR data can not be directly aligned. Here we describe the processing procedure in six steps.
Step 1: For reference dataset selection, two selection criteria are highlighted. First, the acquisition period of the reference dataset (T ref ) should be relatively long among the acquisition periods of the SAR datasets used. Using the data in such a period can guarantee that the selected InSAR measurement points are PS points that remain coherent over a long time. Second, T ref should be located in the middle of the acquisition time of multiple SAR datasets.
Step 2: Ground Control Points (GCPs) are identified and then used to align all SAR datasets to the reference dataset. Specifically, multiple GCPs such as isolated lamp poles or houses from all SAR datasets are identified by using geocoded SAR and optical imagery. Using these identified GCPs and the Helmert transformation (Reit, 2009), the nonreference SAR datasets are aligned to the reference dataset.
Step 3: Point movements are aligned to a public spatial reference point (P ref ). Selection of P ref should obey the following conditions: 1) its scattering characteristics are nearly invariant over time, i.e., its ADI value is low; 2) it is preferred to be an isolated ground target, e.g., a lamp pole or a corner reflector, as these can avoid any multi-path reflection effects; 3) its deformation over time is negligible.
Step 4: Geolocation uncertainty of each PS is expressed as an error ellipsoid (Dheenathayalan et al., 2016;Dheenathayalan et al., 2018). The size of the error ellipsoid is determined by the length of the semiazimuth, the semi-range and the semi-cross-range. These can be parameterized by the variance-covariance matrix of the geopositioning errors, which are due to e.g., azimuth timing delay, solid earth tides and interferometric phase noise. The azimuth and range directions are the satellite flying direction and radar signal direction, respectively, while the cross-range direction is orthogonal to the azimuth-range plane. Estimation of the error ellipsoid for each PS can be realized by either corner reflector experiments or empirical studies (Yang et al., 2019;Chang et al., 2020;Zhang et al., 2021).
Step 5: Monte Carlo methods are used to identify the counterpart(s) from adjacent datasets of a PS from the reference dataset determined in Step 1. If the error ellipsoids of the PS and counterpart(s) have a positive cross volume, they are considered as a TPP. If an error ellipsoid of one PS point from one SAR dataset overlaps with multiple error ellipsoids of the corresponding PS points from another SAR dataset, we use the cross volumes as weights to lump multiple PS points together as a single equivalent PS point.
Step 6: The public TPPs are determined by screening the permanent existence of the TPPs in all SAR datasets that were identified in Step 5 in a chainlike way. For instance, with five SAR datasets, if a PS from the reference dataset (e.g., the third SAR dataset) is part of two TPPs with the two adjacent SAR datasets (e.g., the second and fourth SAR datasets) and the two TPPs link to another two TPPs from following two adjacent SAR datasets (e.g., the first and fifth SAR datasets), then the PS and the corresponding TPPs serve as a public TPP in a long-term time series. If, however, the PS from the reference dataset can not form TPPs from all SAR datasets, then the PS and the corresponding TPPs are not considered as public TPP candidates and are eventually excluded. The public TPPs thus selected inherit the attributes of the PS from the reference dataset.
After running Step 6 for every PS in the reference dataset, the final public TTPs are determined. Heavy computational burdens of the Monte Carlo methods are overcome by narrowing down the search space (e.g., the spatial distance of a potential TPP is less than 300 m) and using parallel computing.

Temporal integration
We next concatenate the deformation time series produced in Section 2.2 and the identified TPPs in Section 2.3 in order to achieve the longterm deformation time series (see Fig. 1). The first epoch in every SAR mission is regarded as the temporal reference, where the best deformation model of each TTP in every T is determined using the MHT method. Using the best model, the deformation time series of each TTP is concatenated. Specifically, if there is a time gap between two adjacent SAR missions, the best deformation time series model of the former SAR data are used to predict the starting point of the corresponding counterpart from the latter SAR data. Besides, if there is a time overlap between two adjacent SAR data, we average the deformation observations from the former and latter SAR data within the time overlap period.

Quality control: precision and accuracy
Following the workflow in Fig. 1, we now move to the quality control part. The Ensemble Coherence (EC) is an indicator to evaluate the consistency between the estimations and the observations, ranging between 0 and 1 (Van Leijen, 2014). When EC approaches 1, it implies that the model can well describe the observations. The EC (γ) is formulated as where m is the number of SAR images and hence m − 1 is the number of interferograms. Further, ϕ k pq is the double-difference phase observation between p from PSC registry and q from PSP registry from the kth interferogram, ϕ k pq,model is the corresponding model phase, e k pq is the double-difference phase residuals, j is the imaginary operator.
As compared with γ depends upon the temporal model, the Spatio-Temporal Consistency (STC) is an indicator that presents the temporal coherence of PS of interest and the spatial coherence between PS of interest and the PS points in its vicinity (Hanssen et al., 2008). The STC of any PS point of interest q is defined by the minimum Root-Mean-Square Error (ρ q ) and is expressed as where q is the PS under consideration, P is the group of PS in the vicinity, λ is the radar wavelength, m is the number of SAR images, k indicates the kth interferometric image and Φ is the unwrapped interferometric phase. The PS with higher spatio-temporal coherence usually has an STC value approaching 0. Regarding accuracy assessment, external geodetic measurements, such as the GNSS measurements from multiple GNSS stations (see Section 3), can be used to validate the MBF-derived deformation estimation and the concatenated long-term deformation time series. The Root-Mean-Square Deviation (RMSD) is used as an indicator for accuracy assessment.

Geological setting
The Groningen gas field reservoir was developed in the Upper Rotliegend Group of the Early Permian age. The total thickness of the Rotliegend Group in the Groningen field ranges from approximately 100 m in the south-southeast to approximately 300 m in the northnorthwest. From east to west, its thickness is fairly uniform (NLOG, 2018). The reservoir itself is located in a sandstone layer at depths between 2600 m and 3200 m (De Jager and Visser, 2017). There are over 1100 extensional faults within the reservoir (displayed in Fig. 2(b)). The main fault orientation is NNW-SSE, while also E-W and N-S oriented faults are also present. The highest density of faults is in the southern sector of the reservoir. Most of the faulting on the reservoir faults occurred during the late Jurassic to early Cretaceous rifting phase. Over the past 27 years, more than 1500 earthquakes have happened in the area. The depth of the event hypocentre is mainly 3 km (Spetzler and Dost, 2017), coinciding with the depth of the Groningen reservoir layer. Besides, a strong correlation was observed between the earthquake locations and the main NNE-SSE oriented faults (see Fig. 2(b)).

Data used
Images from three space-borne SAR missions, ERS-1/2, Envisat and Radarsat-2 were collected for mapping subsidence dynamics over the Groningen gas field. Their spatial coverage and acquisition dates are shown in Figs. 2(a) and 2c) and their basic parameters are displayed in Table 2. GNSS measurements from 15 continuous GPS (CGPS) stations were collected from the Nevada Geodetic Laboratory (Blewitt et al., 2018). Their spatial location is marked in Fig. 7c) and their basic information is summarized in Table 8 in Appendix B.

Data processing
We applied a conventional MT-InSAR technique separately for 32 ERS-1/2, 68 Envisat and 82 Radarsat-2 SAR data. To avoid spatial and temporal decorrelation, we selected the master image for ERS-1/2, Envisat and Radarsat-2 acquired on 24 th August 1997, 18 th February 2007 and 23 rd May 2012, respectively. We kept all perpendicular baselines below the nominal critical baseline (Colesanti et al., 2000), for instance, 1.1 km for ERS-1/2. Fig. 3 displays the temporal and perpendicular baseline configurations of the SAR data. We separately imported the interferograms produced from the three missions into DePSI to carry out time series analysis.
Next, we took the initial deformation time series results using the DePSI tool (shown in Section 5.1) as an input and employed the MHT method with a defined a priori variance of 15 mm 2 for C band SAR data to determine the best deformation time series model and estimate the corresponding posterior variance on a point-by-point basis. Using these detected optimal results, we updated the functional and stochastic model pointwisely using Eqs. (1) and (2) and carried out deformation re-estimation.
We next integrated the data in time and space. Among our three SAR missions, Envisat was selected as the reference data. We then selected eight GCPs to separately conduct the Helmert transformation between the other two SAR missions and the Envisat mission. Next, the public spatial reference point was selected and presented as the red triangle in Fig. 5. As the corner reflector experiment was lacking, we applied the empirical ratio of the semi-axis-length in range, azimuth and cross-range direction for the geopositioning of the error ellipsoid, i.e., 1:5:43 for ERS-1/2, 1:1:15 for Envisat and 1:1:15 for Radarsat-2, respectively. Considering the pixel size of these three SAR missions, we customized the actual size in range, azimuth and cross-range of the error ellipsoids as 10m : 10m : 90m for ERS-1/2, 9m : 9m : 45m for Envisat and 9m : 9m : 45m for Radarsat-2. With this setup, we identified the potential TPPs between every two adjacent SAR data, i.e., ERS-1/2 and Envisat and Envisat and Radarsat-2, within a 300 m search radius. Next, we screened the final public TPPs from these identified TPPs in a chainlike way. For the temporal integration of the final public TPPs, we used the method described in Section 2.4. Specifically, for a public TPP, we estimated the best deformation model of PS from ERS-1/2 by MHT and used the model to predict the starting point of corresponding PS from Envisat at the first epoch. For those PS from Envisat and Radarsat-2, we separately estimated their best deformation models. We then used the two best deformation models to separately predict the deformation values of counterparts for every epoch within the time overlap and averaged them.

Standard SAR data processing
Following the data processing procedures in Section 4, we obtained 31 ERS-1/2 interferograms, 67 Envisat interferograms and 81 Radarsat-2 interferograms. The ADI used for selecting PSC and PSP was separately set equal to 0.35 and 0.45 for ERS-1/2, to 0.3 and 0.45 for Envisat and to 0.3 and 0.45 for Radarsat-2. The initial deformation time series analysis with DePSI produced 380771, 127743 and 43748 PS from ERS-1/2, Envisat and Radarsat-2, respectively, shown in the last row of Table 3.  Fig. 3. The temporal and perpendicular baseline configuration of (a) ERS-1/2, (b) Envisat and (c) Radarsat-2, respectively. The dots represent SAR images while the lines represent the interferometric pairs. The star in every SAR mission indicates the corresponding master image.

Detected best deformation models
The distribution of the detected best deformation models of all PS from the three SAR missions is listed in Table 3. For instance, it shows that 65.8% of the PS from Radarsat-2 retained the default linear model (M 0 ) as the best deformation model, while 0.01%, 2.4%, 20%, 0.3%, 0.49%, 11% of the PS identified M 1 , M 2 , M 3 , M 4 , M 5 and M 6 as the best model, respectively. Table 3 further indicates that the population density for each model is different for the different SAR missions. This is related to the difference in e.g., SAR acquisition time, observation uncertainty, actual changes in PS temporal behavior, ADI threshold selection, wavelength, orbital direction and a priori variance selection.

Deformation map comparison
Prior to the deformation map comparison for the standard MT-InSAR method (Trad) and the MBF method (MBF), we quantitatively evaluated the PS with relatively large deformation velocities. For instance, for PS with a deformation velocity between − 15 and +15 mm yr − 1 , MBF adjusted 44.4%, 53.9% and 14.9% for ERS-1/2, Envisat and Radarsat-2, respectively. For those, the deformation velocity estimation were better improved, falling into the deformation velocity range [ − 15, + 15] mm yr − 1 . As such, the PS in the velocity range [ − 15, + 15] mm yr − 1 , separately increased from 158414 to 228756 for ERS-1/2, from 65495 to 100790 for Envisat and from 36278 to 41687 for Radarsat-2 (see Fig. 4).
When we acknowledge the maximum vertical subsidence with approx. 8 mm yr − 1 from the literature in Section 1, we see that the maximum subsidence in the LOS direction is about 7.4 mm yr − 1 for ERS-1/2 and Envisat and 6.5 mm yr − 1 for Radarsat-2, for the past years. We therefore safely assumed that the deformation velocity range for all PS equals [− 10, 10] mm yr − 1 , given a 3 mm yr − 1 buffer. In addition, by setting identical thresholds on height and STC values for MBF-and Trad-derived PS results, we further reduced the number of plausible PS points. In doing so, we obtained the final PS for MBF and Trad, i.e., 39318 and 39054 PS for ERS-1/2, 47827 and 43633 PS for Envisat and 33902 and 32684 PS for Radarsat-2, respectively. The final number of PS from each SAR mission using MBF increased, as compared with using Trad. For instance, it increased by 0.7%, 9.6% and 3.7% for ERS-1/2, Envisat and Radarsat-2, respectively. This is related to the number of SAR images and the actual temporal behavior of the detected PS in 1995-2000, 2003-2010 and 2009-2015. In essence, more observations offer better parameter estimations. The deformation velocity maps in the LOS direction of these three SAR missions with MBF and Trad are shown in Fig. 5. The expectation, standard deviation and maximum of their average deformation velocity in the LOS direction are provided in Table 4. We only found marginal differences in PS results derived from MBF and Trad.

Deformation time series comparison
From the final PS, we selected 16 PS from these three SAR missions, which are the spatially closest to the GNSS stations and compared their deformation time series as obtained by MBF and Trad. Fig. 6 shows the deformation time series estimation of six of them, labeled PS1-PS6, with Trad in black and MBF in red. The locations of these six PS are marked in Figs. 5 (b), (d) and (f). The comparison for other ten PS are displayed in Fig. S1 and Table S1 and discussed in the supplementary material. Here we converted the 3D displacement of GNSS into the LOS direction using (Hanssen, 2001), where D u , D n and D e separately represent the GNSS displacements in the up-down, north-south and east-west directions, θ inc and θ azi indicate the local incidence angle and the azimuth angle (measured in clockwise direction from the North). We then compared the deformation estimations of GNSS with the PS counterpart derived from Trad and MBF. For PS1, MBF corrected the phase jump on the 31st epoch (17 th September 2020) and subsequent epochs. PS1's EC value increases from 0.2590 with Trad to 0.2798 with MBF. For PS2, MBF identified the cycle slip on the 5th epoch (11 th February 1996) and corrected this phase unwrapping error. Its EC increases from 0.7273 with Trad to 0.8074 with MBF. For PS3, MBF identified and included a temperature-related model to fit the actual deformation trend well. Its EC increases from 0.4460 with Trad to 0.6140 with MBF. In addition, the deformation time series for PS3 fitted well with the DZYL GNSS observations with a RMSD of 2.63 mm, as compared to that with Trad (RMSD= 9.78 mm). For PS4, MBF identified and corrected two cycle slips on the 28th (16 th July 201) and 47th (14 th October 2012) epochs. Its EC increased from 0.6001 with Trad to 0.6761 with MBF. For PS5 and PS6, MBF separately identified and corrected one and two cycle slips. The EC of PS5 increased from 0.8784 with Trad to 0.8845 with MBF and that of PS6 from 0.8048 with Trad to 0.8992 with MBF. The MBF-derived deformation time series estimations of PS5 coincide with the observation at TENP GNSS station, with a RMSD value equal to 1.79 mm. The RMSD values of Trad-and MBF-derived deformation time series estimations for PS6 with respect to the observations at the TJUC GNSS station were equal to 5.89 mm and 3.65 mm, respectively. It shows that the MBF offers a more accurate estimation for PS6.

Spatio-temporal linking of tie-point pairs
Using the defined error ellipsoids defined in Section 4 and Monte Carlo methods defined in Section 2.3, potential TPPs were identified. The corresponding intersection types are displayed in Table 5. For complex intersection types, we used the cross volume as weights to obtain their equivalent deformation values. The spatial distribution of the identified TPPs is illustrated in Fig. 7. It shows that between ERS-1/2 and Envisat, we identified 17693 TPPs that account for 39% of the number of PS from ERS-1/2. Among these, 13096 (74%) TPPs have an Table 3 The detected best deformation models for all PS from three SAR missions with a priori variance 15 mm 2 for ERS-1/2, Envisat and Radarsat-2.  Fig. 4. The histogram of the PS total number for ERS-1/2, Envisat and Radarsat-2 with the MBF method and the standard MT-InSAR method. Trad is the acronym of Traditional and refers to the standard MT-InSAR method.
error ellipsoid of one PS from ERS that intersects with one error ellipsoid of one PS from Envisat (denoted as 1vs1), 3795 (21.4%) TPPs have one error ellipsoid of one PS from ERS that intersects with two error ellipsoids of two PS from Envisat (denoted as 1vs2), while 712 (4%) TPPs and 100 (0.6%) TPPs have one error ellipsoid from ERS that intersects with three and four ellipsoids of PS from Envisat denoted as 1vs3 and 1vs4, respectively. By associating all TPPs with the Radarsat-2 points, 3654 TPPs were eventually identified within the period of ERS, Envisat and Radarsat-2, thus spanning the full period of 20 years.

20-year time series and accuracy assessment
Every SAR mission has distinct incidence angles (θ inc ), i.e., 23.3 • , 22.75 • and 34.1 • for ERS-1/2, Envisat and Radarsat-2, respectively, resulting in different projections of the actual deformation onto the LOS direction. As SAR data are lacking from an opposed viewing geometry, 3D decomposition is not practically possible. Past studies showed that deformation in the vertical direction is predominant (Ketelaar, 2008;Beers, 2018;Van Eijs and Van der Wal, 2017). Hence, to streamline the LOS observations from all missions, we assumed that the deformation in the Groningen gas field mainly occurred in the vertical direction. Specifically, we projected the LOS deformation (d LOS ) onto the vertical direction (d vertical ) by using d vertical = 1 cosθinc d LOS , where θ inc represents the satellite incidence angle. From the 3654 public TPPs, we took tpp1 and tpp2 as examples to validate the accuracy of their deformation time series results with GNSS data, being the nearest PS to the DZYL and VEEN GNSS stations, respectively. For ttp1, each SAR mission offers one PS to constitute the public tie-point pair tpp1 and the intersection circumstances of their error ellipsoids are displayed in Fig. 8(a). Its concatenated 20-year deformation time series is displayed in Fig. 8(b). From 19 th August 1995 to 14 th May 2006, ttp1 had a slow uplift trend equal to 2.4 mm yr − 1 and its maximum uplift value was 28 mm. After that, it experienced continuous subsidence and its accumulative subsidence value was about 60 mm until June 2015. For tpp2, one PS came from ERS-1/2, two PS from Envisat and one PS from Radarsat-2 form the public tie-point pair tpp2. Intersection circumstances of their error ellipsoids are shown in Fig. 8(c). Because one PS came from Radarsat-2 and one PS from ERS-1/2 intersect with two PS from Envisat, we used the number of voxels representing intersection volumes (i.e., 4 and 368) as weights to group the two PS from Envisat into one PS. Fig. 8(d) shows the concatenated 20-year deformation time series: tpp2 experienced continuous subsidence. Its accumulative subsidence is approximately 95 mm and the average deformation velocity is approximately − 4.75 mm yr − 1 in the vertical direction. We found that the PS results fit well with the GNSS observations with a RMSD of 3.67 mm for tpp1 with respect to the DZYL GNSS station and 4.89 mm for tpp2 with respect to the VEEN GNSS station.

Effectiveness of the MBF method with respect to the standard MT-InSAR method
As we showed in Section 5.3, MBF avoids the underestimation of deformation parameters, thereby reclassifying some PS with an alleged large linear trend as based upon Trad into a reasonable deformation velocity range, e.g., [ − 15, + 15] mm yr − 1 , adjusted from 14.9% to 53.9%. It increased the number of reliable PS from 0.7% to 9.6%. In Section 5.4, we demonstrated that MBF can correct anomalies and improve deformation time series and that the corresponding accuracy of deformation time series also is validated with GNSS data.
In addition, MBF mainly introduced optimized nonlinear models on a point-by-point basis. We selected all PS as detected in Section 5.2 with Table 4 The expectation (i.e., mean), standard deviation (std) and maximum (max) of the average deformation velocity in the LOS direction with MBF and Trad. The unit is mm yr − 1 .   nonlinear models, comparing the performance of MBF with respect to Trad and applied the two indicators mentioned in Section 2.5, i.e., EC and STC. We found that the EC values for MBF are higher than those obtained with Trad (see Fig. 9). After using MBF, EC increased by 10% for ERS-1/2, 19% for Envisat and 10% for Radarsat-2, respectively (see Table 6). It demonstrates that the best deformation model based on MBF is more suitable than the default linear deformation model. Besides, MBF can offer more PS with lower STC, which means fewer deformation estimation errors (see Fig. 9). After using MBF, STC decreased by 3% for ERS-1/2, 20% for Envisat and 2% for Radarsat-2, respectively (see Table 6).

Scenarios influencing the effect of the model-backfeed InSAR deformation estimation
We can imagine four scenarios that influence the effect of the modelbackfeed InSAR deformation estimation method. The first scenario is that the best model is not included in our pre-defined deformation model library. For example, a PS deformation time series may follow the periodic function (such as M 5 and M 6 in Eq. (5)) with inconstant amplitude parameters s, c. When such a periodic function is not included in the library, the determined best deformation model of this PS point is the best in the library, but not the actual best. The second scenario is that a poor-quality initial PSI result may lead to an erroneous MHT output. For instance, if the initial PSI result has pronounced deformation estimation errors (e.g., a phase jump), it could lead to the incorrect deformation model being detected by MHT. The third scenario is that the wrapped and unwrapped phase signal is smeared when the phase noise is too high. The fourth scenario is that the unmodeled atmospheric phase cannot be properly separated from the deformation phase, or structured in the stochastic model. As such, the MHT is less sensitive to detect actual deformation patterns and determine the best deformation model per PS.

Factors influencing the total number of tie-point pairs
We observed five factors that would influence the total number of identified TPPs. The first factor is the limited number of PS from every SAR data. Even though many PS are available within a short period of time, e.g., five years for ERS-1/2 or seven years for Envisat in this research, they do not always exist over the entire period of 20 years. The second factor is the poor quality of the reference alignment that negatively affects the 3D spatial distance calculation of potential TPPs. The third factor is the pre-defined size of the error ellipsoid that can be both oversized or too small. It would cause improper intersection types of TPPs or even the lack of any identified TPPs. The fourth factor is the varying geocoding accuracy of different SAR missions, which can be up to the level of several meters. The fifth factor is that only the geometrical matching can be applied. When finding TPPs for PS acquired from different viewing geometries, physical matching is not applicable in a straightforward way in the absence of external knowledge.

Computational efficiency
During the identification of TPPs, our experiments demonstrated that reducing the search space and applying parallel computing dramatically improved computational efficiency by three to six times. For deformation estimation, we tested the computational efficiency of Trad and MBF with Intel Core i7-8700 3.2 GHz × 12 CPU. The results show that the computational efficiency is related to the number of SAR images and selected PS. The deformation estimation calculation time of MBF increased by 71.35% for ERS-1/2, 57.13% for Envisat and 76.99% for Radarsat-2, respectively, compared with Trad, see Table 7. This is because MBF introduced much more complicated nonlinear models.

Conclusion
This study proposed a model-backfeed InSAR deformation estimation method, the so-called MBF method, that iteratively updates and improves deformation estimations by pointwise detecting the best functional and stochastic model of InSAR measurement points, such as point scatterers. The experiment demonstrates that the MBF method alleviates model imperfection, significantly reduces the deformation Fig. 7. The distribution of the identified tie-point pairs appeared in the period of ERS-1/2∼Envisat (from 1995/08 to 2010/10) (a), Envisat∼Radarsat-2 (from 2003/12 to 2009/07) (b) and ERS-1/2∼Envisat∼Radarsat-2 (from 1995/08 to 2015/06) (c). The white boxes present the locations of tpp1 and tpp2 discussed in Section 5.6.    estimation errors, reclassifies the falsely-categorized deformation anomalies and improves the number of reliable point scatterers. Our spatial data integration method using Monte Carlo methods and introducing geolocation uncertainty also shows its potential to precisely pair point scatterers acquired from different SAR data, without degrading the spatial resolution. We employed these methods to retrieve the long surface dynamics over the Groningen gas field for the past 20 years, with ERS-1/2, Envisat and Radarsat-2 SAR data. For this area, we found that the maximum subsidence velocity was up to 7.83 mm yr − 1 in LOS direction, while several point scatterers displayed nonlinear behavior over time.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. acknowledge Nederlandse Aardoli Maatschappij (Jan de Jager and Clemens Visser) for offering the fault data. We also thank Delft radar group for sharing Doris and DePSI toolbox and dr. Freek van Leijen for his valuable input. We thank the Royal Dutch Meteorological Institute (KNMI) for offering borehole data and Nevada Geodetic Laboratory for offering GNSS data.

Appendix A. Potential deformation models
Apart from listing the potential deformation model composition in Table 1, here we detail those models by showing their functional models.
where v is the linear deformation velocity without an instantaneous event and t is the corresponding temporal baseline, t 1 and t 2 indicate the temporal baseline before and after a single-breakpoint date, while v 1 and v 2 indicate the corresponding linear deformation velocities, η is the temperature related parameter, ΔT is the temperature baseline between slave and master image, L represents the deformation jump magnitude, δt indicates the jump knot in time. Note that the product δt⋅L forms a Heaviside step function (Weisstein, 2002). s and c are the amplitude of periodic deformation and A i and x i (i ∈ [0, 1, 2, 3, 4, 5, 6]) represent the design matrix and unknown parameter(s), respectively. Table 8 lists the names, geographic coordinates with longitude and latitude and observation periods of the 15 GNSS stations that were used for the validation of InSAR-derived deformation time series.