# Use of GNSS Tropospheric Delay Measurements for the Parameterization and Validation of WRF High-Resolution Re-Analysis over the Western Gulf of Corinth, Greece: The PaTrop Experiment

Nikolaos Roukounakis, Dimitris Katsanos, Pierre Briole, Panagiotis Elias, Ioannis Kioutsioukis, Athanassios Argiriou, Adrianos Retalis

► To cite this version:

Nikolaos Roukounakis, Dimitris Katsanos, Pierre Briole, Panagiotis Elias, Ioannis Kioutsioukis, et al.. Use of GNSS Tropospheric Delay Measurements for the Parameterization and Validation of WRF High-Resolution Re-Analysis over the Western Gulf of Corinth, Greece: The PaTrop Experiment. Remote Sensing, 2021, 13 (10), pp.1898. 10.3390/rs13101898. hal-03358951

HAL Id: hal-03358951

https://hal.science/hal-03358951

Submitted on 20 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Use of GNSS Tropospheric Delay Measurements for the Parameterization and Validation of WRF High-Resolution Re-Analysis over the Western Gulf of Corinth, Greece: The PaTrop Experiment

Nikolaos Roukounakis<sup>1,2,3,*</sup>,<sup>1 ID</sup> Dimitris Katsanos<sup>1 ID</sup>, Pierre Briole<sup>2 ID</sup>, Panagiotis Elias<sup>4 ID</sup>, Ioannis Kioutsioukis<sup>3 ID</sup>, Athanassios A. Argiriou<sup>3 ID</sup> and Adrianos Retalis<sup>1 ID</sup>

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Roukounakis, N.; Katsanos, D.; Briole, P.; Elias, P.; Kioutsioukis, I.; Argiriou, A.A.; Retalis, A. Use of GNSS Tropospheric Delay Measurements for the Parameterization and Validation of WRF High-Resolution Re-Analysis over the Western Gulf of Corinth, Greece: The PaTrop Experiment.
*Remote Sens.* **2021**, *13*, 1898. https://doi.org/10.3390/rs13101898

Academic Editor: Stefania Bonafoni

Received: 17 March 2021

Accepted: 7 May 2021

Published: 13 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In the last thirty years, Synthetic Aperture Radar interferometry (InSAR) and the Global Navigation Satellite System (GNSS) have become fundamental space geodetic techniques for mapping surface deformations due to tectonic movements. One major limiting factor to those techniques is the effect of the troposphere, as surface velocities are of the order of a few mm yr⁻<fnref n="1" />, and high accuracy (to mm level) is required. The troposphere introduces a path delay in the microwave signal, which, in the case of GNSS Precise Point Positioning (PPP), can nowadays be partly removed with the use of specialized mapping functions. Moreover, tropospheric stratification and short wavelength spatial turbulences produce an additive noise to the low amplitude ground deformations calculated by the (multitemporal) InSAR methodology. InSAR atmospheric phase delay corrections are much more challenging, as opposed to GNSS PPP, due to the single pass geometry and the gridded nature of the acquired data. Thus, the precise knowledge of the tropospheric parameters along the propagation medium is extremely useful for the estimation and correction of the atmospheric phase delay. In this context, the PaTrop experiment aims to maximize the potential of using a high-resolution Limited-Area Model for the calculation and removal of the tropospheric noise from InSAR data, by following a synergistic approach and integrating all the latest advances in the fields of remote sensing meteorology (GNSS and InSAR) and Numerical Weather Forecasting (WRF). In the first phase of the experiment, presented in the current paper, we investigate the extent to which a high-resolution 1 km WRF weather 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 initially operated with varying parameterization, with GNSS measurements providing a benchmark of real atmospheric conditions. Subsequently, the final WRF daily re-analysis run covers an extended period of one year, based on the optimum model parameterization scheme demonstrated by the parametric analysis. The two datasets (predicted and observed) are compared and statistically evaluated, in order to investigate the extent to which meteorological parameters that affect ZTD can be simulated accurately by the model under different weather conditions. Results demonstrate a strong correlation between predicted and observed ZTDs at the 19 GNSS stations throughout the year (R ranges from 0.91 to 0.93), with an average mean bias (MB) of -19.2 mm, indicating that the model tends to slightly underestimate the tropospheric ZTD as compared to the GNSS derived values. With respect to the seasonal component, model performance is better during the autumn period (October–December), followed by spring (April–June). Setting the acceptable bias range at ±23 mm (equal to the amplitude of one Sentinel-1 C-band phase cycle when projected to the zenith distance), it is demonstrated that the model

<footnote n="1">¹ Institute for Environmental Research and Sustainable Development, National Observatory of Athens, 15236 Athens, Greece; katsanos@noa.gr (D.K.); adrianr@noa.gr (A.R.)</footnote>

<footnote n="2">² Laboratoire de Geologie, UMR CNRS ENS PSL 8538, 24 rue Lhomond, CEDEX 05, 75231 Paris, France; briole@ens.fr</footnote>

<footnote n="3">³ Laboratory of Atmospheric Physics, Department of Physics, University of Patras, 26500 Patras, Greece; kioutio@upatras.gr (I.K.); athanarg@upatras.gr (A.A.A.)</footnote>

<footnote n="4">⁴ Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, 15236 Athens, Greece; pelias@noa.gr</footnote>

<footnote>* Correspondence: nrouk@noa.gr</footnote>
produces satisfactory results, with a percentage of ZTD values within the bias margin ranging from 57% in summer to 63% in autumn.

**Keywords:** tropospheric delay; WRF; physics parameterization; high-resolution modeling; GNSS meteorology; InSAR tropospheric corrections; model validation

# 1. Introduction

In recent years, and with the advent of high-performance computer systems (HPCs), high-resolution Limited-Area Model (LAM) applications are being increasingly utilized, both in operational weather forecasting and in atmospheric research as weather re-analysis tools. Atmospheric studies involving dynamical downscaling with high-resolution LAMs focus either on the estimation of local meteorological parameters and climatic patterns (detailed wind fields, precipitation characteristics, local humidity etc.) or on the actual validation of the model used, depending on its parameterization and comparison with observational data. Common LAMs which have been employed for research and climate re-analysis include MM5, RAMS and WRF. In recent years, WRF [1] has become the most widely used model, due to its versatility as an open-source package supported by a large user community and the fact that it can be operated efficiently on parallel computing platforms, making it suitable for use in a wide range of applications at high horizontal resolutions (down to 1 km or less).

The use of high-resolution regional weather models (LAMs) nested within coarser, global, weather models is currently gaining ground in the domain of satellite imagery, as the precise knowledge of meteorological parameters along the lower atmosphere is extremely useful for the correction of remote sensing data, and in particular Synthetic Aperture Radar (SAR) interferograms [2–8]. Tropospheric stratification and short wavelength spatial turbulences produce a phase delay, and consequently an additive noise to the low amplitude ground deformations calculated by the multitemporal InSAR methodology. The removal of the atmospheric phase screen (APS) is a challenging task (as opposed to GNSS, for example, which uses signals in the same microwave bandwidth), due to several reasons, including the small number of satellites, the single pass geometry and the gridded nature of the acquired data. Several methods have so far been pursued to mitigate the tropospheric delay in InSAR data, including local atmospheric data collection [9], GNSS zenith delay estimations [10–12], satellite multispectral imagery analysis [13–15], empirical phase/topography estimations [8,16–18], and zenith delay estimations from Global Atmospheric Models (GAMs) [19–21]. These methods have their limitations, as they rely either on local data assimilation, which is rarely available at the desired spatial resolution, or on empirical estimations which are difficult in situations where deformation and topography are correlated.

