Enhancing GNSS water vapour retrieval via synergistic microwave radiometry: thermodynamic error diagnosis and bias correction

The retrieval of Precipitable Water Vapour (PWV) from Global Navigation Satellite Systems (GNSS) in thermodynamically complex environments is significantly limited by the accuracy of the weighted mean temperature (T_m). This study evaluates the efficacy of static climatological models versus dynamic ground-based microwave radiometry for T_m determination in the Eastern Mediterranean, a region characterized by sharp refractivity gradients. Using the Cyprus GNSS Meteorology Enhancement research project (CYGMEN) infrastructure in Nicosia, the performance of the ERA5-based HGPT2 model and a co-located Microwave Radiometer (MWR) was assessed against radiosonde (RS) profiles during the 2025 warm season (Spring–Summer). Diagnostic analysis reveals that the static HGPT2 model fails to resolve the diurnal thermodynamic decoupling between the boundary layer and the free troposphere, leading to a systematic overestimation of T_m exceeding 6 K during peak solar insolation. Conversely, the MWR captures short-term thermodynamic variability (r=0.98) but exhibits a systematic cold bias of −1.91 K in raw retrievals. It is demonstrated that a site-specific linear bias correction reduces the MWR T_m Root Mean Square Error (RMSE) from 2.32–1.43 K, significantly outperforming the empirical model. Sensitivity analysis confirms that thermodynamic uncertainty dominates the error budget, outweighing uncertainties in refractivity constants by an order of magnitude. Consequently, standard climatological retrievals diverge from the synergistic MWR-GNSS method during extreme hygrometric events, introducing systematic PWV biases exceeding 1.0 mm when moisture levels surpass 45 mm. The synergistic coupling of real-time radiometric T_m with GNSS data is therefore meaningful for generating climate-quality PWV records in semi-arid coastal regions.

1 Introduction

Atmospheric water vapour (WV) is the primary greenhouse gas, contributing approximately 60 % to the natural greenhouse effect and playing a vital role in regulating the Earth's thermodynamic budget (Kiehl and Trenberth, 1997; Trenberth et al., 2005). Furthermore, WV is the main driver of latent heat transport, influencing convective systems and global precipitation patterns. High-frequency variations in Precipitable Water Vapour (PWV) correlate strongly with atmospheric instability and are a key factor in the initiation of severe weather. Specifically, rapid temporal gradients in PWV often precede heavy rainfall and flash floods, acting as a reliable precursor for convective storms (Brenot et al., 2013). Consequently, assimilating high-resolution PWV data into Numerical Weather Prediction (NWP) models significantly improves short-range precipitation “now-casting” (Bennitt and Jupp, 2012). Accurate PWV retrieval is especially crucial for the Eastern Mediterranean, a climate change “hotspot” warming faster than the global average (Giorgi, 2006; Lelieveld et al., 2012; Held and Soden, 2006). This region is characterized by complex topography and land–sea contrasts, which create sharp atmospheric refractivity gradients. The Eastern Mediterranean faces a hydro-climatic paradox: a long-term drying trend ($-\mathrm{0.5}\mathrm{mm}{\mathrm{decade}}^{-\mathrm{1}}$) alongside increasing high-intensity, short-duration flash flood events (Zittis et al., 2019; Ziv et al., 2021b). GNSS-derived PWV in this region exhibits strong diurnal cycles with amplitudes up to 5 mm, which are closely correlated with atmospheric instability (Ziskin Ziv et al., 2021a). Despite this vulnerability, the Eastern Mediterranean currently lacks dense, continuous atmospheric profiling networks. Traditional observation methods, such as radiosondes (RS), fail to resolve these mesoscale events due to low temporal resolution (typically 12 h intervals) and significant spatial gaps (Soden and Lanzante, 1996). While satellite-based passive remote sensing offers global coverage, it is often limited by revisit times, daylight dependence, or data degradation in coastal zones due to land contamination in the microwave footprint (Bennartz and Bauer, 2003).

These limitations underscore the necessity for ground-based remote sensing techniques that offer continuous, all-weather operability. Ground-based Global Navigation Satellite Systems (GNSS) meteorology has emerged as a robust technique for atmospheric sounding since the seminal proposal by Bevis et al. (1992). By estimating the Zenith Total Delay (ZTD) of GNSS signals traversing the neutral atmosphere, the Zenith Wet Delay (ZWD) can be isolated by subtracting the Zenith Hydrostatic Delay (ZHD), which is accurately modeled from surface pressure observations (Saastamoinen, 1972). GNSS-derived PWV offers significant advantages, including high temporal resolution (sub-hourly), all-weather availability, and cost-efficiency by leveraging existing geodetic infrastructure (Guerova et al., 2016; Jones et al., 2020).

The retrieval of PWV from GNSS ZWD relies on a dimensionless conversion factor, Π, which is a function of the atmospheric weighted mean temperature, T_m. Defined physically as $\int (e/T)\mathrm{d}z/\int (e/T^{\mathrm{2}})\mathrm{d}z$ (Askne and Nordius, 1987), T_m encapsulates the thermal state of the atmospheric column. Because the conversion factor (Π) is nearly linearly proportional to T_m, any relative error in the T_m estimation strictly propagates as an equivalent relative error in the final PWV retrieval. During severe moisture events with an PWV of 50 mm, this translates to an absolute error of ∼0.18 mm. Consequently, a 1 % relative error in T_m translates strictly to a 1 % relative error in PWV. Therefore, alongside the substantial errors inherent in ZTD estimation – such as mapping function inaccuracies and surface pressure interpolation for the ZHD (Ning et al., 2016) – the determination of T_m remains one of the primary sources of uncertainty in GNSS meteorology. In the absence of in-situ profiles, T_m is commonly estimated using empirical regression models or global climatological models. However, earlier studies have demonstrated that empirical T_m regressions, such as the Bevis model (Bevis et al., 1992), introduce relative PWV errors of 1 %–2 % due to weak T_s−T_m correlations in coastal and equatorial regions, where annual/semiannual variations are not adequately captured (Yao et al., 2014; Lan et al., 2016). Similarly, global grid-based T_m models like GPT2w achieve RMSE<4 K at ∼80 % of mid-latitude sites but degrade in data-sparse areas like the Middle East and Africa, where reanalysis quality is limited (Böhm et al., 2015; Jiang et al., 2019). The Hourly Global Pressure and Temperature 2 (HGPT2) model represents a major advancement by providing hourly estimates derived from ERA5 reanalysis (Mateus et al., 2021). Despite recent validation of GPT2w and ECMWF models for Precipitable Water Vapour (PWV) retrieval in the Mediterranean (Oikonomou et al., 2018), a critical gap exists: the quantification of vertical interpolation errors in these models, especially over complex coastal terrains. For instance, while recent validation studies in Cyprus demonstrate strong GNSS-PWV correlations (>0.6) with ERA5 during extreme precipitation, persistent reanalysis interpolation errors are highlighted in mountainous areas (Giannadaki et al., 2025). This lack of validation for HGPT2's performance in the complex topography of the Eastern Mediterranean potentially exacerbates PWV biases during extreme events.

An alternative approach to deriving T_m is the use of ground-based Microwave Radiometers (MWR). MWRs measure brightness temperatures at multiple frequencies to retrieve continuous vertical profiles of temperature and humidity. Ground-based MWRs have been shown to retrieve T_m with $\text{RMSE}\sim \mathrm{1}-\mathrm{2}\mathrm{K}$ in mid-latitudes, outperforming empirical models during synoptic anomalies (Cimini et al., 2011; Crewell and Löhnert, 2007; Löhnert and Maier, 2012). While multi-site intercomparisons reveal that MWR retrievals can exhibit upper-tropospheric cold biases (up to 5 K at >2 km altitude) (Van Malderen et al., 2014; Steinke et al., 2015), simple linear corrections can reduce RMSE by 20 %–40 % (Ning and Elgered, 2021). Operational GNSS–MWR synergies for T_m bias correction have documented gains in PWV accuracy (10 %–30 % RMSE reduction) in European networks (Vaquero-Martínez et al., 2018; Li et al., 2020). However, such applications are rare in the semi-arid Eastern Mediterranean, where MWR could critically mitigate reanalysis uncertainties.

