GPS-Based Multi-Temporal Variation in Precipitable Water over the Territory of Poland

: An increase in temperature causes higher evaporation of water from water bodies; con-sequently, the water content in the atmosphere also increases. The precipitable water (PW), as the water content in the atmospheric air column, is therefore an important parameter to consider when studying climate change. The aim of this study was to analyse multi-annual precipitable water data derived from a dense Global Navigational Satellite Systems (GNSS) network. Twelve years of observations from over a hundred ASG-EUPOS stations were used to estimate changes in precipitation water values over Poland. The data were validated by comparison with the available radio-sounding data. The analysis of the GPS-based PW values showed an upward trend in the PW value of 0.078 mm/year. The spatio-temporal distribution of the mean PW values and their ﬂuctuations over the years were studied and visualised in the form of maps. The results are congruent with the fact that Poland lies on the border of inﬂuence of both continental and oceanic climates. Our results are also consistent with other climate research concerning this region. inﬂuence of both continental and oceanic climates. The obtained values of PW, as well as their spatio-temporal distribution, are consistent with the climate research concerning this region. This conﬁrms that GPS data can contribute to meteorological research, including climatic studies. This work also demonstrates that the development of geodetic infrastructure brings wide beneﬁts to earth sciences.


Introduction
Environmental changes visible on Earth, whether natural or caused by human activity, influence climate change on a global scale [1,2]. Therefore, it is necessary to constantly monitor these changes and study the effect of human activity on them. One of the parameters indicating climate change is the systematic increase in temperature for the last 80 years [3,4]. This increase in temperature causes a higher evaporation of water from water bodies, resulting in an increase in the water content in the atmosphere. Since water vapour is one of the main greenhouse gases in Earth's atmosphere, this cycle repeats. The concept of precipitable water (PW) describes the water content (in various states of aggregation) in a column of atmospheric air that has fully condensed to form a layer of a given height (expressed as the height of the layer formed after it is fully condensed) [5]. PW shows a significant correlation with the daily rainfall intensity (i.e., the average yield of a precipitation day), while this pluviometric indicator varies negligibly across the country [6]. For these reasons, PW is considered to be a parameter that is extremely important in the context of studying climate change.
According to [7], the average size of the PW over Poland is 15 mm. Seasonal changes are also visible, with the lowest values recorded in the winter (January), and the highest in the summer (July). The PW fluctuations depend on temperature changes, which, in turn, determine the air moisture capacity and transport closely related to the inflows of air masses with different characteristics. The fall-winter season and early spring are characterised by lower PW values, mainly due to the Atlantic Ocean, which is characterised by relatively low temperatures and limited evaporation. From May to October, a warm half-year in Poland, PW displays significantly higher values [8,9].
Tropospheric water vapour can be obtained through modelling using ground meteorological data as well as through independent measurement sensors, such as radiometers, spectrometers [10], and upper-air radio soundings. In recent years, PW can also be studied using Global Navigational Satellite Systems (GNSS) observations. The GNSS signal passing through the atmosphere is refracted depending on the atmospheric state. The usefulness of GNSS observations in modelling the troposphere comes from the fact that the GNSS signal is delayed (T) when passing through the tropospheric layers [11]. The refractivity of the troposphere (N) depends on the temperature, pressure, and humidity. The nature of the troposphere enables N to be expressed separately for dry gases (hydrostatic component) and water vapour and condensed water in clouds (wet component). The hydrostatic part represents the dominant component and, due to its slow variability in time and space, is easy to model. The wet parts, as well as tropospheric gradients in horizontal directions [12], are estimated during GNSS data processing. The wet part of the troposphere, which is an indicator of air humidity, is then converted to Integrated Water Vapour (IWV) or PW. In the last decade especially, an increase in GNSS observations in meteorology can be observed. The analysis of GNSS-derived atmospheric parameters [13] showed that they are adequately consistent with numerical models or measurements with radiometers. The distribution of GNSS stations makes them a valuable data source for meteorological studies. Due to the high density of measurement data, changes in the water vapour content can be monitored in detail and their characteristics and variability can be thoroughly analysed and visualised. On this basis, it is possible to monitor the passage of atmospheric fronts [14] or follow the route of cyclones [15,16].
In the last decade, there have been several GNSS-based studies of the troposphere over Polish territory. In [17], the authors analysed one year of data to verify the methods of analysis and filtering of zenith delays. Other works have focused on the monitoring of severe weather [18] or near real-time troposphere products [19][20][21]. Long-term analyses were conducted only on the basis of the reanalysis being carried out within EPN-Repro1 or EPN-Repro2 projects [22] on the EUREF Permanent GNSS Network (EPN). In [23] authors estimated the trends in tropospheric delays for five Polish GNSS stations. Depending on the length of the analysed data (16 and 18 years), trends ranged from −0.14 to 0.42 mm/year. For the same five Polish stations, the changes between both data sets were compared in [24]. Two of those stations were investigated in subsequent studies [25]. These studies show the increase in water vapour content in the atmosphere over Poland territory. However, the lack of sufficiently dense data prevents a detailed analysis of their spatial distribution and their changes over time.
The aim of this study was to analyse multi-annual precipitable water data originating from a dense GNSS network established in 2008. More than twelve years of observation is sufficient to estimate changes in precipitation water values over Poland. Most importantly, these data allow the determination of how much the area of Poland is diversified in terms of the spatial and temporal variability of PW. The data, together with cartographic visualisation in the form of maps of analysed variability and changes in PW in the period 2009 to 2020, constitute a valuable resource for studies on climate changes over Poland.

