Monitoring the Water Vapor Content at High Spatio-Temporal Resolution Using a Network of Low-Cost Multi-GNSS Receivers

The global navigation satellite system (GNSS) has the capacity for remote sensing of water vapor content in the atmosphere. Post-processing of GNSS data can provide integrated water vapor (IWV) with accuracies comparable to measurements of traditional sensors, i.e., water vapor radiometers. While GNSS meteorology benefits from thousands of permanent GNSS stations operating worldwide the spatial resolution of GNSS-derived IWV is limited to tens of kilometers. Further densification of GNSS networks is achievable with low-cost GNSS receivers. We investigated the feasibility of low-cost multi-GNSS receivers for monitoring IWV. The post-processing and the real-time (RT) solution are validated against: 1) the results from a geodetic-grade GNSS receiver; 2) colocated water vapor radiometer; and 3) numerical weather model (NWM). Despite the high variability of the IWV during the validation period, the standard deviation of IWV differences with respect to the water vapor radiometer was 1.0 and 1.5 kg/<inline-formula> <tex-math notation="LaTeX">$\text{m}^{2}$ </tex-math></inline-formula> in post-processing and RT, respectively. The city-scale variability of water vapor content in the atmosphere was monitored by a network of 16 low-cost GNSS receivers deployed in the city of Wroclaw, Poland. During rapidly changing weather conditions, the disagreement between the low-cost GNSS-derived IWV field and the NWM reached up to 5.4 kg/<inline-formula> <tex-math notation="LaTeX">$\text{m}^{2}$ </tex-math></inline-formula> and interstation IWV differences exceeded 5 kg/<inline-formula> <tex-math notation="LaTeX">$\text{m}^{2}$ </tex-math></inline-formula>. It has been demonstrated that low-cost GNSS receivers are reliable tools for precise determination of IWV, also in RT. This study is the first to measure the water vapor content with a spatial resolution of single kilometers and to present a significantly diversified city-scale IWV field.


