Mitigation of Atmospheric Artefacts in Multi Temporal InSAR: A Review

The complexity of the atmosphere renders the modelling of the atmospheric delay in multi temporal InSAR difficult. This limits the potential of achieving millimetre accuracy of InSAR-derived deformation measurements. In this paper we review advances in tropospheric delay modelling in InSAR, tropospheric correction methods and integration of the correction methods within existing multi temporal algorithms. Furthermore, we investigate ingestion of the correction techniques by different InSAR applications, accuracy performance metrics and uncertainties of InSAR derived measurements attributed to tropospheric delay. Spatiotemporal modelling of atmospheric delay has evolved and can now be regarded as a spatially correlated turbulent delay with varying degree of anisotropy random in time and topographically correlated seasonal stratified delay. Tropospheric corrections methods performance is restricted to a case by case basis and ingestion of these methods by different applications remains limited due to lack of their integration into existing algorithms. Accuracy and uncertainty assessments remain challenging with most studies adopting simple statistical metrics. While advances have been made in tropospheric modelling, challenges remain for the calibration of atmospheric delay due to lack of data or limited resolution and fusion of multiple techniques for optimal performance.

Numerical models of refractivity (Tarayre and Massonnet 1994) Homogenous atmospheric variation due to topography Stratified delay Calculation of RMS variation in interferogram (Goldstein 1995) The main residual in interferograms is due tropospheric delay Turbulent delay Stacking (Zebker et al. 1997) Atmospheric signals are random in time Turbulent delay Ground-based meteorological data (Delacourt et al. 1998) Homogenous atmospheric variation due to topography Stratified delay GNSS (Williams et al. 1998) GNSS and InSAR suffers same delay Total delay Numerical weather models (Shimada 2000) Stratified atmosphere Stratified delay Time series filtering (Alessandro Ferretti et al. 2000) Atmospheric signals are random in time Turbulent delay Phase-Elevation Relation (Beauducel et al. 2000) Homogenous stratified atmosphere Stratified delay Satellite multispectral imagery (Li et al. 2005) Variation is due to wet delay Wet delay Multi scale phase-based (Lin et al. 2010) Different scale of sensitivities of signals present in interferogram

Stratified delay
Power law method (Bekaert et al. 2015a) Power law dependence of the delays Stratified delay aufgrund fehlender Daten oder begrenzter Auflösung und der Fusion mehrerer Techniken für eine optimale Leistung bestehen.

Introduction
Interferometric Synthetic Aperture Radar (InSAR) has proven to be a powerful geodetic space technique allowing measurements and observation over large area with a sub-centimetre accuracy in applications related to anthropogenic activities (Amelung et al. 1999;Juncu et al. 2018;Xu et al. 2017), landslides (Bianchini et al. 2013;Ciampalini et al. 2015;Rott and Nagler 2006) and geophysical processes (Béjar-Pizarro et al. 2013;Biggs et al. 2007;Hooper et al. 2004). This success was augmented by improved data availability with missions like Sentinel-1 from the European Space Agency (ESA) granting free access to data. Similarly introduction and subsequent improvement of multi-temporal InSAR (MTI) algorithms by Ferretti et al. (2000), Berardino et al. (2002) and Hooper et al. (2012) increased the capacity of InSAR to achieve mm accuracy. Advances have been made in mitigating other error sources in InSAR but atmospheric-phase screen (APS) remains a significant challenge in optimal utilization of InSAR (Foster et al. 2013;Cao et al. 2019). Though earlier studies on use of InSAR as geodetic deformation tool reported success, some fringes were reported that could not be associated with any deformation signals, which was later interpreted as atmospheric delay (Tarayre and Massonnet 1996).
Since then research has been geared towards understanding the nature of atmospheric delay in InSAR (Goldstein 1995;Jolivet et al. 2011;Wei et al. 2019) and errors associated with tropospheric delay (Cao et al. 2018; González and Fernández 2011;Hanssen 2001). The atmospheric errors in deformation measurements could blur the deformation signal. Zebker et al. (1997) postulated that a change in relative humidity of 20% could lead to an error of 10-14 cm in deformation measurements while Fattahi and Amelung (2015) reported that a 10 cm amplitude due to seasonality of atmospheric delay could lead to a bias of up to 23.8 cm. To mitigate the effects of tropospheric delay, earlier methods relied on models that had earlier been applied to space-based techniques (Delacourt et al. 1998;Rosen et al. 1996;Tarayre and Massonnet 1996).
On the same account, the delay estimated from GPS was also applied to mitigate delay in InSAR (Williams et al. 1998).
Stacking was later introduced for the differential InSAR techniques (Wright et al. 2001;Zebker et al. 1997). MTI techniques addressed delay through filtering (Ferretti et al. 2000) and later on methods that rely on the phase-elevation relationship (Cavalié et al. 2007). Due to limitations and the complexity of modelling the troposphere, other methods that rely on external data sets, such as numerical weather models (Doin et al. 2009) and spectrometers (Li et al. 2005), were introduced. Table 1 shows the temporal evolution of correction techniques.
Mitigation of atmospheric delay resulted in InSARderived measurements producing consistent results with other co-located geodetic techniques. For instance, Delacourt et al. (1998) reported a decrease of up to 35% in magnitude of deformation consistent with tiltmeter estimates after mitigation of atmospheric delay to 1991-1993 Mt Etna eruption compared to what had been reported by Massonnet et al. (1995). In some cases, previous findings have been invalidated after proper mitigation of delay. Yip et al. (2019) showed that 5.7 cm/year deformation at Agung volcano that had been reported by Chaussard and Amelung (2012) was due to tropospheric delay and not deflation. Similarly, Rémy et al. (2015) reported no deformation had occurred at Llaima volcano as earlier reported by Bathke (2011), as any deformation less than ± 7 cm could not be detectable as this was equivalent to tropospheric delay at the volcano. Moreover, due to limitations associated with each of the mitigation methods, proper selection should be done to get reliable measurements (Fournier et al. 2011). Albino et al. (2020), for instance, report on the September 2017 Agung interferogram that had been previously treated as deformation, yet it was found to be a tropospheric delay. The information shared was misleading to the evacuation process. Reliability of the InSAR-derived measurements is tied to the error sources, tropospheric delay being the major error source, result in bias and uncertainty in the measurements (Fattahi and Amelung 2015). We, therefore, seek to review the modelling of tropospheric delay, discuss existing methods of mitigation of tropospheric delay and how MTI algorithms have treated tropospheric delay together with how different InSAR applications have handled tropospheric delay. We also discuss the performance metrics and the impact of the accuracy of the methods on MTI-derived measurements and finally provide a workflow that can guide in selecting the appropriate technique.
The paper is organized as follows: Sect. 2 outlines atmospheric delay in InSAR, its components, and their spatiotemporal characteristics. In Sect. 3, we discuss tropospheric correction techniques from the literature. Section 4 outlines the integration of correction techniques within the existing MTI algorithms. Section 5 provides an overview of the ingestion of the correction techniques by different MTI applications. In Sect. 6, we discuss the accuracy metrics and uncertainty attributed to tropospheric delay. Section 7 discusses the selection of the appropriate technique(s). Section 8 concludes the study.

