5 The tropospheric effect is a significant error source in radar observations. To accurately 6 calculate displacement fields, it is essential to apply advanced techniques to mitigate the 7 tropospheric effect in the interferograms. This study examines land subsidence in Tehran 8 province, Iran, using Sentinel -1A data from 2014 to 2021. Methods such as t he Generic 9 Atmospheric Correction Online Service (GACOS), ray tracing, and a new integration method 10 were utilized to estimate the tropospheric effect. The integration method and ray tracing 11 technique employed ERA5 reanalysis data from the European Centre f or Medium -Range 12 Weather Forecasts (ECMWF). In the final stage, the InSAR interferograms were refined using 13 tropospheric corrections derived from these methods. Comparison of the subsidence data 14 obtained from the three tropospheric correction techniques wit h Global Positioning System 15 (GPS) observations showed that the results were closely aligned. However, the subsidence 16 velocity derived from the new integration approach proved to be the most accurate. 17 Keywords: InSAR , Tropospheric Correction , Ray Tracing , Subsidence 18 19 20 21 2 1. Introduction 22 Natural and human activities such as tectonic plate movements, volcanic activity, earthquakes, and changes in groundwater levels can lead to vertical deformation of the Earth's surface, known as subsidence. Efforts have been made to forecast and monitor groundwater level depletion using machine learning (Haji -Aghajany et al., 2023) and to predict subsidence using similar methods (Tasan et al., 2023). Regular and optimal monitoring of subsidence can help prevent hazards in both residential and natural areas. In recent years, factors such as industrial activities and population growth have led to excessive groundwate r extraction, resulting in aquifer depletion and subsequent land subsidence in various regions. These depletions are significant contributors to displacement and can cause subsidence of up to 17 cm in some areas (Gumilar et al., 2015). The primary cause of global subsidence is the change in groundwater levels, influenced by natural parameters and human activities, including precipitation, groundwater extraction, and soil absorption. Subsidence, the gradual or sudden downward movement of the Earth's surface, can pose serious risks to the stability of aquifers and human - made structures, making the study and monitoring of subsidence crucial worldwide. Iran has faced significant spatiotemporal limitations in groundwater resources in recent years, primarily due t o extensive extraction for industrial activities and urban development (Motagh et al., 2008). This paper focuses on an area south of Tehran in Tehran province, Iran. Land subsidence in this region can be studied using hydrogeological methods and geodetic t ools such as GPS and InSAR (Abidin, 2008). InSAR is a precise and powerful technique with high spatial resolution, capable of measuring Earth's surface displacement with millimeter accuracy (Ferretti et al., 2001). However, the 3 tropospheric effect is a ma jor limiting factor in InSAR outputs (Zebker et al., 1997; Jolivet et al., 2014), necessitating its correction or reduction. Various tools are available for calculating the tropospheric effect, including meteorological data and radiosondes. Numerous previo us studies have focused on mitigating tropospheric effects using various approaches. One common method involves the interferometric combination of a series of radar images, assuming that the atmospheric effect is random over time while the displacement signal correlates with time (Ferretti et al., 2001; Tang et al., 2016). This approach employs temporal filtering of large radar image time series to eliminate atmospheric influences. Another approach relies on connected stacks of interferograms to characteriz e phase delay patterns (Foster et al., 2006; Lauknes, 2011). Additional methods analyze the relationship between elevation and tropospheric delay (Remy et al., 2003; Elliott et al., 2008) or utilize external meteorological data (Fournier et al., 2011; Dee et al., 2011; Jolivet et al., 2014). Recently, spectral diversity has been applied to reduce tropospheric effects (Jung and Lee, 2015; Fattahi et al., 2017; Mastro et al., 2020). Moreover, some researchers have employed three - dimensional (3D) and various t wo-dimensional (2D) ray tracing approaches (Haji -Aghajany et al., 2018; Haji -Aghajany et al., 2019; Haji -Aghajany and Amerian, 2020a), as well as outputs from numerical weather models (Haji -Aghajany and Amerian, 2020b). Troposphere tomography, a powerful t echnique for reconstructing water vapor, has also shown promise and has been applied in fields such as precipitation downscaling (Izanlou et al., 2024; Haji -Aghajany, 2021; Haji -Aghajany et al., 2022). In a recent study aimed at enhancing the accuracy of I nSAR tropospheric corrections, Maddahi et al. (2024) utilized tomography to reduce tropospheric effects in InSAR measurements (Maddahi et al., 2024). This study reviews the impact of the tropospheric layer on InSAR measurements, using GACOS products, a ray tracing algorithm, and 4 a new integration method to calculate tropospheric delay. These computed delays are then applied to the radar interferograms. Finally, land subsidence estimates derived from the three tropospheric correction methods are compared wit h GPS measurements to assess their accuracy. 2. InSAR process ing The InSAR technique relies on phase extraction from radar images. By calculating the phase difference recorded in corresponding pixels of two radar images taken from the same area at different times, an interferogram is generated, which serves as the primary measurement in InSAR. The resulting interferogram contains a mixture of signals: deformation, topographical effects, orbital errors, tropospheric effects, and noise due to diffusio n variations. The phase differences observed in the interferogram represent both the deformation signal and contributions from topographical changes, orbital errors, atmospheric effects, and diffusion changes (Hanssen, 2001) : def atm topo orbit noise           (1) Where def , atm , topo , orbit and noise are phase differences caused by deformation, tropospheric delay , topography, orbit al error , and thermal noise, respectively . Since the goal of InSAR analysis is to compute deformation, other contributions to the InSAR phase must be removed. In this paper, the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) with a 10 -meter resolution is used to mitigate the topographical effect, and orbital files are used to correct orbital errors. Intrusive parameters in InSAR measurements can significantly affect the calculated displacement fields caused by phenomena such as subsidence (Khalili et 5 al., 2023) or earthquakes (Haji -Aghajany et al., 2020). This study aims to explore methods for reducing these effects. In the traditional InSAR method, only two radar images are used to calculate the deformation that occurred between the acquisition times of these images. To addre ss the limitations of this approach, Persistent Scatterer Interferometry (PSI) and Small Baseline Subset (SBAS) algorithms have been developed to compute deformation time series using multi -temporal radar images (Berardino et al., 2002). These algorithms r ely on a large number of images and persistent scatterers (PS) in the study area. The PSI method allows for the calculation of linear displacement velocity and the reconstruction of nonlinear time series for each PS. The SBAS algorithm, on the other hand, uses PSs that maintain coherence above a selected threshold across all interferograms (Samsonov, 2010). Both algorithms offer high accuracy in calculating deformation over time, represented as time series. Recently, InSAR processing methods have evolved to combine the strengths of the PSI and SBAS algorithms (Hooper et al., 2008). 3. Tropospheric effect mitigation Tropospheric delay in radar signals arises from variations in tropospheric parameters such as pressure, temperature, and water vapor. As a radar s ignal traverses the troposphere, its velocity decreases, causing variable delays that affect observations (Jolivet et al., 2014). These variations can induce localized phase gradients in InSAR measurements, primarily due to changes in tropospheric paramete rs. Generally, temporal variations in pressure and temperature are not significant enough to create localized phase gradients in InSAR data (Hanssen, 1998). The primary factor influencing InSAR data is temporal changes in water vapor, which is typically pr esent in the lower troposphere (Hanssen, 1998; Hanssen, 2001). Therefore, it is essential to compute these phase gradients caused by tropospheric variations. 6 3.1. Tropospheric correction using new integration technique The tropospheric path delay in the zenith direction, known as the zenith tropospheric delay, is computed using numerical integration, which can be expressed as follows: 1 1 3 2 1 1 2 30.776 0.71 , 6 3.7 , 5 100 0z k R R ee -6 1 d dZTD = 10 P(z )+ ((k - k ) +k ) kdz2 1 3 2 g R TT z m k KPa KPa Pv k K a       (2) where z0 is the height, z is the highest height of the troposphere, Rv (461.495 Jkg-1K-1) and Rd (287.05 Jkg-1K-1) are the specific gas constant s for water v apor and dry air, respectively , P(z 0) is the surface pressure , e is the water vapor pressure, T is the temperature in K , and gm is the averaged gravitational acceleration. To mitigate the tropospheric effect from the radar signal, the zenith tropospheric delay must be converted to the line -of-sight direction (Cao et al., 2021). However, previous res earch has shown that this conversion, typically achieved using mapping functions, can introduce errors into the computation and final outputs. These errors may even result in misinterpreting the deformation signal in InSAR analysis. Therefore, this paper introduces a new, simplified integration method along the line -of-sight direction to circumvent errors associated with mapping functions in the results. Similar to the classic integration method along the zenith direction, this approach requires meteorologi cal data in 3D space. The meteorological data along the line of sight are estimated using interpolation techniques. Figure 1 illustrates the comparison between the classic and new integration methods. In this study, spline interpolation is applied in the v ertical direction, while kriging interpolation is used in the horizontal direction (Dubrule, 1984). For more detailed information about this method, please refer to Cao et al. (2021). 7 Figure 1. Comparison between classic integration and new mode (Khalili et al., 2023) 3.3. Tropospheric correction using GACOS products The GACOS system, developed at Newcastle University in the United Kingdom, relies on meteorological data to address tropospheric delay in radar signals. It utilizes the 90 m resolut ion SRTM DEM (60° S - 60° N) and the Terra ASTER Global Digital Elevation Model (GDEM) in its processing. Additionally, GACOS integrates High -Resolution (HRES) datasets with spatial resolutions of up to 0.1° x 0.1° lat/long grid and 137 vertical pressure l evels every 6 hours, derived from the European Centre for Medium -Range Weather Forecasts (ECMWF) (Xiao et al., 2021). Through the implementation of an iterative tropospheric decomposition model developed by Yu et al. (2018), GACOS separates stratified and turbulent components from tropospheric delays, producing high spatial resolution Zenith Total Delay (ZTD) maps. Yu et al. (2018) evaluated GACOS performance using various sites worldwide. Beyond theoretical and technological validations, GACOS offers a sig nificant innovation in its operational mode: it provides a free online service, reducing barriers to InSAR tropospheric correction. Users only 8 need to submit a request with the location of their study region and radar data, and GACOS tropospheric products will be generated for them. Unlike GPS techniques that rely on ground measurements or meteorological data with potential delays of weeks, GACOS offers global coverage and near real -time availability. 3.2. Tropospheric correction using ray tracing method The ray tracing method stands as a potent technique for accurately reconstructing the path of wave propagation using tropospheric indices, commonly employed in estimating tropospheric effects. These approaches are typically categorized into two classes: 2D and 3D algorithms (Hofmeister, 2016). In the 2D approach, ray tracing occurs within a plane with a fixed azimuth, while in the 3D method, the signal's path can be reconstructed in 3D space. One notable 2D ray tracing method is the piecewise -linear (PWL) a lgorithm, known for its adequate accuracy and high processing speed. This algorithm employs the refractive index at various tropospheric levels to apply Snell’s law of refraction and compute the tropospheric effect, utilizing recursive relations for tropos pheric delay computation (Hofmeister, 2016). Another refinement of the 2D ray tracing technique is the refined piecewise -linear (RPWL) algorithm, which enhances the accuracy of the PWL method by incorporating intermediate compression levels. This approach aims to improve the precision of ray tracing processes. Further insights into these methodologies and their application can be found in previous studies (Hofmeister, 2016). For this study, the RPWL algorithm is employed to compute the tropospheric delay. 4. Study area and data set Since the 1960s, Tehran has experienced land subsidence, a phenomenon exacerbated by the city's burgeoning population, which now exceeds 15 million inhabitants. Situated in northern Iran, the Tehran plain spans 2,250 km² and is bounded by the Alborz Mounta ins to the north 9 and the Arad and Fashapouye mountains to the south. Land subsidence, primarily attributed to excessive groundwater extraction, predominantly affects the southwestern portion of the basin. Precise leveling operations initially revealed subs idence in this region. Tehran, Iran's capital, lies within a semiarid zone particularly susceptible to rapid land subsidence. Over recent years, subsidence in the Tehran Plain has encroached upon densely populated urban areas, particularly impacting the so uthern and southwestern sectors of the city (Esmaeili and Motagh, 2016; Azadnejad et al., 2019). This area has endured significant adverse effects from land subsidence, primarily driven by excessive groundwater withdrawal. Numerous studies have investigate d the extent of subsidence in Tehran (Ajourlou et al., 2019; Haghighi and Motagh, 2019; Azadnejad et al., 2019, 2020). The first attempt to estimate deformation in the Tehran Plain using the InSAR method was undertaken by Shemshaki et al. (2005). Figure 2 provides the geographical context of the study area. Figure 2. Location of the study area To conduct the InSAR technique, a series of radar images captured by the Sentinel -1A satellite over two areas have been utilized. The specifications of the employe d images are detailed in Table 1, while the spatiotemporal distribution is depicted in Figure 3. In this study, ERA5 reanalysis data from ECMWF has been employed for both ray tracing and new integration techniques. ERA5 provides reanalysis data 10 covering th e period from 1950 to the present, furnishing meteorological data on 137 pressure levels from the surface up to 1 Pa (approximately 80 km). The spatial resolution of ERA5 data is approximately 31 km (Hersbach and Dee, 2016). For each radar image, the close st ERA5 data in time is selected. Temperature, pressure, and water vapor values at each pressure level are then interpolated. Kriging and spline interpolation methods are applied along the horizontal and vertical directions, respectively. Additionally, GPS observations are utilized to assess the derived displacement rates. Table.1. Information of the used radar images IW Beam Mode 28 Path 112 Frame Ascending Flight Direction VV+VH Polarization 44900 Absolute Orbit 16 Jan 2016 – 1 Jan 2021 Date Figure 3. Spatiotemporal distribution of images 5. Processes A series of Sentinel -1A radar images have been employed in the InSAR process utilizing the SBAS method. To estimate the displacement velocity field, intrusive effects must be mitigated or minimized from the calculated interferograms. Following the removal of topographical and orbital effects using DEM and orbital data, the remaining interferograms c ontain only tropospheric effects and deformation phase. 11 In this study, the tropospheric effect has been computed using the aforementioned methods, including the new integration algorithm, GACOS products, and ray tracing. Below, samples of computed troposp heric delay from these approaches can be observed in Figure 4. Figure 4. Sample s of obtained tropospheric delay from new integration algorithm, GACOS products and ray tracing technique (left to right) . After applying the computed tropospheric delays to the obtained interferograms, the displacement rate of the study area has been determined (Figure 5). To facilitate comparison of the results, two separate profiles have been selected from the derived displacement velocity fields. Figure 6 illustrates that when employing the ray tracing method, numerous uplift and subsidence signals are observed across the study area. However, the utilization of GACOS products and the new integration technique has notably reduced the number of artificial signals. The presenc e of high displacement rates and a considerable number of signals in the displacement rate map obtained via the spatiotemporal filter method may be attributed to the lingering tropospheric effects in the obtained interferograms. Conversely, the use of the new integration technique and GACOS products could potentially mitigate both subsidence and uplift signals. 12 Figure 5. Obtained displacement velocity fields 5. Figure 6. Considered profiles on the obtained displacement velocity maps 6. Validation with GPS observation To assess the efficiency of the results, it is essential to compare the computed displacement velocities with the displacement rates derived from accurate measurements, such as GPS data. The validation using GPS data provides more acc urate and reliable results due to the time interval overlap between the GPS data and radar images. To better evaluate the obtained results, the InSAR time series have been compared with the GPS time series. This comparison, conducted at the position of GPS stations, is depicted in Figure 7. 13 Figure 7. Samples of InSAR time -series and water level decline of four piezometric wells from up to down, respectively . In order to comprehensively evaluate the mentioned methods, various statistical parameters have been computed. These statistical parameters, calculated to compare the obtained results with GPS observations, are presented in Table 2. Table.2. Statistical parameters of different methods Root mean square error (mm) Mean absolute difference (mm) Method 12.71 10.60 New integration 16.89 21.22 GACOS 20.14 18.96 Ray tracing Based on the obtained results, it can be concluded that the new integration method demonstrates superior effectiveness in tropospheric calculations, leading to a more accurate displacement field compared to the other two approaches. Considering the expecte d accuracy of InSAR, which typically falls within millimeters, the differences observed in the results are statistically significant. However, it's important to note that implementing the new integration method and the ray tracing technique necessitates pr ecise calculations and programming 14 efforts. In contrast, the use of GACOS products is notably simpler and faster, requiring minimal user calculations. 7. Conclusions Tropospheric delay represents a significant constraint in radar measurements, capable of introducing errors of over one hundred millimeters to InSAR outputs. In this study, three methods —GACOS products, ray tracing technique, and the new integration approach —were compared to estimate tropospheric delay and enhance InSAR displacement accuracy. To investigate the impact of these approaches on InSAR displacement, a region was selected in the central part of Iran, within Tehran province. Following the computation of tropospheric correction using the aforementioned methods and subsequent correction of obtained interferograms, time series and displacement velocities were calculated. Ultimately, comparison of the outputs with GPS data revealed that the new integration technique yielded a more accurate displacement field compared to the other two approa ches.