I. INTRODUCTION
T HE water vapor is an environmentally significant constituent of the atmosphere, i.e., it is a key element of the hydrological cycle, a prominent greenhouse gas, one of key source of clouds and precipitation. Most of its quantity is confined to 4 km above the surface, but its distribution is highly variable over time and space [1], [2]. Complex lifecycle of water vapor makes its observation highly demanding [3]. Conversely, rather than the 3-D distribution of water vapor, an integrated water vapor (IWV), i.e., a liquid equivalent of the total water vapor contained in the air column, is measured. Although different devices and techniques can measure IWV, among which the most commonly used are microwave radiometers [4], sun photometers [5], and Fourier transform infrared spectrometers [6], only the global navigation satellite system (GNSS) operates continuously and in all-weather conditions [7]. A network of GNSS ground stations turns into a powerful tool for remote sensing of the troposphere, which is also called GNSS meteorology [8], [9]. The role of IWV monitoring with GNSS can hardly be overestimated. Water vapor is the subject of the climate change debate [10], [11] and GNSS reprocessed data provide homogenized decadal time series of IWV [12], [13]. Therefore, variations of IWV are continuously the subject of analysis ranging from global studies [14], [15], through regional scale [12], [16], [17], to individual sites [18], [19]. IWV allows for reconnaissance and evolution monitoring of atmospheric rivers, which are responsible for heavy precipitation [20], thus they need to be represented in ensemble forecast systems [21], [22]. IWV variability is also important for forecasting heavy rainfalls [23], monitoring and nowcasting severe storm events [17], [24], thunderstorms classification [25], and hailstorm suppression [26]. Last but not least, GNSS troposphere products, i.e., signal delays, horizontal gradients, and refractivity profiles, are independent data sources for numerical weather models (NWMs). Therefore, multiple meteorological agencies worldwide operationally assimilate GNSS-derived troposphere products [27], [28], [29], [30]. In this way, the revised initial state of a NWM leads to improvements for the forecast of geopotential heights and precipitation amount [31], more accurate water vapor field representations [32], and for the overcast conditions [27]. The capacity of GNSS as an accurate sensor of atmospheric water vapor was already demonstrated decades ago [8], [33]. The main troposphere product of GNSS, i.e., zenith total delay (ZTD) divides into hydrostatic and wet parts. The latter can be used to quantify the IWV using other well-determined meteorological data i.e., pressure and temperature [34], [35] of GNSS observations provides IWV with the accuracy of about 1-2 kg · m −2 [36], [37], [38]. As the GNSS technology and products have improved, the potential of GNSS meteorology as a valuable source of IWV has been noticed. Nowadays deriving water vapor from near real-time (RT) GNSS troposphere products [39], [40], [41] exists as an operational service provided by the European Meteorological Network (EUMETNET) GNSS water vapor program (E-GVAP, http://egvap.dmi.dk). The reported quality of ZTD estimates from near RT global positioning system (GPS) data processing are 3-10 mm [40], [42], which corresponds to 1-3 kg · m −2 of IWV roughly. Over the last decade, major efforts toward improved troposphere products were made, mostly within European Cooperation in Science and Technology (COST) Action ES1206 "Advanced GNSS Tropospheric Products for Monitoring Severe Weather Events and Climate" [43]. The development of new GNSS constellations led to the increasing number of operational GNSS satellites, thus allowing to sense the atmosphere more evenly. Therefore, the capacity of various GNSS was investigated [44], [45] and operational multi-GNSS services are being developed [46]. Moreover, efforts toward RT GNSS meteorology were made [47], [48], [49]. Although the accuracy of RT satellite orbits, which is not as good as the one of post-process orbits [50], [51], [52], limits the accuracy of RT troposphere products, recent studies reported the RT ZTD accuracy of 5-10 mm, IWV equivalent of c.a. 1.5-3 kg · m −2 [47]. GNSS meteorology benefits from thousands of GNSS permanent stations operating worldwide, mostly for geodetic purposes. With over 2000 station in Europe, troposphere is sensed with a spatial resolution of c.a. 30 to 100 km. GEONET in Japan (https://www.gsi.go.jp), which is considered as the most densified country-wide GNSS network with over 1300 stations, improves the spatial resolution to c.a. 20 km. Yet, none of these networks allows to monitor the troposphere with a spatial resolution of single kilometers, which is already considered in regional climate and weather models [53], [54]. Neither these networks, nor any other sensing technique, are feasible to monitor regional variability of water vapor distribution. However, deployment of new GNSS stations is limited by economic reasons. Therefore, further densification of GNSS networks is achievable with mass-market chipsets capable of logging carrier phase measurements, which are available for a few years. It has been demonstrated that dual-frequency low-cost GNSS receivers are suitable for precise positioning [55], [56], [57] and hazard monitoring [58], [59], [60] under favorable conditions, but their position determination performance is worse than the geodetic-grade receivers [61]. Despite this, post-processed ZTD estimates from low-cost GNSS receivers are of similar accuracy as from geodetic-grade receivers [62], [63], [64], thus showing a great potential of low-cost devices for meteorological applications. So far, limited research on low-cost GNSS meteorology focuses on post-process ZTD results and singlestation experiments. In this study, the city-scale variability of water vapor content in the atmosphere is monitored by a dense network of low-cost GNSS receivers for the first time.
Moreover, in addition to post-processing regime, RT solution is demonstrated, and both solutions are validated against the colocated water vapor radiometer and NWM of highest spatiotemporal resolution.

A. Low-Cost Receivers and Antennas
In all experiments, the in-house developed low-cost receivers were used. The devices were equipped with a ublox ZED-F9P high-precision GNSS module, microcomputer Raspberry Pi 3B+, and the global system for mobile communications (GSM) modem (Fig. 1). The receiver module tracked the following signals: GPS (G) L1 C/A and L2C, Galileo (E) E1 B/C and E5b, Global'naya Navigatsionnaya Sputnikovaya Sistema (GLONASS) (R) L1OF and L2OF. BeiDou signals were not recorded due to the limited number of channels of the GNSS receiver module, i.e., it could trace up to 32 satellites simultaneously. A self-made Python script was developed to store carrier-phase and code observations on both frequencies in receiver independent exchange system (RINEX) 3.04 format.
The u-blox C099-F9P was delivered with a u-blox ANN-MB-00-00 patch antenna and a circular ground plane. Although the maximum phase center offset (PCO) and phase center variations (PCV) information was provided by the manufacturer, they did not agree with the GPS-only PCO/PCV developed by [65]. Therefore, we used an in-house developed quad-system PCO model for u-blox antennas [66]. One device (BX02) was connected to a survey-grade ArduSimple AS-ANT2B-SUR-L1L2-25SMA-00 antenna, for which the PCO/PCV information remained unavailable.

B. Study Area
From December 5, 2020, the low-cost devices were successively deployed in the urban area of Wrocław city, Poland, and its suburbs (Fig. 2). Between February 27 till March 28, 2021, 16 devices were active. The test area of 243 km 2 (extension of 14 km North-South and 24 km East-West) was covered.

C. GNSS Processing Strategy
Two strategies were used to process GNSS observations, namely RT and final (FIN). Both strategies applied pseudorange P and carrier phase L dual-frequency undifferenced where i is the number of the frequency f with wavelength λ i and a corresponding ambiguity N s i ; ρ s 0 is the geometric distance between the position of the satellite s and antenna phase center (X s , Y s , Z s ) and the a priori position of the receiver r antenna phase center (X r , Y r , Z r ); c is the speed of the light; δt s is satellite clock offset; s r represents satellite, receiver, and site displacement effect corrections according to International Earth Rotation and Reference Systems Service (IERS) Convention 2010 [69]; m s h and m s w are hydrostatic and wet-mapping functions for the a priori zenith hydrostatic delay Z h and estimated zenith wet delay (ZWD) Z w , respectively; e is the direction vector and δX r = δ X r δY r δ Z r T is the position correction vector in the Earth centered Earth fixed (ECEF) frame; δt r is the receiver clock offset; δt sys G is the receiver intersystem bias between GPS and another system sys; I s is the slant ionosphere delay. Corrections δX r to receiver static coordinates were estimated together with Z w , I s and N s , whereas the Vienna mapping functions 1 (VMF-1, [70]) were the source of Z h , m s h and m s w . For the RT strategy, an in-house developed software GNSS-WARP [71] was used to process GPS and Galileo observations. The selection of GNSS in RT strategy was motivated by previous research by Hadas and Hobiger [72], as it results in increased precision and suppresses orbit-related artificial signals. For FIN strategy, GPS and GLONASS observations were processed with the Canadian spatial reference system (CSRS)-PPP service (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php), because the CSRS-PPP service supports only these two constellations. Such limitations result from the availability of the satellite orbit and clock corrections of the highest accuracy [51], [73], [74]. Therefore, in the RT strategy, the calibration model from the International GNSS Service (IGS) repro3 campaign was used, because it contained Galileo PCO/PCV, whereas for the FIN strategy the official IGS model was applied. The processing strategies are summarized in Table I in more detail.

D. Determination of IWV
ZWDs can be converted into a time series of precipitable water (PW [m]) using the following formula [33], [75]: where ρ w = 997 kg · m −3 is the density of liquid water, R w = 461.525 J·kg −1 ·K −1 is the gas constant of water vapor, T m [K] is the mean temperature along the GNSS signal propagation path in the troposphere and is calculated from the surface temperature T s [K] at the location of a GNSS station, M d = 28.9644 g · mol −1 and M w = 15.9994 g · mol −1 are a molar mass of dry and wet air, respectively, k 1 = 0.776890 K · Pa −1 , k 2 = 0.712952 K · Pa −1 , and k 3 = 3754.63 K 2 · Pa −1 . The PW is a measure of the equivalent height of the column formed from the condensed water. Therefore, a PW value of 1 mm is equivalent to an IWV value of 1 kg · m −2 . The estimated ZWD ( Z w ) absorbs the residual hydrostatic delay δ Z h due to inaccurate a priori value of the zenith hydrostatic delay Z' h Therefore, it is the zenith total delay ZTD G that remains the only unbiased troposphere parameter in GNSS processing The reconstruction of the Z w from the ZTD G requires an unbiased Z h from a reliable source, such as an NWM. However, the ambiguous transformation of equipotential surfaces, i.e., pressure-level data, in NWM to the vertical reference frame of GNSS, i.e., the orthogonal height system, offsets NWM-derived hydrostatic delay Z M h and wet delay Z M w , and so the NWM-derived total delay ZTD M Since Z w is one order of magnitude smaller than Z h , the Z w is one order of magnitude smaller than Z h , so Z w is negligible. Hence, the misalignment between ZTD G and ZTD M results directly from the Z h . Therefore, the commonly accepted procedure to determine the Z h is to calculate the mean difference between the two ZTDs over a certain past period of time where · is the notation for mean. Combining (6)-(9), the following formula for obtaining the unbiased Z w is achieved:

E. Reference Data
Final troposphere estimates for station WROC were obtained from the Warsaw University of Technology (WUT), Poland, which is the analysis center of the Regional Reference Frame Sub-Commission for Europe (EUREF) permanent network (EPN). Estimated uncertainties of ZTD and horizontal gradients were at the level of 0.3 mm and below 0.04 mm, respectively. The surface temperature was obtained from the high-resolution NWM weather research and forecasting Fig. 4. Probability histograms of ZTD differences between the two processing strategies and the reference product for station BX02 (left) and WROC (right), respectively. model (WRF; [76], [77]. Moreover, Z M h , Z M w , and ZTD M were determined from the WRF model using a ray-tracing method [78]. All atmosphere parameters were obtained for locations of low-cost devices as well as for their orthogonal projections on the horizontal layer at the height of station WROC. Last but not least, discrete IWV measurements from the water vapor radiometer RPG-HATPRO-G5 (Humidity and Temperature PROfiler) were used. The radiometer was colocated with GNSS station WROC, i.e., the distance between the two instruments was smaller than 10 m. RPG-HATPRO-G5 is a series of microwave radiometers based on passive observations of atmospheric components' thermal emissions in the microwave spectral range [79]. The observed radiances are converted into brightness temperatures; further, the meteorological variables (e.g., IWV) are retrieved using artificial Fig. 5. Time series of IWV from the two GNSS receivers obtained with RT and FIN strategies, and from the colocated water vapor radiometer; σ and μ indicate standard deviation and mean difference between IWV from GNSS and IWV from radiometer. IWV [kg/m 2 ] scatterplots with correlation coefficient r for the two colocated stations BX02 and WROC, and both processing strategies, i.e., RT and FIN; color represents normalized density; r indicates Pearson's correlation between solutions. neural network algorithms [80]. The estimated uncertainty of the radiometer's ZTD measurements equals 7.4 mm in terms of the root mean square error (RMSE) when compared to the radiosonde observations [81].

III. RESULTS
Two research periods were distinguished, namely the validation period and the campaign period. The validation period A. Quality Control 1) ZTD Validation: Due to a short distance of 7.2 m between stations WROC and BX02, as well as a small height difference between both stations, i.e., 0.17 m, the direct comparison of ZTD from the low-cost device and the geodetic grade receiver was possible. Troposphere estimates for the two colocated stations were obtained using both processing strategies, i.e., RT and FIN (Fig. 3). A ZTD of estimated uncertainty δZTD larger than three times the RMSE of δZTD was identified as an outlier and rejected from further analysis. The percentage of obtained solutions during the test period was 100% for both FIN solutions, whereas for RT solutions was 97.3% and 97.8% for station BX02 and station WROC, respectively. The major gap in RT solutions around midnight of December 25, was due to the failure in the reception of RT satellite orbit and clock correction products. During the validation period, ZTD varied from 2270 mm to almost 2420 mm, with several major and rapid ZTD changes. Nevertheless, a high agreement was found among the obtained time series, as well as with the reference product from the EPN. In general, δZTDs were smaller for station WROC than for station BX02. For both FIN solutions, regular peaks in δZTD were noticed, which were due to midnight discontinuities in the satellite orbit and clock products, but did not affect estimated ZTDs noticeably. Unexpectedly, for station WROC larger δZTDs were obtained with the FIN strategy than with the RT strategy, which indicated differences in the PPP filter implemented in both processing engines. Although underestimated uncertainties of products are typical for an epoch-wise RT PPP processing, the centimeter-level δZTD resulted from tight coordinate constraining, which might lead to a smoothed ZTD product.
The ZTD differences (ZTD) between estimated values and the reference product from EPN (Fig. 4) proved the high reliability of obtained solutions. The mean ZTD of 0.3, −0.1, 0.2, and −0.4 mm for RT BX02, FIN BX02, RT WROC, and FIN WROC, respectively, indicated that the obtained results were free of systematic errors. Standard deviations of ZTD at station WROC varied from 1.6 to 3.6 mm for the FIN and RT solution, respectively, and at station BX02 from 6.2 to 8.8 mm, respectively. Therefore, all results met the target requirements of the E GVAP program both for subhourly and post process ZTD products, set to 10 and 7 mm, respectively, [82].
For BX02, a limited improvement in accuracy and precision of ZTD was noticed with the FIN strategy when compared to the RT strategy (Fig. 4 left), whereas for station WROC a significant improvement in precision was achieved with the FIN strategy (Fig. 4 right), i.e., the percentage of ZTD within the range of +/−5 mm increased from 80% to 98%.
2) IWV Validation: The unbiased Z w were obtained with (10). The corresponding IWVs were calculated with (4) and compared with IWV from the colocated water vapor radiometer (Fig. 5). Irrespective of the type of the GNSS station (i.e., geodetic-grade or low-cost) and processing strategy (postprocessing or RT), GNSS-derived IWV was underestimated (less moist) with respect to the radiometer by 0.3 kg/m 2 . Even more significant biases between GNSS and radiometer data were identified in several studies [81], [83] and could be explained by the low long-term stability of water vapor radiometers and the requirement of the calibration procedure [84], [85]. Despite the high variability of the IWV during the validation period, i.e., from 5 kg/m 2 (December 27, 2020) to 25 kg/m 2 (December 23, 2020), the dynamic changes of IWV were captured well by the GNSS and fit the water vapor radiometer data. For the FIN strategy, standard deviations σ of differences between GNSS and water vapor radiometer ranged from 0.5 to 1.0 kg/m 2 for station WROC and station BX02, respectively. For the RT strategy, the corresponding σ were 0.8 and 1.5 kg/m 2 . Therefore, the IWV precision obtained from a low-cost receiver compared to IWV obtained with a geodetic grade receiver was worse by a factor of 2, whereas the change of the strategy from FIN to RT reduced the IWV precision by a factor of 1.5.
Strong correlations between IWV time series were found (Fig. 6), revealing a very good agreement between all solutions and for the entire range of calculated IWV values. The largest correlation coefficient r equaled 0.99 and was found for the IWV time series at station WROC obtained with RT and FIN strategies. Therefore, for a geodetic grade receiver, the lowest accuracy of RT satellite orbits and clocks was of minor importance and did not significantly affect ZTD which was used to calculate IWV. Slightly weaker correlations were found for the IWV time series obtained with two different stations; r equaled 0.97 and 0.98 for the RT and FIN strategies, respectively. This comparison revealed the differences in the quality of observations acquired with a geodetic grade and a low-cost receiver. The lowest r equaled 0.95 and was found for the IWV time series at the low-cost station BX02. This indicated that IWV products obtained with the low-cost receiver benefited from the higher accuracy of post-process orbit and clock products.

