Tropospheric Correction of Sentinel-1 Synthetic Aperture Radar Interferograms Using a High-Resolution Weather Model Validated by GNSS Measurements

Synthetic Aperture Radar Interferometry (InSAR) is a space geodetic technique used for mapping deformations of the Earth’s surface. It has been developed and used increasingly during the last thirty years to measure displacements produced by earthquakes, volcanic activity and other crustal deformations. A limiting factor to this technique is the effect of the troposphere, as spatial and temporal variations in temperature, pressure, and relative humidity introduce significant phase delays in the microwave imagery, thus “masking” surface displacements due to tectonic or other geophysical processes. The use of Numerical Weather Prediction (NWP) models as a tropospheric correction method in InSAR can tackle several of the problems faced with other correction techniques (such as timing, spatial coverage and data availability issues). High-resolution tropospheric modelling is particularly useful in the case of single interferograms, where the removal of the atmospheric phase screen (and especially the highly variable turbulent component) can reveal large-amplitude deformation signals (as in the case of an earthquake). In the western Gulf of Corinth, prominent topography makes the removal of both the stratified and turbulent atmospheric phase screens a challenging task. Here, we investigate the extent to which a high-resolution WRF 1-km re-analysis can produce detailed tropospheric delay maps of the required accuracy by coupling its output (in terms of Zenith Total Delay or ZTD) with the vertical delay component in GNSS measurements. The model is operated with varying physical parameterization in order to identify the best configuration, and validated with GNSS zenithal tropospheric delays, providing a benchmark of real atmospheric conditions. We correct sixteen Sentinel-1A interferograms with differential delay maps at the line-of-sight (LOS) produced by WRF re-analysis. In most cases, corrections lead to a decrease in the phase gradient, with average root-mean-square (RMS) and standard deviation (SD) reductions in the wrapped phase of 6.0% and 19.3%, respectively. Results suggest a high potential of the model to reproduce both the long-wavelength stratified atmospheric signal and the short-wave turbulent atmospheric component which are evident in the interferograms.