This study leverages the infrastructure of the CYGMEN (Cyprus GNSS Meteorology Enhancement) project, which is establishing a dense, multi-sensor meteorological network in Cyprus. The network, termed CyMETEO, integrates a dense array of continuous GNSS stations distributed across the island. Due to the high cost and operational complexity of radiometric profiling, the network currently features a single, centralized thermodynamic “supersite” at the Athalassa observatory, where a GNSS receiver is strictly co-located with a MWR and a RS launch facility. This unique instrumental setup provides an ideal testbed for inter-comparing atmospheric retrieval techniques in a coastal, semi-arid environment. The primary objective of this manuscript is to evaluate the accuracy of GNSS-derived PWV over the Eastern Mediterranean by assessing the performance of different T_m determination strategies. We specifically investigate the efficacy of the HGPT2 model compared to MWR-derived T_m and RS benchmarks. The study aims to quantify the error budget of GNSS-PWV and determine whether the inclusion of MWR data provides statistically significant improvements over the state-of-the-art HGPT2 model. The manuscript is organized as follows: Section 2 describes the study area and the instrumentation of the CyMETEO network; Sect. 3 details the methodology for GNSS processing, ZTD estimation, and the mathematical derivation of T_m from different sources; Sect. 4 presents the validation results and statistical analysis against RS reference data; and Sect. 5 concludes with recommendations for operational PWV monitoring strategies in the region.

2 Data and Methodology

2.1 Observational Site and CYGMEN Infrastructure

The observational campaign was conducted at the Athalassa atmospheric observatory in Nicosia, Cyprus (35.15° N, 33.40° E, 160 m a.s.l. (meter above sea level)), situated in the central Mesaoria plain. The site is characterized by complex topography, bounded by the Troodos Mountain to the southwest and the Pentadaktylos Mountain to the north, as shown in Fig. 1a. This study presents the first comprehensive analysis of radiometric data acquired under the CYGMEN infrastructure project, established to monitor the thermodynamic state of the Eastern Mediterranean atmosphere. To ensure robust thermodynamic profiling and validation, three primary datasets were collated, as shown in Table 1.

Figure 1 Location and instrumentation at Athalassa, Cyprus. (a) Site location on the island's elevation map. (b) GNSS reference station. (c) RPG-HATPRO radiometer. (d) Radiosonde balloon launching.

2.2 Instrumentation and Data Processing

2.2.1 Microwave Radiometry (MWR)

The RPG-HATPRO radiometer observes downwelling atmospheric brightness temperatures (T_B) across 14 channels: seven in the K-band (22–31 GHz) sensitive to water vapour, and seven in the V-band (51–58 GHz) sensitive to oxygen for temperature profiling. This instrument enables the continuous retrieval of temperature (T) and absolute humidity (ρ_v) profiles on a standardized grid of 94 vertical levels from the surface up to 10 km. The vertical resolution is optimized for the planetary boundary layer (PBL), ranging from 10–30 m up to 500 m, and decreasing to 100–500 m in the free troposphere. For this study, high-frequency MWR observations were resampled to 15 min intervals to align with GNSS processing epochs. It is a well-documented limitation of passive microwave radiometry that retrieval accuracy degrades significantly during precipitation, as liquid water on the instrument's radome heavily contaminates the measured brightness temperatures (Foth et al., 2024; Parde et al., 2025; Pakkattil et al., 2025; Ware et al., 2004). Because this study focused on the warm, dry season in the Eastern Mediterranean (March–October 2025), rainfall events were naturally sparse. Nevertheless, to ensure data integrity, real-time precipitation flags generated by the co-located Vaisala WXT536 surface weather transmitter were utilized as a strict quality-control filter. Any MWR profiles retrieved during active precipitation events were excluded from the dataset to prevent wet-radome anomalies from artificially skewing the thermodynamic bias analysis. To diagnose potential biases in the MWR-derived T_m, the dataset was split into a training Set (April–June 2025) for regression modeling and a validation Set (July–October 2025) for independent testing. In addition to thermodynamic profiling, the MWR's native retrieval algorithm possesses the capacity to directly estimate PWV from its K-band brightness temperatures.

Any MWR profiles retrieved during active precipitation events were excluded from the dataset to prevent wet-radome anomalies from artificially skewing the thermodynamic bias analysis. To mitigate T_m errors in MWR, a supervised linear regression model was developed to calibrate the MWR observations. For robust independent validation, the collocated dataset was separated into two distinct temporal subsets: the training Set (April–June 2025), which was used to derive the regression coefficients, and the validation Set (July–October 2025), which was used exclusively to test the correction's performance on unseen data. A simple linear correction model was fitted to the training data using Ordinary Least Squares (OLS) minimization. The relationship is defined in Eq. (1) as:

$$\begin{array}{}\text{(1)}& T_{\text{m,corr}}=\mathit{\alpha }\cdot T_{\text{m,MWR}}+\mathit{\beta }\end{array}$$

where T_m,corr is the corrected MWR temperature, T_m,MWR is the raw T_m derived from the radiometer and α (slope) and β (intercept) are the learned coefficients minimizing the residual sum of squares between the MWR and RS values. Based on our training Set, the derived coefficients applied to the validation Set were α=1.0623 and $\mathit{\beta }=-\mathrm{15.6062}\mathrm{K}$.

2.2.2 Radiosonde Data Processing

To establish a rigorous validation dataset, PWV was derived from high-resolution vertical profiles obtained from collocated radiosonde launches. A strict collocation window was applied, where MWR profiles were averaged within ±30 minutes of the balloon launch time. The raw telemetry data, comprising pressure (P), temperature, and dew point temperature (T_d), were processed to derive the total columnar water vapour content (in kg m⁻²) through the vertical integration of specific humidity, assuming the atmosphere is in hydrostatic equilibrium. The determination of the necessary moisture variables relied on the Magnus-Tetens approximation, which provides a widely accepted empirical relationship for saturation vapour pressure. First, the actual vapour pressure (e, in hPa) was computed directly from the dew point temperature (T_d, in °C). This calculation utilized the coefficients defined by Bolton (1980), which are optimized for saturation vapour pressure over liquid water in the meteorological temperature range, as shown in Eq. (2):

$$\begin{array}{}\text{(2)}& e=\mathrm{6.112}\cdot \exp \left({\displaystyle \frac{\mathrm{17.67}\cdot T_{\mathrm{d}}}{T_{\mathrm{d}}+\mathrm{243.5}}}\right)\end{array}$$

Subsequently, the specific humidity (q, in kg kg⁻¹) was derived via Eq. (3), representing the mass mixing ratio of water vapour to the total moist air parcel:

$$\begin{array}{}\text{(3)}& q={\displaystyle \frac{\mathit{ϵ}\cdot e}{P-(\mathrm{1}-\mathit{ϵ})\cdot e}}\end{array}$$

where P is the static pressure (hPa) and ϵ≈0.622 represents the ratio of the molecular weight of water vapour to that of dry air. Once the specific humidity profile was established, the PWV was calculated by integrating q with respect to pressure. The retrieval algorithm employed the trapezoidal rule for numerical integration, which approximates the integral as the sum of discrete atmospheric layers (Eq. 4):

$$\begin{array}{}\text{(4)}& \text{IWV}={\displaystyle \frac{\mathrm{1}}{g}}\sum _{i=\mathrm{0}}^{N-\mathrm{1}}{\displaystyle \frac{q_i+q_{i+\mathrm{1}}}{\mathrm{2}}}\cdot |P_{i+\mathrm{1}}-P_i|\end{array}$$