B. Field Campaign
The field campaign started on February 27, 2021, with three active low-cost devices and the geodetic-grade station WROC. Further low-cost stations were successively deployed, but on March 10 the BX02 station, colocated with station WROC, became malfunctioning. A maximum of 15 stations operating simultaneously was reached on March 16 and, with minor outages, lasted for the following 13 days, until March 28, 2021. The average availability of GNSS observations during the operation of receivers was 97% (Fig. 7). 1) ZTD: ZTD estimates for all test stations were obtained using both processing strategies, i.e., RT and FIN. For each receiver, the first 24 h of results were removed, thus estimates from the initialization period of the PPP were disregarded. With the RT strategy, no solution was obtained for a few periods due to the interrupted transmission of RT satellite orbit and clock corrections. The direct comparison of obtained ZTD time series (Fig. 8) revealed the following. Varying conditions of the atmosphere were captured during the test campaign, with ZTD change amplitude exceeding 140 mm and a few dynamic events, e.g., on March 11 and 27, 2021. ZTD time series obtained with the FIN strategy were more smooth than those obtained with the RT strategy. This was due to: 1) the use of precise satellite orbits and clock corrections in the FIN strategy, rather than RT products and 2) the tighter constraining of the ZTD in the PPP filter of the FIN strategy, justified by the primary goal of the CSRS-PPP service, i.e., precise coordinate determination. Moreover, offsets between stations were noticed, which were justified by different heights of the test receivers, ranging from 155.38 m for station BX10 to 190.09 m for station BX22. Therefore, before IWV determination, height correction to ZTD (Z H ) was applied for each station independently, so that the subsequent spatial analyses of IWV were conducted on a common height level, i.e., at the height of station WROC (180.81 m). Z H was determined using ZTD difference in the WRF at station location but at two levels, i.e., the estimated height and the height of station WROC. Z H ranged from −3 mm (station BX22), through 0 mm (station BX02), to +8 mm (station BX10).
2) IWV: The Z h was determined for each station independently with (9), using the height-corrected ZTD estimated with GNSS and ZTD from WRF at the height of station WROC. Z h varied from −2 mm (stations BX17 and BX13) to +22 mm (station BX10). Z h differed not only among stations but also between the two strategies by several mm, which we justified by: 1) the missing PCO model in the FIN solution (Table I) and 2) the limited accuracy of the receiver height determination, especially with the RT strategy. Nevertheless, the application of Z h led to the unbiased determination of the Z w with (10), and consequently-to the unbiased IWV using (4). The obtained IWV time series were quite consistent among test stations. The correlation coefficient r for IWV estimated at two stations varied from 0.69 to 0.98 for the RT strategy and from 0.86 to 1.00 for the FIN strategy (Fig. 9). Several missing r between BX02 and other stations were due to overlapping periods of operation shorter than three days. Relatively lower r was obtained with the RT strategy for stations BX05 and BX17. The worse performance of both stations resulted from partially obstructed rather than open-sky measuring conditions, revealing that low-cost receivers were sensitive to multipath interference.
Although the highest r was found for the two colocated stations, i.e., WROC and BX02 (r = 0.98 with the RT strategy and r = 1.0 with the FIN strategy), there was no clear evidence that the correlation coefficients drop with the increasing interdistance between stations (Fig. 10). Therefore, either the size of the network was too small to capture such an expected phenomenon, or the precision of the estimated IWV was insufficient. Another reason could be that for most of the time IWV field remained relatively uniform for the entire test area, whereas a diversified IWV field occurred during short events, thus not affecting r significantly. This was confirmed by investigating the time evolution of the standard deviation of IWVs at test stations. Typically, for the RT strategy, it remained below 1.0 kg/m 2 but exceeded 2.0 kg/m 2 during four periods, with a maximum value of 5.2 kg/m 2 on March 23, 2021. For the FIN strategy, the standard deviation of IWV usually oscillated between 0.4 and 0.7 kg/m 2 , with only three periods with the standard deviations larger than 1 kg/m 2 , and the maximum value of 1.2 kg/m 2 on March 27, 2021.
All stations captured varying conditions, including dry and moist events, with IWV ranging from c.a. 0.8 to 25.0 kg/m 2 . IWV at stations BX05 and BX17 often varied from IWV estimated for other stations, hence these two stations introduced the significant noise on the lasagna plot of the IWV time series (Fig. 11). Most of the time, IWV estimates were similar among stations. However, a few periods with significantly varying estimates were noticed, with the most prominent example on March 23, 2021, for the RT strategy, with IWV ranging from 1.9 to 21.6 kg/m 2 . A maximum IWV difference between two stations reached 16.6 kg/m 2 (stations BX10 and BX14) and 3.6 kg/m 2 (stations BX15 and BX09), for RT and FIN strategy, respectively, whereas typically, the interstation IWV differences remained below 2.4 kg/m 2 and 1.6 kg/m 2 for the RT and FIN strategy, respectively. Therefore, IWV estimates on March 23 indicated significant variability of the water vapor on the regional scale.
GNSS-derived IWVs usually agreed with the WRF model (Fig. 12). The absolute IWV differences between GNSS and WRF were smaller than 1 kg/m 2 for 60% and 75% for RT and FIN strategies, respectively. However, several periods were noticed with larger discrepancies at all test stations. The WRF model underestimated IWV with respect to the RT strategy by up to 12.2 kg/m 2 (station BX14, March 10, 2022) Station-specific IWV differences (GNSS minus WRF) were zero mean, with standard deviations smaller than 1.5 and 1.0 kg/m 2 , for RT and FIN strategies, respectively (Fig. 13). Exceptions were found only for stations BX05 and BX17, for which standard deviations reached 2.0 km/m 2 for the RT strategy and 1.4 kg/m 2 for the FIN strategy. Therefore, a good agreement between GNSS and WRF was found. Outliers were identified as IWV differences greater than three times the station-specific standard deviation of IWV differences. There were 1.4% and 0.7% of outliers in RT and FIN strategy, respectively, thus indicating the intermittent occurrence of periods with significant inconsistencies between GNSS and WRF during the test campaign.
In aiming to analyze the spatio-temporal dynamics of the IWV, IWV fields for the test area were investigated epoch by epoch. Discrete IWV values at locations of test stations were complemented with IWV values determined from the WRF model at those model nodes, which were outside the area covered by operating stations. Therefore, GNSS-derived IWV fields were observed on top of the background IWV fields from the WRF model. Representative epochs were selected, as shown in Figs. 14 and 15 for RT and FIN strategies, respectively.
With the RT strategy, both dry and moist periods were observed which matched well the NWM WRF [ Fig. 14(a) and (b)], thus reflecting a unified and unbiased IWV field, despite changing air humidity. A unified state of the IWV field, also during moderate wet periods [ Fig. 14(c)], was typical during the test campaign, which proved the high  RT strategy, the differences between the GNSS-derived IWV and the WRF model were significant. Therefore, it was the WRF model which derived an inaccurate IWV field and RT GNSS estimates could potentially calibrate the WRF model, thus improving its accuracy. During some extreme cases, not only did the interstation IWV differences fairly exceed 5 kg/m 2 [ Fig. 14(f)], with a maximum interstation difference of 16 kg/m 2 , but varied significantly from the IWV field of the WRF model [ Fig. 14(g) and (h)]. Whereas the WRF model returned an inaccurate and unified IWV field, with GNSS a spatially varying IWV field was observed. Such IWV differences proved, that water vapor content in the atmosphere may vary significantly on a regional scale, and dense regional GNSS networks may add extra information to a NWM. Moreover, during highly dynamic changes of the IWV, e.g., from less than 5 kg/m 2 to over 20 kg/m 2 within 6 h [ Fig. 14(g) and (h)], the GNSS captured even more extreme IWV states than the WRF model. While the WRF model flattered the water vapor content and provided less dynamic conditions, the loosely constrained GNSS estimates reflected a significantly diversified IWV field, with a clear dry island in the eastern part of the city in the morning [ Fig. 14(g)] and two wet islands in the suburbs at noon [Fig. 14(h)].
With the FIN strategy, similar observations were made as with the RT strategy. Typically, the agreement with WRF and GNSS was good under changing but spatially unified IWV field [ Fig. 15(a)-(c)]. Offsets between the IWV fields from WRF and GNSS were usually present during the same epochs as they were noticed with the RT strategy [ Fig. 15(c) and (d)], however the differences, in the absolute sense, were smaller, but still significant. On the other hand, for just a single period of c.a. half a day, it happened that the IWV from the FIN strategy did not reflected the varying state of the water vapor content as captured with the RT strategy [ Fig. 15(f)]. Moreover, the FIN strategy resulted with the IWV field that was in a very good agreement with the WRF model. Such an inconsistency between the GNSS results obtained with both strategies were justified by the missing RT orbit and clock corrections on March 22, shortly before the period of interest, that degraded the RT results. The FIN strategy successfully Fig. 15. IWV field for the test area obtained with the FIN strategy at representative epochs (a-h); GNSS-derived IWV is interpolated between the low-cost GNSS stations, whereas the surrounding IWV comes from the WRF model. captured the dry island in the eastern part of the city in the morning of March 27 [ Fig. 15(g)]. However, the tighter constraining of the ZTD in the PPP filter of the FIN strategy smoothed the dynamic changes of the IWV and prevented from capturing extremely wet conditions, i.e., the estimated IWVs were even smaller than in the WRF model [ Fig. 15(h)].