Atmospheric Delay in InSAR
Geodetic-based methods that rely on electromagnetic propagation for measurement, such as very long baseline interferometry (VLBI), Global Navigation Satellite Systems (GNSS) and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), suffer from atmospheric delay (Teke et al. 2013). InSAR is subject to the same interference as it is based on the same observation principles. The variability of pressure, temperature, water vapour, clouds and total electrons content induces change in atmospheric refractivity index resulting in signal delay or advance of the electromagnetic signal. The air refractivity N, according to Hanssen (2001), Kursinski et al. (1997), and Smith and Weintraub (1953) can be described as: Fig. 1 a Shows atmospheric delay in InSAR, depicting the two components where θ is the satellites incidence angle. b Shows temporal variability of relative humidity and temperature estimated from ERA-5 reanalysis weather model data for a point in Nairobi P d represents dry air partial pressure, e is water vapour partial pressure, T is Temperature, w cloud water vapour content, TEC is the total electron content while f is the signal frequency. k 1, k 2 and k 3 are constants as determined by Bean & Dutton (1966), Smith &Weintraub (1953), andThayer (1974). The first four terms are due to non-dispersive troposphere while the last term is due to the dispersive nature of ionosphere. The atmospheric delay can thus be treated as a composite of the two distinct ionospheric and tropospheric components as shown in (Fig. 1).

Ionospheric Delay
Ionospheric delay is caused by the dispersive nature of the ionosphere resulting in the advancement of the propagated signal. The signal is advanced proportional to the total electrons in the ionosphere and inversely proportional to the square of the signal frequency, hence more pronounced in low-frequency signals, such as L-band (Gomba et al. 2016a). The advancement of the signal does not depend on the electrons along the path but the total number of electrons, which is dependent on the solar activity, latitude, hour and geomagnetic conditions in the atmosphere (Meyer 2010). The impact of ionospheric advancement on InSAR includes distortion of image geometry and focus, Faraday rotation and azimuthal offsets plus phase screens, which can be mitigated by either range or azimuth spectrum techniques (Fattahi et al. 2017). They stated that a unit ionospheric delay expected on L-band (1.27 GHz) will be of magnitude of 1/3.9 for S-band (2.5 GHz), 1/16.1 for C-band (5.405 GHz) and 1/57.7 for X-band (9.65 GHz), assuming the same ionospheric conditions and observation geometry. Gomba et al (2016b) reported minimal effects of ionospheric interference on high frequencies like Sentinel-1 C-band. Interferences manifested as azimuthal strikes, which are usually overshadowed by tropospheric signals and can be sufficiently mitigated by split spectrum methods. Liang et al. (2019) reported an increase in errors when ionospheric corrections are applied to C-band interferograms at mid-latitudes which they opined that corrections are not necessary. Moreover, registration algorithms, such as enhanced spectral diversity (Prats-Iraola et al. 2012), improve the accuracy of registration, limiting the errors associated with ionospheric delay. Accordingly ionospheric influence, unlike tropospheric delay, does not pose a major challenge to InSAR processing, particularly on high-frequency data which are the most common. We, thus, limit our discussions to tropospheric delay.

Tropospheric Delay
Tropospheric delay is the dominant error source in InSAR as different tropospheric conditions during SAR acquisitions introduce varying delays. Among the variables, water vapour exhibits high temporal variability as shown in Fig. 1. Tropospheric delay is composed of the stratified delay, which is correlated with topography, and turbulent delay attributed to random processes in the atmosphere, which is independent of topography (Hanssen 2001). However, there is an overlap between the two: the so-called wet delay which is a component of stratified delay is due to water vapour variation, which is also the dominant variable in turbulent delay. The two components exhibit different characteristics in space and time as will be discussed in the next sub-sections.