Introduction
Synthetic Aperture Radar Interferometry (InSAR) is a technique successfully deployed in recent years for mapping accurate fields (at the millimetre level) of ground deformations itable water vapour (PWV) accurately at high spatial resolution (250-1000 m) makes them a highly efficient technique under suitable conditions. The Medium Resolution Imaging Spectrometer (MERIS) PWV accuracy has been estimated close to 1 mm, equivalent to 6 mm of Zenith Wet Delay (ZWD) for each epoch, or 9 mm between two epochs [35]. This is equivalent to approximately 1 cm in radar line-of-sight for ENVISAT data with an incidence angle of 23 • . With respect to the Moderate Resolution Imaging Spectroradiometer (MODIS), PWV accuracy has been estimated at best equal to that of MERIS, and at worst twice that of MERIS [22].
With regards to NWP models, recent studies have focused both on the use of Global Atmospheric Models (GAMs) and Limited-Area Models (LAMs) to predict delays at the time of SAR acquisitions and correct for the stratified tropospheric delays. In a study where data from three GAMs (ERA40, OPERA and NARR) were used [29], it was shown that tropospheric artefacts were better removed compared with InSAR derived delay/elevation ratios in cases where the correlation between elevation and displacement is large. Further attempts to exploit output from GAM products, such as ERA-Interim, NARR or MERRA [19], for corrections in different geographical and tectonic environments, demonstrated a better estimation of stratified tropospheric delay and rather poor results for estimating turbulent patterns on single interferograms. Global Atmospheric Models suffer from low temporal and spatial resolution, and output data need to be interpolated in space and time in order to match the resolution of an interferogram and the exact acquisition times. Therefore, the technique of using high-resolution regional weather models (LAMs) nested within coarser, global weather models to estimate the atmospheric delay is gaining ground [22,[30][31][32][33][36][37][38]. Foster et al. [31] employed the MM5 regional model at high horizontal resolution (3 km) to obtain tropospheric delay fields over the Island of Hawaii and Mount St Helens in the United States (US) with mixed success, as the model configuration failed to accurately predict tropospheric delays at shorter wavelengths (under 8 km) in most cases. Wadge et al. [30] performed tropospheric delay corrections, using the Unified Regional Mesoscale Model at high-resolution over Mount Etna in Sicily. The study concluded that LAMs show promising potential in calculating the path delay due to tropospheric water vapour in regions of high relief and high water vapour variance, such as Mount Etna, and that the model was able, under certain conditions, to simulate local modifications of water vapour content by mechanisms such as land-sea breezes. Nico et al. [38] produced synthetic phase delay maps by deploying WRF at high-resolution (1 km) over Lisbon and the Azores in Portugal, with results confirming a good agreement between the interferograms and WRF-derived phase delays and minor statistically relevant differences between them. More recently, Bekaert et al. [22] produced tropospheric delay fields with the WRF regional model over three test regions with complex topography. The model is nested at 7 km horizontal resolution and is initiated with data from the Global Forecast System (GFS). Results are compared with tropospheric delays obtained with other state-of-the-art methods, including MERIS and MODIS spectroscopy, ERA-Interim, and both the conventional linear and power-law empirical methods, demonstrating that weather models offered better performance when atmospheric turbulence and dynamic local weather is present. On the whole, recent studies which employ LAMs for the removal of tropospheric artifacts from InSAR data demonstrate the ability of weather models to calculate detailed tropospheric delay fields under any atmospheric conditions and at the exact times of SAR acquisitions, thus offering a reliable tool for tropospheric corrections. Indeed, it is evident that in cases of atmospheric turbulence and dynamic local weather, weather models can offer better performance [22,37]. However, the generic configuration and parameterization of the LAMs used in these studies have prohibited, so far, the full exploitation of the method.
Here, our work aims to maximise the potential of using LAMs for the calculation and removal of the tropospheric component from SAR interferograms by evaluating a novel methodology that integrates recent advances in the fields of remote sensing meteorology (GNSS and InSAR) and high-resolution numerical weather forecasting (WRF). We inves-tigate the extent to which a high-resolution WRF 1-km re-analysis can produce detailed tropospheric delay maps of the required accuracy by coupling its output, in terms of Zenith Total Delay (ZTD), with the vertical delay component in GNSS measurements. One of the main limitations in the estimation of precise ZTD fields is the presence of highly variable water vapour signals, both in space and in time, which are exhibited in the differential atmosphere as densely distributed short-wavelength phase gradients [19]. It is expected that the high horizontal resolution and dense vertical layering of the model will be capable of capturing near-surface atmospheric processes, such as sea breezes, orographic flows, turbulent boundary layer interactions, etc., which influence the distribution of water vapour in settings of complex topography, such as the western Gulf of Corinth in Greece, thus overcoming this important limitation. The model is operated with varying physical parameterization in order to identify the best possible configuration, with GNSS measurements used as reference data for fine-tuning and validating the model. We then compare sixteen Sentinel-1 interferograms with differential delay maps at the line-of-sight (LOS) produced by WRF re-analysis and perform tropospheric corrections of the atmospheric phase screen (APS) in the wrapped interferograms.

Description of the Study Area and Experimental Setup
The experiment, with the code name PaTrop (Patras-Troposphere), was implemented to provide the data needed for this study. The PaTrop test site covers an area of approximately 130 × 90 km in the region of the Western Gulf of Corinth (GoC). It is one of the most active intra-continental rifts in the world, and therefore a region with intense seismic activity. Geodetic studies conducted over the past 20 years with continuous GNSS and InSAR observations have revealed north-south extension rates up to 1.5 cm yr −1 [39,40], one of the highest worldwide. A network of GNSS stations used to monitor surface displacements in the area continuously provided the in situ zenithal tropospheric delay measurements for our study. These include nineteen permanent Topcon GB1000 and Topcon Net G3A receivers fitted with Topcon PG-A1 antennas, the locations of which are shown in Figure 1b. The stations cover a wide geographical extent while capturing a variety of different topographical and meteorological conditions (i.e., coastal, inland, or mountainous terrain) at elevations between 0 and 1020 m above sea level (ASL). This was necessary in order to account for water vapour variations resulting from orographic, coastal, and frontal gradients that could be present. The first objective of the PaTrop experiment was to couple the zenithal tropospheric delay (ZTD) derived from GNSS data with the ZTDs derived from the output of a highresolution meteorological model (WRF) in order to investigate the model's capability to reproduce the tropospheric conditions that contribute to the InSAR noise signal (in par-  The first objective of the PaTrop experiment was to couple the zenithal tropospheric delay (ZTD) derived from GNSS data with the ZTDs derived from the output of a highresolution meteorological model (WRF) in order to investigate the model's capability to reproduce the tropospheric conditions that contribute to the InSAR noise signal (in particular the highly variable water vapour distribution) and provide a benchmark of real observational data for validating the model output. A series of re-analysis runs were produced to determine the best possible configuration for our study, keeping the spatial resolution of the inner domain at 1 km, which is an approach supported theoretically and practically to tackle uncertainties in high-resolution modelling [41,42]. A parametric analysis was performed for a two-week period (17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29) June 2016), during which the output of five different model configurations was tested against GNSS tropospheric measurements from the network of permanent stations in the study area [43]. The optimum model configuration which resulted from the analysis was subsequently employed to provide the data for the main part of the PaTrop experiment.