IV. CONCLUSION
In this study, we investigated the feasibility of low-cost multi-GNSS receivers for monitoring regional dynamics of water vapor content in the atmosphere. We used in-house developed hardware based on the u-blox ZED-F9P highprecision GNSS module and deployed 16 stations in and around the city of Wrocław, Poland. Using multi-GNSS observations, ZWDs were estimated in RT using the original software WARP and in the post-processing using the online tool CSRS-PPP. The IWV was calculated using estimated ZWD and surface temperature determined from the high-resolution NWM WRF. Since one of the low-cost stations was colocated with the IGS station WROC and the water vapor radiometer, it was possible to validate estimated ZTDs and IWVs against external GNSS product and independent technique, respectively. During the validation period, i.e., December 9-28, 2020, ZTD ranged from 2270 to 2420 mm with several major and rapid ZTD changes. Standard deviation of ZTD differences between low-cost GNSS and IGS Final product was 6.2 and 8.8 mm in post-processing and RT, respectively. Therefore, the target requirements of the E-GVAP program (7 and 10 mm, respectively) were met. Although IWVs obtained with low-cost GNSS were underestimated with respect to the measurement from the water vapor radiometer by 0.3 kg/m 2 , occurrence of such a bias is well explained in the literature [see Section III-A2)]. Despite the high variability of the IWV during the validation period, i.e., from c.a. 5 kg/m2 to almost 25 kg/m2, the standard deviation of IWV differences was 1.0 and 1.5 kg/m 2 in post-processing and RT, respectively. Therefore, it has been demonstrated that low-cost GNSS receivers are reliable tools for precise determination of IWV, also in the RT. For 30 consecutive days, i.e., from the beginning of February 27 till the end of March 28, 2021, 16 low-cost GNSS stations operated in the field, covering an area of c.a. 14 × 24 km. GNSS observations were processed both in post-processing and RT, and the corresponding IWVs were determined. We noticed that low-cost GNSS receivers suffered from operating in partially obstructed sky conditions, indicating their sensitivity to multipath interference. Nevertheless, high agreement between IWV time-series among stations was noticed, irrespective of the processing regime. No evidence was found that the correlation coefficient for IWV estimated at two stations drops with the increasing distance. It was justified by a small size of the test area, over which the IWV field remained uniform for most of the time. Several periods with significant disagreement between GNSS-derived IWV field and WRF model were noticed, up to 5.4 kg/m 2 , especially during rapidly changing weather conditions. While GNSS captured more extreme dry and wet conditions, the WRF model, operating at 6 h interval, flattered the water vapor field over time and space. Moreover, periods with significant interstation IWV differences were noticed, i.e., exceeding 5 and 3 kg/m 2 in RT and post-processing, respectively. To our knowledge, it was the first time when water vapor content was measured with the spatial resolution of single kilometers and when a diversified regional-scale IWV field was observed.

DATA AVAILABILITY
The dataset containing troposphere products for the campaign period are available for download from https://zenodo. org/record/5196311. Raw GNSS observations from low-cost receivers are available at https://zenodo.org/record/6488497.

DECLARATION OF COMPETING INTEREST
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.
ACKNOWLEDGMENT Software for managing the network of low-cost receivers was provided by U+Geo Ltd. (www.uplusgeo.pl). The authors acknowledge colleagues who provided their facilities for the location of low-cost receivers during the campaign period.