Stratified Delay
The elevation-dependent delay leads to points at lower elevation experiencing stronger delays compared to points at higher elevation due to longer interaction time with the troposphere as shown by tropospheric delay map in Fig. 2. The delay can be estimated through integration of Eq. 1 between ground elevation and top of atmosphere, typically a height is chosen above which no change in refractivity is expected. ( z o is the surface elevation while z ref is the elevation above which no change in refractivity is expected. Liquid water constitutes an insignificant delay and due to the dispersive nature of cloud droplets it can be ignored for short-wavelength data (Hanssen 2001;Jung et al. 2014;Saastamoinen 1972). The first two terms of Eq. 2 can be expressed using the equation of state and the relationship between dry and water vapour partial pressures.
where p d is the dry air partial pressure, R d is the specific gas constant for dry air (R d = 287.05 J/kg/K), R v is the specific gas constant for water vapour (R v = 461.495 J/kg/K), is the density of moist air. Under the condition that hydrostatic equilibrium is satisfied the hydrostatic term can be integrated separately as shown in Eq. (5).
P zo is the surface pressure and g m is the average gravitational acceleration within the integration range. Reducing the delay to Eq. (6) as expressed by Davis et al. (1985), Doin et al. (2009), andHanssen (2001) The first term represents hydrostatic delay, which is dependent on surface pressure and constitutes a significant portion of the stratified delay. The remaining terms describe wet delay, which, though of small magnitude, varies rapidly due to spatiotemporal variation of water vapour in the lower section of the troposphere. Most studies have focused on mitigating wet delay as it is thought to be the major contributor to tropospheric delay in InSAR. Though hydrostatic delay is considered invariant within SAR image extends, Elliott et al. (2008) and Jolivet et al. (2014) showed that it should not be neglected particularly in regions with significant topographic variation. Similarly, Doin et al. (2009) argue that it should be considered where the temperature range is 10 degrees Celsius or greater. Moreover, the availability of wide scene images like Sentinel-1 and the ability to concatenate various scenes, which might stretch into different climatic conditions, cause hydrostatic delay to no longer be considered homogenous with a scene as earlier postulated by Hanssen (2001).

Mapping Function
The integration of refractivity is computed vertically and, due to the side-looking nature of SAR, a cosine function is usually used to convert vertically integrated delay into slant range, utilizing the satellite incidence angle ( ).
The slant range conversion assumes homogenous refractivity within the horizontal distance between the signal path and the ground cell as integrated vertically. For example using existing weather models, satellite data acquired with an incidence angle of 40 o results in horizontal distance of 17 km between the ground cell as integrated vertically and the satellite path ( Fig. 1). At this distance lateral tropospheric variation is expected. To account for lateral variation, Hobiger et al. (2010) proposed ray tracing based on the Eikonal equation (Paris and Hurd 1969), which was, however, hindered by computational complexity and the low resolution of weather models. The potential of ray tracing for the correction of other space geodetic techniques has been documented (Haji-Aghajany and Amerian 2018; Hofmeister 2016). Similarly Kinoshita et al. (2013) concluded that the availability of meso-scale weather models necessitate accurate mapping functions to take full advantage of their resolution. Recent studies report ray tracing offering better performance in mitigation of tropospheric delay in InSAR (Haji-Aghajany et al. 2019; Haji-Aghajany and Amerian 2020). Other mapping functions, such as vector-based techniques, based on establishing a common coordinate reference system between the ground cells and satellite have been put forward (Hu and Mallorquí 2019). The slant range conversion is, however, still common due to the lack of high-resolution data able to capture atmospheric variability and due to computational complexity of ray tracing techniques.

Seasonality of Stratified Delay
Tropospheric delay estimation as shown in Eq. (6) relies on tropospheric variables that vary with seasons, which impact the temporal variation of stratified delay. Earlier techniques treated atmospheric delay as random in time and could be mitigated with methods, such as stacking or filtering. However, residuals correlating with topography may remain after the application of these methods. Those residuals could be misinterpreted as deformation despite not being associated with any actual physical process (Jung et al. 2014;Zhu et al. 2017). Global weather models revealed seasonality of stratified delay characterized by periodicities that take the form of inter-annual or annual sinusoidal variation (Doin et al. 2009;Jolivet et al. 2011Jolivet et al. , 2014. This seasonality influenced by both wet and hydrostatic delays as shown Fig. 3 which were estimated from ERA-5 data, the latest atmospheric reanalysis by the European Centre for Medium-Range Weather Forecasts (ECMWF). Wet delay maximum oscillations occurs in summer while minimum oscillations occur in winter (Dong et al. 2019;Samsonov et al. 2014). However, the relationship can be more complex for tropical regions (Doin et al. 2009). The seasonality for the hydrostatic delay is smooth in time, the wet delay though seasonal is still highly variable. Ignoring seasonality of stratified delay could introduce bias in InSAR-derived measurements. Fattahi and Amelung (2015) argue that bias associated with seasonality of stratified delay though linear in amplitude is non-linear in acquisitions time and inverse of the total time coverage. A small number of images cannot capture seasonality; this might be an issue for earlier missions but is usually not a problem for current satellites with short revisit times.

Turbulent Delay
Turbulence refers to delay associated with random processes in the atmosphere introducing variations of the refractivity index that cannot be deterministically modelled by equation of refractivity. Random processes in the atmosphere can be attributed to local weather conditions (Tarayre and Massonnet 1996), strong convective effects (Anber et al. 2014), variations of local land covers and different ecosystems (Mahmood et al. 2014). It affects both flat and mountainous regions as shown in interferogram of Fig. 2 and for proper Fig. 3 Temporal difference of two points along the Kenyan Rift with elevation difference of 1241 m depicting seasonality of stratified delay estimated from ERA-5 reanalysis data fitted with sinusoidal functions mitigation, the turbulent delay needs to be considered (Yu et al. 2018a). Although turbulent delay cannot be deterministically solved, stochastic methods can be used for qualitative data description and adjustment through the use of power spectrum, fractural dimension and structure functions (González and Fernández 2011;Hanssen 2001). Turbulent delay, according to the Kolmogorov turbulence theory by Tatarski (1961), can be estimated by spatial structure functions. The structure function between two points in space for a variable x, can be described as where r o is the reference location, h is the distance from the reference location and the angle brackets show the ensemble average. For a stochastic process the structure function may be described by a power law dependence (Williams et al. 1998) as where C is the scale characteristics of the process, is the spectral index and is the power index. Treuhaft and Lanyi (1987) showed through VLBI measurements that structure functions of turbulent water vapour exhibit power law dependence, several other studies have reported similar pattern in InSAR attributed to turbulent delay (Goldstein 1995;Hanssen 1998;Jonsson 2002;Li et al. 2007). The power law index exhibits different magnitudes at different scales. Hanssen (2001) reported a power index of − 8/3 for wavelengths greater than 2 km, − 5/3 for scale between 0.5-2 km and − 2/3 for scale less than 0.5 km which he attributes to noise. Similar observations were noted by Jonsson (2002). The power spectrum is important in determining the filtering parameter needed to separate various signals (Ferretti et al. 2000), determining the interpolation regime (Williams et al. 1998) and evaluating the effective scales at which a particular method has effectively accounted for atmospheric delay (Hobiger et al. 2010). Structure function can be approximated with a variogram where an experimental variogram is defined as where n is the number of sample points, x i is the location and h is the lag. Comparison with the structure function reduces the relationship to D x (h) =2 (h) . This relationship provides a means of quantifying the variability of the turbulent delay in space as well as providing statistical means of data adjustment together with expected uncertainties in the results.

Spatial Autocorrelation of Turbulent Delay
As it will be discussed in Sect. 3, external data sets utilized for atmospheric calibration are sparse in nature necessitating spatial interpolation. Therefore, understanding the spatial dependence of tropospheric delay is important to designing an interpolation strategy. Turbulent delay manifests at shorter wavelengths and is the component affected most by the assumption of data independence as near points exhibit some form of similarity. Spatial autocorrelation can be modelled by covariance functions estimated from variograms under the assumption of second order stationarity to show spatial relationships between two points in space (Jonsson 2002;Knospe and Jónsson 2007) as follows: This reduces the relationship between the variogram and the covariance (σ 2 ) to Equation 11 only describes the relation in space. Onn and Zebker (2006) extended this relation to include the temporal aspect. They assumed temporal autocorrelation to be guided by the frozen law hypothesis, which can be expressed as is the spatial position, t is the acquisition time for the two images and Δh is the delay as function of height. The importance of including temporal auto-covariance was later demonstrated by Xu et al.(2011) who showed the method by Onn and Zebker (2006) to be superior for the interpolation of water vapour delay.

Isotropy of Atmospheric Signals in InSAR
Spatial structures modelling variation of turbulent delay in space may consider distance only, assuming variation to be independent of direction, or use simplified models. However, anisotropic characteristics in interferograms associated with turbulent delay have previously been reported by Hanssen (1998) and Jonsson (2002). They extended the isotropic variogram to capture the anisotropy as follows: With additional terms representing the change in x and y directions, models that only consider variation as isotropic to some extent have successfully captured variability, for instance, Xu et al. (2011) and Cao et al. (2018) modelled atmospheric delay as isotropic using the spherical theoretical variogram. Their focus was not on characterization of anisotropy and they probably chose a simple model that could still describe spatial variability. However, Knospe and Jonsson (2010) argue that stochastic process characteristics dictate the type of theoretical variograms used to fit the experimental variogram and therefore should capture the statistical and physical properties of the stochastic process. Some variograms that take into account anisotropy have been put forward, e.g., by Wei et al. (2019) who applied exponential theoretical variogram to model anisotropy associated with turbulent delay and reported a mixture of both isotropy and anisotropy at different ranges.
Their choice was motivated by the works of Knospe and Jonsson (2010) who had used the Marten model that represented well the physical properties of turbulent delay and is equivalent to the exponential model when applied with a smoothing factor of 0.5. The exponential model can be fitted with few samples and does not require many terms like the Marten model. We generated variogram maps to highlight isotropic characterstics due to tropospheric delay over a region with minimal change in elevation. As can be seen in Fig. 4 variogram maps exhibit different isotropic characteristics. Anisotropy can also be visualized using the radon transform, which captures intensity variations in a particular direction (Hanssen 1998;Jonsson 2002;Li et al. 2007). These studies also reported a mixture of isotropic and anisotropic characteristics on interferograms suggesting that anisotropy is dependent on temporal atmospheric conditions. Understanding anisotropy is useful for the interpolation of external data used for the correction of the tropospheric delay. Kriging accuracy, for instance, is dependent on proper modelling of variogram. Improved results when anisotropy is considered have been reported (Knospe and Jónsson 2007;Refice et al. 2010;Wei et al. 2019). However, Wei et al. (2019) argue that when density sampling is low like sparse GNSS data sets, anisotropy can be ignored as it has no impact on variogram modelling.

MTI Tropospheric Correction Techniques
Tropospheric delay mitigation within MTI can be classified as empirical or auxiliary methods. Empirical methods calculate the delay directly from the interferogram based on the spatiotemporal properties of tropospheric delay signals. Auxiliary methods estimate the delay using external data. The data products can either be the delay or tropospheric variables that can be used to estimate the delay. Due to limitations associated with each of the techniques, these methods have also been used together. Table 2 summarizes the mitigation methods and the target delay within MTI.

Time Series Filtering
Time series Filtering is the conventional technique in MTI processing. MTI algorithms strive to separate the deformation signal with other artefacts based on a statistical analysis of the spectral behaviour of different contributing phases within the interferogram. APS is considered to be highly correlated in space and temporally uncorrelated, hence can be eliminated through a low-pass filter in space and high-pass filter in time (Berardino et al. 2002;Colesanti et al. 2003;Ferretti et al. 2000;Hooper et al. 2004).
For the unwrapped interferogram i , the differential phase is the sum of multiple contributing signals.
where c elevation error, b i is the baseline, c v velocity of the target, a x i t i tropospheric error and is the processing errors.
Where the first term in Eq. 16 is the contribution from the slave images and the second term represents contribution from the master image. HPT is high pass filter in time and LPS in is low pass filter in space. After accounting for the residual topographic error, the remaining component is majorly composed of the deformation signal, the tropospheric delay and random noise. The method by Berardino et al. (2002) is based on the same principle but filtering is done on the single value decomposition output and relies on estimated overall APS not limited to the scatters. Improper choice of spatial and temporal filters might mask out important deformation signals reducing the reliability of the InSAR-derived measurements (Cao et al. 2018;Hu et al. 2018). However, a priori information on the behaviour of the deformation signal or atmospheric conditions can be of help in choosing the optimal window. Gong et al. (2015) argued that numerical weather predictions models could be used to provide a priori information on atmospheric variability needed for optimal selection of filters. The assumption of temporal randomness negates the seasonality of the stratified delay leading to a bias in the derived deformation estimates (Doin et al. 2009;Jolivet et al. 2011). Hence when applied filtering can mitigate the turbulent delay and might leave seasonal delay which might be interpreted as deformation particularly for areas with high topographic variation. With prior information on spatial variability of tropospheric delay, time series filtering can mitigate the turbulent delay after seasonal delay has been properly accounted for. Similarly, it can be applied in regions of minimal topographical range where stratified delay is negligible. While it has the advantage of being independent of external data its major downside is, it mitigates mostly the turbulent delay in regions with limited topographic variation subject to choosing optimal filtering parameters. There is also the risk of filtering out nonlinear deformation which has a similar variation in time like the turbulent delay.

Phase-based Method
Models that define the relationship between delay and elevation have been put forward including simple global linear relationship (Cavalié et al. 2007), third-degree polynomial (Béjar-Pizarro et al. 2013), piecewise polynomial form of cubic splines (Remy et al. 2003), quadratic (Jung et al. 2014 and power law (Bekaert et al. 2015a). Jolivet et al. (2011) argued that models should be constrained to either linear or quadratic relationships to prevent introduction of high parameters in inversion. Additionally, complex relationships might depict false impression due to overfitting. Linear models assume homogenous stratified atmosphere and thus tropospheric delay is due to elevation change.
where trop is the delay, k is the slope factor and c is the constant shift estimated from non-deforming region (Cavalié et al. 2007) or removed prior estimation of the delay (Elliott et al. 2008). The linear methods assume delay defined by one relation k, which might not be the case for wide scenes. Similarly, non-deforming area or a priori information on deformation is required to estimate c . Lin et al. (2010) proposed a multiscale approach through band pass filtering to overcome global uncertainties in the slope function on assumption of band range insensitive to deformation. Bekaert et al. (2015a) extended the multiscale approach and argued that the delay could be approximated by a power law function and could be estimated on overlapping windows to account for spatial variability.
where k is the window slope factor, h o is the reference height above which change in atmospheric delay is negligible and is the exponential decay factor. Bekaert et al. (2015b) modified the original power law method to allow for disconnected blocks. The success of the power law method is dependent on the delays of master and slave images exhibiting the power law which might not be the case. Shen et al. (2019) demonstrated the limitations of the power law method and instead proposed to estimate the relationship between the delay and elevation from high-resolution weather models. Once a relationship has been defined, it could be used for scaling the delay from the phase information. Similarly, Hu et al. (2018) proposed a linear model applied to phase differences weighted by a covariance matrix which could mitigate stratified delay and partly mitigate turbulent delay. This method, due its empirical nature, is not affected by propagation errors associated with processing external data. Its use is, however, limited to regions with significant topographic variation.

Spectrometers
Precipitable water vapour can be derived at near-infrared or infrared bands by algorithms that rely on the observation of attenuation of water vapour. The Moderate Resolution Imaging Spectroradiometer (MODIS) on board of the National Aeronautics and Space Administration (NASA) Aqua and Terra satellites was launched in 1999 and 2002, respectively. MODIS has a resolution of 1 km with a revisit time of 2-3 days. Envisat carried the Medium Resolution Imaging Spectroradiometer (MERIS) with spatial resolutions of 300 m and 1200 m. The precipitable water vapour (pwv) derived from spectrometer is converted to zenith wet delay as where Π is as dimensionless factor as determined by Bevis et al. (1994) is the density of water, w ratio of molecular masses and T m is the weighted temperature along the integration path. T m can be estimated from sounding or numerical weather models. Li et al. (2003a) compared MODIS-derived zenith wet delay with GNSS and Radiosonde data, though MODIS overestimated the delay, they argued that it could be calibrated with a linear model. Similarly, Li et al. (2005) showed the potential MODIS-derived Pwv to mitigate wet delay in InSAR. Moreover, Li et al. (2011) noted an improved accuracy due from APS mitigation by MERIS due to its ability to correct both stratified and turbulent delay. Puysségur et al. (2007) had argued that the MERIS algorithm could either overestimate or underestimates pwv on cloud shadows and on undetected clouds. A major limitation of the spectrometer near-infrared Pwv is that it can only be applied to daytime cloud-free acquisition. However, MODIS has the potential to estimate infrared Pwv, as it has both day and night capabilities for tropospheric correction (Chang et al. 2014). The temporal difference in acquisition between SAR images and MODIS may also introduce errors given the high temporal variability of water vapour. Taking the Kenyan Rift as an example, Sentinel-1 images from ascending orbit track 130 are acquired at 15:56 and for descending track 152 at 03:19 while MODIS are imaged between 8:00 and 8:30. Hence when used as correction method it is only applicable to the ascending track which is further hindered by the 8 hour difference. Also, spectrometer data may require de-stripping due to the radiometric calibration for MODIS and smile effects for MERIS which may limit its use (Li et al. 2009). Limitations associated with spectrometers may introduce significant errors, for instance, Caijun et al. (2011) reported inconsistent results of up to 4.5 cm between MODIS-corrected InSAR measurements and GNSS. Moreover, it only corrects for wet delay and therefore other methods may need to be incorporated with spectrometer products for an optimal mitigation of tropospheric delay.

Global Navigation Satellite Systems
GNSS and InSAR suffer from similar atmospheric interference though the GNSS's delay differs from delay in InSAR. GNSS's delay represents the average of path delays estimated from more than four satellites in space constrained by the elevation angle used to resolve the GNSS position. Tropospheric delay is regarded as an error term in GNSS processing, but it can be used for InSAR corrections or water vapour mapping with a pointwise accuracy of 1-3 mm compared to a radiosonde (Li et al. 2003b). GNSS can provide point representative zenith wet delay at a highest temporal resolution of up to 15 s. The high temporal resolution was exploited by Onn and Zebker (2006)

Numerical Weather Models
Predictive weather models, reanalysis products and operational forecasting models have been utilized in mitigating tropospheric corrections in MTI (Table 3). The models provide atmospheric variables at a regular grid with fixed time intervals, which are used to estimate the delay based on Eq. (7). The different parametrization of the weather models may result in significant discrepancies among the tropospheric estimates characterized by different spatial patterns (example in Fig. 5). Unless there is a reference technique, it becomes difficult to choose which model to apply as it may degrade the results.
The performances of Numerical Weather Models (NWM) have been characterized by a mixed rate of success. Some studies have reported unsatisfactory results (Fournier et al. 2011;Liu et al. 2009;Stephens et al. 2020). Poor performance has been attributed to coarse spatial resolution of the models, particularly the earlier versions, limited ground data used for constraining the models and mismatch with SAR acquisitions (Gong et al. 2010). Other studies have reported better mitigation of stratified delay (Albino et al. 2020;Kinoshita et al. 2013;Rosell et al. 2019). Improved performance in some cases is due to improved modelling and resolution of the reanalysis products as well as improved parametrization of predictive models. The complexity of the atmospheric conditions limits the NWM to better mitigate  Albino et al. (2020) total delay subject to how strong is the turbulence (Jolivet et al. 2014). Gong et al. (2013) opine that NWM incorrectly estimates water vapour, limiting its use to deformation applications with a magnitude of more than 20 mm. A further modelling example by Yu et al. (2018b) shows promising results to mitigate total delay. The benefit of NWM is their potential to provide estimates which can be used to support or determine a priori estimates of other techniques (Gong et al. 2015;Zhu et al. 2017). Their performance is dependent on the scale, the complexity of the atmosphere and the parametization of the model. However, the performance of high-resolution NWM remains unknown in regions, such as Africa, western China and South America, where limited data are available to constrain the models (Shen et al. 2019).

Combination
Each of the tropospheric correction techniques has strengths and limitations due to the nature of tropospheric delay signals, the complexity of the troposphere, and the external calibration data capability and availability. Furthermore, each of the techniques may result in different estimates with different spatial distribution within interferograms as shown in Fig. 6. A similar observation was also reported by Welch and Schmidt (2017). Applying each of the methods may introduce errors. Dogru (2020) reports on overestimations of delay by NWM and phase-based method delay depending on the magnitude of the signal. However, a combination reduced the amount of overestimated corrections. Tang et al. (2016) report on the inability of filtering to account for tropospheric delay in mountainous region which was solved by incorporating a weather model leading to deformation estimates consistent with GPS measurements.
For optimal correction, studies have proposed to combine methods which act complimentary to each other. For instance, Puysségur et al. (2007)  and Forecasting (WRF) to correct for stratified delay and the time series filtering to account for turbulent delay. Liao et al. (2013) adopted phase-based and filtering to correct for stratified delay and turbulent delay, respectively. Combinations have also been used to fill the data gaps. Li et al. (2009) proposed the combination of MERIS and MODIS to overcome drawbacks, such as the presence of clouds or unavailability of one of the datasets. Similarly, Chang and He (2011) reported an increase in accuracy when sparse GNSS data gaps are filled with numerical weather data instead of GNSS only which increased the error bounds. The outputs of methods, such as weather models, have also been utilized to determine the modelling parameters for phase-based methods (Shen et al. 2019;Zhu et al. 2017). Also due to different resolutions of weather models, Dong et al. (2019) proposed a fusion of multiple NWM models for optimal performance through a weighted approach. Moreover, it provides redundant data for the effective modelling of tropospheric into stratified and turbulent delay (Yu et al. 2018b). In cases where combination has been adopted, improved accuracy in performance has been reported. For instance, Albino et al. (2020) report an increase in accuracy through the integration of a phase-based method and weather model as shown in Fig. 7. Shamshiri et al. (2020) reported improved accuracy of up to 83% though combining small baselines interferograms and GNNS compared to 50% for the ERA-I weather model. Similarly, Yu et al. (2020) reported 61% improvement in accuracy by combining GACOS to remove the long wavelength and filtering to remove the short-wavelength components compared to conventional filtering method.

Tropospheric Corrections within MTI Processing Chain
Conventional MTI processing chains considered APS random in time and correlated in space and therefore it could be mitigated through filtering as described in Sect. 3.1. Though some studies still treat tropospheric delay as random in time,  (Hooper et al. 2012) Filtering, phase-based, spectrometers, weather models Open source SAPROZ (Perissin 2014) Inverted residuals with options of stratified Commercial LiCSBAS (Morishita et al. 2020) Weather model (GACOS) Open source GAMMA (Werner et al. 2000) Spatial-temporal filtering Commercial GIAnT (Agram et al. 2016) Phase-based and weather models from pyaps Open source OSARIS (Loibl et al. 2019) Weather model (GACOS) Open source MintPy (Yunjun et al. 2019) Phase-based and numerical weather models Open source Pyrate (Garthwaite et al. 2017) Spatio-temporal filtering Open source others have incorporated other techniques in their processing chain. It is common to apply corrections to interferograms in the MTI processing chain.
Corrections are applied to the unwrapped interferograms though applying correction to the wrapped interferograms may improve the unwrapping accuracy (Li et al. 2006a;Pinel et al. 2008;Jolivet et al. 2011;Jung et al. 2014). Except for filtering that considers the spatio-temporal properties for the entire stack in estimation of tropospheric delay, other methods are based on determining the tropospheric delay differences between the two dates of an interferogram.
where Φ t 1 t 2 is the unwrapped interferogram that has been corrected for residual topographic error, Φ � t 1 t 2 is the corrected interferogram while Φ t2 trop and Φ t1 trop are the tropospheric estimates from the external data sets. For the phasebased method, the tropospheric component is estimated as described in Sect. 3.2 for each interferometric pair and Φ � t 1 t 2 in obtained in a similar way. To estimate tropospheric delay for the entire stack, Eq. (21) is applied to a parallel network similar to that of the MTI technique. Support for various correction techniques is dependent on the MTI software and open source might support other techniques. The conventional procedures are listed in Table 4 that do not require further modification by the user. There are also stand-alone software products that provide estimates that can be integrated in MTI processing like pyaps (Jolivet et al. 2011) and TRAIN (Bekaert et al. 2015b). The emergence of online service, such as GACOS, described by Yu et al. (2018a) provides a faster means of evaluating and correcting tropospheric artefacts in interferograms. Current MTI platforms that are concerned with large-scale processing of InSAR, such as Loibl et al. (2019) and Morishita et al. (2020), have adopted GACOS as their main correction tool, probably due to its coverage and ready-to-use delay maps that can mitigate total delay.

Tropospheric Corrections in InSAR Applications
We also reviewed studies that utilised MTI in their applications related to anthropogenic activities, landslides, earthquakes and volcanic activity, whose focus was not on modelling atmospheric delay. This is meant to give the reader an overview of how tropospheric delay is treated in InSAR application and the reliability of such measurements. Rémy et al. (2015) and Yip et al. (2019) have demonstrated that improper modelling of atmospheric delay can lead to wrong results. We discussed the limitations of each method in Sect. 3. The search criteria were that the study should not have tropospheric delay mitigation or quantification as part of its objective and should at least describe the tropospheric correction technique(s) applied in their methodology. A total of 120 peer reviewed articles that met the criteria were reviewed, 30 for each of the four applications. The reference list of the articles is in the appendix. Table 5 presents a summary of the applied correction methods. Filtering is the most frequently utilised technique as it is the conventional method and most MTI software have not integrated other techniques in their processing chain. The use of other correction techniques is either not supported or the users have to perform additional post-processing steps themselves, which is subject to their own capabilities. A lack of popularity of the NWM methods for applications related to fluid extraction and landslides is probably related to the spatial resolution of the NWM. Most studies in these categories focus on small regions at scales lower than the resolution of most numerical weather products. It is, however, important to point out that the lack of studies that utilize a particular technique does not mean that those studies do not exist. They could exist but did not meet the search criteria or were not accessible. The use of spectrometers, as discussed in Sect. 3.3, is well documented but most studies either focus on the potential of spectrometers for water vapour mapping or use it as a reference to determine the performance of other techniques. These studies therefore do not meet the search criteria. The sparse GNSS network seems to limit its use as a correction technique, hence its poor representation compared to other techniques. Studies that employ phase-based methods are often located in regions with significant topography where the elevation dependence can be easily determined and applied by fitting polynomials. The results highlight that though multiple techniques have been proposed, their implementation remains limited due to the lack of data, the lack of integration of the methods within the MTI processing chain and probably the user capabilities.

Accuracy Assessment and Uncertainty in Measurements due to Tropospheric Delay
Tropospheric noise in interferograms is usually identified through visual examination before and after correction. Though subjective, this can provide qualitative assessment of the correction method. Another approach is to introduce noise with the same properties as tropospheric delay using simulated data. Performance gauged on the ability to retrieve a priori error Tymofyeyeva and Fialko 2015). Similarly, statistical diagnostic outputs can be used to assess accuracy, which can be determined by comparing the results with the values from other measurements (Li et al. 2003a;Gong et al. 2013) or estimated from mean calculated from a window or the whole of the interferogram (Fu et al. 2018). Reduced standard deviation or increased correlation after the correction indicates an improved performance. However, these are global estimates and may not provide more information on spatial variability. Moreover, diagnostic outputs, such as root-mean-square error, are not suitable for spatially correlated data, such as atmospheric delay, as they may not be normally distributed (Zandbergen 2008). Variograms discussed in Sect. 2.3.1 provide robust accuracy metrics with useful information on the variability at different wavelength scales (Foster et al. 2013;Kinoshita et al. 2013). Similarly, variance spectra comparisons of the corrected interferograms against original interferograms not only show the variability but also the wavelength scales at which the applied technique mitigated atmospheric artefacts (Foster et al. 2006;Hobiger et al. 2010). While no standard metrics exist, models that capture spatial variability provide accuracy indicators that respect the underlying tropospheric characteristics. Performance indicators can be extended to estimate the uncertainty in measurements, which are directly correlated to the quality of the mitigation of tropospheric delay. The potential of MTI to achieve mm accuracy is documented but most of the MTI algorithms or studies do not provide uncertainties (Cao et al. 2018). The uncertainties in MTI measurements, due to other artefacts superimposed over the deformation signal, can largely be attributed to tropospheric delay (Fattahi and Amelung 2015). Uncertainty inclusion in MTI modelling is important to determine error bounds (Barnhart and Lohman 2013), filtering parameters (Gong et al. 2015) or the necessary duration to achieve the required deformation accuracy (Emardson et al. 2003). The phase in the interferogram can be considered to be composed of functional and stochastic parts (González & Fernández, 2011).
where j i,def is the phase delay due to deformation, j i,topo is due to topographic error, j i,orb is due to the orbital error, j i,atm is the atmospheric-phase delay while j i, represents other noise due to decorrelation, processing and the satellite's own systems. Orbital errors can be mitigated by applying polynomial function, though for modern satellites, it is minimal (Fattahi and Amelung 2014). The use of coherent pixels leads to decorrelation error being treated as minimal (Cao et al. 2018). The deterministic part also includes the stratified delay component unless it has been accounted for, which introduces bias in measurements. The bias can be estimated by comparing elevation change with delay (Delacourt et al. 1998;Taylor and Peltzer 2006) or inferred from modelling of geophysical phenomena parameters. The error bounds before and after tropospheric correction indicating the bias (Scott and Lohman 2016). Alternatively, it can estimate the bias due to data source used in the mitigation of stratified delay (Mateus et al. 2015). The stochastic part can thus be attributed to turbulent delay with the following variance covariance matrix (Hanssen 2001): σ 2 is the weighted a priori variance usually taken as 1 and Q e i e i is the cofactor matrix. Variance covariance can be estimated from interferograms, like in the case of González and Fernández (2011) though they ignored the interchanging role of SAR as master or slave. Cao et al. (2018) advanced the proposed network-based model to estimate both variance covariance of InSAR measurements and atmospheric delay. Alternatively, Parker et al. (2015) suggested the use of variance covariance estimated from NWM time series to provide a priori uncertainties and to determine the amount of images required to detect low-magnitude deformation in volcanoes. Similarly, Fattahi and Amelung (2015) estimated the uncertainity from differences between MODIS and ERA-I data.
The magnitude of uncertainty is dependent on the scale, period and availability of data (Table 6). 1. The nature of terrain also determines the type of the correction. Low topographic variation imply a weak stratified delay, while both the hydrostatic delay and the wet delay must be considered in areas with high topographic variation (Elliott et al. 2008). Auxiliary data availability may restrict the choice a correction technique. GNSS for instance is limited to areas with good network while data like spectrometers though globally available are only applicable to images acquired during daytime. On the other hand, redundant data provide the opportunity to apply different methods and, based on accuracy assessment, the method that performs better is adopted. Moreover, data like numerical weather models can be used to provide information on the nature of spatial temporal variation like the cases of Gong et al. (2015) and Shen et al. (2019). In instances where the total delay is present, methods can be integrated to act complimentary to each other. For example, the phase-based method can be used to mitigate the stratified delay while filtering is used to mitigate the turbulent delay. Not all methods are supported by the current MTI processing software packages. In Table 4, we have provided a summary of the methods that are supported by the different algorithms. This can limit the the available options. The course resolutions of the numerical weather models prevent its applicability in areas with limited acreage. Similarly, vast areas require methods that can handle both short-wavelength and long-wavelength tropospheric  Figure 8 provides a schematic workflow that could guide in choosing the appropriate technique. It begins with the identification of the tropospheric signals present in the interferograms. Based on the type of signal, topographic variation, image acquisition time, data availability and the supported methods within the processing chain appropriate technique can be adopted. MTI applications though not directly related to a mitigation method determine the expected tropospheric signal superimposed with the deformation signal. For instance, volcanoes are characterized by high topographic variation leading to a strongly stratified delay (Pinel et al. 2011), landslides too are influenced by a change in slope and will contain a component of the stratified delay. Methods that target stratified delay have to be incorporated for landslides and volcano monitoring. Moreover, methods that capture seasonal variation should be applied for landslides studies to prevent misidentification of seasonal variation as landslides (Dong et al. 2019). For applications, such as earthquake and crustal deformation monitoring covering large areas, the tropospheric delay will be characterized by short-and long-wavelength components. In such instances, methods that target both components should be applied . Additionally, if it stretches to regions of different climatic conditions, the estimation of the hydrostatic delay must be included due to the expected change in pressure and temperature (Doin et al. 2009). For applications such as fluid extraction or urban monitoring that cover a small extend characterized with limited topographic variation and the expected deformation signal is linear, methods that target turbulent delay are appropriate. However, for geophysical applications that exhibit nonlinear deformation, if subjected to filtering some of the deformation could be masked out. The location of the study also has an influence. For instance, tropical regions are characterized by a complex troposphere cannot be modelled easily (Ebmeier et al. 2013;Seimon and Phillipps 2011), and cannot be captured well by the NWM. This limits the use of NWM as mitigation technique in such regions (Stephens et al. 2020;Webb et al. 2020). For applications that are concerned with handling huge amount of data, such as global volcano monitoring, use of already determined delay maps, such as GACOS, can be appropriate. Similarly operational applications can benefit from GACOS due to its near real-time availability (Albino et al. 2020). The expected deformation magnitude is also a factor, for applications that are in the scale of mm the mitigation technique should have lower uncertainty than the deformation signal which can be estimated as described in Sect. 6.

Conclusion and Outlook
While success has been achieved in modelling tropospheric delay characteristics in InSAR, their mitigation remains a major challenge due to the complexity of the atmosphere, the lack of external data or limited availability of data needed for the calibration of tropospheric delay. Though conventional MTI algorithms treat tropospheric delay as random in time, some component of the delay is seasonal. All components of the tropospheric delay need to be accounted for, as ignoring the hydrostatic delay may lead to unreliable results particularly for large-scale processing that stretches through different climatic zones. The specific correction methods have reported varying rates of success. It is difficult to say which method is the best, as each method performs differently depending on the region, topographic variation, spatial coverage and period. However, the combination of multiple methods takes advantage of the strengths of each individual method and seems to provide the best overall results. This is, however, subject to proper understanding of the components of tropospheric artefacts present in the interferograms. Moreover, improvements in the accuracy of numerical weather models provide exciting opportunities for the correction of tropospheric effects in InSAR. This is supported by optimized techniques, such as GACOS, which have posted success and are now being integrated into large-scale multitemporal processing. Most studies use the conventional time series filtering as a mitigation method due to lack of integration of other tropospheric correction methods to the existing multi-temporal algorithms leading to limited uptake by InSAR applications studies. Simple accuracy performance metrics, such as the root-mean-square error, adopted by major studies, are not satisfactory, as this does not account for spatial variability. Spatial covariance models, some accounting for anisotropy, not only show the performance metrics but also provide an avenue for statistical adjustments to evaluate uncertainty and to determine the error bounds. The choice of the appropriate technique is a multi-faced procedure that is influenced by the nature of the tropospheric signal present in the interferograms, by the magnitude and temporal behaviour of the deformation signal, topographic variation, by the SAR acquisition time, by the climatic conditions and by the external data availability.
The future of InSAR is promising with the expected launch of NASA-ISRO SAR Mission in 2022 adding to the already huge historical data sets. This may enable new applications, such as global deformation monitoring, but it also comes with new challenges to processing and the handling of errors, such as tropospheric delay. Machine learning has been proposed but only so far limited to classifying observation as deformation or tropospheric delay (Anantrasirichai et al. 2018). Hence, future mitigation methods will focus on dealing with challenges associated with automatic processing of such data. Already the processed weather models, such as GACOS, allow near realtime tropospheric estimates (Albino et al. 2020) with a lag of two days. The large data sets also provide an opportunity to estimate tropospheric delay empirically or through the intergration with other methods. Given the challenges associated with tropospheric delay, the estimation and the discrepancy in estimates betwen the different correction methods stochastic modelling is inevitable. This highlights the need for the inclusion of uncertainties and error bounds in InSAR-derived measurements.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.