where g is the gravity dependent on altitude, g(ϕ,h), where ϕ represents the Geodetic latitude and h is the orthometric height. P is converted to Pascals (Pa) prior to integration and N represents the total number of vertical levels in the RS profile. It should be noted that while the trapezoidal rule can theoretically overestimate the integral of an exponentially decaying profile, the Vaisala RS41-SGP provides high-frequency 1 s telemetry (yielding a vertical spatial resolution of approximately 5–8 m). At this exceptionally fine resolution, the linear approximation between measurement levels effectively converges with the true atmospheric profile, rendering any systematic integration bias mathematically negligible. Also, it should be noted that IWV, representing the mass column integral in kg m⁻², is physically and numerically equivalent to PWV expressed as a depth in millimeters (mm), assuming the standard density of liquid water (1000 kg m⁻³). A strict quality assurance protocol was implemented to ensure vertical completeness; only radiosonde flights that successfully maintained continuous telemetry up to the 10 km AGL integration ceiling were included in the final comparative dataset. While the term PWV is frequently used when discussing direct profile integration from the MWR and RS, this study uses PWV (mm) as the standardized final retrieval metric to align with operational meteorological and forecasting conventions. It is important to note that while the nominal manufacturer uncertainty for the Vaisala RS41 humidity sensor is stated as 4 % for individual profile measurements, the uncertainty of the resulting PWV is significantly lower. Because PWV is computed by integrating hundreds of discrete measurements across the vertical column (Eq. 3), uncorrelated random sensor noise is largely suppressed through statistical cancellation. Consequently, the integrated variables derived from the radiosonde, such as PWV and the T_m, possess the requisite precision to serve as a robust “ground truth” standard for evaluating the finer relative uncertainties (1 %–2 %) associated with the GNSS and MWR retrievals. To ensure a rigorous and direct intercomparison with the active MWR, the radiosonde integration was strictly confined to a maximum altitude of 10 km Above Ground Level (AGL). This vertical cutoff was deliberately chosen to exactly match the 10 km ceiling of the standard RPG-HATPRO retrieval grid. While GNSS integrates delays through the entire atmosphere, bounding the in-situ reference data is mathematically necessary to isolate profiling performance. It is well established that this 10 km threshold does not introduce a systematic dry bias when comparing against total-column GNSS (Van Baelen et al., 2005). Furthermore, ambient temperatures at this altitude range from −40 to −50 °C, strictly limiting the saturation vapour pressure. Consequently, the residual water vapour between 10 km and the tropopause is thermodynamically constrained to fractions of a millimeter. Omitting this minute residual mass is functionally negligible, as it falls well within the overall baseline uncertainty (typically 1–2 mm) of the total-column radiosonde PWV retrieval.

Table 1 Summary of Instrumentation and Datasets.

2.2.3 GNSS Data Processing

ZTD estimates were derived from the collocated Leica GR50 receiver (station NICO) using the Tefnut PP software (Douša and Václavovic, 2014). The processing employed a Precise Point Positioning (PPP) strategy with an elevation cutoff angle of 10°. To account for tropospheric mapping errors, the Vienna Mapping Function 1 (VMF1) was applied. Station coordinates were constrained to the IGS14 reference frame, and satellite orbits and clock corrections were utilized from IGS Ultra-Rapid products. While IGS Final products are the gold standard for historical climate reprocessing due to their minimal orbital uncertainty, this study deliberately utilized IGS Ultra-Rapid products to evaluate the proposed synergistic retrieval architecture under near real-time operational constraints. Because a primary application of continuous GNSS-PWV is its assimilation into short-range NWP for severe weather “nowcasting”, it is crucial to assess system performance using the satellite orbits and clocks actually available during active forecasting. Although Ultra-Rapid products introduce a slight degradation in ZTD precision compared to final products, this uncertainty (typically fractions of a millimeter in PWV) remains negligible compared to the massive, multi-millimeter systematic errors introduced by static thermodynamic modeling, which is the primary focus of this investigation. While modern Numerical Weather Prediction systems frequently assimilate ZTD directly to avoid conversion uncertainties, deriving an accurate physical PWV product remains essential. PWV serves as an intuitive, absolute moisture metric heavily utilized by operational forecasters for severe weather nowcasting, and is fundamentally necessary for building long-term, cross-instrument climatological records. To isolate the ZWD, the ZHD was precisely calculated using continuous, co-located surface pressure observations obtained directly from the Vaisala WXT536 weather transmitter installed at the site, rather than relying on interpolated pressure fields. ZTD values were estimated at 15 min intervals, directly aligning with the temporal resolution of the MWR. It must also be noted that the computation of ZHD is significantly dependent on the assumed value of the dry refractivity constant, k₁. As established by Bevis et al. (1994) and further evaluated by Healy (2011), while k₁ is known to a high degree of relative accuracy, its residual fractional uncertainty introduces a persistent systematic bias into the ZHD estimation. Because ZWD is isolated by subtracting ZHD from the total delay, this k₁-induced bias directly propagates into the final PWV error budget, acting alongside the conversion uncertainties analyzed later in this study.

2.3 Thermodynamic Modeling and Synergistic Retrieval Strategy

The conversion of GNSS-derived ZWD to PWV is governed by a proportionality factor, Π, whose accuracy is largely dictated by the T_m. To assess the fidelity of thermodynamic inputs for GNSS meteorology, we evaluated three distinct T_m derivation strategies. For profile-resolving instruments (MWR and RS), T_m values were computed by integrating the vertical profiles of physical temperature, T(z) (K), and absolute humidity, ρ_v(z) (kg m⁻³). Consistent with Bevis et al. (1992), T_m is defined as the mean temperature of the atmosphere weighted by the water vapour partial pressure, which can be expressed in terms of vapour density as shown in Eq. (5):

$$\begin{array}{}\text{(5)}& T_{\mathrm{m}}={\displaystyle \frac{{\int }_{z_{\text{surf}}}^{z_{\text{top}}}{\mathit{\rho }}_{\mathrm{v}}\left(z\right)\mathrm{d}z}{{\int }_{z_{\text{surf}}}^{z_{\text{top}}}\frac{{\mathit{\rho }}_{\mathrm{v}}\left(z\right)}{T\left(z\right)}\mathrm{d}z}}\end{array}$$

In practice, the continuous integrals were discretized using the trapezoidal rule from the surface (z_surf) to the highest available profile level (z_top). This approach assumes linear variation of T and ρ_v between measurement levels. For standalone GNSS retrieval (where no dynamic profiles are available), T_m was derived from the HGPT2 (Hourly Global Pressure and Temperature 2) model (Mateus et al., 2021). HGPT2 is an advanced “blind” empirical model, meaning its outputs are independent of the specific observational year. While dynamic NWP models provide superior real-time meteorological data, “blind” models like HGPT2 remain heavily utilized in standard geodetic GNSS processing where real-time meteorological or NWP data streams are unavailable. It is constructed from a comprehensive 20 year historical baseline of atmospheric data from the ERA5 global reanalysis. Unlike standard static climatologies, HGPT2 leverages the full ERA5 spatial resolution (0.25°×0.25°) and provides temporal resolution at 1 h intervals for any given Day of Year (DOY). It achieves this by employing a time-segmentation concept, modeling thermodynamic variables via long-term mean values combined with annual, semi-annual, and quarterly periodic functions.