Materials and Methods
In the presented research, the observations from reference stations belonging to the Polish network ASG-EUPOS were used [26]. The network was launched in 2008 as the official densification of the European Terrestrial Reference System 1989 in Poland [27]. From the beginning, it has been collecting GNSS observations from c.a. 100 stations. All collected data stored in Receiver Independent Exchange System format v.2.11 were processed according to the Guidelines for the International Association of Geodesy Reference Frame Sub-commission for Europe (EUREF) Densifications in the GAMIT software [28]. Forty-seven stations belonging to the EPN [29] were included in the analysis as a reference and to improve the geometry of the network. Additionally, cross-border partner stations from the Czech Republic and Slovakia were added to the calculation. Finally, the analysis covered 157 stations ( Figure 1). Forty-seven stations belonging to the EPN [29] were included in the analysis as a reference and to improve the geometry of the network. Additionally, cross-border partner stations from the Czech Republic and Slovakia were added to the calculation. Finally, the analysis covered 157 stations ( Figure 1). The processing was based on GPS observations only. All archival data from the period 2 June 2008 to 28 January 2017 were reprocessed using the International GNSS Service products, IGS08, the same as in the study [30]. Observations from 29 January 2017 till 31 December 2020 were processed using the IGS14. Tropospheric delay was modelled as the functions of elevation angle (e) and azimuth (a): where ZHD denotes the Zenith Hydrostatic Delay, ZWD denotes the Zenith Wet Delay, is the north-gradient component, and is the east-gradient component. The three mapping functions were used for transition from the slant direction to the zenithnamely: the hydrostatic component ( ( )), wet component ( ( )), and gradients ( ( )). Vienna Mapping Function grids [31] were used as a priori values for ZHD as well as for coefficients for mapping functions (mfh and mfg). ZWD were estimated in hourly intervals. Additionally, the values of the horizontal gradients were estimated once a day. A summary of the processing parameters is presented in Table 1.  The processing was based on GPS observations only. All archival data from the period 2 June 2008 to 28 January 2017 were reprocessed using the International GNSS Service products, IGS08, the same as in the study [30]. Observations from 29 January 2017 till 31 December 2020 were processed using the IGS14. Tropospheric delay was modelled as the functions of elevation angle (e) and azimuth (a): T(e, a) = m f h (e)·ZHD + m f w (e)·ZWD + m f g (e)·[cos(a)·G N + sin(a)·G E ], (1) where ZHD denotes the Zenith Hydrostatic Delay, ZWD denotes the Zenith Wet Delay, G N is the north-gradient component, and G E is the east-gradient component. The three mapping functions were used for transition from the slant direction to the zenith-namely: the hydrostatic component (m f h (e)), wet component (m f w (e)), and gradients (m f g (e)). Vienna Mapping Function grids [31] were used as a priori values for ZHD as well as for coefficients for mapping functions (mf h and mf g ). ZWD were estimated in hourly intervals. Additionally, the values of the horizontal gradients were estimated once a day. A summary of the processing parameters is presented in Table 1.
The value of the mean temperature (T m ) was determined according to the formula (3) based on the temperature from the GNSS station (T s ). Only fifteen stations are equipped with meteo sensors. Therefore, for all stations T s was taken from the ERA5 [32].
The obtained GPS PW values were averaged to daily and monthly intervals. The study focused on monthly, annual, and long-term variations in PW. Although they were relatively high with rapidly passing atmospheric fronts, the daily GPS PW fluctuations were not considered in the context of long-term analyses. Values determined as outliers according to the 3σ criterion were removed. Prepared in this way, two sets of time series were used to further analyse data, preceded by homogenisation. If the data gaps were too large, the period was excluded from the analysis. Analysis was limited to a full year from 1 January 2009 to 31 December 2020. Finally, the GPS PW mean values for each day and each station were determined. A total of 437,842 data points were received. For the analysis of long-term changes, only stations with at least 85% nominal observations were selected, and the data gaps were not longer than half a year ( Based on the hourly values of the estimated ZWD, the PW values were calculated. Meteorological parameters were taken from the work [11], which are, respectively, The value of the mean temperature (Tm) was determined according to the formula (3) based on the temperature from the GNSS station (Ts). Only fifteen stations are equipped with meteo sensors. Therefore, for all stations Ts was taken from the ERA5 [32].
The obtained GPS PW values were averaged to daily and monthly intervals. The study focused on monthly, annual, and long-term variations in PW. Although they were relatively high with rapidly passing atmospheric fronts, the daily GPS PW fluctuations were not considered in the context of long-term analyses. Values determined as outliers according to the 3σ criterion were removed. Prepared in this way, two sets of time series were used to further analyse data, preceded by homogenisation. If the data gaps were too large, the period was excluded from the analysis. Analysis was limited to a full year from 1 January 2009 to 31 December 2020. Finally, the GPS PW mean values for each day and each station were determined. A total of 437,842 data points were received. For the analysis of long-term changes, only stations with at least 85% nominal observations were selected, and the data gaps were not longer than half a year ( Figure 2a). The values of the linear trend, as well as the annual and semi-annual signal amplitude, were modelled (

Results
This section presents the results of the analyses, which were divided into three stages. In first subsection, we present the overall changes in the estimated values of GPS precipitable water. We then focused on the results from the LS analysis. Finally, we compared the estimated GPS PW with the available radio sounding (RS) data.

Results
This section presents the results of the analyses, which were divided into three stages. In first subsection, we present the overall changes in the estimated values of GPS precipitable water. We then focused on the results from the LS analysis. Finally, we compared the estimated GPS PW with the available radio sounding (RS) data.

General Statistics on the GPS PW Changes
The mean value of ZWD formal error is 7.58 mm, which translates into a mean GPS PW error of 0.91 mm. The estimated hourly GPS PW values exceed 40 mm in the summer period. For full hydrological annual cycles (2009-2020), the mean value of precipitable water for each station was calculated. The total mean equals 15.05 mm and ranges for individual stations in Poland from 14.27 mm (Koscierzyna, Pomerania) to 16.04 mm (Wroclaw, Lower Silesia). The distribution of the multi-annual mean values of GPS PW is presented in Figure 3a.

General Statistics on the GPS PW Changes
The mean value of ZWD formal error is 7.58 mm, which translates into a mean GPS PW error of 0.91 mm. The estimated hourly GPS PW values exceed 40 mm in the summer period. For full hydrological annual cycles (2009-2020), the mean value of precipitable water for each station was calculated. The total mean equals 15.05 mm and ranges for individual stations in Poland from 14.27 mm (Koscierzyna, Pomerania) to 16.04 mm (Wroclaw, Lower Silesia). The distribution of the multi-annual mean values of GPS PW is presented in Figure 3a.  (Table 2). Furthermore, almost all extreme GPS PW values (over 40 mm) were recorded in July (e.g., 44.11 mm on 19 July 2018 or 43.81 mm on 29 July 2010). Maps of the monthly PW means (Figure 4) also show the climatic influence of the continent; in winter, the eastern part of the country shows lower values of GPS PW, while in the summer the eastern part shows higher values than the western one.   (Table 2). Furthermore, almost all extreme GPS PW values (over 40 mm) were recorded in July (e.g., 44.11 mm on 19 July 2018 or 43.81 mm on 29 July 2010). Maps of the monthly PW means (Figure 4) also show the climatic influence of the continent; in winter, the eastern part of the country shows lower values of GPS PW, while in the summer the eastern part shows higher values than the western one.

Longterm Variation of PW
LS analysis was conducted by identifying shifts in the GPS PW time series occurring upon changing antennas. In the analysed set of cases, the changes caused a shift in the GPS PW series by 0.15 mm on average. Only in 3 out of 72 affected stations did it exceed 0.5 mm. The low value of visible shifts may be due to the individually calibrated antennas mounted at the ASG-EUPOS stations. The shifts' impact translates to a maximum of 0.01 mm/year, which could be considered insignificant. The noise level of the residual time series is about 4-5 mm (Figures 2c and 3b). Even with the modelling of higher harmonics for 1/3 and 1/4 of the year, the standard deviation of the residues is still over 4 mm.
It should be noted here that any erroneous estimate of the discontinuity may affect the value of the estimated trend. In the analysed 12-year period, the maximum error of trend determination is about 8% of the mis-modelled shifts' value. However, this applies to exceptional cases where the shifts occur in the middle of the analysed period. In our case, the error of trend determination caused by discontinuities does not exceed 0.01 mm/year. The formal error of the estimated trends is usually between 0.01 and 0.02 mm/year, which only confirms the fact that, in the analysed case, the antenna changes did not cause disturbances that could significantly affect the obtained results. However, this issue should be kept in mind when interpreting the results.
Throughout the entire country, the trend is positive and ranges from 0.03 to 0.13 mm/year. Given the above considerations relating to their error and shift estimation, one should acknowledge their values as reliable. The apparent systematic increase in the value of GPS PW (Figure 5c) is confirmed by the well-known fact that the PW increases with an increase in temperature, which has been reported in Poland since 1988 [33].   (Figures 2c and 3b). Even with the modelling of higher harmonics for 1/3 and 1/4 of the year, the standard deviation of the residues is still over 4 mm.

Longterm Variation of PW
It should be noted here that any erroneous estimate of the discontinuity may affect the value of the estimated trend. In the analysed 12-year period, the maximum error of trend determination is about 8% of the mis-modelled shifts' value. However, this applies to exceptional cases where the shifts occur in the middle of the analysed period. In our case, the error of trend determination caused by discontinuities does not exceed 0.01 mm/year. The formal error of the estimated trends is usually between 0.01 and 0.02 mm/year, which only confirms the fact that, in the analysed case, the antenna changes did not cause disturbances that could significantly affect the obtained results. However, this issue should be kept in mind when interpreting the results.
Throughout the entire country, the trend is positive and ranges from 0.03 to 0.13 mm/year. Given the above considerations relating to their error and shift estimation, one should acknowledge their values as reliable. The apparent systematic increase in the value of GPS PW (Figure 5c) is confirmed by the well-known fact that the PW increases with an increase in temperature, which has been reported in Poland since 1988 [33]. The illustrations showing the amplitudes of annual and semi-annual oscillations (Figure 5a,b), reflect the influences of the continent and the ocean on the climate of Poland. The latitudinal system of mountains and lowland areas, as well as the dominant zonal direction of the movement of air masses from west to east, facilitates the advection of humid air from the Atlantic Ocean through Western Europe over Poland, especially in the western part. Towards the east, the influence of the oceanic air masses gradually diminishes. The nature of isoamplitudes corresponds to the course of the average annual amplitudes of air temperature and the spatial distribution of thermal continentalism indicators (according to Chromow, Ewert, Conrad, and Johansson-Ringleb) [34]. Higher values of the amplitude of annual and semi-annual oscillations indicate that eastern Poland is characterised by greater annual contrasts in precipitation as well as in other climate elements (e.g., air temperature, snow cover duration, etc.) [8,9,34,35].

Comparison with the RS Data
The obtained results were validated by comparison with the available radio sounding (RS) data. Within the study area, we obtained data from three sounding points located in Poland (Figure 1), for which sufficiently long observational data were available for three sounding points from the National Oceanic and Atmospheric Administration, Earth System Research Laboratories (ESRL), and Radiosonde Database [36]. The values of RS PW were calculated according to [37] as: where |g| = 9.8 ms −2 is the magnitude of gravitational acceleration; ρ = 1 is the liquid water density; PB and PT are the air pressures in the bottom and top layer, respectively; rT is a column-average of the total water mixing ratio. These results were compared with PW data from the nearest GPS station and are summarised in Table 3. The obtained differences (GPS minus RS) for the analysed four pairs range from −1.21 mm to −0.08 mm, accounting for 0.5% to 8% of the GPS PW values. The highest differences were recorded for the station in Wroclaw. For all pairs, higher amplitudes of annual and semi-annual signals in GPS PW were observed, which may be due to additional artefacts related to GPS observations. The illustrations showing the amplitudes of annual and semi-annual oscillations (Figure 5a,b), reflect the influences of the continent and the ocean on the climate of Poland. The latitudinal system of mountains and lowland areas, as well as the dominant zonal direction of the movement of air masses from west to east, facilitates the advection of humid air from the Atlantic Ocean through Western Europe over Poland, especially in the western part. Towards the east, the influence of the oceanic air masses gradually diminishes. The nature of isoamplitudes corresponds to the course of the average annual amplitudes of air temperature and the spatial distribution of thermal continentalism indicators (according to Chromow, Ewert, Conrad, and Johansson-Ringleb) [34]. Higher values of the amplitude of annual and semi-annual oscillations indicate that eastern Poland is characterised by greater annual contrasts in precipitation as well as in other climate elements (e.g., air temperature, snow cover duration, etc.) [8,9,34,35].

Comparison with the RS Data
The obtained results were validated by comparison with the available radio sounding (RS) data. Within the study area, we obtained data from three sounding points located in Poland (Figure 1), for which sufficiently long observational data were available for three sounding points from the National Oceanic and Atmospheric Administration, Earth System Research Laboratories (ESRL), and Radiosonde Database [36]. The values of RS PW were calculated according to [37] as: where |g| = 9.8 ms −2 is the magnitude of gravitational acceleration; ρ = 1 is the liquid water density; P B and P T are the air pressures in the bottom and top layer, respectively; r T is a column-average of the total water mixing ratio. These results were compared with PW data from the nearest GPS station and are summarised in Table 3. The obtained differences (GPS minus RS) for the analysed four pairs range from −1.21 mm to −0.08 mm, accounting for 0.5% to 8% of the GPS PW values. The highest differences were recorded for the station in Wroclaw. For all pairs, higher amplitudes of annual and semi-annual signals in GPS PW were observed, which may be due to additional artefacts related to GPS observations.

Discussion
Poland is situated in moderate latitudes between 49 and 55 north parallels (Figure 1). According to the Köppen-Geiger classification [38], it is located on the border of two climatic zones. The classification is based on average monthly temperatures and the amount and distribution of annual precipitation in relation to latitude. The first one, Dfb, belongs to the continental zone. It has snowy climates and covers the eastern part of the country. The second one is Cfb-oceanic climate, mild, with no dry season. Both zones are quite humid and characterised by warm summers. The most important factors shaping the climate in Poland include: latitude, the influence of the Baltic Sea and the Atlantic Ocean, and the layout of lowlands and mountainous areas [39] (Figure 6a). The impacts of oceanic and continental climate factors are visible in the spatial distribution of climate elements ( Figure 6). The precipitable water over the territory of Poland is diversified both spatially and temporally. The developed GPS PW distributions (Figures 3 and 4) confirm the quoted regularities.

Discussion
Poland is situated in moderate latitudes between 49 and 55 north parallels (Figure 1). According to the Köppen-Geiger classification [38], it is located on the border of two climatic zones. The classification is based on average monthly temperatures and the amount and distribution of annual precipitation in relation to latitude. The first one, Dfb, belongs to the continental zone. It has snowy climates and covers the eastern part of the country. The second one is Cfb-oceanic climate, mild, with no dry season. Both zones are quite humid and characterised by warm summers. The most important factors shaping the climate in Poland include: latitude, the influence of the Baltic Sea and the Atlantic Ocean, and the layout of lowlands and mountainous areas [39] (Figure 6a). The impacts of oceanic and continental climate factors are visible in the spatial distribution of climate elements ( Figure 6). The precipitable water over the territory of Poland is diversified both spatially and temporally. The developed GPS PW distributions (Figures 3 and 4) confirm the quoted regularities. The presented mean values of the GPS PW illustrate that the eastern part of Poland is under great influence from the continental climate. In the cool season (Figure 6a), it is characterised by a lower GPS PW than the western regions under the greater influence of the Atlantic Ocean [41]. Higher values of PW in the summer season (Figure 6b) in the south-eastern part of the country are due to more intense convective processes-typical for areas with continental features (e.g., this region has one of the largest numbers of days The presented mean values of the GPS PW illustrate that the eastern part of Poland is under great influence from the continental climate. In the cool season (Figure 6a), it is characterised by a lower GPS PW than the western regions under the greater influence of the Atlantic Ocean [41]. Higher values of PW in the summer season (Figure 6b) in the south-eastern part of the country are due to more intense convective processes-typical for areas with continental features (e.g., this region has one of the largest numbers of days with storms during the year at over 30) and the influence of seasonal pressure changes (over the area of Eastern Europe) [9,35,41]. This is even more visible on the maps of the annual cycle of the GPS PW (Figure 5a) or its overall variability in the analysed 12 years (Figure 3b), where the isolines coincide with the boundary of the continent's influence [34]. The above regularities refer to the course of latitudinal gradients-i.e., water vapour pressure and air temperature. The former is higher by about 0.1 hPa/1 • longitude in the western part of Poland. The latter increases in summer by 0.2 • C for each degree of longitude, while in winter it decreases by 0.3 • C [8].
The analysis of the GPS PW value in the studied period showed a clear upward trend in the PW value. This is also confirmed by the LS analysis conducted (see Section 3.2). A positive trend was obtained for the entire region (Figure 5c). The average trend for the Polish region is 0.078 mm/year. The obtained regularity confirms the studies conducted earlier [25,42], where positive trends were obtained for selected Polish EPN stations. A direct comparison with those studies is not possible due to the different observation periods. The length of the analysed period is also important here, which, according to [23], affects the value of the estimated trends. The lack of full consistency in the obtained results may be affected by the fact that, out of the 12 years analysed in this work, nine are among the warmest in the period 1880 to2020 [43]. The spatial distribution of trends itself is difficult to link directly with climatic zones in Poland. Only the larger trends in the coastal belt are characteristic (potentially the influence of the Baltic and the Atlantic Ocean). Higher values of trends were also obtained for stations near Warsaw (BOGO = 0.138 mm/year, BOGO = 0.124 mm/year, JOZE = 0.123 mm/year, JOZ2 = 0.107 mm/year) as well as for Krakow (KRA1 = 0.132 mm/year; KRAW = 0.109 mm/year), which may suggest the influence of an urban heat island [44] on the distribution and temporal changes in PW (UHI). Smooth spatial changes in the annual and semi-annual signal amplitude confirm the influence of the continental climate. The values of the obtained amplitudes are generally higher than the data from the radio soundings. Nevertheless, they are consistent with the values found in other works [25,42]. The annual signal reaches a maximum between July 18 (south-eastern Poland) and July 25 (north-western Poland). This suggests the greater influence of the ocean than dependence on latitude.
Long-term GPS observations also allow the identification of the extreme and anomalous seasons. Let us take an example of the year 2010. According to the bulletins of the Polish Institute of Meteorology and the Water Management National Research Institute, it was wet and cold [33]. This can also be seen in the graph of GPS PW annual mean (Figure 7a). The year 2010, similarly to 2014, shows a positive anomaly. The maximum PW values recorded in 2010 particularly show the scale of the phenomenon (Figure 7b). At the same time, October 2010 was extremely dry and cool [45]. The mean value of GPS PW in October 2010 was 3.87 mm (Figure 8), which is lower than the 12-year average for this period.
Remote Sens. 2021, 13, 2960 9 of 13 with storms during the year at over 30) and the influence of seasonal pressure changes (over the area of Eastern Europe) [9,35,41]. This is even more visible on the maps of the annual cycle of the GPS PW (Figure 5a) or its overall variability in the analysed 12 years (Figure 3b), where the isolines coincide with the boundary of the continent's influence [34]. The above regularities refer to the course of latitudinal gradients-i.e., water vapour pressure and air temperature. The former is higher by about 0.1 hPa/1° longitude in the western part of Poland. The latter increases in summer by 0.2 °C for each degree of longitude, while in winter it decreases by 0.3 °C [8].
The analysis of the GPS PW value in the studied period showed a clear upward trend in the PW value. This is also confirmed by the LS analysis conducted (see Section 3.2). A positive trend was obtained for the entire region (Figure 5c). The average trend for the Polish region is 0.078 mm/year. The obtained regularity confirms the studies conducted earlier [25,42], where positive trends were obtained for selected Polish EPN stations. A direct comparison with those studies is not possible due to the different observation periods. The length of the analysed period is also important here, which, according to [23], affects the value of the estimated trends. The lack of full consistency in the obtained results may be affected by the fact that, out of the 12 years analysed in this work, nine are among the warmest in the period 1880 to2020 [43]. The spatial distribution of trends itself is difficult to link directly with climatic zones in Poland. Only the larger trends in the coastal belt are characteristic (potentially the influence of the Baltic and the Atlantic Ocean). Higher values of trends were also obtained for stations near Warsaw (BOGO = 0.138 mm/year, BOGO = 0.124 mm/year, JOZE = 0.123 mm/year, JOZ2 = 0.107 mm/year) as well as for Krakow (KRA1 = 0.132 mm/year; KRAW = 0.109 mm/year), which may suggest the influence of an urban heat island [44] on the distribution and temporal changes in PW (UHI). Smooth spatial changes in the annual and semi-annual signal amplitude confirm the influence of the continental climate. The values of the obtained amplitudes are generally higher than the data from the radio soundings. Nevertheless, they are consistent with the values found in other works [25,42]. The annual signal reaches a maximum between July 18 (south-eastern Poland) and July 25 (north-western Poland). This suggests the greater influence of the ocean than dependence on latitude.
Long-term GPS observations also allow the identification of the extreme and anomalous seasons. Let us take an example of the year 2010. According to the bulletins of the Polish Institute of Meteorology and the Water Management National Research Institute, it was wet and cold [33]. This can also be seen in the graph of GPS PW annual mean ( Figure  7a). The year 2010, similarly to 2014, shows a positive anomaly. The maximum PW values recorded in 2010 particularly show the scale of the phenomenon (Figure 7b). At the same time, October 2010 was extremely dry and cool [45]. The mean value of GPS PW in October 2010 was 3.87 mm (Figure 8), which is lower than the 12-year average for this period.   Our study shows the distribution of precipitable water in the atmosphere over Poland and its changes in recent years. The maps already presented show a spatial variation, which cannot be obtained without sufficiently dense data. Although it is possible to use a long observation period using EPN stations, there are only several such stations [29] in Poland, and, as shown in works [23][24][25]42], only a few of them have observations for long periods. The data from either or both sets of stations-i.e., full ( Figure 9a) and EPN only (Figure 9b)-indicate a minimum value of PW in the Suwalki Lake District (the northeastern edge of Poland) and an increase towards the south-east. However, denser data provide more detailed information about regional extremes. The reduced number of stations resulted in a different distribution of PW in the country, including no clear influence of the sea and ocean in the north-western part and a much less marked continental influence in the central-eastern part of Poland (Figure 9).

Conclusions
Determining the spatial and temporal variability of precipitable water is challenging. Data obtained through meteorological modelling depend on the model, while radiosonde measurement data are sparse and not evenly distributed. A convenient alternative is GNSS-based precipitable water data. The detailed spatial distribution of the PW is particularly beneficial in areas with many different climatic influences (e.g., continental, maritime, mountainous, oceanic), such as Poland. This can be supplied with the use of GNSS data from EPN and ASG-EUPOS stations. Studies have shown that by reducing the period of observation we can use extra data and significantly increase the spatial resolution. Our study shows the distribution of precipitable water in the atmosphere over Poland and its changes in recent years. The maps already presented show a spatial variation, which cannot be obtained without sufficiently dense data. Although it is possible to use a long observation period using EPN stations, there are only several such stations [29] in Poland, and, as shown in works [23][24][25]42], only a few of them have observations for long periods. The data from either or both sets of stations-i.e., full ( Figure 9a) and EPN only (Figure 9b)indicate a minimum value of PW in the Suwalki Lake District (the north-eastern edge of Poland) and an increase towards the south-east. However, denser data provide more detailed information about regional extremes. The reduced number of stations resulted in a different distribution of PW in the country, including no clear influence of the sea and ocean in the north-western part and a much less marked continental influence in the central-eastern part of Poland (Figure 9). Our study shows the distribution of precipitable water in the atmosphere over Poland and its changes in recent years. The maps already presented show a spatial variation, which cannot be obtained without sufficiently dense data. Although it is possible to use a long observation period using EPN stations, there are only several such stations [29] in Poland, and, as shown in works [23][24][25]42], only a few of them have observations for long periods. The data from either or both sets of stations-i.e., full ( Figure 9a) and EPN only (Figure 9b)-indicate a minimum value of PW in the Suwalki Lake District (the northeastern edge of Poland) and an increase towards the south-east. However, denser data provide more detailed information about regional extremes. The reduced number of stations resulted in a different distribution of PW in the country, including no clear influence of the sea and ocean in the north-western part and a much less marked continental influence in the central-eastern part of Poland ( Figure 9).