The advantage of using LAMs over other techniques is the ability to produce detailed delay maps at high spatial resolution (down to 1 km), by estimating atmospheric parameters, such as water vapour and temperature, in the entire vertical column which cannot be easily estimated otherwise. Foster et al. [2,3], 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 U.S. 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. [5], performed tropospheric delay corrections, using the Unified Regional Mesoscale Model at high-resolution over Mount Etna in Sicily. The study concluded that LAMs show a promising potential in calculating the path delay due to tropospheric water vapor 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. More recently, Bekaert et al. [8] 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. The statistical analysis demonstrated that spectrometers provided the largest RMSE reduction, but only under daylight and cloud-free conditions. Phase-based methods (linear and power-law) outperformed the weather models in regions where tropospheric delays were correlated with topography, but in regions where this was less apparent, due to atmospheric turbulence and dynamic local weather, weather models offered better performance. 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 [8,22]. However, the generic configuration and parameterization of the LAMs used in these studies have prohibited, so far, the full exploitation of the method.

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 which integrates the latest advances in the fields of remote sensing meteorology (GNSS and InSAR) and high-resolution Numerical Weather Forecasting (WRF). Our research focuses on the interaction of the three aforementioned techniques in the area of the western Gulf of Corinth in Greece and aims to optimise their synergic output by identifying strengths, weaknesses and uncertainties with respect to the measurement of zenith tropospheric delays. 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 (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. 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, thus overcoming this important limitation.

The model is initially operated with varying parameterization in order to demonstrate the best possible configuration, with GNSS measurements providing a benchmark of real atmospheric conditions. The two datasets (predicted and observed) are compared and statistically evaluated for a period of one year, in order to investigate the extent to which meteorological parameters that affect ZTD, can be simulated accurately by the model under different weather conditions. In the next phase, we compare twenty Sentinel-1A interferograms with differential delay maps at the LOS produced by WRF re-analysis. Results of the second phase of the experiment are presented in a separate paper.

## 2. Materials and Methods

### 2.1. *Description of the Study Area and Experimental Setup*

The experiment, with the code name PaTrop (Patras–Troposphere), was implemented for providing the data needed for this study. The PaTrop test site covers an area of approximately 150 × 90 km in the region of the Western Gulf of Corinth (GoC). It is identified as a site of intense geophysical activity and forms one of the most active intra-continental rifts in the world. Geodetic studies conducted over the past 20 years based on GNSS and InSAR observations have revealed North-South extension rates up to 1.5 cm yr⁻¹ [23,24], one of the highest worldwide. In this context, the existing network of permanent GNSS stations used by the Corinth Rift Laboratory (CRL) to continuously
monitor surface displacements in the area could be exploited for the needs of our study. A network of nineteen permanent Topcon GB1000 and Topcon Net G3A GNSS receivers fitted with Topcon PG-A1 antennas provided the GNSS data and the subsequent in-situ zenithal tropospheric delay measurements, as shown in Figure 1b. Ten of those receivers were installed during the PaTrop campaign in order to expand the existing CRL network. The stations' locations were selected to cover the whole geographic extent of the study area, 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 deemed necessary, in order to account mainly for water vapour variations resulting from orographic, coastal, and frontal gradients that could be present.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) Left, map showing the four nested domains (d01-d04) used for WRF weather re-analysis over the Western GoC. (b) Right, map of the PaTrop study area, Western GoC. Previously installed permanent GNSS stations are marked in blue colour, newly installed GNSS stations marked in red colour. The black box indicates the location of the 1 × 1 km inner domain of the WRF re-analysis run (Domain d04).</figcaption></figure>

Concurrently with the field campaign, the modelling setup and configuration was also implemented. The primary objective of PaTrop is to couple the zenithal tropospheric delay (ZTD) derived from GNSS data with the ZTDs derived from the output of a high-resolution meteorological model (WRF), in order to investigate the model's capability to reproduce the tropospheric conditions that contribute to the noise signal (in particular the highly variable water vapour distribution), and provide a benchmark of real observational data for validating the model output. The Weather Research and Forecasting (WRF) model [1] is a widely used open-source weather forecasting and re-analysis model that can be configured by the user according to the specific needs of each study. WRF has a proven record of producing high-resolution meteorological simulations down to a scale of hundreds of meters. It can be installed and operated in parallel mode using multi-processor computing resources, and therefore computational time can be greatly reduced for high-resolution simulations covering large geographical areas, such as PaTrop. Keeping the spatial resolution of the inner domain at 1 km, a series of re-analysis runs were produced to demonstrate the best possible configuration for our study, which is an approach supported theoretically and practically to tackle uncertainties in high-resolution modelling [25,26]. The parametric analysis was performed for a two-week period (17–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. The optimum
model configuration which resulted from the analysis was subsequently employed to provide the data for the main part of the PaTrop experiment. In order to establish the correlation between the observed and predicted time series, the following metrics are used, which give us a quantitative indication of the prediction skill of different model schemes with respect to the GNSS network dataset:

$$ \text{Mean Bias (MB): } \bar{\sigma} = \frac{\sum(f_i - o_i)}{N} \qquad (1) $$

$$ \text{Mean Absolute Bias (MAB): } |\bar{\sigma}| = \frac{\sum(|f_i - o_i|)}{N} \qquad (2) $$

$$ \text{Root Mean Square Error (RMSE): } RMSE = \sqrt{\frac{\sum (f_i - o_i)^2}{N}} \qquad (3) $$

$$ \text{Pearson Correlation Coefficient (R): } r = \frac{\frac{1}{N} \sum_i (f_i - \bar{f})(o_i - \bar{o})}{\sigma_f \sigma_o} \qquad (4) $$

where $f_i$ denotes the model value, $o_i$ the observational value, $N$ is the number of pairs in the examined time series, $\sigma_f$ is the standard deviation of the model values and $\sigma_o$ the standard deviation of the observations.

## 2.2. Data Processing

Processing of the GNSS data started after the end of the experimental campaign and included calculations of zenith tropospheric delays (ZTDs) every 30 min using the JPL NASA Precise Point Positioning (PPP) GIPSY-OASIS 6.4 software. High-precision International GNSS Service (IGS) final orbit data were used, while the quality of the sites ensured that multipath effects would not bias the estimated ZTDs. Static tropospheric processing with no mapping function was used, and a priori ZHD was estimated from an elevation-dependent function. It was intended to use a simple set of parameters in order to compare the WRF model output with a geometrical GNSS calculation which does not assimilate data from external sources (i.e., weather models). The exact settings of the processing protocol used are listed in Table 1.

<table><caption>**Table 1.** CRL tropospheric solution settings—all stations.</caption><thead><tr><th>Software</th><th>Gipsy 6.4</th></tr></thead><tbody><tr><td>GNSS</td><td>GPS and GLONASS</td></tr><tr><td>Troposphere estimated parameters</td><td>ZTD (5 min)</td></tr><tr><td>Ionosphere</td><td>HOI included</td></tr><tr><td>Ocean tides</td><td>FES2004</td></tr><tr><td>ZTD timestamp</td><td>hh:00 and hh:30</td></tr><tr><td>Elevation cut-off</td><td>5</td></tr></tbody></table>

With regards to the processing of WRF tropospheric data, values of surface pressure ($P_s$), as well as air temperature and water vapour pressure ($T, p_o$) in the vertical column (for each of the 45 vertical pressure levels) were derived every 30 min, at hh:00 and hh:30 h, at the nearest 1 km grid point from each GNSS station, to coincide with the observational time series. An automatic routine calculates the “dry” and “wet” delay terms separately, from the three parameters ($P_s, T, p_v$), together with the point elevation and layer heights, finally adding up the two to calculate the ZTD value [27–29]:

$$ \text{ZTD (mm)} = \text{ZHD} + \text{ZWD} \qquad (5) $$

With:

$$ \text{ZHD} = (2.2779 \pm 0.0024) P_s / f(\lambda, \text{H}) \qquad (6) $$
where $P_s$ is the total pressure (mbar) at the Earth's surface, and:

$$f(\lambda, H) = (1 - 0.00266 \cos2\lambda - 0.00028H) \quad (7)$$

accounts for the variation in gravitational acceleration with latitude $\lambda$ and the height $H$ of the surface above the ellipsoid (in km).

And:

$$ZWD = (0.382 \pm 0.004)K^2 mbar^{-1} \int \left( \frac{p_v}{T^2} \right) dz \quad (8)$$

where $p_v$ is the water vapour pressure (mbar), and $T$ the air temperature (K), integrated along the zenith path $z$.

WRF derived ZTDs were calculated at the exact elevation of the GNSS receiver, by vertically interpolating these parameters, thus minimising errors due to vertical height differences between the two datasets. Figure 2 illustrates graphically the GNSS-WRF geometry used in the calculation of the respective ZTD values.

<figure><img src="image_6.png" /><figcaption>Figure 2. Example of GNSS-WRF geometry for ANOC Station (h = 1020 m). Dark blue dotted lines depict the ZTD path calculated by WRF at the four nearest grid points, light blue lines depict the GNSS slant paths measured by the receiver at the time of acquisition and light blue dotted line depicts the ZTD path calculated by the GIPSY tropospheric processing.</figcaption></figure>

### 2.3. Model Configuration and Parameterization of Physical Components

For the high-resolution dynamical downscaling simulation performed with WRF v 3.7.1 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; two-way nested, i.e., feedback from nest to its parent domain. The vertical layer distribution consists of 45 sigma levels up to a height of about 20 km (0.1 hPa), with bottom layers being more densely populated. Boundary conditions for the model initialization were taken from the ERA-Interim global climate re-analysis database, with a 75 km horizontal resolution, 35 vertical layers and 6 h temporal resolution. The model was initiated every day from the ERA-Interim input data at 18:00 local time, producing 30 h simulations with the first 6 h being spin-up time. Model output was recorded every 30 min, from which Zenith Hydrostatic Delays (ZHD) and Zenith Wet Delays (ZWD) were calculated as previously described. 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.

For the model physical and dynamical parameterization, five different schemes were tested, in order to evaluate each scheme for its forecasting skill. There have been numerous studies validating the output of different model configurations with observations under specific conditions [25,30–32], showing that globally there is no optimal scheme, but rather different schemes produce better results with respect to application, domain, season, variable, etc. Therefore, a model parameterization test was performed for a two-week period (17–29 June 2016), during which the output of the five different model configurations was tested against GNSS tropospheric measurements from 16 permanent stations in the study area (d04).

The schemes were selected based on existing studies where similar high-resolution WRF simulations were used. All five schemes use the same parameterization for radiation physics (shortwave and longwave) and cumulus convection. In the first three schemes (MOD1, MOD2 and MOD3), cumulus convection is modelled in the 27 km domain according to the Kain–Fritsch scheme [33,34]. The cumulus scheme is not activated for the 9, 3 and 1 km domains, because at higher resolution the model can theoretically resolve convection explicitly [34]. In schemes MOD4 and MOD5, cumulus convection was turned on for d02 (9 km) in order to test the effect of cumulus parameterization in a smaller domain. Convection plays an important role for cloud formation and is controlled by micro-scale processes such as mixtures of updrafts and downdrafts. These simulated convective features are less distinguishable as model resolution becomes coarser; therefore, parameterization becomes necessary, although computationally demanding. Furthermore, in locations such as the Gulf of Corinth, where cloud formation is strongly influenced by the intense topography (land-sea contrasts and mountainous features), it is expected that cumulus parameterization in the 9 km domain will better represent the effects of sub-grid scale processes on the grid variables, particularly in the case of squall line formation, thunderstorms and other strong convection events [35].

Longwave radiation was simulated by the RRTM scheme [36] with the default diffusion scheme selected. The shortwave radiation was simulated by the Dudhia scheme [37]. Both are typical schemes for high-resolution WRF simulations. With respect to microphysics, land surface and planetary boundary layer options, the model was originally configured with the basic options for a less computationally demanding ERA-Interim dynamical downscaling at 1 km horizontal resolution. In the second scheme (MOD2), the microphysics scheme was changed to Morrison double-moment [38], in order to investigate the effect of a more complex model which uses double-moment ice, snow, rain and graupel for cloud-resolving simulations [25,39]. In the third scheme (MOD3), the same configuration as MOD2 was used, with a change in land surface parameters, in which the Pleim–Xiu land surface model and Pleim–Xiu surface layer physics was used [40,41]. In the fourth scheme (MOD4), the same configuration as MOD3 was used, this time with the cumulus scheme activated in the second domain (9 km), as explained above. Finally, the fifth scheme (MOD5) used the parameterization of Zhang et al. [42], as it refers to a study with similar dynamical downscaling characteristics (input dataset and horizontal resolution), performed over a complex terrain with land-sea contrasts (Hong Kong island). The parameters changed compared with the previous schemes were microphysics (SBU-YLin model was used, a sophisticated scheme that has ice, snow and graupel processes, suitable for real-data high-resolution simulations) [43], PBL scheme (YSU used in conjunction with SBU-YLin), and surface layer scheme (MM5 similarity). The NOAH model was used for land surface. Again, cumulus convection is activated in the second domain (9 km). Table 2 lists the parameterization used in each scheme.
<table><caption>**Table 2.** WRF parameterization options used for the PaTrop sensitivity analysis.</caption><thead><tr><th></th><th>MOD1</th><th>MOD2</th><th>MOD3</th><th>MOD4</th><th>MOD5</th></tr></thead><tbody><tr><td>Microphysics (mp)</td><td>WSM3</td><td>Morrison</td><td>Morrison</td><td>Morrison</td><td>SBU-YLin</td></tr><tr><td>Land surface (sf)</td><td>NOAH</td><td>NOAH</td><td>Pleim-Xiu</td><td>Pleim-Xiu</td><td>NOAH</td></tr><tr><td>Surface layer physics (sfclay)</td><td>Monin-Obukhov</td><td>Monin-Obukhov</td><td>Pleim-Xiu</td><td>Pleim-Xiu</td><td>MM5 similarity</td></tr><tr><td>Radiation physics (sw)</td><td>Dudhia</td><td>Dudhia</td><td>Dudhia</td><td>Dudhia</td><td>Dudhia</td></tr><tr><td>Radiation physics (lw)</td><td>RRTM</td><td>RRTM</td><td>RRTM</td><td>RRTM</td><td>RRTM</td></tr><tr><td>Planetary boundary layer physics (pbl)</td><td>MYJ</td><td>MYJ</td><td>ACM2</td><td>ACM2</td><td>YSU</td></tr><tr><td>Cloud physics (cu)</td><td>Kain-Fritsch at 27 km</td><td>Kain-Fritsch at 27 km</td><td>Kain-Fritsch at 27 km</td><td>Kain-Fritsch at 27 and 9 km</td><td>Kain-Fritsch at 27 and 9 km</td></tr></tbody></table>

# 3. Results

## 3.1. *Parametric Analysis and Evaluation 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 forecasting skill of different WRF parameterization schemes. 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. The selection of the test period for the parametric analysis was based on weather conditions (i.e., high temperatures and the occurrence of a convective storm event during the last days of June), which theoretically would produce a high degree of ZTD variability. Out of the 19 GNSS stations of the PaTrop network, 16 were fully operational during this period, while three (GALA, PSAT and AIGI) were not included in the analysis. Time series of WRF MOD1-MOD5 vs. GNSS ZTDs at 4 PaTrop stations (KALA, MESA, EYPA, PSAR) for the test period (17–29 June 2016) are presented in Figure 3, and corresponding bias plots of predicted (WRF MOD5) minus observed (GNSS) ZTDs are shown in Figure 4. Temporal resolution is 30 min and all times are in UTC.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Observed and simulated time series of WRF MOD1-MOD5 vs. GNSS ZTDs at 4 PaTrop stations, KALA (a), MESA (b), EYPA (c) and PSAR (d) over the test period (17–29 June 2016).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Bias plots of predicted (WRF MOD5) minus observed (GNSS) ZTDs at 4 PaTrop stations, KALA (a), MESA (b), EYPA (c) and PSAR (d) over the test period (17–29 June 2016). The red line indicates the mean bias.</figcaption></figure>

Results at most stations exhibit a strong correlation between WRF derived ZTDs and GNSS values, with the exception of the last two days of the experiment (28–29/6) where the model shows significant deviation from the observed data. All schemes seem to follow the overall trend (including the storm event of 25–26 June). The Correlation Coefficient R at the 16 locations, for the test period, ranges from 0.57 at LIDO to 0.85 at PAT0, both with the MOD5 scheme (Table 3). Overall, MOD3, MOD4 and MOD5 exhibit the strongest correlation, with MOD5 having the highest average R score (0.74).

<table><caption>**Table 3.** Pearson Correlation Co-efficient and RMSE results for all schemes (17/6–29/6/2016).</caption><thead><tr><th rowspan="2">Station</th><th rowspan="2">Elevation (m)</th><th colspan="6">R</th><th colspan="6">RMSE (mm)</th></tr><tr><th>MOD1</th><th>MOD2</th><th>MOD3</th><th>MOD4</th><th>MOD5</th><th>MOD1</th><th>MOD2</th><th>MOD3</th><th>MOD4</th><th>MOD5</th></tr></thead><tbody><tr><td>MESO</td><td>2</td><td>0.73</td><td>0.74</td><td>0.78</td><td>0.77</td><td>0.81</td><td>30.5</td><td>30.9</td><td>28.3</td><td>26.1</td><td>25.1</td></tr><tr><td>XILI</td><td>4</td><td>0.70</td><td>0.71</td><td>0.76</td><td>0.78</td><td>0.77</td><td>32.4</td><td>32.6</td><td>30.3</td><td>26.6</td><td>26.7</td></tr><tr><td>VALI</td><td>9</td><td>0.70</td><td>0.68</td><td>0.70</td><td>0.71</td><td>0.73</td><td>25.2</td><td>26.4</td><td>25.4</td><td>24.4</td><td>24.4</td></tr><tr><td>LAMB</td><td>10</td><td>0.70</td><td>0.71</td><td>0.74</td><td>0.75</td><td>0.76</td><td>25.2</td><td>25.4</td><td>23.1</td><td>22.8</td><td>22.9</td></tr><tr><td>TRIZ</td><td>25</td><td>0.63</td><td>0.61</td><td>0.70</td><td>0.70</td><td>0.70</td><td>27.0</td><td>28.4</td><td>26.6</td><td>24.8</td><td>24.9</td></tr><tr><td>PSAR</td><td>55</td><td>0.69</td><td>0.68</td><td>0.75</td><td>0.73</td><td>0.73</td><td>26.7</td><td>28.0</td><td>26.4</td><td>24.1</td><td>23.1</td></tr><tr><td>PAT0</td><td>91</td><td>0.77</td><td>0.78</td><td>0.81</td><td>0.84</td><td>0.85</td><td>29.9</td><td>30.5</td><td>25.0</td><td>24.6</td><td>27.6</td></tr><tr><td>ARSA</td><td>115</td><td>0.74</td><td>0.75</td><td>0.79</td><td>0.82</td><td>0.83</td><td>27.0</td><td>27.5</td><td>23.0</td><td>22.1</td><td>23.3</td></tr><tr><td>EYPA</td><td>166</td><td>0.67</td><td>0.68</td><td>0.74</td><td>0.75</td><td>0.75</td><td>27.3</td><td>27.6</td><td>24.7</td><td>22.8</td><td>22.8</td></tr><tr><td>ROD3</td><td>452</td><td>0.71</td><td>0.71</td><td>0.74</td><td>0.76</td><td>0.77</td><td>28.4</td><td>28.9</td><td>26.7</td><td>24.8</td><td>26.7</td></tr><tr><td>MESA</td><td>477</td><td>0.67</td><td>0.69</td><td>0.74</td><td>0.74</td><td>0.74</td><td>28.5</td><td>28.9</td><td>27.3</td><td>25.4</td><td>25.1</td></tr><tr><td>LIDO</td><td>550</td><td>0.60</td><td>0.60</td><td>0.67</td><td>0.60</td><td>0.57</td><td>29.8</td><td>30.0</td><td>29.5</td><td>28.0</td><td>26.6</td></tr><tr><td>KOUN</td><td>564</td><td>0.71</td><td>0.69</td><td>0.71</td><td>0.73</td><td>0.75</td><td>25.7</td><td>27.0</td><td>25.0</td><td>23.5</td><td>23.1</td></tr><tr><td>KALA</td><td>716</td><td>0.66</td><td>0.66</td><td>0.69</td><td>0.70</td><td>0.73</td><td>22.3</td><td>22.6</td><td>20.1</td><td>20.6</td><td>21.4</td></tr><tr><td>KRIN</td><td>758</td><td>0.68</td><td>0.66</td><td>0.67</td><td>0.70</td><td>0.72</td><td>22.8</td><td>24.5</td><td>23.0</td><td>21.7</td><td>22.2</td></tr><tr><td>ANOC</td><td>1020</td><td>0.66</td><td>0.67</td><td>0.72</td><td>0.72</td><td>0.70</td><td>35.5</td><td>36.5</td><td>35.1</td><td>32.2</td><td>31.4</td></tr><tr><td>AVG16</td><td></td><td>0.69</td><td>0.69</td><td>0.73</td><td>0.73</td><td>0.74</td><td>27.8</td><td>28.5</td><td>26.2</td><td>24.6</td><td>24.4</td></tr></tbody></table>

The model in general seems to slightly underestimate the observed data, especially for high ZTD values. Mean bias (MB) values indicate that, on average, this offset ranges from -14.9 mm (MOD1) to -11.6 mm (MOD5). Mean absolute bias (MAB) is a measure of the average absolute error between the two time series. Mean absolute bias values at
the 16 locations range from about 15 mm at KALA to 26 mm at ANOC. Schemes with more sophisticated physical parameterization, better suited for high-resolution re-analysis (MOD4 and MOD5) exhibit the lowest overall MAB values (19.5 mm). Similarly, root mean square error (RMSE) values range from 21 mm at KALA to 31 mm at ANOC. Figure 5 illustrates the RMSE spread between the 16 points, for each model configuration. It is shown that schemes with a more complex physical parameterization (MOD4 and MOD5) exhibit lower RMSE values, as well as narrower spread than initial WRF configurations. Results also indicate that the RMSE distribution is more homogeneous among coastal and inland stations (blue and green colours), than stations in mountainous locations (orange and red). This finding could be related to the specific meteorological conditions which characterize regions with intense topography, as well as the fact that input data (e.g., measurements from surface stations) injected into the model as initial conditions, may be sparse in certain remote locations, leading to a poorer model prediction skill in these areas.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** RMSE distribution at the 16 stations, for each model configuration. Stations are colour-coded, according to their location: blue—coastal; green—inland; orange—upland north; red—upland south.</figcaption></figure>

A graphical summary of validation metrics for MOD5 is found in Figure 6. We observe that although R is decreasing with increasing station elevation (as previously discussed), MAB and RMSE are not significantly increasing, which is a result of lower absolute ZTD values at higher altitudes where the tropospheric layer is thinner.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Distribution of model validation metrics (MAB, RMSE and PPC) per station elevation h—MOD5.</figcaption></figure>
Under turbulent and humid atmospheric conditions (e.g., afternoon of 25 June), different schemes behave in different manners. Two separate groups can be distinguished (MOD1, MOD2 and MOD3 vs. MOD4 and MOD5), with the second group exhibiting a better prediction skill, possibly as a result of the Kain–Fritsch cumulus convection scheme being turned on at the 9 km domain. A significant deviation event can also be distinguished during the second half of 28 June, where the model produces a sudden “drop” in all stations (of the order of 60–120 mm) in relation to the observed ZTDs.

## 3.2. Validation of WRF Derived Tropospheric Delay Maps with GNSS ZTD Measurements for the PaTrop Period (January–December 2016)

Following the initial configuration of the WRF model and parameterization based on the short-scale analysis, as described in the previous section, the main part of the PaTrop experiment extends into a whole year of validation of model re-analysis output with the use of observational tropospheric data from the CRL GNSS network in the Western Gulf of Corinth. The WRF scheme which provided the best simulation results based on the parametric test was selected (MOD5), having a more complex physical parameterization, and hence being better suited for high-resolution re-analysis simulations. This includes: 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, MOD5 uses 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. The MOD5 output, represented as Zenith Total Delay (ZTD) values calculated from specific atmospheric parameters (surface pressure, air temperature and water vapour profiles) over the entire 1 × 1 km grid is validated against a dataset of GNSS derived ZTD values, providing point measurements at the 19 points where the stations are located.

Macroscopically, the two ZTD datasets are closely correlated in all 19 locations (Figure 7), with values peaking during the warm months (July–September) and subduing during the cold period (December–March), as expected (ZTD is proportional to surface pressure and air temperature which both increase during summer). It is also visible from the oscillating pattern of the time series that in absolute terms the tropospheric delay signal shows a higher degree of variability during the second half of the year than during the first half, as a result of a combination of more intense temperature and water vapour fluctuations. Bias plots of predicted (WRF) minus observed (GNSS) ZTDs (Figure 8) indicate a slightly stronger bias during the warm period. For a more detailed analysis of the correlation between the two datasets, the same statistical indices are used as in the physical parameterization study (R, MB, MAB, RMSE), which give us a quantitative indication of the accuracy and variability of the model prediction with respect to the dataset of the GNSS monitoring network. A summary of the regression statistics for the whole period of study is presented in Table 4.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Observed and simulated time series of WRF MOD5 vs. GNSS ZTDs at 4 PaTrop stations, KALA (a), MESA (b), EYPA (c) and PSAR (d) over January–December 2016.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Bias plots of predicted (WRF MOD5) minus observed (GNSS) ZTDs at 4 PaTrop stations, KALA (a), MESA (b), EYPA (c) and PSAR (d) over January–December 2016. The yellow line indicates the mean bias (MB).</figcaption></figure>
<table><caption>**Table 4.** Statistical indices of complete WRF ZTS vs. GNSS ZTD time series—January–December 2016.</caption><thead><tr><th>Station</th><th>Elevation ASL (m)</th><th>MB</th><th>MAB</th><th>RMSE</th><th>R</th></tr></thead><tbody><tr><td>ANOC</td><td>1020</td><td>-23.8</td><td>24.5</td><td>28.3</td><td>0.92</td></tr><tr><td>ARSA</td><td>115</td><td>-16.4</td><td>19.2</td><td>23.9</td><td>0.92</td></tr><tr><td>AIGI</td><td>142</td><td>-18.1</td><td>20.7</td><td>25.3</td><td>0.92</td></tr><tr><td>GALA</td><td>33</td><td>-19.6</td><td>22.1</td><td>27.3</td><td>0.93</td></tr><tr><td>EYPA</td><td>166</td><td>-19.5</td><td>21.3</td><td>25.7</td><td>0.93</td></tr><tr><td>KALA</td><td>716</td><td>-11.1</td><td>14.9</td><td>19.1</td><td>0.92</td></tr><tr><td>KOUN</td><td>564</td><td>-16.8</td><td>19.1</td><td>23.4</td><td>0.92</td></tr><tr><td>KRIN</td><td>758</td><td>-19.4</td><td>21.0</td><td>26.4</td><td>0.91</td></tr><tr><td>LAMB</td><td>10</td><td>-14.0</td><td>19.0</td><td>22.4</td><td>0.91</td></tr><tr><td>LIDO</td><td>550</td><td>-12.1</td><td>16.1</td><td>20.3</td><td>0.92</td></tr><tr><td>MESA</td><td>477</td><td>-22.5</td><td>23.7</td><td>27.9</td><td>0.93</td></tr><tr><td>MESO</td><td>2</td><td>-18.3</td><td>21.4</td><td>26.7</td><td>0.91</td></tr><tr><td>PAT0</td><td>91</td><td>-25.5</td><td>26.5</td><td>30.9</td><td>0.93</td></tr><tr><td>PSAR</td><td>55</td><td>-15.1</td><td>18.1</td><td>22.7</td><td>0.93</td></tr><tr><td>PSAT</td><td>19</td><td>-28.2</td><td>29.0</td><td>32.9</td><td>0.92</td></tr><tr><td>ROD3</td><td>452</td><td>-20.8</td><td>22.2</td><td>26.5</td><td>0.92</td></tr><tr><td>TRIZ</td><td>25</td><td>-21.4</td><td>22.9</td><td>27.4</td><td>0.93</td></tr><tr><td>VALI</td><td>9</td><td>-14.5</td><td>18.1</td><td>22.9</td><td>0.93</td></tr><tr><td>XILI</td><td>4</td><td>-28.0</td><td>28.6</td><td>33.1</td><td>0.93</td></tr><tr><td>AVG19</td><td></td><td>-19.2</td><td>21.5</td><td>25.9</td><td>0.92</td></tr></tbody></table>

* The Pearson correlation co-efficient R measures the extent to which two variables are linearly related. An R value of 1 means that the two variables are perfectly positively linearly related and that the points in the scatter plot lie exactly on a straight line (y = x). The R at the 19 locations, for the entire annual time series, ranges from 0.91 to 0.93, i.e., it is fairly uniform, indicating that the model's variability matches the variability of the observed tropospheric delay about 90% of the time.

* Mean bias (MB) is a measure of the accuracy of the model's ZTD output with respect to the observational dataset. Mean bias values for ZTD (GNSS-WRF), range from -11.1 mm (KALA station) to -28.2 mm (PSAT station), and indicate that the model tends to slightly underestimate the tropospheric ZTD as compared to the GNSS derived values. This finding is in-line with similar WRF evaluation studies [27,33] reporting consistently negative differences in relative humidity (a primary physical parameter in calculating the ZTD) with respect to ground observations in high-resolution WRF re-analysis scenarios, which are attributed to differences in vertical mixing strength and entrainment.

* Mean absolute bias (MAB) and root mean square error (RMSE) are both a measure of the absolute error between the two time series and are particularly useful, as the correction of the tropospheric component in InSAR interferograms is dependent on the model's capability to produce high-resolution differential meteograms of tropospheric delay with the minimum absolute error (of the order of magnitude of one interferometric phase cycle $\pi$). Mean absolute bias values at the 19 locations range from 14.9 mm (KALA station) to 29.0 mm (PSAT station), with RMSE values covering a similar range.

A further look into the model validation metrics for the entire PaTrop period reveals some additional trends:

* As shown in Figure 9a, the RMSE seems to be independent of the horizontal distance $s$ between the GNSS station and the nearest WRF grid point where the calculation of the predicted ZTD is performed. Therefore, we can conclude that the horizontal resolution of 1 km used for the WRF simulation is adequate.

* With respect to station elevation $h$, a small reduction of MB in terms of its absolute value is evident with increasing $h$, as expected, due to smaller ZTD values (Figure 9b).
Out of the 19 stations, the three highest mean negative biases are in XILI, PSAT and PAT0 (elevations ASL 4, 19 and 91 m), while the two lowest are in LIDO and KALA (550 and 716 m). The graphical summary of validation metrics for the entire period (Figure 10) also reveals that while R is fairly constant with increasing station elevation, MAB and RMSE exhibit a reduction.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Plots of: RMSE vs. horizontal distances between WRF grid point and GNSS station (a); MB vs. GNSS station elevation h (b).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Distribution of model validation metrics (MAB, RMSE and R) per station elevation h.</figcaption></figure>

### 3.3. *Seasonal Characteristics of WRF vs. GNSS ZTD and Evaluation of Model Performance*

Beyond the annual time series analysis and evaluation, it is useful to perform a more in-depth investigation of seasonal trends, which can provide an important insight into the model's forecast skill under different meteorological conditions. The year is therefore divided into four seasons (S1-S4) based on their distinct climatological characteristics as follows:

* Winter S1: January–March.

* Spring S2: April–June.

* Summer S3: July–September.

* Autumn S4: October–December.

Table 5 lists the corresponding model validation metrics per station and for each season. In general, model performance is better during autumn (S4), as correlation with the observed ZTD time series is high at all stations (R = 0.93 with range 0.91–0.94) and MAB and RMSE are lower (average MAB = 20.5; RMSE = 24.6). In spring (S2), WRF reproduces
88% of the ZTD variability ($R = 0.88$ with range 0.85–0.89) and exhibits a slightly higher bias than autumn (average $MAB = 21.1$; RMSE = 25.7). In winter (S1), the model reproduces 85% of the ZTD variability; however, with a wider correlation range between stations ($R = 0.85$ with range 0.75–0.88), and bias is similar to S2 (average $MAB = 21.2$; RMSE = 25.2). Model forecasting skill deteriorates during summer (S3), where correlation is weaker (average $R = 0.83$) and bias indicators are higher ($MAB = 22.8$, RMSE = 28.0). These results are in-line with seasonal climate characteristics in the region of Western Greece, as stable weather conditions and frontal precipitation patterns which prevail during autumn are well simulated by the model, whereas the combination of high temperatures and turbulent conditions during some summer days, resulting in convective storms are more difficult to predict. A graphical illustration of the RMSE variability for each season, at the location of each station, is seen in Figure 11.

When we look at the monthly variability of the mean absolute bias (Figure 12), a distinct “peak” in February and March is evident in most coastal and inland stations, possibly due to the winter storms which are common during this period. A second “peak” during the summer period (June–September) is explained by a combination of high temperatures and the occurrence of intense convective storm events, as discussed earlier. Stations which are located inland exhibit a more uniform monthly variability of the $MAB$ (co-variance is high), while coastal stations exhibit a more diverse variability, particularly during the “hot” season S3. In stations located upland, two separate groups are distinct: the first one including the Peloponnese stations (south) plus LIDO exhibit a smaller degree of seasonal $MAB$ variability, while ANOC and MESA (both located in the northern Aitoloakarnania region) exhibit higher seasonal variability (distinct “peaks” during S1 and S3 and higher $MAB$). We can further investigate local inhomogeneities between stations by calculating the correlation coefficient $R$ of their difference, i.e., between the GNSS $\Delta$ZTD and WRF $\Delta$ZTD of stations MESA and KALA, for seasons S1–S4. Thus, the annual and semi-annual component of ZTD variability is removed from the time series and the ability of the model to resolve small-scale inhomogeneities is demonstrated. Results indicate that model performance is better during the second half of the year, as WRF reproduces 54% of the observed $\Delta$ZTD variability between the two stations in summer (S3) and 52% in autumn (S4) ($R$ is 0.54 and 0.52 respectively). The correlation is slightly weaker during the first half, with $R = 0.48$ in winter (S1) and $R = 0.49$ in spring (S2).
<table><caption>**Table 5.** Seasonal statistical indices of WRF ZTS vs. GNSS ZTD time series—2016.</caption><thead><tr><th rowspan="2">Season</th><th colspan="3">S1</th><th colspan="3">S2</th><th colspan="3">S3</th><th colspan="3">S4</th></tr><tr><th>|σ|</th><th>RMSE</th><th>R</th><th>|σ|</th><th>RMSE</th><th>R</th><th>|σ|</th><th>RMSE</th><th>R</th><th>|σ|</th><th>RMSE</th><th>R</th></tr></thead><tbody><tr><td>ANOC</td><td>24.9</td><td>28.0</td><td>0.85</td><td>25.1</td><td>28.9</td><td>0.89</td><td>25.8</td><td>30.5</td><td>0.82</td><td>22.3</td><td>25.6</td><td>0.93</td></tr><tr><td>ARSA</td><td>19.8</td><td>24.2</td><td>0.86</td><td>18.3</td><td>22.8</td><td>0.89</td><td>20.6</td><td>25.8</td><td>0.84</td><td>18.1</td><td>22.4</td><td>0.93</td></tr><tr><td>AIGI</td><td>19.6</td><td>23.5</td><td>0.81</td><td>20.4</td><td>25.2</td><td>0.89</td><td>22.2</td><td>27.6</td><td>0.84</td><td>20.8</td><td>24.6</td><td>0.93</td></tr><tr><td>GALA</td><td>22.1</td><td>27.1</td><td>0.84</td><td>22.3</td><td>27.4</td><td>0.87</td><td>23.1</td><td>28.6</td><td>0.86</td><td>-</td><td>-</td><td>-</td></tr><tr><td>EYPA</td><td>20.9</td><td>24.6</td><td>0.87</td><td>20.3</td><td>24.9</td><td>0.89</td><td>22.9</td><td>27.9</td><td>0.84</td><td>21.1</td><td>25.0</td><td>0.94</td></tr><tr><td>KALA</td><td>14.2</td><td>18.3</td><td>0.82</td><td>15.0</td><td>19.0</td><td>0.89</td><td>16.7</td><td>21.5</td><td>0.83</td><td>13.4</td><td>17.0</td><td>0.93</td></tr><tr><td>KOUN</td><td>18.9</td><td>22.6</td><td>0.83</td><td>17.7</td><td>22.0</td><td>0.89</td><td>20.2</td><td>25.4</td><td>0.84</td><td>19.6</td><td>23.4</td><td>0.93</td></tr><tr><td>KRIN</td><td>22.7</td><td>27.9</td><td>0.75</td><td>20.7</td><td>26.6</td><td>0.87</td><td>20.4</td><td>25.7</td><td>0.82</td><td>18.3</td><td>23.3</td><td>0.91</td></tr><tr><td>LAMB</td><td>18.5</td><td>22.8</td><td>0.84</td><td>17.5</td><td>22.4</td><td>0.87</td><td>18.4</td><td>23.1</td><td>0.84</td><td>16.6</td><td>21.4</td><td>0.93</td></tr><tr><td>LIDO</td><td>16.8</td><td>20.6</td><td>0.88</td><td>16.5</td><td>20.9</td><td>0.85</td><td>16.2</td><td>21.1</td><td>0.85</td><td>15.5</td><td>19.2</td><td>0.93</td></tr><tr><td>MESA</td><td>22.5</td><td>26.0</td><td>0.87</td><td>22.6</td><td>26.8</td><td>0.88</td><td>26.3</td><td>31.0</td><td>0.82</td><td>23.1</td><td>27.1</td><td>0.93</td></tr><tr><td>MESO</td><td>20.0</td><td>23.9</td><td>0.88</td><td>20.9</td><td>25.2</td><td>0.89</td><td>25.3</td><td>32.2</td><td>0.71</td><td>20.3</td><td>24.9</td><td>0.94</td></tr><tr><td>PAT0</td><td>26.0</td><td>29.4</td><td>0.88</td><td>26.1</td><td>30.4</td><td>0.88</td><td>28.6</td><td>33.1</td><td>0.82</td><td>26.4</td><td>30.4</td><td>0.93</td></tr><tr><td>PSAR</td><td>16.9</td><td>21.0</td><td>0.85</td><td>17.8</td><td>22.4</td><td>0.89</td><td>20.3</td><td>25.5</td><td>0.85</td><td>17.3</td><td>21.6</td><td>0.94</td></tr><tr><td>PSAT</td><td>30.1</td><td>33.2</td><td>0.86</td><td>28.2</td><td>32.6</td><td>0.87</td><td>28.0</td><td>33.1</td><td>0.84</td><td>29.1</td><td>32.8</td><td>0.94</td></tr><tr><td>ROD3</td><td>19.3</td><td>23.8</td><td>0.84</td><td>23.3</td><td>27.4</td><td>0.89</td><td>24.6</td><td>29.3</td><td>0.85</td><td>21.6</td><td>26.4</td><td>0.92</td></tr><tr><td>TRIZ</td><td>22.5</td><td>26.4</td><td>0.87</td><td>22.8</td><td>27.2</td><td>0.89</td><td>24.9</td><td>30.2</td><td>0.84</td><td>22.5</td><td>26.5</td><td>0.94</td></tr><tr><td>VALI</td><td>17.5</td><td>21.7</td><td>0.84</td><td>17.8</td><td>22.8</td><td>0.89</td><td>20.0</td><td>25.5</td><td>0.84</td><td>17.2</td><td>21.5</td><td>0.94</td></tr><tr><td>XILI</td><td>30.5</td><td>33.6</td><td>0.86</td><td>27.5</td><td>32.7</td><td>0.88</td><td>29.8</td><td>35.0</td><td>0.84</td><td>26.3</td><td>30.6</td><td>0.94</td></tr><tr><td>AVG19</td><td>21.2</td><td>25.2</td><td>0.85</td><td>21.1</td><td>25.7</td><td>0.88</td><td>22.8</td><td>28.0</td><td>0.83</td><td>20.5</td><td>24.6</td><td>0.93</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** RMSE per station and season (S1 top left, S2 top right, S3 bottom left, S4 bottom right).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Monthly variability of MAB at the 19 PaTrop stations, classified per terrain type: coastal (a); inland (b); upland (red = north, orange = south) (c).</figcaption></figure>

The overall model performance with respect to estimating the ZTD parameter for the successful tropospheric correction of InSAR interferograms is evaluated by setting the acceptable bias range as the amplitude of one Sentinel-1 C-band cycle phase $\pi$ (equal to $\lambda/2$, where $\lambda$ is the SAR signal wavelength), i.e., $\pm 23$ mm of tropospheric delay, when projected to the zenithal distance (considering a mean incidence angle of 35°). This allows us to determine the percentage of $\sigma$ lying within the accepted range per season (Table 6), which in theory represents the probability that the model successfully detects and removes the tropospheric delay.

<table><caption>**Table 6.** $\sigma$ values within error range—WRF vs. GNSS ZTDs average of all stations.</caption><thead><tr><th>Season</th><th>Ratio of σ above π (23 mm)</th><th>Ratio of σ below -π (-23 mm)</th><th>Ratio of σ between π and -π</th></tr></thead><tbody><tr><td>S1</td><td>&lt;0.01</td><td>0.39</td><td>0.61</td></tr><tr><td>S2</td><td>0.01</td><td>0.37</td><td>0.62</td></tr><tr><td>S3</td><td>0.02</td><td>0.41</td><td>0.57</td></tr><tr><td>S4</td><td>&lt;0.01</td><td>0.37</td><td>0.63</td></tr></tbody></table>

## 4. Discussion

Results of the physical parameterization analysis indicate that WRF configurations MOD4 and MOD5 overall exhibit a better prediction skill during the test period, with small differences between them. Turning on the Kain-Fritsch cumulus convection scheme on the 9 km domain improves model output, particularly during intense frontal events, as demonstrated by validation metrics of MOD4 in comparison with MOD3 (all physical parameters are the same except the cumulus scheme). Moreover, the use of more complex microphysics schemes (such as Morrison and SBU-YLin), surface layer (Pleim-Xiu and
MM5) and planetary boundary layer (ACM2 and YSU) models, although computationally demanding, further enhances the model prediction skill. Indeed, model setup and physical parameterization can affect the intrinsic water balance [44], resulting in a more accurate representation of precipitation and relative humidity fields [25,45]. Based on the results of the parametric analysis, MOD5 is selected as the optimum configuration for producing tropospheric delay maps for the entire PaTrop period (January–December 2016), having a slightly higher R, lower MB and better clustering of RMSE than MOD4.

On an annual basis, the statistical analysis of the results demonstrated that the correlation between predicted and observed ZTDs at the 19 stations is good throughout the year (correlation co-efficient R ranges from 0.91 to 0.93), with mean bias (MB) ranging from –11.1 mm (KALA station) to –28.2 mm (PSAT station), indicating that the model tends to slightly underestimate the tropospheric ZTD as compared to the GNSS derived values. Results are in-line with previous WRF high-resolution studies [25,46,47], which report consistent underestimation of RH at the ground level (ranging from 9 to 15%). With respect to the seasonal component, model performance is better during the autumn period (October–December), followed by the spring period (April–June). Correlation with the observed ZTD time series is high at all stations (correlation co-efficient R is 0.93 and 0.88, respectively) and MAB and RMSE values are low. On the contrary, model forecasting skill seems to deteriorate during summer (July–September), where correlation is weaker (R = 0.83) and MAB and RMSE values are higher (average MAB = 22.8, RMSE = 28.0). Again, results re-confirm research findings from similar comparisons of GNSS tropospheric products and LAMs in Europe [48–50], reporting a pronounced seasonal cycle of GNSS-NWP model bias and standard deviation with the largest errors in summer and smallest in winter. Setting the acceptable bias range at ±23 mm (equal to the amplitude of one Sentinel-1 C-band phase cycle when projected to the zenith distance), it is demonstrated that the model produces satisfactory results, with a percentage of ZTD values within the bias margin ranging from 57% in summer to 63% in autumn.

## 5. Conclusions

The experimental approach followed aims to maximise the potential of using LAMs for the calculation and removal of the tropospheric noise from InSAR data. Five different WRF parameterization schemes were tested, ranging from schemes with relatively simple physical and dynamical parameterization to more complex schemes which require longer computational times. A parametric analysis was performed for a two-week period (17–29 June 2016), during which the output of different model configurations, in terms of ZTD, was compared with GNSS tropospheric measurements from 16 permanent stations in the study area (d04). Results were tested for their statistical significance and demonstrated that the optimum WRF configuration to be used for the entire period of the PaTrop experiment was MOD5. Consequently, daily runs of the optimal scheme were performed for a one-year period (January–December 2016), and model output was again validated with the use of the observational GNSS dataset from the CRL network in the Western Gulf of Corinth.

Overall, the current study confirms the high potential and effectiveness of using high-resolution atmospheric modelling for calculating the ZTD component, and thus correcting atmospheric phase gradients in interferograms. The innovative elements that our methodology brings in comparison with previous studies are: (a) High-resolution re-analysis with WRF: our model is downscaled at 1 km horizontal resolution and therefore offers a detailed reconstruction of the 3-D tropospheric field over the study area. The model was locally configured and parameterized in the area of the western Gulf of Corinth, and numerous complex schemes were tested in order to demonstrate the optimal configuration at the specific location; (b) Validation of the model: WRF output was validated with the use of GNSS tropospheric data retrieved from a dense array of stations covering the study area. Model validation with vertical column data (GNSS zenith 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 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 are exhibited in the differential atmosphere as densely distributed short-wavelength phase gradients.

**Author Contributions:** Conceptualization, N.R. and P.B.; methodology, N.R. and P.B.; software, N.R., D.K. and P.B.; validation, N.R.; formal analysis, N.R.; investigation, N.R.; resources, N.R., D.K., P.E. and P.B.; data curation, N.R. and P.B.; writing—original draft preparation, N.R.; writing—review and editing, N.R., P.B., I.K., D.K., P.E., A.A.A. and A.R.; visualization, N.R.; supervision, P.B., I.K., A.A.A. and A.R.; project administration, P.B.; funding acquisition, P.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Not applicable.

**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.

**Acknowledgments:** This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility—ARIS—under project IDs PaTrop, PaTrop2 and Atmo-InSAR.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Skamarock, W.; Klemp, J.B.; Dudhia, J.; Gill, D.O.; Barker, D.; Duda, M.G.; Huang, X.-Y.; Wang, W. A description of the advanced research WRF version 3. In *NCAR Technical Note NCAR/TN-475 + STR*; University Corporation for Atmospheric Research: Boulder, CO, USA, 2008. [CrossRef]

2. Foster, J.; Brooks, B.; Cherubini, T.; Shacat, C.; Businger, S.; Werner, C.L. Mitigating atmospheric noise for InSAR using a high resolution weather model. *Geophys. Res. Lett.* **2006**, *33*, L16304. [CrossRef]

3. Foster, J.; Kealy, J.; Cherubini, T.; Businger, S.; Lu, Z.; Murphy, M. The utility of atmospheric analyses for the mitigation of artifacts in InSAR. *J. Geophys. Res. Solid Earth* **2013**, *118*, 748–758. [CrossRef]

4. Puysségur, B.; Michel, R.; Avouac, J.P. Tropospheric phase delay in interferometric synthetic aperture radar estimated from meteorological model and multispectral imagery. *J. Geophys. Res.* **2007**, *112*, B05419. [CrossRef]

5. Wadge, G.; Zhu, M.; Holley, R.; James, I.; Clark, P.; Wang, C.; Woodage, M. Correction of Atmospheric Delay Effects in Radar Interferometry Using a Nested Mesoscale Atmospheric Model. *J. Appl. Geophys.* **2010**, *72*, 141–149. [CrossRef]

6. Eff-Darwich, A.; Perez, J.; Fernandez, J.; Garcia-Lorenzo, B. Using a mesoscale meteorological model to reduce the effect of tropospheric water vapour from DInSAR data: A case study for the island of Tenerife. *Pure Appl. Geophys.* **2012**, *169*, 1425–1441. [CrossRef]

7. Kinoshita, Y.; Furuya, M.; Hobiger, T.; Ichikawa, R. Are numerical weather model outputs helpful to reduce tropospheric delay signals in InSAR data? *J. Geod.* **2013**, *87*, 267–277. [CrossRef]

8. Bekaert, D.P.S.; Walters, R.J.; Wright, T.J.; Hooper, A.; Parker, D.J. Statistical comparison of InSAR tropospheric correction techniques. *Remote Sens. Environ.* **2015**, *170*, 40–47. [CrossRef]

9. Delacourt, C.; Briole, P.; Achache, J.A. Tropospheric corrections of SAR interferograms with strong topography. *Appl. Etna. Geophys. Res. Lett.* **1998**, *25*, 2849–2852. [CrossRef]

10. Webley, P.W.; Bingley, R.M.; Dodson, A.H.; Wadge, G.; Waugh, S.J.; James, I.N. Atmospheric water vapour correction to InSAR surface motion measurements on mountains: Results from a dense GPS network on Mount Etna. *Phys. Chem. Earth* **2002**, *27*, 363–370. [CrossRef]

11. Li, Z.; Muller, J.; Cross, P.; Fielding, E. Interferometric synthetic aperture radar (InSAR) atmospheric correction: GPS, Moderate Resolution Imaging Spectroradiometer (MODIS), and InSAR integration. *J. Geophys. Res.* **2005**, *110*. [CrossRef]

12. Löfgren, J.; Björndahl, F.; Moore, A.; Webb, F.; Fielding, E.; Fishbein, E. Tropospheric correction for InSAR using interpolated ECMWF data and GPS Zenith Total Delay from the Southern California Integrated GPS Network. In Proceedings of the 2010 IEEE International Geoscience and Remote Sensing Symposium, Honolulu, HI, USA, 25–30 July 2010; pp. 4503–4506. [CrossRef]

13. Li, Z.; Fielding, E.; Cross, P.; Preusker, R. Advanced InSAR atmospheric correction: MERIS/MODIS combination and stacked water vapour models. *Int. J. Remote Sens.* **2009**, *30*, 3343–3363. [CrossRef]

14. Walters, R.; Elliott, J.; Li, Z.; Parsons, B. Rapid strain accumulation on the Ashkabad fault (Turkmenistan) from atmosphere-corrected InSAR. *J. Geophys. Res. Atmos.* **2013**, *118*, 1–17. [CrossRef]

15. Xu, W.; Li, Z.; Ding, X.; Feng, G.; Hu, D.; Long, J.; Yin, H.; Yang, E. Correcting atmospheric effects in ASAR interferogram with MERIS integrated water vapor data. *Chin. J. Geophys.* **2010**, *53*, 1073–1084.

16. Elliott, J.; Biggs, J.; Parsons, B.; Wright, T. InSAR slip rate determination on the Altyn Tagh Fault, northen Tibet, in the presence of topographically correlated atmopsheric delays. *Geophys. Res. Lett.* **2008**, *35*, L12309. [CrossRef]
17. Cavalié, O.; Lasserre, C.; Doin, M.; Peltzer, G.; Sun, J.; Xu, X.; Shen, Z. Measurement of interseismic strain across the Haiyuan fault (Gansu, China) by InSAR. Earth and Plan. Sci. Lett. **2008**, *275*, 246–257. [CrossRef]

18. Lin, Y.N.; Simons, M.; Hetland, E.A.; Muse, P.; DiCaprio, C.A. Multiscale approach to estimating topographically correlated propagation delays in radar interferograms. *Geochem. Geophys. Geosyst.* **2010**, *11*, Q09002. [CrossRef]

19. Doin, M.P.; Lasserre, C.; Peltzer, G.; Cavali, O.; Doubre, C. Corrections of stratified tropospheric delays in SAR interferometry: Validation with global atmospheric models. *J. Appl. Geophys.* **2009**, *69*, 35–50. [CrossRef]

20. Jolivet, R.; Grandin, R.; Lasserre, C.; Doin, M.P.; Peltzer, G. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. *Geophys. Res. Lett.* **2011**, *38*, L17311. [CrossRef]

21. Jolivet, R.; Agram, P.S.; Lin, N.Y.; Simons, M.; Doin, M.P.; Peltzer, G.; Li, Z. Improving InSAR geodesy using Global Atmospheric Models. *J. Geophys. Res. Solid Earth* **2014**, *119*, 2324–2341. [CrossRef]

22. Yun, Y.; Zeng, Q.; Green, B.W.; Zhang, F. Mitigating atmospheric effects in InSAR measurements through high-resolution data assimilation and numerical simulations with a weather prediction model. *Int. J. Remote Sens.* **2015**, *36*, 2129–2147. [CrossRef]

23. Briole, P.; Rigo, A.; Lyon-Caen, H.; Ruegg, J.C.; Papazissi, K.; Mitsakaki, C.; Balodimou, A.; Veis, G.; Hatzfeld, D.; Deschamps, A. Results from repeated Global Positioning System surveys between 1990 and 1995. *J. Geophys. Res.* **2000**, *105*, 25605–25625. [CrossRef]

24. Avallone, A.; Briole, P.; Agatza-Balodimou, A.M.; Billiris, H.; Charade, O.; Mitsakaki, C.; Nercessian, A.; Papazissi, K.; Paradissis, D.; Veis, G. Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area. *C.R. Geosci.* **2004**, *336*, 301–312. [CrossRef]

25. Kioutsioukis, I.; De Meij, A.; Jakobs, H.; Katragkou, E.; Vinuesa, J.F.; Kazantzidis, A. High resolution WRF ensemble forecasting for irrigation: Multi-variable evaluation. *Atmos. Res.* **2015**, *167*, 156–174. [CrossRef]

26. Garcia-Diez, M.; Fernandez, J.; Vautard, R. An RCM multi-physics ensemble over Europe: Multi-variable evaluation to avoid error compensation. *Clim. Dyn.* **2015**, *45*, 3141–3156. [CrossRef]

27. Bevis, M.; Businger, S.; Herring, T.; Rocken, C.; Anthes, R.; Ware, R. GPS meteorology: Remote sensing of atmospheric water vapor using the global positioning system. *J. Geophys. Res.* **1992**, *97*, 15787–15801. [CrossRef]

28. Elgered, G.; Davis, J.L.; Herring, T.A.; Shapiro, I.I. Geodesy by radio interferometry: Water vapor radiometry for estimation of the wet delay. *J. Geophys. Res.* **1991**, *96*, 6541–6555. [CrossRef]

29. Davis, J.L.; Herring, T.A.; Shapiro, I.I.; Rogers, A.E.; Elgered, G. Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. *Radio Sci.* **1985**, *20*, 1593–1607. [CrossRef]

30. Mooney, P.A.; Mulligan, F.J.; Fealy, R. Evaluation of the sensitivity of the weather research and forecasting model to parameterization schemes for regional climates of Europe over the period 1990–95. *J. Clim.* **2013**, *26*, 1002–1017. [CrossRef]

31. Kotlarski, S.; Keuler, K.; Christensen, O.B.; Colette, A.; Déqué, M.; Gobiet, A.; Goergen, K.; Jacob, D.; Lüthi, D.; van Meijgaard, E.; et al. Regional climate modeling on European scales: A joint standard evaluation of the EURO-CORDEX RCM ensemble. *Geosci. Model Dev.* **2014**, *7*, 1297–1333. [CrossRef]

32. Garcia-Diez, M.; Fernandez, J.; Fita, L.; Yague, C. Seasonal dependence of WRF model biases and sensitivity to PBL schemes over Europe. *Q. J. R. Meteorol. Soc.* **2013**, *139*, 501–514. [CrossRef]

33. Kain, J.S.; Fritsch, M. A one-dimensional entraining/detraining plume model and its application in convective parameterization. *J. Atmos. Sci.* **1990**, *47*, 2784–2802. [CrossRef]

34. Kain, J.S. The Kain–Fritsch convective parameterization: An update. *J. Appl. Meteorol.** *2004*, *43*, 170–181. [CrossRef]

35. Zheng, Y.; Alapaty, K.; Herwehe, J.; Delgenio, A.; Niyogi, D. Improving high-resolution weather forecasts using the weather research and forecasting (WRF) model with an updated Kain–Fritsch scheme. *Mon. Weather Rev.* **2015**, *144*, 150528113541009. [CrossRef]

36. Mlawer, E.J.; Taubman, S.J.; Brown, P.D.; Iacono, M.J.; Clough, S.A. Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *J. Geophys. Res.* **1997**, *102*, 16663–16682. [CrossRef]

37. Dudhia, J. Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model. *J. Atmos. Sci.* **1989**, *46*, 3077–3107. [CrossRef]

38. Morrison, H.; Thompson, G.; Tatarskii, V. Impact of cloud microphysics on the development of trailing stratiform precipitation in a simulated squall line: Comparison of one- and two-moment schemes. *Mon. Weather Rev.* **2009**, *137*, 991–1007. [CrossRef]

39. Warrach-Sagi, K.; Schwitalla, T.; Wulfmeyer, V.; Bauer, H. Evaluation of a climate simulation in Europe based on the WRF–NOAH model system: Precipitation in Germany. *Clim. Dyn.* **2013**, *41*, 755–774. [CrossRef]

40. Pleim, J.E.; Xiu, A. Development and testing of a surface flux and planetary boundary layer model for application in mesoscale models. *J. Appl. Meteor.* **1995**, *34*, 16–32. [CrossRef]

41. Xiu, A.; Pleim, J.E. Development of a land surface model part I: Application in a mesoscale meteorology model. *J. Appl. Meteor.* **2001**, *40*, 192–209. [CrossRef]

42. Zhang, C.; Lin, H.; Chen, M.; Yang, L. Scale matching of multiscale digital elevation model (DEM) data and the Weather Research and Forecasting (WRF) model: A case study of meteorological simulation in Hong Kong. *Arab J. Geosci.* **2014**, *7*, 2215–2223. [CrossRef]

43. Lin, Y.; Colle, B.A. A new bulk microphysical scheme that includes riming intensity and temperature-dependent ice characteristics. *Mon. Weather Rev.* **2011**, *139*, 1013–1035. [CrossRef]
44. Fersch, B.; Kunstmann, H. Atmospheric and terrestrial water budgets: Sensitivity and performance of configurations and global driving data for long term continental scale WRF simulations. *Clim. Dyn.* **2014**, *42*, 2367–2396. [CrossRef]

45. Madala, S.; Satyanarayana, A.; Srinivas, C.; Bhishma, T. Performance evaluation of PBL schemes of ARW model in simulating thermo-dynamical structure of pre-monsoon convective episodes over Kharagpur using STORM data sets. *Pure Appl. Geophys.* **2016**, *173*, 1803–1827. [CrossRef]

46. Schwitalla, T.; Bauer, H.-S.; Wulfmeyer, V.; Aoshima, F. High resolution simulation over central Europe: Assimilation experiments during COPS IOP 9c. *Q. J. Roy. Meteor. Soc.* **2011**, *137*, 156–175. [CrossRef]

47. Politi, N.; Nastos, P.T.; Sfetsos, A.; Vlachogiannis, D.; Dalezios, N.R. Evaluation of the AWR-WRF model configuration at high resolution over the domain of Greece. *Atmos. Res.* **2018**, *208*, 229–245. [CrossRef]

48. Guerova, G.; Brockmann, E.; Quiby, J.; Schubiger, F.; Matzler, C. Validation of NWP mesoscale models with Swiss GPS Network AGNES. *J. App. Meteorol. 2003*, *42*, 141–150. [CrossRef]

49. Guerova, G.; Brockmann, E.; Schubiger, F.; Morland, J.; Matzler, C. An integrated assessment of measured and modeled IWV in Switzerland for the period 2001–2003. *J. App. Meteorol. 2005*, *44*, 1033–1044. [CrossRef]

50. Haase, J.; Ge, M.; Vedel, H.; Calais, E. Accuracy and variability of GPS tropospheric delay measurements of water vapor in the western Mediterranean. *J. App. Meteorol. 2003*, *42*, 1547–1568. [CrossRef]