Applying the linear correction model (as formulated in Sect. 2.2.1) successfully re-centers the error distribution. To quantify the benefits of sensor synergy in integrated water vapour estimation, this study defines and contrasts two distinct GNSS PWV retrieval architectures. The first, “Standard Retrieval” which is a control method utilizes the ZTD_GNSS combined with the T_m derived empirically from the HGPT2 climatological model (Böhm et al., 2015). Second “Synergistic Retrieval” which proposed method couples ZTD_GNSS with a physical T_m derived directly from a collocated MWR. For the synergistic approach, the dimensionless conversion factor (Π) was calculated dynamically using the MWR-derived T_m following Eq. (6) and Eq. (7).

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\text{PWV}=\mathrm{\Pi }\cdot \text{ZWD}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{\Pi }={\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{{\mathit{\rho }}_{\mathrm{w}}{R}_{\mathrm{v}}\left[k_{\mathrm{2}}^{\prime }+(k_{\mathrm{3}}/T_{\mathrm{m}})\right]}}\end{array}$$

where ρ_w represents the density of liquid water (1000 kg m⁻³) and R_v is the specific gas constant for water vapour (461.52 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). To assess the sensitivity of the Π to the choice of thermodynamic coefficients, three widely used formulations were employed in this study, following Davis et al. (1985)/Thayer (1974), Bevis et al. (1994), and Rüeger (2002), as shown in Table 2.

Table 2 Refractivity constants used in the sensitivity analysis of the Π factor, based on three commonly adopted formulations.

To rigorously quantify the uncertainty in the final PWV retrieval and avoid fragmented error attributions, standard error propagation must be applied to the fundamental conversion Eq. (6). Assuming the uncertainties in the wet delay and the conversion factor are uncorrelated, the variance of the final PWV (${\mathit{\sigma }}_{\text{PWV}}^{\mathrm{2}}$) is expressed using partial derivatives as in Eq. (8):

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{\text{PWV}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial \text{PWV}}{\partial \text{ZWD}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\text{ZWD}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \text{PWV}}{\partial \mathrm{\Pi }}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Pi }}^{\mathrm{2}}\end{array}$$

Evaluating these primary partial derivatives yields the proportional contributions of the geodetic and thermodynamic components as shown in Eq. (9):

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{\text{PWV}}^{\mathrm{2}}={\mathrm{\Pi }}^{\mathrm{2}}{\mathit{\sigma }}_{\text{ZWD}}^{\mathrm{2}}+{\text{ZWD}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{\Pi }}^{\mathrm{2}}\end{array}$$