Conclusions
Determining the spatial and temporal variability of precipitable water is challenging. Data obtained through meteorological modelling depend on the model, while radiosonde measurement data are sparse and not evenly distributed. A convenient alternative is GNSS-based precipitable water data. The detailed spatial distribution of the PW is particularly beneficial in areas with many different climatic influences (e.g., continental, maritime, mountainous, oceanic), such as Poland. This can be supplied with the use of GNSS data from EPN and ASG-EUPOS stations. Studies have shown that by reducing the period of observation we can use extra data and significantly increase the spatial resolution.

Conclusions
Determining the spatial and temporal variability of precipitable water is challenging. Data obtained through meteorological modelling depend on the model, while radiosonde measurement data are sparse and not evenly distributed. A convenient alternative is GNSSbased precipitable water data. The detailed spatial distribution of the PW is particularly beneficial in areas with many different climatic influences (e.g., continental, maritime, mountainous, oceanic), such as Poland. This can be supplied with the use of GNSS data from EPN and ASG-EUPOS stations. Studies have shown that by reducing the period of observation we can use extra data and significantly increase the spatial resolution.
The use of GPS data confirmed that PW values are systematically increasing for the entire territory of Poland. The spatial distribution of the mean PW values and their fluctuations over the years is congruent with the fact that Poland lies on the border of the influence of both continental and oceanic climates. The obtained values of PW, as well as their spatio-temporal distribution, are consistent with the climate research concerning this region. This confirms that GPS data can contribute to meteorological research, including climatic studies. This work also demonstrates that the development of geodetic infrastructure brings wide benefits to earth sciences.