WRF Configuration and Parameterization of Physical Components
For the high-resolution dynamical downscaling simulation performed with WRF v 3.7.1 [44] over the PaTrop area of the western Gulf of Corinth, four nested domains were used (d01-d04), with a horizontal resolution of 27, 9, 3 and 1 km, respectively, as shown in Figure 1a where P s is the total pressure (mbar) at the Earth's surface, and: accounts for the variation in gravitational acceleration with latitude λ and the height H of the surface above the ellipsoid (in km). And: where p v is the water vapour pressure (mbar), and T the air temperature (K), integrated along the zenith path z. In practice, we calculated one ZWD per vertical layer, and we added the 45 values to obtain the total ZWD. Land-use categories were taken from USGS 24-category data, which are available for different horizontal resolutions (10 , 5 , 2 , 30"). The initial land topography dataset used was the Global 30 Arc-Second Elevation Model (GTOPO30) provided by United States Geological Survey (USGS), with a 30" resolution for the smaller domain (d04), and coarser resolutions (10 , 5 , 2 ) for domains d01-d03, respectively. However, in order to test the impact of a more detailed topography on the re-analysis output, a high-resolution terrestrial dataset of d04 was later introduced (ASTER 1" global GDEM v2), with a horizontal grid of 30 m.
Five different physical parameterization schemes were tested, as listed in Table 1, in order to evaluate each scheme for its forecasting skill. The schemes were selected based Remote Sens. 2021, 13, 2258 6 of 17 on existing studies where similar high-resolution WRF simulations were used [41,[45][46][47], ranging from relatively simple to more complex and computationally demanding physical parameterizations suited for high-resolution runs. There have been numerous studies validating the output of different model configurations with observations under specific conditions [41,[45][46][47], showing that globally there is no optimal scheme, but rather, different schemes produce better results with respect to application, domain, season, variable, etc. Similarly, we performed a parametric analysis for a two-week period (17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29) June 2016), during which the output of the five different model configurations was validated against GNSS tropospheric measurements from 16 permanent stations in the study area (d04).

Tropospheric Correction of SAR Interferograms
Based on the results of the parametric analysis, output data from the optimal WRF high-resolution 1-km re-analysis configuration were used to calculate precise ZTD fields over the PaTrop study area, at the exact times of Sentinel-1 SAR acquisitions for the ascending track 175 and descending track 80 (i.e., 1630 and 0430 UTM, respectively). Resulting Zenith Total Delay (ZTD) values calculated from specific model output parameters (surface pressure, air temperature and water vapour profiles) were validated against a dataset of GNSS derived ZTD values, providing point measurements at the 19 points where the stations were located.
Sentinel-1 SAR data were used in this study for the generation of InSAR interferograms for 2016. The two ESA Sentinel-1 satellites (S1A and S1B) have a 6-day repeat time and carry a C-band synthetic aperture radar with a 56 mm wavelength and four operating modes. Acquisitions with a 5 × 20 m resolution in the Interferometric Wide (IW) mode were used. In total, 17 acquisitions were combined to produce 8 interferograms for the ascending (S-N) track 175 and 8 interferograms for the descending (N-S) track 80. Temporal baselines ranged from 6 days to 42 days, and perpendicular baselines were in the order of 150 m. Multilook was 6 in range and 2 in azimuth. Processing of InSAR data was done with the European Space Agency's Sentinel Application Platform (SNAP) version 5.0 software [48], following several steps: SAR image formation, co-registration, interferogram formation, flattening (using precise orbits from ESA), and topography removing using a three arc-second (about 90 m) Shuttle Radar Topography Mission (SRTM) DEM. The final georeferenced product was resampled at 25 m using bilinear interpolation. The tropospheric correction process took place before unwrapping the interferogram. This was done in order to minimise phase ambiguities and improve the reliability of interferogram unwrapping in a region such as the western GoC, where coherence is low due to the vegetation and rough topography, which results in geometric decorrelation [49].
The InSAR observations examined are listed in Table 3. Average values of WRF vs. GNSS ZTD bias are also listed in the Table. ∆bias was calculated by averaging the absolute bias (ZTD WRF -ZTD GNSS ) differences between epoch 1 and epoch 2, at the 19 PaTrop stations, and is an indication of the model's performance with respect to the observational data.
where f i1 and f i2 denote the model value at epochs 1 and 2, respectively, o i1 and o i2 the observational value at epochs 1 and 2, respectively, and N is the number of observations. For the tropospheric correction, the first step was to generate differential delay maps of the total single-path tropospheric delay at the line-of-sight (LOS) at the times of the two SAR acquisitions. This was done by subtracting the 1 × 1 km single ZTD map produced by the d04 WRF output of epoch 2 from the corresponding ZTD map of epoch 1 (Figure 2). In the resulting differential delay map, LOS total delay values were calculated at each 1-km grid cell by multiplying the corresponding ZTD value with cosθ, where θ is the average incidence angle of the S1 swath in IW mode (35 • ). These values were then horizontally and vertically interpolated, using the weighted average inverse distance to a power gridding method, to a new 25 × 25 m grid corresponding to the pixel resolution of the interferogram. The resulting differential delay map was then wrapped (LOS total differences are converted into 2π interferometric phase fringes) and subtracted from the wrapped interferogram to produce a phase map of residuals, as illustrated in Figures 7 and 8 in the Results Section. Before the phase subtraction, the map of differential delay was "shifted" by minimizing the RMS between the two GeoTIFF maps so that their average zero phases were aligned.

Parametric Analysis and Validation of WRF Schemes with GNSS Data
The parametric analysis was performed for a two-week period (17-29 June 2016), with the aim to test the effect of different WRF parameterization schemes on the model Figure 2. Example of producing a differential ZTD map from the subtraction of two single ZTD maps (epoch 1-epoch 2), produced with WRF at the times of InSAR acquisitions. Pair corresponds to dates 18-30 September 2016. The size of the area is 130 × 90 km, corresponding to the inner WRF domain d04. The resulting differential ZTD map is then converted to LOS total delay map and wrapped (arithmetic values were transformed into 2π phase values) to correct the corresponding interferogram.
All Sentinel-1 interferograms examined have temporal baselines between 6 and 42 days, ensuring the highest possible coherence and negligible amounts of crustal deformation between acquisitions. The coherence maps of both the ascending and descending tracks are shown in Figure 3. There were no major tectonic events or anthropogenic processes within the selected time series, and therefore the phase signal should mostly be attributed to tropospheric effects. Unlike historic SAR sensors, orbital errors for Sentinel-1A/B are expected to be small due to highly accurate orbit information [15].
high-resolution WRF-derived delay fields led to a decrease in the phase gradient, as demonstrated by the corresponding RMS and SD reductions ( Table 3). The RMS of the corrected interferogram was improved in 15 out of 16 cases, with reductions ranging from 2.3% to 14.1% (6.0% on average), while SD was improved in all 16 cases, with reductions ranging from 13.4% to 32.7% (19.3% on average). Furthermore, the degree of tropospheric delay correction was correlated with WRF-GNSS average bias differences (Δbias) at the times of acquisitions. This is demonstrated in Figure 4, where the slope of the reductions in both indicators with respect to Δbias was inversely proportional. Cases with medium coherence (in green) had a slightly wider error margin than cases with high coherence (in red). When the correlation between the two indexes was plotted (Figure 6), it was confirmed that their variability was not random and that coherence was a determining factor (interferograms with high coherence are more likely to produce better tropospheric corrections). Apart from the quantitative assessment, corrections were also assessed qualitatively, case by case, by identifying visible improvements in fringe continuity. In wrapped interferograms, the use of numerical indicators was not always adequate for assessing the degree of tropospheric phase gradient improvement due to potential problems with pixel

Parametric Analysis and Validation of WRF Schemes with GNSS Data
The parametric analysis was performed for a two-week period (17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29) June 2016), with the aim to test the effect of different WRF parameterization schemes on the model output. Tropospheric data (i.e., ZTD values) from the CRL GNSS network provided the benchmark for the analysis for comparison with ZTDs derived from WRF at the nearest grid point. Five different physical parameterization schemes were evaluated, as outlined in Table 2. Results were tested for their statistical significance in terms of correlation and bias with respect to the observational dataset and indicated that WRF schemes MOD4 and MOD5 exhibited a better prediction skill during the test period, with small differences between them. A summary of all validation metrics (Pearson co-efficient, mean bias, mean absolute bias and root-mean-square error) for the five parameterization schemes is presented in Table 2. MOD5 overall exhibited the strongest correlation (Pearson correlation co-efficient R) with observations, together with a lower mean bias (MB), lower root-mean-square error (RMSE) than other parameterizations, and smaller spread of RMSE between stations than MOD4 ( Figure 4). It was therefore selected as the optimum configuration for producing tropospheric delay maps for the entire PaTrop period (January-December 2016) [43].
The parametric analysis demonstrated that WRF schemes that included more complex physical parameterization (MOD4 and MOD5), although more computationally demanding, are better suited for high-resolution re-analysis simulations. The physical parameterization of MOD5 scheme included the SBU-YLin microphysics model, with a more sophisticated scheme for ice, snow, rain and graupel processes in the lower troposphere, the NOAH land surface model, the MM5 similarity scheme for surface layer physics and the YSU planetary boundary layer model. In addition, MOD4 and MOD5 used a cumulus convection scheme in both the 27-km and 9-km domains, thus simulating processes such as convective fluxes and the associated evaporation or condensation of water more coherently over a complex terrain with land-sea contrasts. MOD5 overall exhibited the strongest correlation (Pearson correlation co-efficient R with observations, together with a lower mean bias (MB), lower root-mean-square erro (RMSE) than other parameterizations, and smaller spread of RMSE between stations tha MOD4 ( Figure 3). It was therefore selected as the optimum configuration for producin tropospheric delay maps for the entire PaTrop period (January-December 2016) [43].
The parametric analysis demonstrated that WRF schemes that included more com plex physical parameterization (MOD4 and MOD5), although more computationally de manding, are better suited for high-resolution re-analysis simulations. The physical pa rameterization of MOD5 scheme included the SBU-YLin microphysics model, with a mor sophisticated scheme for ice, snow, rain and graupel processes in the lower tropospher the NOAH land surface model, the MM5 similarity scheme for surface layer physics an the YSU planetary boundary layer model. In addition, MOD4 and MOD5 used a cumulu convection scheme in both the 27-km and 9-km domains, thus simulating processes suc as convective fluxes and the associated evaporation or condensation of water more cohe ently over a complex terrain with land-sea contrasts.

InSAR Tropospheric Correction with the Use of WRF Derived Delay Maps
In order to provide a quantitative assessment of the corrections applied in every cas the root-mean-square (RMS) and standard deviation (SD) of the phase distribution of bot the original and the corrected interferograms were calculated and their differences rec orded. A reduction in the RMS or SD of the wrapped interferogram after the correctio wasapplied is a clear indication that there was a decrease in the phase gradient and fring continuity was smoother. Table 3 lists the RMS and SD results together with the corre sponding Δbias value for all 16 cases examined, while Figure 4 graphically illustrates th correlation of RMS and SD with Δbias.

InSAR Tropospheric Correction with the Use of WRF Derived Delay Maps
In order to provide a quantitative assessment of the corrections applied in every case, the root-mean-square (RMS) and standard deviation (SD) of the phase distribution of both the original and the corrected interferograms were calculated and their differences recorded. A reduction in the RMS or SD of the wrapped interferogram after the correction was applied is a clear indication that there was a decrease in the phase gradient and fringe continuity was smoother. Table 3 lists the RMS and SD results together with the corresponding ∆bias value for all 16 cases examined, while Figure 5 graphically illustrates the correlation of RMS and SD with ∆ bias .    In most cases, corrections applied to the wrapped interferograms with the use of high-resolution WRF-derived delay fields led to a decrease in the phase gradient, as demonstrated by the corresponding RMS and SD reductions ( Table 3). The RMS of the corrected interferogram was improved in 15 out of 16 cases, with reductions ranging from 2.3% to 14.1% (6.0% on average), while SD was improved in all 16 cases, with reductions ranging from 13.4% to 32.7% (19.3% on average). Furthermore, the degree of tropospheric delay correction was correlated with WRF-GNSS average bias differences (∆ bias ) at the times of acquisitions. This is demonstrated in Figure 5, where the slope of the reductions in both indicators with respect to ∆ bias was inversely proportional. Cases with medium coherence (in green) had a slightly wider error margin than cases with high coherence (in red). When the correlation between the two indexes was plotted (Figure 6), it was confirmed that their variability was not random and that coherence was a determining factor (interferograms with high coherence are more likely to produce better tropospheric corrections). In most cases, corrections applied to the wrapped interferograms with the use o high-resolution WRF-derived delay fields led to a decrease in the phase gradient, a demonstrated by the corresponding RMS and SD reductions ( Table 3). The RMS of th corrected interferogram was improved in 15 out of 16 cases, with reductions ranging from 2.3% to 14.1% (6.0% on average), while SD was improved in all 16 cases, with reduction ranging from 13.4% to 32.7% (19.3% on average). Furthermore, the degree of tropospheri delay correction was correlated with WRF-GNSS average bias differences (Δbias) at th times of acquisitions. This is demonstrated in Figure 4, where the slope of the reduction in both indicators with respect to Δbias was inversely proportional. Cases with medium coherence (in green) had a slightly wider error margin than cases with high coherence (i red). When the correlation between the two indexes was plotted (Figure 6), it was con firmed that their variability was not random and that coherence was a determining facto (interferograms with high coherence are more likely to produce better tropospheric cor rections). Apart from the quantitative assessment, corrections were also assessed qualitatively case by case, by identifying visible improvements in fringe continuity. In wrapped inter ferograms, the use of numerical indicators was not always adequate for assessing the de Apart from the quantitative assessment, corrections were also assessed qualitatively, case by case, by identifying visible improvements in fringe continuity. In wrapped interferograms, the use of numerical indicators was not always adequate for assessing the degree of tropospheric phase gradient improvement due to potential problems with pixel decorrelation (low coherence) in parts of the interferogram or the existence of other components which contribute to the phase gradient. Here, two examples (out of the 16 examined) are shown, one for ascending track 175 ( Figure 7) and one for descending track 80 (Figure 8).
Four additional examples (Figures S1-S4) can be found in the Supplementary Materials (addendum). Case-by-case comments can be found in the legend of each figure. Summarising, we see that as a general rule, in examples where interferogram coherence is high and the forecasting skill of the WRF simulation is good, as predicted by GNSS measurements, the differential troposphere is significantly removed, and the residual phase map exhibits smoother fringe continuity. However, corrections are not always visible across the whole interferogram, first of all, because of low coherence in the western part of the image (and other parts as well), second, because of other errors (geometrical, etc.), third, because the model does not always recreate the differential atmosphere properly. More specifically, out of the sixteen interferograms that we examined: (a) in cases where coherence was good (temporal baselines usually <30 days) and ∆ bias was small (between 0 and 20 mm), the degree of tropospheric correction was high, resulting in a significant reduction in the density of tropospheric fringes in large sections of the interferogram, as illustrated in example S10 ( Figure 8); (b) in cases where coherence was low (temporal baselines usually >30 days) and ∆ bias was small (between 0 and 20 mm), the degree of tropospheric correction was high only in areas with high coherence, as illustrated in example S2 ( Figure 7); (c) in cases where WRF-GNSS average bias differences were high (∆ bias > 20 mm), the density of tropospheric fringes was reduced at a lesser degree, and the correction was localised in smaller areas of the interferogram. ferogram, as illustrated in example S10 ( Figure 8); (b) in cases where coherence was lo (temporal baselines usually > 30 days) and Δbias was small (between 0 and 20 mm), t degree of tropospheric correction was high only in areas with high coherence, as illu trated in example S2 (Figure 7); (c) in cases where WRF-GNSS average bias differenc were high (Δbias > 20 mm), the density of tropospheric fringes was reduced at a lesser d gree, and the correction was localised in smaller areas of the interferogram. Corresponding WRF-derived wrapped differential LOS delay map (top right). The correlation between interferogram and meteogram is directly visible, for example, in the Mornos valley (red box) or around the Panachaikon mountain (black box) as the model's forecasting skill was high in both acquisition epochs (∆bias low). The wrapped differential LOS delay map generated from WRF output data was subtracted from the corresponding wrapped interferogram to produce a map of residual 2π phase cycles (bottom). The resulting residual map, in this case, demonstrates tropospheric corrections in several areas of the interferogram where coherence was high, leading to a decrease in the phase gradient.
geometry (top left). Corresponding WRF-derived wrapped differential LOS delay map (top right). The correlation between interferogram and meteogram is directly visible, for example, in the Mornos valley (red box) or around the Panachaikon mountain (black box) as the model's forecasting skill was high in both acquisition epochs (Δbias low). The wrapped differential LOS delay map generated from WRF output data was subtracted from the corresponding wrapped interferogram to produce a map of residual 2π phase cycles (bottom). The resulting residual map, in this case, demonstrates tropospheric corrections in several areas of the interferogram where coherence was high, leading to a decrease in the phase gradient. Corresponding WRF-derived wrapped differential LOS delay map (top right). An example of good consistency between WRF and GNSS is also reflected in the good correlation between the interferogram and the delay map, with short and long wavelengths being observed. In the corresponding residual map after subtraction of WRF-derived wrapped differential LOS delay map from the interferogram (bottom), good tropospheric corrections and an overall reduction in the phase gradient were observed across the whole extent of the interferogram.

Discussion
High-resolution tropospheric modeling is particularly useful in the case of single i terferograms, as the precise knowledge of meteorological parameters can remove t main source of error from the InSAR data, which is the delay due to the atmospheric r fraction of the signal, revealing large-amplitude deformation signals (as in the case of earthquake). In fact, it is the wet component of the delay present in the interferogram which can be used to generate accurate maps of the atmosphere's precipitable water vap Corresponding WRF-derived wrapped differential LOS delay map (top right). An example of good consistency between WRF and GNSS is also reflected in the good correlation between the interferogram and the delay map, with short and long wavelengths being observed. In the corresponding residual map after subtraction of WRF-derived wrapped differential LOS delay map from the interferogram (bottom), good tropospheric corrections and an overall reduction in the phase gradient were observed across the whole extent of the interferogram.

Discussion
High-resolution tropospheric modeling is particularly useful in the case of single interferograms, as the precise knowledge of meteorological parameters can remove the main source of error from the InSAR data, which is the delay due to the atmospheric refraction of the signal, revealing large-amplitude deformation signals (as in the case of an earthquake). In fact, it is the wet component of the delay present in the interferograms, which can be used to generate accurate maps of the atmosphere's precipitable water vapor (PWV) over large areas [50,51], thus providing a useful tool in the domain of InSAR meteorology and vice-versa for the assimilation of interferometric SAR data into numerical weather models.
Although computationally demanding, downscaling WRF at 1 km horizontal resolution offers a detailed reconstruction of the 3-D tropospheric field over the study area, necessary for capturing small-scale atmospheric processes (such as sea breezes, orographic flows, turbulent boundary layer interactions, etc.) which cannot be captured by coarser LAMS or GAMs. Furthermore, model output was validated with the use of GNSS tropospheric data retrieved from a dense array of stations (nineteen receivers) covering the study area. Model validation with vertical column data (GNSS zenithal delays) instead of ground measurements offers the capability of evaluating the model's forecasting skill over the entire 3-D field, thus enabling fine-tuning of its physical parameterization. It was demonstrated by the parametric analysis that schemes with a more complex physical parameterization (MOD4 and MOD5) exhibited lower bias and better correlation with respect to the GNSS observations. Additionally, turning the cumulus convection scheme on for both the 27-km and 9-km domains simulated processes, such as convective fluxes and the associated evaporation or condensation of water, more coherently over complex terrain with land-sea contrasts.
The proposed methodology was tested in a location of complex topography (western Gulf of Corinth), where the presence of highly variable water vapour signals made the removal of short-wavelength phase gradients even more challenging. Our results suggested a reduction in both long-wavelength (5-50 km) and short-wavelength (<5 km) phase delays of the wrapped phase. Residual maps exhibited a reduction in the stratified topographycorrelated atmospheric signal, but most importantly, a reduction in the difficult to detect turbulent atmospheric signal (i.e., "wet" delay) in complex topographical structures of the scale of a few km, such as river incisions, valleys, etc.
Interferograms with a small temporal baseline (less than 12 days) had a higher coherence (i.e., pixel correlation between the two acquisitions) and the reduction in tropospheric "fringes" was more prominent, both visually and quantitatively. In interferograms with a larger temporal baseline, low coherence (i.e., pixel decorrelation due to system noise or terrain changes) masked many of the tropospheric artefacts. Therefore, a combination of high interferometric coherence and low ∆bias (average model bias difference) between the two epochs is a good indication that the tropospheric correction will lead to a reduction in the phase gradient.

Conclusions
Overall, our study demonstrated the high potential and effectiveness of using highresolution atmospheric modelling (WRF in this instance), for correcting the effects of tropospheric delay on InSAR observations and correcting atmospheric phase gradients in interferograms. The proposed methodology augments the model's ability to predict zenithal delays in the western Gulf of Corinth, by fine-tuning its physical parameterization with the use of in situ ZTD measurements from a network of permanent GNSS stations. Furthermore, by introducing a high-resolution topography (ASTER 1s DEM), the calculation of ZTD delay fields became more accurate, and bias was minimised. The use of high-resolution Limited-Area Models (LAMs), validated by GNSS measurements, has a number of advantages over other methods which are currently used for removing the atmospheric phase screen in InSAR observations. First of all, it can be used day and night and under any weather conditions. The method can be applied in any geographical location, as long as the LAM is locally configured and parameterized. Model output data can be retrieved at the exact times of InSAR acquisitions, and the high spatial resolution (1-km) and dense vertical layering are capable of capturing near-surface atmospheric processes where complex topography is present. Finally, model validation with vertical column data (GNSS zenithal delays) instead of ground measurements offered the capability of evaluating the model's forecasting skill over the entire 3-D field, thus enabling fine-tuning of its physical parameterization with the use of a high-accuracy representative observational dataset. This is particularly useful when it comes to estimating the highly variable water vapour signals which were exhibited in the differential atmosphere of the interferogram as densely distributed short-wavelength phase gradients.
Tropospheric corrections performed over a set of 16 wrapped Sentinel-1 interferograms with the use of high-resolution WRF-derived delay fields led to significant reductions in atmospherically-related phase gradients, with average root-mean-square and standard deviation reductions in the wrapped phase of 6.0% and 19.3%, respectively. The actual degree of correction was related to the WRF-GNSS ZTD average bias difference (∆ bias ) between the two acquisition epochs, and this can be a useful indicator for determining the effectiveness of the approach based on the model's forecasting skill. Results suggest that both the stratified and the turbulent atmospheric signal can be reduced from wrapped interferograms. This is a fair improvement compared with predictive methods based on coarser GAMs, which are effective in reducing only lateral variations in stratification.
The removal of the differential tropospheric signal before the unwrapping process is beneficial for the correct estimation of the remaining noise in the interferogram. The phase ambiguity due to the aliasing of atmospheric gradients in regions of rough topography is reduced, and this minimises unwrapping errors. This will eventually lead to more reliable final products, thus enabling the detection of ground deformation signals in single interferograms (i.e., in the case of an earthquake) and improving velocity field estimates by resolving lateral variations in stratification in InSAR time series analysis.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to extensive size.