The uncertainty in the conversion factor (${\mathit{\sigma }}_{\mathrm{\Pi }}^{\mathrm{2}}$) is itself a compound term driven by the T_m and the static atmospheric refractivity constants (k₂′ and k₃). Its variance is defined via partial derivatives as shown in Eq. (10):

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{\Pi }}^{\mathrm{2}}={\left({\displaystyle \frac{\partial \mathrm{\Pi }}{\partial T_{\mathrm{m}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{T_{\mathrm{m}}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{\Pi }}{\partial k_{\mathrm{2}}^{\prime }}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{k_{\mathrm{2}}^{\prime }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathrm{\Pi }}{\partial k_{\mathrm{3}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{k_{\mathrm{3}}}^{\mathrm{2}}\end{array}$$

The sensitivity of the conversion factor strictly to T_m (the dynamic thermodynamic variable evaluated in this study) is quantified by its partial derivative and it represented as Eq. (11):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathrm{\Pi }}{\partial T_{\mathrm{m}}}}=\mathrm{\Pi }\left[{\displaystyle \frac{k_{\mathrm{3}}}{T_{\mathrm{m}}^{\mathrm{2}}\left(k_{\mathrm{2}}^{\prime }+\frac{k_{\mathrm{3}}}{T_{\mathrm{m}}}\right)}}\right]\end{array}$$

This consolidated formulation establishes the exact mathematical limits of thermodynamic error propagation. As demonstrated in the sensitivity analysis (Sect. 3.4), this framework accurately isolates the dynamic uncertainties driven by T_m from the baseline static biases introduced by the chosen refractivity constants.

2.4 Diagnostic Parameters and Error Analysis

The vertical structure of the atmosphere was analyzed by segregating the dataset into two regimes: the PBL (0–2 km), where water vapour is concentrated, and the Free Troposphere (>2 km). Additionally, the water vapour scale height (H_v) was calculated to parameterize the vertical distribution of moisture. H_v was derived for both RS and MWR by fitting an exponential decay function (Eq. 12) to the absolute humidity profile (ρ_v).

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_{\mathrm{v}}\left(z\right)={\mathit{\rho }}_{\mathrm{v},\mathrm{0}}\cdot \exp \left(-{\displaystyle \frac{z}{H_{\mathrm{v}}}}\right)\end{array}$$

where ρ_v(z) is the absolute humidity at height z, and ρ_v,0 is the surface humidity. This curve fitting was deliberately restricted to the lowest 4 km of the atmosphere. Because this layer contains the vast majority (>90 %) of the tropospheric water vapour mass, bounding the fit prevents the algorithm from heavily weighting near-zero, noisy upper-tropospheric values that mathematically degrade the fit for the boundary layer. Furthermore, the scale height metric fundamentally assumes the atmosphere conforms to a well-behaved exponential decay. Profiles yielding H_v values outside the physically realistic range of 0.1–4.0 km were excluded from the statistical analysis to prevent artificial statistical skewing during complex meteorological states (e.g., deep convective mixing) where the underlying exponential model is invalid. Forcing a mathematical fit onto these non-exponential profiles yields physically meaningless artifacts. Therefore, a Quality Assurance filter was applied, bounding the analysis to the physically realistic range of $\mathrm{0.1}\mathrm{km}<H_{\mathrm{v}}<\mathrm{4.0}\mathrm{km}$. Profiles yielding values outside this range were discarded because they indicate the underlying exponential model itself is invalid for that specific atmospheric profile, preventing artificial statistical skewing in the instrument intercomparison. Profiles yielding H_v values outside the physically realistic range of 0.1–4.0 km were excluded from the statistical analysis. To evaluate the performance limitations of standard climatological models under varying hygrometric conditions, the systematic error (ΔPWV) was defined as the residual between the synergistic and standard approaches (Eq. 13):

$$\begin{array}{}\text{(13)}& \mathrm{\Delta }\text{PWV}={\text{PWV}}_{\text{Synergistic}}-{\text{PWV}}_{\text{Standard}}\end{array}$$

The dataset was stratified into discrete bins of 5 mm PWV to isolate regimes of moisture abundance. Within each bin, the mean bias and ±1σ uncertainty were computed. These statistics were utilized to determine the “Systematic Bias Threshold,” defined herein as the specific hygrometric threshold where the systematic model error exceeds 1 mm. Finally, the propagation of thermodynamic uncertainty into the moisture retrieval was quantified via linear regression analysis. This compared the relative error in T_m (HGPT2 vs. MWR) against the resulting relative error in PWV, serving as an empirical verification of the theoretical sensitivity approximation given in Eq. (14):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{\Delta }\text{PWV}}{\text{PWV}}}\approx {\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{m}}}{T_{\mathrm{m}}}}\end{array}$$

3 Results

The following evaluation follows a top-down diagnostic approach. First, the macroscopic baseline performance of the final derived moisture products is established. Subsequently, the underlying thermodynamic variables driving these discrepancies are isolated, culminating in the development of a targeted calibration scheme to mitigate the identified biases.

3.1 Temperature and Humidity Profile Validation

MWR-retrieved temperature T and ρ_v profiles were validated against collocated RS observations at 00:00 UTC and 12:00 UTC during March–October 2025. Profiles were stratified into the planetary boundary layer (PBL; 0–2 km) and free troposphere (>2 km), as shown in Figs. 2 and 3. Mean vertical temperature profiles show agreement between MWR and RS (Fig. 2a and b). In the boundary layer (0–2 km), MWR retrieves temperature with high precision (r>0.98, RMSE<1.5 K). Above 2 km, a cold bias is observed in the MWR retrieval, reaching −5.16 K at 12 UTC (Fig. 2f). Despite this bias, the linearity remains strong (r≈0.97), indicating the sensor captures relative thermal variations aloft despite the absolute offset. This confirms the trend observed in the mean profiles, where the MWR underestimates temperatures in the mid-to-upper troposphere. Furthermore, horizontal balloon drift driven by prevailing winds inevitably causes the radiosonde to sample a different atmospheric volume than the MWR's strict zenith view. While this spatiotemporal mismatch introduces random scatter into the upper-level comparisons, it does not artificially skew the systematic biases identified in this study. Consequently, the RMSE increases substantially to approximately 6.4–6.7 °C. The stark contrast in accuracy between the lower and upper troposphere is a known characteristic of ground-based microwave radiometry (Parde et al., 2025; Pakkattil et al., 2025). The high accuracy below 2 km is attributed to the high information content of the opaque V-band channels (51–58 GHz), whose weighting functions peak near the surface. Above 2 km, these weighting functions broaden significantly, reducing vertical resolution and causing a “smearing” effect where the instrument provides a volume-averaged temperature rather than a precise point measurement. The observed cold bias is likely a result of the retrieval algorithm (e.g., neural network) relying heavily on a climatological a priori dataset that does not perfectly represent the thermal conditions of the transition season observed, or systematic offsets in the radiative transfer model (absorption coefficients) used for training.

Figure 2 Comparison of radiosonde and microwave radiometer (MWR) temperature profiles: (a, b) Mean vertical temperature (T) profiles at 00 and 12 UTC with variability shading; (c–f) Scatter comparisons for the lower (0–2 km) and upper (>2 km) atmosphere at both times.

Figure 3 Comparison of radiosonde and microwave radiometer (MWR) absolute humidity (ρ_v) profiles: (a, b) Mean vertical ρ_v profiles at 00 and 12 UTC with variability shading; (c–f) Scatter comparisons for the lower (0–2 km) and upper (>2 km) atmosphere at both times.

The mean ρ_v profiles (Fig. 3a and b) show the expected exponential decrease of moisture with height. At 00 UTC, the profiles align reasonably well. However, at 12 UTC, the MWR profile exhibits a structural deviation between 1–2 km, failing to capture the smooth moisture gradient recorded by the RS. This discrepancy may be attributed to the MWR's limited vertical resolution during periods of active daytime mixing or complex humidity layering. The retrieval of humidity in the lower atmosphere shows moderate agreement but is less accurate than the temperature retrievals. Performance is notably better at night (00 UTC) with r of 0.878 and RMSE of 1.98 g m⁻³. At 12 UTC, the correlation drops to 0.744, and the scatter increases ($\text{RMSE}=\mathrm{2.31}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$). Because the integrated mass of the water vapour column is physically and numerically equivalent to its depth (assuming the standard density of liquid water), this physical quantity is exclusively referred to as PWV expressed in millimeters (mm) throughout this study to align with operational meteorological conventions. A negative bias persists at both times ($-\mathrm{0.51}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$ at 00 UTC and $-\mathrm{0.91}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$ at 12 UTC), indicating a tendency for the MWR to underestimate moisture content in the boundary layer, particularly during the day. Surprisingly, the statistical linearity for ρ_v improves slightly or remains stable above 2 km, likely due to the lower overall magnitude of humidity at these heights. The correlation coefficients remain stable (∼0.87). In contrast to the lower levels, the bias shifts to slightly positive values (0.23 g m⁻³ at 00 UTC and 0.46 g m⁻³ at 12 UTC), suggesting a slight moist bias in the MWR retrievals aloft. The linear fits (Fig. 3e and f) align closely with the 1:1 line, with slopes near unity (0.90 and 1.00), indicating that the MWR effectively captures the free tropospheric humidity trends despite the lower absolute values. The difficulty in retrieving accurate ρ_v profiles, particularly at 12 UTC, stems from the limited vertical resolution of the K-band channels (22–31 GHz). Unlike temperature profiling, humidity profiling offers very few independent degrees of freedom (typically<3), making it difficult for the MWR to resolve sharp vertical gradients often present at the top of the convective boundary layer during the daytime. The structural deviation and underestimation are common issues linked to the “smoothing” error inherent in passive radiometry, where sharp moisture inversions are averaged out. Furthermore, the persistent bias suggests potential uncertainties in the water vapour absorption models (spectroscopic parameters) or non-representative training data used in the retrieval algorithm.

3.2 Precipitable Water Vapour (PWV) and Scale Height (H_v) Validation

Unlike vertical profiling, the MWR excels in measuring total column quantities. The comparison with RS yields an excellent correlation (r=0.971) and a low RMSE of 1.72 kg m⁻² (Fig. 4). This performance disparity – superior PWV versus degraded profiles – confirms that while the sensor cannot resolve vertical structural details due to smoothing error, the radiometric brightness temperature in the K-band remains strictly proportional to the total precipitable water mass. The GNSS-derived PWV shows a slight negative bias relative to RS ($-\mathrm{0.45}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$), whereas it relative to the MWR exhibits a positive bias ($+\mathrm{0.45}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$). The cumulative offset observed in the MWR-GNSS intercomparison ($+\mathrm{0.89}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$) highlights the systematic differences in calibration and retrieval assumptions between active (GNSS) and passive (MWR) techniques. The GNSS underestimation is likely driven by errors in the T_m derived from the static HGPT2 model, a hypothesis further explored in Sect. 3.4. To further diagnose the structural limitations of the retrievals, we evaluated the water vapour H_v. While H_v is admittedly a single-parameter representation of the complex atmospheric moisture profile, it is a crucial parameter that provides a representative value for the rate at which water vapour decreases with altitude – a key factor in understanding atmospheric stability, cloud formation, and radiative transfer processes. In this study, it is utilized specifically as a diagnostic metric to quantify the vertical structural limitations of passive microwave remote sensing. The comparison of H_v calculated from RS and MWR profiles is shown in Fig. 5. Unlike the high-fidelity PWV retrievals, the MWR-derived scale height shows negligible correlation with RS observations (r=0.25, $R^{\mathrm{2}}=-\mathrm{2.87}$) and a massive systematic positive bias of 0.62 km. The histograms (Fig. 5b) further elucidate this discrepancy: while the RS scale heights follow a narrow, physically realistic distribution centered around a mean (μ) of 1.51 km, the MWR distribution is artificially broad and shifted to significantly higher values (μ=2.13 km).

Figure 4 Intercomparison of precipitable water vapour (PWV) retrieved from Microwave Radiometer (MWR), GNSS, and Radiosonde observations. (a) MWR PWV versus radiosonde PWV, (b) GNSS PWV (derived using HGPT2 T_m) versus radiosonde PWV, and (c) MWR PWV versus GNSS PWV (derived using HGPT2 T_m).

The large scatter and ambiguity in the MWR estimates – which completely dwarf the individual least-squares fit uncertainties of the exponential regression – are a direct consequence of the instrument's physical limitations. H_v is highly sensitive to the sharp vertical gradient of humidity at the top of the planetary boundary layer. However, the K-band channels (22–31 GHz) utilized for humidity profiling possess broad weighting functions, restricting the vertical degrees of freedom to typically fewer than three. Because the MWR lacks the vertical resolution to capture sharp moisture inversions, the retrieval algorithm mathematically smears the moisture mass upward. This inherent “smoothing error” artificially elongates the vertical moisture profile, effectively inflating the calculated e-folding depth. Therefore, the inclusion of this H_v analysis serves to transparently demonstrate a critical operational boundary: while the MWR is an excellent standard for PWV, it is significantly unreliable and mathematically unsuited for characterizing vertical moisture compactness. The significant deviations observed in these macroscopic retrieval products necessitate a deeper investigation into the intermediate thermodynamic variables driving the conversion process. Consequently, the isolated performance of the T_m is evaluated in Sect. 3.3, followed by the introduction of a post-retrieval MWR calibration scheme in Sect. 3.4 designed to mitigate these native biases.

Figure 5 Comparison of scale height from radiosonde (RS) and microwave radiometer (MWR): (a) Scatter plot with 1:1 line and linear fit, including summary statistics; (b) Frequency distributions showing mean scale heights for RS and MWR.

3.3 Weighted Mean Temperature (T_m) Validation

The accurate estimation of the T_m is critical for converting GNSS-derived ZWD into PWV. The performance of T_m derived from the MWR and the empirical GPT2w model (HGPT2) was evaluated against RS measurements, which serve as the “ground truth.” The results are presented in Fig. 6. The time series (Fig. 6a) illustrates the seasonal evolution of T_m from April–October 2025. The Radiosonde observations (black dots) show significant variability, capturing synoptic-scale weather fluctuations. The MWR-derived T_m (orange dots) tracks these fluctuations with remarkable precision, overlaying the RS points almost perfectly. In stark contrast, the HGPT2 model (blue dots) provides a smooth, climatological curve. While it captures the general seasonal trend, it completely misses the day-to-day thermodynamic variability, often overestimating T_m during cooler transient events and underestimating it during warmer anomalies. The empirical model shows only moderate performance (r=0.800) with a substantial spread (RMSE=4.54 K). A systematic positive bias of 1.67 K indicates that HGPT2 generally overestimates the atmospheric temperature profile in this region. The scatter plot reveals a diffuse, “cloud-like” distribution, confirming its inability to capture real-time atmospheric dynamics. The MWR demonstrates superior performance, achieving a near-perfect correlation (r=0.981). The RMSE is significantly reduced to 2.33 K, which is nearly half the error of the empirical model. Interestingly, the MWR exhibits a negative bias of −1.91 K, suggesting a systematic underestimation of T_m. Crucially, this bias does not originate in the free troposphere, but rather in the planetary boundary layer (0–3 km). Since T_m is weighted by water vapour pressure, this “cold bias” indicates the MWR is underestimating the intense near-surface heating or the sharp lapse rates characteristic of the Nicosia environment. Despite this offset, the tight linearity indicates that MWR is an excellent source for capturing real-time T_m variations. Comparing the large dataset of MWR against HGPT2 (N=12 493) confirms the discrepancy between dynamic and static modeling. The correlation is lower (r=0.770) and the scatter is large (RMSE=4.89 K), further proving that static empirical models are insufficient for high-precision GNSS meteorology compared to dynamic radiometer measurements. While errors in ZTD estimation contribute significantly to the overall uncertainty budget, the specific error introduced during the conversion from delay to water vapour is linearly dependent on the accuracy of T_m. Assuming a given ZTD, a standard rule of thumb states that a 1 % relative error in T_m translates to roughly a 1 % relative error in the resulting PWV. By switching from a static model (HGPT2, ∼4.5 K error) to a dynamic sensor (MWR, ∼2.3 K error), the uncertainty in the GNSS water vapour product is effectively halved. This validates the “synergistic” approach of using collocated MWR thermal data to process GNSS signals.

Figure 6 Comparison of weighted mean temperature (T_m) derived from HGPT2, MWR, and Radiosonde (RS) during March–November 2025. (a) Time series of T_m estimates from all three sources. (b–d) Scatter plots showing statistical comparisons between: (b) HGPT2 vs. RS, (c) MWR vs. RS, and (d) HGPT2 vs. MWR.

3.4 Diagnostic Analysis of Thermodynamic Conversion Uncertainty

3.4.1 Diurnal Bias Amplification in Static Models

To pinpoint the physical origin of the HGPT2 model's deficiency, a diurnal cycle analysis was performed (Fig. 7). While the previous statistical metrics indicated a general positive bias, the temporal breakdown in Fig. 7a reveals that this error is not uniform, but is driven by a fundamental misrepresentation of atmospheric thermodynamics. The MWR-derived T_m (orange line) exhibits a physically realistic, dampened diurnal amplitude of approximately 1.5 K. This stability reflects the high thermal inertia of the tropospheric column, which does not heat rapidly in response to surface insolation. In stark contrast, the HGPT2 model (blue line) displays an exaggerated diurnal wave with an amplitude exceeding 8.5 K, peaking synchronously with solar noon (12:00 UTC). As previously documented in the literature (Wang, 2005; Bock et al., 2021), deriving T_m via empirical regression on surface temperature (T_s) is known to introduce spurious diurnal cycles. Our observations confirm this intrinsic limitation: because the empirical model's periodic functions are overly sensitive to T_s, it assumes intense surfcae-level heating propagates uniformly through the column, failing to capture the true thermodynamic decoupling between the turbulent planetary boundary layer and the stable free troposphere. During the hours of peak solar insolation (11:00–14:00 UTC), the coastal environment experiences active convective mixing and the onset of the sea breeze, which dramatically alters the vertical distribution of water vapour. If the underlying reanalysis climatology fails to adequately resolve the sharp moisture capping inversion at the top of the daytime planetary boundary layer (PBL), it will misrepresent the T_m weighting function. Specifically, if the model traps too much moisture near the intensely heated surface – or fails to capture the thermodynamic decoupling between the turbulent PBL and the stable free troposphere – the integral will disproportionately weight the hottest atmospheric layers. This coupled temperature-humidity mechanism physically manifests as the severe diurnal bias peak effect observed in Fig. 7b, where the systematic bias surges to over +6 K. This demonstrates that high-precision GNSS meteorology requires synergistic MWR data to capture both the true thermal stability and the dynamic vertical moisture weighting of the atmosphere.

Figure 7 Diurnal variation of weighted mean atmospheric temperature T_m derived from microwave radiometer (MWR) observations and HGPT2 model simulations (top panel). The bottom panel shows the corresponding hourly mean bias (T_m,HGPT2−T_m,MWR), with shaded envelopes indicating variability. The yellow shaded region highlights the period of maximum daytime heating.

3.4.2 Calibration and Bias Correction of MWR T_m

Figure 8 presents a statistical validation of the MWR derived T_m against co-located RS observations. The analysis highlights the necessity and efficacy of a linear bias correction scheme to improve GNSS-PWV conversion accuracy. The scatter plot (Fig. 2a) reveals a distinct systematic deviation in the original MWR retrieval relative to the RS reference. The data points consistently fall below the 1:1 identity line, indicating a negative bias in the raw MWR T_m product. The original RMSE is 2.32 K. This error is largely driven by the systematic offset rather than random scatter, as evidenced by the high linearity (R²) of the relationship. The thermodynamic profiles were retrieved using the manufacturer's standard Neural Network (NN) algorithm, trained on Region historical RS data.

The Probability Density Function (PDF) of the errors (T_m,MWR−T_m,RS) in Fig. 8b clearly visualizes the bias shift. The pre-correction distribution is non-Gaussian and shifted significantly to the negative domain, with a mean bias (μ) of −1.90 K. In the context of GNSS meteorology, a T_m error of ≈2 K translates to a relative PWV error of approximately 0.7 %–1.0 %. For climate monitoring, this represents a significant systematic dry bias. Applying the linear correction model (as formulated in Sect. 2.2.1) successfully re-centers the error distribution. The post-correction bias is reduced to 0.50 K, and the histogram aligns symmetrically around the zero-error line. The correction reduces the RMSE to 1.43 K, which is consistent with the theoretical accuracy limit of ground-based radiometric profiling (typically 1–2 K). The remaining spread (width of the green histogram) represents the random error component, likely attributable to instrumental noise and the imperfect spatiotemporal matching between the instantaneous MWR zenith view and the drifting radiosonde balloon. The correction methodology effectively removes the systematic instrumental bias without artificially compressing the natural variability of the atmosphere. The reduction of RMSE by ∼38 % (from 2.32–1.43 K) confirms that site-specific calibration of T_m is a mandatory processing step for generating climate-quality GNSS-PWV datasets.

Figure 8 Evaluation of weighted mean temperature T_m correction against Radiosonde (RS) observations. (a) scatter plots of original and bias-corrected MWR-derived T_m versus RS T_m with the dashed line indicating perfect agreement. (b) presents the probability density of errors (T_m,MWR−T_m,RS) before and after correction, demonstrating a substantial reduction in cold bias and RMSE.

3.4.3 Uncertainty Budget Analysis

In standard GNSS network processing, the largest source of PWV uncertainty is often the interpolation or modeling of surface pressure required to calculate the ZHD (Van Malderen et al., 2022). However, the CYGMEN observatory setup mitigates this spatial interpolation error by utilizing the co-located Vaisala WXT536 sensor, which has a stated pressure accuracy of ±0.5 hPa. A 0.5 hPa pressure uncertainty propagates to approximately 1.15 mm of error in the ZHD. After applying the Π conversion factor, this restricts the pressure-induced PWV uncertainty to roughly ±0.17 mm. Because this high-precision localized pressure data effectively minimizes ZHD uncertainty, the accuracy of the T_m parameterization emerges as the dominant remaining variable in the PWV error budget for this site.

Figure 9 (a) PWV uncertainty attributed to thermodynamic assumptions and to the choice of refractivity constants. (b) Difference in GNSS-derived PWV resulting from the use of alternative refractivity constant formulations relative to Bevis et al. (1994).

It is important to note that the complete error budget for GNSS-derived PWV encompasses significant uncertainties originating from the ZTD estimation phase itself. These include geodetic errors such as satellite orbit and clock uncertainties, mapping function inaccuracies, and site-dependent electromagnetic effects like signal scattering and multipath. While these geodetic factors are critical, the following component-wise uncertainty analysis (Fig. 9) specifically isolates the errors introduced during the subsequent conversion step (Π). To decouple these retrieval contributions, two primary sources of uncertainty were isolated: the thermodynamic parameterization of T_m and the selection of atmospheric refractivity constants ($k_{\mathrm{2}}^{\prime },k_{\mathrm{3}}$). When decoupling these retrieval contributions, it is critical to distinguish between the statistical nature of the underlying error sources. As demonstrated by Healy (2011), uncertainties in the atmospheric refractivity constants ($k_{\mathrm{2}}^{\prime },k_{\mathrm{3}}$) act strictly as static systematic biases; selecting a different set of published constants permanently shifts the baseline of the Π by a fixed margin. Conversely, the uncertainty originating from the T_m parameterization is a dynamic, compound error. As highlighted by Wang et al. (2005) and Bock et al. (2021), empirical T_m models derived from surface temperatures often fail to capture the true profile variance, introducing both a systematic bias (the model's mean regional offset) and a substantial random error component (the statistical scatter, or RMSE, driven by real-time thermodynamic variability and diurnal decoupling). While Fig. 9 juxtaposes these two distinct sources to illustrate their relative bounding magnitude on the final PWV product, their significantly different statistical behaviors – static bias versus dynamic scatter – must be acknowledged. As illustrated in Fig. 9a, and explicitly evaluating the components of the conversion uncertainty framework established in Eq. (10), the variance introduced by the T_m estimation strategy (${\mathit{\sigma }}_{T_{\mathrm{m}}}^{\mathrm{2}}$) significantly outweighs the influence of the physical constants (${\mathit{\sigma }}_{k_{\mathrm{2}}^{\prime }}^{\mathrm{2}},{\mathit{\sigma }}_{k_{\mathrm{3}}}^{\mathrm{2}}$). Feeding our empirically derived thermodynamic uncertainties into the partial derivative formulation defined in Eq. (11) specifically, substituting the HGPT2 RMSE of 4.54 K versus the corrected MWR RMSE of 1.43 K as our ${\mathit{\sigma }}_{T_{\mathrm{m}}}$values – yields an isolated PWV retrieval error of approximately 1–2 mm due to stochastic thermodynamic variability. In contrast, evaluating the exact mathematical limits of the refractivity coefficients (${\mathit{\sigma }}_{k_{\mathrm{2}}^{\prime }},{\mathit{\sigma }}_{k_{\mathrm{3}}}$) defined here as the maximum divergence between the historical Davis et al. (1985), the standard Bevis et al. (1994), and the updated Rüeger (2002) formulations – results in an uncertainty an order of magnitude smaller. Figure 9b further resolves the impact of the refractivity constants, showing the differential bias between the oldest (Davis) and newest (Rüeger) standards. The relationship is linear and proportional to the total water vapour content, consistent with a scaling of the Π factor. While the transition to the Rüeger (2002) constants introduces a systematic positive shift, the magnitude of this correction (typically<0.2 mm for standard loading) is negligible for synoptic meteorological applications compared to the noise induced by T_m errors. However, for long-term climatological trend analysis where stability is paramount, consistent adherence to the Rüeger (2002) standard is recommended to eliminate this small, but persistent systematic bias. Overall, the correction of the T_m is 2.5 times more important than selection of the constant.

3.5 Error Propagation and Synergistic Retrieval Assessment

In this section, the PWV was derived using bias-corrected mean temperature (T_m) and constant values based on the study by Rüeger (2002), as mentioned in the Sect. 3.4. The impact of T_m errors on the final PWV product was analyzed to quantify the benefits of the synergistic retrieval method. Figure 10 visualizes the direct relationship between the relative error in T_m and the resulting relative error in PWV. The plot reveals a strictly linear relationship (R²=0.984) with a slope of 0.981. This confirms the theoretical approximation that $(\mathrm{\Delta }\text{PWV}/\text{PWV})\approx (\mathrm{\Delta }{T}_{\mathrm{m}}/T_{\mathrm{m}})$. The color gradient indicates that this linear error propagation holds true across all PWV magnitudes (from<10 mm to >45 mm). This implies that temperature errors propagate directly into moisture errors regardless of the humidity level, making accurate T_m crucial at all times. Figure 11 investigates the systematic difference (ΔPWV) between the synergistic retrieval (using MWR T_m) and a standard retrieval (using empirical T_m) as a function of moisture abundance (PWV magnitude). For drier conditions (PWV<25 mm), the difference is minimal (near zero), and the uncertainty (shaded region) is low. This suggests that for low humidity, the choice of T_m source is less critical. As atmospheric moisture increases (>25 mm), a significant negative bias emerges. The curve dips sharply, reaching nearly −1.0 mm at extreme humidity (45+ mm). The “Systematic Bias Threshold” marker indicates that beyond 45 mm, the discrepancy exceeds 1.0 mm. The fact that the bias magnitude scales directly with total PWV provides physical confirmation that the error source is located in the boundary layer, where the bulk of the water vapour resides. The growing negative bias demonstrates that standard GNSS processing (using static models like HGPT2) systematically overestimates water vapour during extreme events compared to the more accurate synergistic method. Rather than extrapolating these localized errors to regional hydrological impacts, we emphasize the primary empirical observation: the systemic deviation of the standard empirical model scales proportionally with the magnitude of the PWV regime. Crucially, this systematic overestimation of moisture during extreme events is deeply intertwined with the diurnal cycle of the local atmosphere. This analysis quantifies the specific operational penalty of utilizing static climatological models in this region, demonstrating that HGPT2 incurs an PWV error exceeding 1.0 mm during severe thermodynamic events. As previously established (Fig. 7), the static HGPT2 model displays an exaggerated diurnal wave with an amplitude exceeding 8.5 K. Because the static model fails to account for the thermodynamic decoupling between the heated boundary layer and the cooler free troposphere during the day, this T_m error artificially inflates the amplitude of the GNSS-derived PWV diurnal cycle during peak solar insolation. By utilizing the synergistic retrieval approach, this spurious daytime moisture amplification is effectively mitigated. While further multi-site, long-term studies are required to assess the broader impacts on regional operational forecasting, our localized dataset clearly indicates that integrating real-time MWR thermal data successfully removes diurnal artifacts and reduces systematic measurement biases at this site.

Figure 10 Driver of model failure: Impact of weighted mean temperature (T_m) accuracy on PWV retrieval.

Figure 11 Systematic breakdown and instability of the Standard GNSS model under extreme thermodynamic conditions.

4 Discussions

The results of this study necessitate a fundamental re-evaluation of how T_m parameterization errors are parameterized in GNSS meteorology, particularly within thermodynamically complex, semi-arid coastal environments like the Eastern Mediterranean. The pronounced failure of the static HGPT2 model to capture the diurnal T_m cycle reveals a structural limitation inherent to empirical modeling. The observed “diurnal bias peak” effect is not merely a statistical anomaly; it represents a physical disconnect. Static empirical models rely heavily on T_s, effectively assuming that intense surface-level heating propagates uniformly through the atmospheric column. This assumption critically breaks down during the daytime in the EM, where the turbulent planetary boundary layer (PBL) aggressively decouples from the stable free troposphere. Evidence for this severe decoupling is explicitly documented in the high-vertical-resolution RS profiles collected during the campaign. Because the passive MWR struggles to effectively capture this sharp boundary – a direct result of the broad weighting functions and degraded vertical resolution inherent to its K-band observations – the instrument exhibits a “smoothing error” across the inversion layer. This structural limitation highlights exactly why applying a site-specific bias correction to the MWR's native output is a necessary prerequisite for precision GNSS meteorology. Furthermore, the failure of the reanalysis climatology to properly resolve the sharp moisture capping inversion during the onset of the daytime sea-breeze significantly corrupts the moisture-weighted T_m integral. Ground-based microwave radiometry overcomes this structural blindness by directly measuring the integrated thermal emissions of the column.

However, the performance of the MWR in this study highlights the duality of passive microwave remote sensing: it is highly proficient at retrieving integral quantities but degrades severely when resolving differential or gradient-based parameters. The successful reduction of the T_m RMSE via site-specific linear correction confirms that the MWR's K-band and V-band channels effectively capture the true thermal inertia of the troposphere. The initial systematic cold bias observed aloft is a known artifact of ill-posed neural network retrievals (Cimini et al., 2006; Löhnert and Maier, 2012). Because the vertical resolution of passive microwave observations degrades rapidly with height, the retrievals become heavily constrained by historical training datasets (the climatological prior), which often fail to capture localized, transition-season lapse rates in the free troposphere. Conversely, the complete failure of the MWR to derive a physically realistic water vapour scale height (H_v) exposes the “smoothing error” inherent to passive radiometry. Because the broad weighting functions of the K-band channels cannot resolve sharp boundary layer moisture inversions, the retrieval algorithm mathematically smears the moisture mass upward. This confirms that while MWR serves as a robust standard for total column mass, researchers must exercise extreme caution when utilizing its smoothed profiles to characterize vertical moisture compactness.

While this study relies on a single-site, multi-month dataset, the physical mechanisms identified have broad relevance beyond the Nicosia region. The Eastern Mediterranean serves as a highly representative climatic hotspot for semi-arid coastal environments experiencing enhanced warming and intensified hydrological cycles. It is important to note that the specific threshold of >45 mm identified here is characteristic of the climatological moisture capacity of the Eastern Mediterranean during extreme summer anomalies. While the exact numerical value of this “Systematic Bias Threshold” will vary geographically depending on local atmospheric dynamics and latitude, the underlying physical principle remains universal: empirical T_m models systematically degrade proportionally to the total atmospheric moisture mass during severe local extremes. The core vulnerability exposed in this research – that static global models are structurally blind to sharp boundary layer thermodynamic decoupling during peak insolation – is a fundamental physics problem, not a local anomaly. Therefore, the proposed synergistic MWR-GNSS retrieval architecture provides a universally applicable solution for mitigating systematic dry biases in any complex terrain or coastal environment globally. While the simple linear regression applied in this study proved highly effective at correcting systematic T_m biases for operational GNSS conversions, there remains room for algorithmic improvement. As the CYGMEN infrastructure accumulates a multi-year climatological database of high-resolution radiosonde profiles, future work should focus on complementary Neural Network (NN) training. By retraining the MWR retrieval algorithms using site-specific radiative transfer modeling rather than relying on the manufacturer's regional historical priors, the native temperature and humidity profiles can be further optimized at the retrieval level.

Finally, our component-wise uncertainty analysis clarifies the error propagation chain in the GNSS-PWV conversion process, shifting the paradigm of where optimization efforts should be focused. Historically, significant effort within the geodetic community has been expended on refining atmospheric refractivity constants. However, we demonstrate that the error induced by transitioning from the historical Davis et al. (1985) formulations to the modern Rüeger (2002) constants is practically negligible (<0.2 mm) for synoptic meteorological applications. The true “weak link” in the retrieval chain is unequivocally the thermodynamic parameterization, which introduces errors an order of magnitude larger.

5 Conclusion

This study demonstrated that the accuracy of GNSS-derived Precipitable Water Vapour (PWV) in the Eastern Mediterranean region, is significantly affected by the thermodynamic rigidity of static climatological models. By implementing a synergistic retrieval strategy that couples GNSS delays with real-time ground-based microwave radiometry (MWR), we successfully quantified and mitigated these limitations. The investigation yielded three primary methodological conclusions. First, we established that standard empirical models (e.g., HGPT2) are structurally incapable of resolving the diurnal thermodynamic decoupling between the boundary layer and free troposphere. This deficiency leads to severe systematic errors (the “diurnal bias peak” effect) exceeding 6 K in weighted mean temperature (T_m) during peak solar insolation, which directly propagates into a PWV bias>1.0 mm during extreme hygrometric events. Second, the MWR proved to be a superior source for T_m parameterization errors, provided that site-specific calibration is applied. The development of a linear bias correction scheme reduced the MWR T_m root-mean-square error from 2.32–1.43 K. This correction substantially reduces the conversion-related uncertainty in the GNSS water vapour product compared to standard climatological approaches. Third, the component-wise sensitivity analysis confirmed that thermodynamic parameterization is a highly significant source of uncertainty that exacerbates existing geodetic ZTD errors, outweighing uncertainties in refractive index constants by an order of magnitude. Consequently, the proposed combined retrieval represents a highly valuable architectural upgrade for monitoring severe weather in complex coastal environments like the Eastern Mediterranean. However, it must be acknowledged that there are many sites worldwide where the deployment of microwave radiometers may not be justified. Given the high capital and operational costs of radiometric hardware, the presence of other unmitigated geodetic uncertainties, and the adequate performance of static T_m models in less thermodynamically complex regions, this synergistic approach is best reserved for targeted deployments in highly vulnerable climatic hotspots.

For the climate-sensitive Eastern Mediterranean region, relying on static models for GNSS processing risks systematically masking moisture trends during heatwaves and deep convection. We therefore recommend the operational integration of collocated MWR observations into national GNSS processing chains. Where collocation is not feasible, future work should focus on assimilating MWR-derived diurnal shape functions into static models to bridge the gap between climatology and reality. This study establishes the “Corrected Synergistic Method” as a robust benchmark for generation of climate-quality water vapour datasets in complex thermodynamic environments. From an operational perspective, relying exclusively on MWRs for regional moisture monitoring is constrained by high capital costs, maintenance complexity, and signal degradation during precipitation events. Conversely, GNSS networks provide highly cost-effective, dense, and all-weather monitoring capabilities. The primary operational interest of this proposed methodology is the “supersite” calibration strategy: utilizing a centralized MWR to capture the true, real-time thermodynamic diurnal variations that static models like HGPT2 miss, and subsequently assimilating these dynamic T_m corrections over a much wider, regional network of standard GNSS receivers. This synergy allows forecasting centers to leverage the superior thermodynamic accuracy of a single MWR to drastically improve the high-resolution, continuous PWV datasets generated by dense, low-cost GNSS infrastructure.