Interannual Variability of the GNSS Precipitable Water Vapor in the Global Tropics

This paper addresses the subject of inter-annual variability of the tropical precipitable water vapor (PWV) derived from 18 years of global navigation satellite system (GNSS) observations. Non-linear trends of retrieved GNSS PWV were investigated using the singular spectrum analysis (SSA) along with various climate indices. For most of the analyzed stations (~49%) the GNSS PWV anomaly was related to the El Niño Southern Oscillation (ENSO), although its influence on the PWV variability was not homogeneous. The cross-correlations coefficient values estimated between the Multivariate ENSO Index (MEI) and PWV were up to 0.78. A strong cross-correlation was also found for regional climate pattern expressed through CAR, DMI, HAW, NPGO, TNA and TSA indices. A distinct agreement was also found when instead of climate indices, the local sea surface temperature was examined (average correlation 0.60). The SSA method made it also possible to distinguish small-scale phenomena that affect PWV, such as local droughts or wetter rainy seasons. The overall nature of the investigated changes was also verified through linear trend analysis. In general, not a single station was characterized by a negative trend and its weighted mean value, calculated for all stations was equal to 0.08 ± 0.01 mm/year.


Introduction
The tropical region is a major center of atmospheric convection, accompanied by significant latent heat release [1]. The Earth's climate system's global circulation depends on latent heat flux from the surface to the atmosphere through the chain of evaporationclouds-precipitation processes [2], while the transport of heat and momentum by Hadley Cells is vital to the maintenance of the global heat balance. Since precipitation rates are the largest in the global tropics, so is the amount of latent heat transferred to the atmosphere. Thus, variability in the tropical water vapor impacts global circulation patterns. As a result, tropical weather patterns across numerous timescales have a great repercussions on the weather and climate in mid-latitudes through atmospheric teleconnections [3,4]. The equatorial humidity changes resulting from, e.g., El Niño-Southern Oscillation (ENSO), which is the most prominent year-to-year climate fluctuation on Earth [4], has a strong impact on tropical weather and climate [5][6][7]. However, footprints of ENSO are also observed nearly world-wide, including both poles, in the Earth's climate system [8][9][10]. Therefore, long-term, reliable monitoring of water vapor variability in the tropics is an important task in understanding global climate, especially considering the fact that water vapor is characterized by the single highest positive feedback on the surface temperature [11][12][13][14] and dominates the effect of the Earth's surface temperature increase [15].
The amount of water vapor in the atmospheric column, usually expressed as integrated water vapor (IWV) or precipitable water vapor (PWV), can be measured using various techniques. Traditionally, meteorological balloons that perform high-resolution profiles of atmospheric temperature, humidity, pressure and horizontal winds up to 30 km height (radio sounding, RS) are used for this purpose. Although this is one of the most commonly used method for PWV measurements, RS data suffer from certain limitations. For one, RS profiles are not purely vertical but may suffer from a drift due to horizontal winds [16,17], resulting in ascent through air masses characterized by different moisture content. Furthermore, the spatial distribution of the RS measurement network is rather coarse (from~200 km up to~700 km, depending on the location), highly inhomogeneous and land-locked, with only a handful of research vessels providing data over open ocean. Operational temporal resolution of RS launches is daily or twice-daily [18], which is more than sufficient for monitoring of climate phenomena, but limits the ability to investigate scale interactions, especially in the tropics where the primary mode of variability is diurnal cycle [19,20]. Balloon launch frequency can be increased over certain areas for a limited period for a targeted observations, but coordinated enhancement over large regions (e.g., the entire maritime continent) requires tremendous coordinated effort on an international level [21]. Finally, RS operations require ground staff during balloon launch and their cost itself poses strong limitations for increased spatial and temporal resolution of these data. Such limitations are even more important for the research of unpredictable climate patterns, such as ENSO. For example, the Tropical Ocean-Global Atmosphere (TOGA) Coupled Ocean Atmosphere Response Experiment (COARE) intensive observation period (IOP) from November 1992 to March 1993 has been designed to study atmosphere-ocean interactions, organization of atmospheric convection, ocean variability and multi-scale interactions within the western Pacific warm pool in order to gain better understanding of ENSO, among others [22]. Although a tremendous amount of research alongside the project's goals has been delivered as a result of this endeavor, the IOP period has not encountered clear ENSO conditions, which means that the resources deployed as part of TOGA-COARE provided little insight into ENSO dynamics [23,24]. Recently, the National Oceanic and Atmospheric Administration (NOAA) El Niño Rapid Response (ENRR) [25] project has been carried out during a major El Niño event in 2016, with a ship, aircrafts and land stations deployed for targeted observations. This project aimed to improve the understanding of the atmospheric response to El Niño conditions and its implications on heavy weather conditions in the US. The success of NOAA ENRR confirms the fundamental need for continuous research aimed at an in-depth understanding of the dynamics of tropical climate patterns and their teleconnections. However, its relatively short period of 12 weeks did not cover many aspects of the 2016 event, such as the initial evolution of the atmospheric moisture pattern due to shifting SSTs. Therefore, a reliable and continuous monitoring of atmospheric moisture content across the Pacific Ocean would greatly augment capacity of targeted efforts such as NOAA ENRR.
In recent years, a new technique used for the purpose of short-and long-term PWV variability analysis has been developed. It is based on the global navigation satellite systems (GNSS) ground and satellite infrastructure. Since water vapor affects the propagation of the GNSS signal through changing atmospheric refractivity, each of the GNSS permanent stations can provide long term series of unique data of atmospheric refraction. Firstly, Bevis et al. [26] proposed the concept of GNSS meteorology, by calculating the PWV based on GNSS observations together with surface temperature and pressure. Since then, many studies have proved the great utility of this technique in monitoring troposphere humidity changes [27][28][29][30][31].
The great potential of the GNSS technique in the monitoring of PWV variability emerges from the high temporal and spatial resolution of the observations, which are collected regardless of weather conditions. In addition, stations belonging to the GNSS permanent network, such as the International GNSS Service (IGS), have been performing measurements even since 1994, thus providing an over 20-year time series. Another advantage is the automatic data collection and cost-effectiveness of GNSS data acquisition systems (e.g., single unit measures for many months as opposed to disposable RS Atmosphere 2021, 12, 1698 3 of 23 units). Although it is not a direct method of PWV measurements, the application of an appropriate processing strategy for GNSS observations enables the use of GNSS PWV time series as a valuable source for PWV long term analysis, with accuracy similar to the PWV measured by radio soundings [32,33]. Studies conducted so far, on the global [34][35][36][37][38] and regional [32,[39][40][41] scales have clearly proven the efficacy of GNSS PWV data in climate monitoring.
This paper attempts to establish a link between GNSS PWV decadal variability and global, regional, and local climate phenomenon, which are main drivers of moisture interannual changes. It focuses on long-term variability in PWV over GNSS stations located in the tropical regions, both from a linear and non-linear point of view. Special attention has been paid to the ENSO to demonstrate the possibilities of using the GNSS infrastructure to monitor this phenomenon as well as its impact on weather conditions. The structure of the paper is as follows: section two contains information about the dataset used for this study and describes the methodology; section three presents results and discussion. The last section is devoted to conclusions.

Materials and Methods
GNSS signals passing through the Earth's troposphere are disturbed by the temperature, air pressure, and water vapor, which cause their delay. In general, this delay is expressed towards the zenith and is called zenith tropospheric delay (ZTD). It is estimated through GNSS observation processing and consists of a delay caused by the hydrostatic (zenith hydrostatic delay, ZHD) and wet (zenith wet delay, ZWD) components of the atmosphere. In practice, a priori ZHD is used in the process of ZTD estimation, and correction calculated to its value is often equated to the ZWD. However, the ZWD obtained in such a way includes both the ZWD itself and correction for formerly modelled ZHD. Therefore, in GNSS meteorology, ZHD is also calculated a posteriori, using, e.g., the Saastamoinen formula [42], and in such a form is subtracted from GNSS ZTD. The resulting ZWD is strictly related to the humidity in the atmosphere and can serve IWV estimation using the Bevis formula [26]. The water vapor weighted mean temperature, which is necessary for this procedure, can be estimated through, e.g., its linear relationship with the surface temperature [43,44].
In this study, we focus on the inter-seasonal variability of the PWV, estimated based on the GNSS observations derived from the selected stations belonging to the IGS network [45]. It should be noted that station selection aimed data consistency among analyzed locations. Therefore, selected stations belong to only one of many global, pan-regional and regional networks and should be considered a subset of all GNSS data. The selected stations are located between 30 • S and 30 • N and have been operating between January 2001 and December 2018, and were selected to ensure long-term coherent data sets available for analyses. The collected 30-s GPS observations were processed in Bernese GNSS Software ver. 5.2 [46], using the precise point positioning (PPP) approach. We decided to use this method, since it prevents error propagation between stations and provides one of the most reliable results from the point of view of GNSS climate monitoring [33]. Precise ephemeris and clocks from second CODE reprocessing [47] were used. In addition, observations below 5 • elevation mask were excluded, to limit the negative effects of phenomena such as multipath. Hourly ZTD values were estimated together with tropospheric gradients using Vienna Mapping Function 1 [48] and Chen and Herring mapping function [49]. Estimation of the ZWD were possible by way of the Saastamoinen model and air pressure derived from the ERA5 model [50]. Conversion between ZWD and PWV was performed using water vapor weighted mean temperature estimated based on surface temperature data from ERA5. This enables the estimation of long-term PWV with mean bias below 2 mm and linear trend difference below 0.01 mm/year in comparison to the RS data [33]. Based on the so-processed PWV we estimated daily solutions. Stations with less than 70% of the available data were excluded from further analysis. As a result, the average completeness of the solutions at selected stations was at a level of 85%. The location of the Atmosphere 2021, 12, 1698 4 of 23 chosen 43 stations, together with information about their height above sea level, is given in Figure 1 and Table A1 (Appendix A). mm/year in comparison to the RS data [33]. Based on the so-processed PWV we estimated daily solutions. Stations with less than 70% of the available data were excluded from further analysis. As a result, the average completeness of the solutions at selected stations was at a level of 85%. The location of the chosen 43 stations, together with information about their height above sea level, is given in Figure 1 and Table A1 (Appendix  A). A PWV analysis was supported by SST data from the Fleet Numerical Meteorology and Oceanography Center (FNMOC) High-Resolution SST/Sea Ice Analysis for the Global High-Resolution SST project [51]. These data are available every 6 h at a near global Mercator projection grid, with 12 km resolution at the equator. Daily anomalies of SST were computed by removing the seasonal cycle, defined as the mean plus the first three annual harmonics fitted to an average annual cycle calculated between 1 January 2006-31 December 2018.
The singular spectrum analysis (SSA) method [52], [53] was used to extract the nonlinear trend, which we assume is the slowly-varying version of a time series that comprises long-term variations but without seasonal periodicities. SSA is a non-parametric and model-free method which can be applied to analysis and prediction of any series. Its advantage over classical spectral methods and its importance in non-linear dynamics results from the adaptive character of eigenvalues. In non-linear dynamics, the SSA method was introduced by Broomhead and King [54], and is based on principal component analysis (PCA). Its particular importance in the time series analysis results from the fact that the power spectrum of the analyzed signal is the sum of spectra of individual components. In other words, the time series can be represented as a sum of different components such as trend (defined as a slowly varying series), periodicities, and noise. Therefore, the SSA enables the decomposition of the signal and obtainment of useful information about a specific component. Thus, SSA has applications in many areas, including smoothing, noise reduction, non-linear trend estimation, extraction of periodicities, estimation of volatility, etc. Because the SSA method requires continuous time series, we filled the gaps considering the temporal correlation present in the data and the iterative SSA described by Kondrashov and Ghil [55]. Then, homogenization of PWV time series was performed using an approach presented by Hoseini et al. [56]. However, instead of using ERA-Interim as in the original, we used higher-resolution data from the newest ERA5 reanalysis [50]. After finding all the change points and correcting the mean shifts, we extracted the nonlinear trend as one of the reconstructed components of the original time series. Figure 2 presents the prepared PWV time series (top) for GNSS station GUAM (red) and the first five reconstructed components (RC) derived from the SSA A PWV analysis was supported by SST data from the Fleet Numerical Meteorology and Oceanography Center (FNMOC) High-Resolution SST/Sea Ice Analysis for the Global High-Resolution SST project [51]. These data are available every 6 h at a near global Mercator projection grid, with 12 km resolution at the equator. Daily anomalies of SST were computed by removing the seasonal cycle, defined as the mean plus the first three annual harmonics fitted to an average annual cycle calculated between 1 January 2006-31 December 2018.
The singular spectrum analysis (SSA) method [52,53] was used to extract the nonlinear trend, which we assume is the slowly-varying version of a time series that comprises longterm variations but without seasonal periodicities. SSA is a non-parametric and model-free method which can be applied to analysis and prediction of any series. Its advantage over classical spectral methods and its importance in non-linear dynamics results from the adaptive character of eigenvalues. In non-linear dynamics, the SSA method was introduced by Broomhead and King [54], and is based on principal component analysis (PCA). Its particular importance in the time series analysis results from the fact that the power spectrum of the analyzed signal is the sum of spectra of individual components. In other words, the time series can be represented as a sum of different components such as trend (defined as a slowly varying series), periodicities, and noise. Therefore, the SSA enables the decomposition of the signal and obtainment of useful information about a specific component. Thus, SSA has applications in many areas, including smoothing, noise reduction, non-linear trend estimation, extraction of periodicities, estimation of volatility, etc. Because the SSA method requires continuous time series, we filled the gaps considering the temporal correlation present in the data and the iterative SSA described by Kondrashov and Ghil [55]. Then, homogenization of PWV time series was performed using an approach presented by Hoseini et al. [56]. However, instead of using ERA-Interim as in the original, we used higher-resolution data from the newest ERA5 reanalysis [50]. After finding all the change points and correcting the mean shifts, we extracted the nonlinear trend as one of the reconstructed components of the original time series. Figure 2 presents the prepared PWV time series (top) for GNSS station GUAM (red) and the first five reconstructed components (RC) derived from the SSA method. RC1 and RC2 include information about annual oscillation, so we show the aggregate value. Similarly, in the case of RC4 and RC5, which present components of semi-annual oscillation. RC3 can be interpreted as a non-linear trend. For comparison purposes, the grey color represents PWV signal and its decomposition based on data from a neighboring RS station. Also, correlation values are presented, which are very high for the RCs shown (above 0.9). This shows that GNSS data are comparable with RS signals and confirms the ability of the GNSS data to provide reliable monitoring of integrated water vapor properties. method. RC1 and RC2 include information about annual oscillation, so we show the aggregate value. Similarly, in the case of RC4 and RC5, which present components of semi-annual oscillation. RC3 can be interpreted as a non-linear trend. For comparison purposes, the grey color represents PWV signal and its decomposition based on data from a neighboring RS station. Also, correlation values are presented, which are very high for the RCs shown (above 0.9). This shows that GNSS data are comparable with RS signals and confirms the ability of the GNSS data to provide reliable monitoring of integrated water vapor properties.

Figure 2.
Example of precipitable water vapor (PWV) time series for GNSS station GUAM (red) and radio sounding station 91212 (red) (top) and first five reconstructed components (RC) obtained by the singular spectrum analysis (SSA) method. Annual and semi-annual oscillations are represented by RC1-2 and RC4-5, respectively. Nonlinear trend is represented by RC3.

Results
This section presents the results and a discussion of the investigation of the PWV's long-term variability. We have divided it into analysis of non-linear and linear changes. In case of non-linear variability, the results are presented and discussed with a division based on main climate indicators that affect inter-annual variability. The similarity between a given index and non-linear PWV variability was assessed through a cross-correlation coefficient. We used this method since it allowed us to obtain the strength of correlation regardless of the time difference between the occurrence of a phenomenon and its potential impact on the local environment. We adopted an absolute value of the cross-correlation coefficient equal to at least 0.4, as an indicator of a statistically significant correlation between a given phenomenon and changes in PWV. The lag window was adopted as 12 months. Some stations show a significant correlation with more than one climate index, which is related to the global teleconnections. As a consequence, a detailed analysis is conducted with comparison to the climate index, that, according to the cross-correlation process, has the biggest impact on the PWV variability. Since global or regional climate indices may not be fully representative for each of the analyzed stations, we also conducted a detailed analysis of the influence of local SST on the PWV variability.

Pacific Related Climate Indices
One of the commonly used ENSO indices, that combines both oceanic and atmospheric variability is the Multivariate ENSO Index Version 2 (hereby called MEI) [57]. Next to the SST, this index includes information about the sea-level pressure (SLP), winds

Results
This section presents the results and a discussion of the investigation of the PWV's long-term variability. We have divided it into analysis of non-linear and linear changes. In case of non-linear variability, the results are presented and discussed with a division based on main climate indicators that affect inter-annual variability. The similarity between a given index and non-linear PWV variability was assessed through a cross-correlation coefficient. We used this method since it allowed us to obtain the strength of correlation regardless of the time difference between the occurrence of a phenomenon and its potential impact on the local environment. We adopted an absolute value of the cross-correlation coefficient equal to at least 0.4, as an indicator of a statistically significant correlation between a given phenomenon and changes in PWV. The lag window was adopted as 12 months. Some stations show a significant correlation with more than one climate index, which is related to the global teleconnections. As a consequence, a detailed analysis is conducted with comparison to the climate index, that, according to the cross-correlation process, has the biggest impact on the PWV variability. Since global or regional climate indices may not be fully representative for each of the analyzed stations, we also conducted a detailed analysis of the influence of local SST on the PWV variability.

Pacific Related Climate Indices
One of the commonly used ENSO indices, that combines both oceanic and atmospheric variability is the Multivariate ENSO Index Version 2 (hereby called MEI) [57]. Next to the SST, this index includes information about the sea-level pressure (SLP), winds and Outgoing Longwave Radiation for the Pacific region between 30 • S-30 • N and 100 • E-70 • W. In principle, positive MEI events (El Niño) are characterized by warmer than usual SSTs accompanied by lower than usual SLP in the central-east tropical Pacific, anonymously high SLP over Indonesia and western equatorial Pacific, reduction of tropical convection over Indonesia and the western Pacific with its reinforcement over the central tropical Pacific, and strong changes of trade winds through their reduction or even reversal. ENSO influences global atmospheric and oceanic circulation. Therefore, one could expect to see its clear impact on each of the analyzed stations located in the tropical belt. However, various PWV variabilities were found and some stations had higher agreement with the ENSO than the others.
Generally, the statistically significant impact (cross-correlation exceeding 0.4) of ENSO on PWV variability was found across the whole of the analyzed area, although its strength and character differed. 21 stations (~49%) were characterized by absolute coefficient of cross-correlation higher than 0.4 (12 stations positive, 9 stations negative) ( Table 1). The highest agreement was found within the Pacific Ocean basin. In this area stations like TUVA (Funafuti, Tuvalu), THTI (Papeete, French Polynesia), TONG (Nuku'alofa, Tonga), PIMO (Quezon City, Philippines) and KOUC (Koumac, New Caledonia) were characterized by a cross-correlation coefficient equal to 0.78, 0.70, −0.66, −0.66 and −0.63, respectively. No lags between ENSO and PWV changes were found at TONG and KOUC stations, which indicate that they varied at the same time as occurrence of El Niño and La Niña episodes. For the stations with a cross-correlation coefficient exceeding 0.5 it can be noticed that estimated lag is within 0 to +5 months. This imply that PWV changes, that have a similar character to changes in the MEI index, occurred together with or after the phenomenon indicated by this index. There were also stations for which statistically significant crosscorrelations were found, but estimated lags were −8 months (KOKB) or even −10 months (KOKV). Such lags would suggest that changes in PWV preceded the ENSO phenomenon, which in these cases is rather a statistical artifact. In addition, PWV over these stations is mostly modulated by other climate patterns (presented later in the article). Positive lags were mostly concentrated in the western Pacific and central Africa, while negative in the Maritime Continent and northern Australia ( Figure 3). Detailed information about the lags and cross-correlation is given in Table 1.
Atmosphere 2021, 12, x FOR PEER REVIEW and Outgoing Longwave Radiation for the Pacific region between 30° S-30° N an E-70° W. In principle, positive MEI events (El Niño) are characterized by warme usual SSTs accompanied by lower than usual SLP in the central-east tropical P anonymously high SLP over Indonesia and western equatorial Pacific, reduct tropical convection over Indonesia and the western Pacific with its reinforcemen the central tropical Pacific, and strong changes of trade winds through their reduc even reversal. ENSO influences global atmospheric and oceanic circulation. The one could expect to see its clear impact on each of the analyzed stations located tropical belt. However, various PWV variabilities were found and some station higher agreement with the ENSO than the others. Generally, the statistically significant impact (cross-correlation exceeding ENSO on PWV variability was found across the whole of the analyzed area, altho strength and character differed. 21 stations (~49%) were characterized by absolute cient of cross-correlation higher than 0.4 (12 stations positive, 9 stations negative) 1). The highest agreement was found within the Pacific Ocean basin. In this area s like TUVA (Funafuti, Tuvalu), THTI (Papeete, French Polynesia), TONG (Nuku Tonga), PIMO (Quezon City, Philippines) and KOUC (Koumac, New Caledonia characterized by a cross-correlation coefficient equal to 0.78, 0.70, −0.66, −0.66 and respectively. No lags between ENSO and PWV changes were found at TONG and stations, which indicate that they varied at the same time as occurrence of El Niño Niña episodes. For the stations with a cross-correlation coefficient exceeding 0.5 it noticed that estimated lag is within 0 to +5 months. This imply that PWV change have a similar character to changes in the MEI index, occurred together with or af phenomenon indicated by this index. There were also stations for which statistica nificant cross-correlations were found, but estimated lags were −8 months (KO even −10 months (KOKV). Such lags would suggest that changes in PWV preced ENSO phenomenon, which in these cases is rather a statistical artifact. In addition over these stations is mostly modulated by other climate patterns (presented later article). Positive lags were mostly concentrated in the western Pacific and central while negative in the Maritime Continent and northern Australia ( Figure 3). D information about the lags and cross-correlation is given in Table 1.   The fact that a PWV non-linear trend was correlated with ENSO did not mean that this was the major phenomenon modulating air humidity. As a consequence, in this subsection we present results for stations, for which ENSO had the biggest impact out of all the analyzed indices. Figure 4a presents stations for which a PWV non-linear trend is positively correlated with the ENSO. These stations are located in the southern-west Pacific, namely TUVA, THTI, FALE (Faleolo, Samoa), ASPA, and in the central-west Pacific, namely PNGM (Lombrum, Papua New Guinea). In their case it can be observed how El Niño is associated with an increase in moisture, and La Niña with a decrease, especially considering the strong La Niña episode of 2010-2012. From the beginning of 2010 (the end of El Niño) to the middle of 2011 (the strongest phase of La Niña), the mean monthly PWV dropped by 4.0 mm, 4.2 mm, 4.4 mm, 4.9 mm and 6.5 mm for FALE, ASPA, THTI, PNGM and TUVA, respectively. Likewise, during strong La Niña in 2007-2009, all stations indicate a decrease in the mean PWV, with the ASPA station being the only exception. Although this station is located about 150 km from the FALE station, it only shows distinct PWV variability during the 2010-2012 La Niña. Similar results for this station were, however, obtained by Wang et al. [58], who also used the SSA method to verify long-term nonlinear changes. As in our case, they did not find a negative PWV anomaly in 2008 during another strong La Niña event, but rather constantly decreasing values from a positive anomaly in 2006 (indicated in both studies) to a negative anomaly in 2011/2012 (also in both studies). A higher disagreement between PWV anomalies estimated using the SSA method can be found in case of the THTI station. While in previously mentioned study [58], the value of PWV is rather constantly dropping since approximately 2005 to 2011, in the case of our results it can be noticed that the highest PWV value, from the same period of time, occurred in January 2007 (during El Niño), then was clearly negative in August 2008 (La Niña), positive in January 2010 (El Niño) and again negative in March 2011 (La Niña). These differences result probably from the calculation strategy adopted for processing GNSS observations, which has a great impact on GNSS PWV reliability. Nevertheless, a comparative analysis of local SST and PWV anomalies (discussed later in Section 3.1.4) for this station showed that the periods highlighted by us, for which the highest and lowest PWV values were observed, are in strong agreement with the highest and lowest SST anomalies. The SST analysis also confirmed the correctness of the non-linear trend for the ASPA station, for which low PWV variability was related to the small SST anomalies during the analyzed period.  Within the Pacific Ocean basin, we also paid special attention to the Hawaii region, whose climate is under a major influence from the Pacific Ocean surrounding the islands. Figure 5 presents long-term PWV anomaly for five stations located in this region, namely HNLC (Oahu island), KOKB and KOKV (both located near each other at Kauai island), MAUI (Maui island) and MKEA (Hawaii island). As can be noticed, all stations present similar long-term variability and are characterized by occurrence of the highest and the lowest anomalies at the same time, regardless of altitude differences. As an example, the most negative anomaly for the HNLC station (22 m) was found in April 2010, whereas in A higher discrepancy between ENSO and PWV non-linear variability can also be observed for the PNGM station. According to Smith et al. [59] the ENSO impact on Papua New Guinea's rainfall is not obvious and homogeneous due to the high topographical complexity. This region is also under strong influence from the Maritime Continent/Australian monsoon that drives the sub-seasonal atmospheric variability and shows two-way interactions with ENSO [60,61]. As can be noticed in Figure 4 (left), PWV anomalies tend to be negative during La Niña events, while their response to the El Niño events is not obvious. These results are in line with those obtained by Murphy et al. [62], who have analyzed various responses of the Pacific islands' weather and climate to the ENSO state. For the station located in Kavieng (New Ireland province, Papua New Guinea), nearby PNGM, they found the occurrence of dryer than normal conditions during La Niña, but wetter than normal conditions during only specific types of El Niño, which support our results and general concept that reliable estimated GNSS PWV can be successfully used for long-term inter-annual moisture variability studies. Figure 4 (right), in turn, presents stations that, in contrast to the previous cases, are characterized by negative PWV anomaly during El Niño episodes, and positive during La Niña. Since these stations are located at much greater distances from each other than in the case of positively correlated stations, the PWV response to the ENSO was not so homogeneous. Nevertheless, the highest PWV anomalies were, again, observed during the 2010-2012 La Niña. For this period PWV anomaly reached 5.3 mm, 3.7 mm, 3.5 mm, 2.7 mm, 1.7 mm and 1.5 mm for TOW2 (Cape Ferguson, Australia), TONG, KOUC, PIMO (Quezon City, Philippines), LAUT and NTUS (Singapore, Singapore), respectively. Strong, negative correlation to the ENSO phases is especially worth emphasizing in the case of KOUC (−0.63), LAUT (−0.48) and TONG (−0.66) stations, that are located close to the previously discussed ASPA, FALE, THTI and TUVA stations. Although they are also in the southern-west Pacific, PWV non-linear trends represent an opposite response to the ENSO. Similar results were found by Murphy et al. [62] who confirmed analogous differences in November-April rainfall for the stations located on the islands of Tuvalu (in our case TUVA station), Fiji (in our case LAUT station) and Tonga (in our case TONG station). In their study, Tuvalu was characterized by negative SST and rainfall anomalies during La Niña events and positive anomalies during El Niño events, which is in close agreement with the results presented here. Furthermore, stations located on the islands of Fiji and Tonga were characterized by negative rainfall anomalies during El Niño events and positive anomalies during La Niña, which is also in line with our results. According to Murphy et al. similar PWV variations occur at Apia island (FALE station in our case), while according to results obtained by us, the results are opposite. In this case, however, the cross-correlation coefficient obtained by us was the smallest among all the stations presented here (0.44) and may not be representative for seasonal rainfall, which was mainly discussed in their study. Additionally, a lack of GNSS data for specific periods may have negatively affected the accuracy of SSA trend estimation, thereby reducing overall consistency between PWV and MEI (more detail in Section 3.1.4).
Within the Pacific Ocean basin, we also paid special attention to the Hawaii region, whose climate is under a major influence from the Pacific Ocean surrounding the islands. Figure 5 presents long-term PWV anomaly for five stations located in this region, namely HNLC (Oahu island), KOKB and KOKV (both located near each other at Kauai island), MAUI (Maui island) and MKEA (Hawaii island). As can be noticed, all stations present similar long-term variability and are characterized by occurrence of the highest and the lowest anomalies at the same time, regardless of altitude differences. As an example, the most negative anomaly for the HNLC station (22 m) was found in April 2010, whereas in the case of the MKEA station (3754 m) it was May 2010. The highest positive anomaly, in turn, was found in August 2010 and June 2010 for HNLC and MKEA, respectively. Altitude, however, affects the range of PWV variability. The widest differences were observed at HNLC station, and varied between −2.4 mm and 2.4 mm, whereas the lowest were found at MAUI (3062 m) and varied between −1.4 mm and 1.4 mm (MKEA was up to 2.0 mm). Anomalies for KOKB and KOKV stations (1168 m) usually fell between results for HNLC and MAUI. According to Chu and Chen [63], this archipelago tends to be dry during El Niño events and wet during La Niña, although such conditions may not be observed for each ENSO phase. Even though results obtained for these stations show correlation with the ENSO (cross-correlation coefficients given in Table 1, with the exception of HNLC station for which it amounted to 0.38), such a dependency can not be observed easily. Consequently, we used the HAW index (https://psl.noaa.gov/forecasts/sstlim/timeseries/, accessed on 1 March 2021), which describes SST anomalies in the vicinity of Hawaii (22 • N-18 • N, 150 • W-140 • W). As can be noticed in Figure 5, rather negative SST anomalies between 2008 and 2013 occurred with mostly negative PWV anomalies. A positive SST anomaly since 2013, in turn, was accompanied by mostly positive anomalies in PWV signal. A distinct positive anomaly in 2018 most likely reflects the effects of hurricanes and cool season storms which contribute to the above average wetter conditions over most of Hawaii [64], while distinct negative in 2010 reflect moderate up to exceptional drought conditions which afflicted nearly half or more of Hawaii [65]. In general, cross-correlation coefficients calculated using PWV non-linear trend and HAW index were 0.65 (0 months lag), 0.47 (0 months lag), 0.54 (0 months lag), 0.42 (−1 months lag) and 0.56 (with 0 months lag), which is higher than in the case of MEI. Although ENSO is a global climate phenomenon, it is not the main driver of SST variability within Subtropical Pacific. In the case of the Hawaii area, SST is also strongly modulated by North Pacific Subtropical Gyre (NPSG). NPSG is a system of ocean currents rotating clockwise across the North Pacific, which transports water to the Hawaii region from the northeast Pacific. Depend on the speed of these currents, changes the time needed to transport this originally cool water from higher to lower latitudes. As a result, there is more (slower speed) or less (higher speed) time to warm water up before it reached Hawaii, and in consequence higher or lower PWV anomalies. The dependency between North Pacific Gyre Oscillation (NPGO) index [66], that represent changes in NPSG speed via SST and sea surface height anomaly, and PWV variability over Hawaii is presented at Figure 5  sults for HNLC and MAUI. According to Chu and Chen [63], this archipelago tends to be dry during El Niño events and wet during La Niña, although such conditions may not be observed for each ENSO phase. Even though results obtained for these stations show correlation with the ENSO (cross-correlation coefficients given in Table 1, with the exception of HNLC station for which it amounted to 0.38), such a dependency can not be observed easily. Consequently, we used the HAW index (https://psl.noaa.gov/forecasts/sstlim/timeseries/, accessed on 1 March 2021), which describes SST anomalies in the vicinity of Hawaii (22° N-18° N, 150° W-140° W). As can be noticed in Figure 5, rather negative SST anomalies between 2008 and 2013 occurred with mostly negative PWV anomalies. A positive SST anomaly since 2013, in turn, was accompanied by mostly positive anomalies in PWV signal. A distinct positive anomaly in 2018 most likely reflects the effects of hurricanes and cool season storms which contribute to the above average wetter conditions over most of Hawaii [64], while distinct negative in 2010 reflect moderate up to exceptional drought conditions which afflicted nearly half or more of Hawaii [65]. In general, cross-correlation coefficients calculated using PWV non-linear trend and HAW index were 0.65 (0 months lag), 0.47 (0 months lag), 0.54 (0 months lag), 0.42 (−1 months lag) and 0.56 (with 0 months lag), which is higher than in the case of MEI. Although ENSO is a global climate phenomenon, it is not the main driver of SST variability within Subtropical Pacific. In the case of the Hawaii area, SST is also strongly modulated by North Pacific Subtropical Gyre (NPSG). NPSG is a system of ocean currents rotating clockwise across the North Pacific, which transports water to the Hawaii region from the northeast Pacific. Depend on the speed of these currents, changes the time needed to transport this originally cool water from higher to lower latitudes. As a result, there is more (slower speed) or less (higher speed) time to warm water up before it reached Hawaii, and in consequence higher or lower PWV anomalies. The dependency between North Pacific Gyre Oscillation (NPGO) index [66], that represent changes in NPSG speed via SST and sea surface height anomaly, and PWV variability over Hawaii is presented at Figure 5

Indian Ocean Related Climate Indices
Another, large-scale climate phenomenon that has been examined for its impact on the PWV long-term variability is related to the Indian Ocean. The Indian Ocean Dipole (IOD) shows interaction with the ENSO [67][68][69] and presents similar mechanism to the El Niño/La Niña events. A positive phase of IOD takes place when SST in the south-eastern tropical Indian Ocean is cooler than usual, while SST across the Indian Ocean, in the western region, is warmer than usual. The negative IOD phase is characterized by the opposite situation. The IOD phases can be expressed by the dipole mode index (DMI) [70], which refers to the SST anomalies.
Since in contrast to the ENSO, the IOD influence is rather limited to the Indian Ocean basin, occurrence of the patterns associated with it was only investigated for stations in this region. Depending on whether the stations were located on the eastern or western side of the Indian Ocean, we used DMIeast or DMIwest indexes, that are more suitable for these regions. Results show that among all the considered stations, five were clearly affected by IOD. These were DGAR (Diego Garcia Island, United Kingdom), MBAR (Mbarara, Uganda) stations (Figure 6a), BAKO (Cibinong, Indonesia), COCO (Cocos (Keeling) Island, Australia) and DARW (Darwin, Australia) (Figure 6b).
during positive SST anomalies in the western Indian Ocean basin, the westerly winds weaken along the equator, thus allowing warm and humid air masses to shift towards central Africa [71]. A distinct consistency between PWV anomalies and DMI can also be found for DGAR station. Despite a smaller cross-correlation coefficient between the PWV and DMIwest index (0.58 with 0 months lag), IOD is a major factor in modulating the climatic variations of PWV. The highest positive PWV anomaly (2.7 mm) occurred in December 2015, when the SST in the western Indian Ocean basin was also the highest during the analyzed time span. The highest negative (−1.7 mm), in turn, took place in 2011 when the DMIwest index was also negative. Although during the analyzed time span there were also similar or even greater negative variations in this index, the changes in PWV caused by them were less significant. However, since this station is also positively correlated with ENSO (0.56 with +2 moths lag), such an effect may be partly explained by a modulating effect of the extreme 2010-2012 La Niña.  From westerly located stations (Figure 6a), the highest value of the cross-correlation coefficient was found for the MBAR station and DMIwest index (0.67 with 0 months lag). Despite the fact, that this station is located in eastern-central Africa at more than 1340 m.a.s.l, it is clearly affected by the IOD and noted variance simultaneously with its phases. The highest PWV anomaly (1.3 mm) was observed in January 2016, while the second highest (0.9 mm) in March 2016, both together with the highest positive anomalies of the DMIwest index. Such a distinct increase in humidity results from the fact that during positive SST anomalies in the western Indian Ocean basin, the westerly winds weaken along the equator, thus allowing warm and humid air masses to shift towards central Africa [71]. A distinct consistency between PWV anomalies and DMI can also be found for DGAR station. Despite a smaller cross-correlation coefficient between the PWV and DMIwest index (0.58 with 0 months lag), IOD is a major factor in modulating the climatic variations of PWV. The highest positive PWV anomaly (2.7 mm) occurred in December 2015, when the SST in the western Indian Ocean basin was also the highest during the analyzed time span. The highest negative (−1.7 mm), in turn, took place in 2011 when the DMIwest index was also negative. Although during the analyzed time span there were also similar or even greater negative variations in this index, the changes in PWV caused by them were less significant. However, since this station is also positively correlated with ENSO (0.56 with +2 moths lag), such an effect may be partly explained by a modulating effect of the extreme 2010-2012 La Niña.
From easterly located stations, the highest influence of the IOD on air humidity was found for the DARW station. In this case, the cross-correlation coefficient was equal to 0.59, although the lag time was also distinct and amounted to +9 months. This is visible on Figure 6b where the highest (4.4 mm) and second highest  [72] and the general rule that La Niña increases humidity and rainfall mostly over northeast Australia, while El Niño brings droughts. They also partly explained the lack of full agreement between DARW PWV non-linear trend variability and the DMIeast index. Similar variations in PWV were observed at the BAKO station. In this case, IOD also affects the air humidity the most. The cross-correlation coefficient value calculated based on the DMIeast index and PWV anomaly was equal to 0.52, and, similarly to the DARW station, was characterized by a distinct lag, equal to +8 months. The highest PWV anomaly was observed not in 2010, at the time of the record-breaking La Niña, but in September 2016 (4.3 mm), simultaneously with the highest DMIeast index anomalies, most likely due to the fact that this station is less correlated to the ENSO (−0.42 with −5-month lag). The lowest, but still significant, similarity between IOD phases and PWV non-linear trend variability was found for the COCO station. The value of a cross-correlation coefficient between PWV anomalies and DMIeast index at this station was 0.48, with a +6-month lag. As in the case of the DARW and BAKO stations, the two highest PWV values were observed together with the highest SST anomalies expressed through the DMIeast index and amounted to 5.2 mm and 3.1 mm for September 2010 and April 2016, respectively. On the other hand, the highest anomaly was observed during the record-breaking La Niña (2010-2012), although this station is not correlated with the ENSO (0.26 with +8 months lag). This indicates that despite a general lack of dependency between ENSO and PWV long-term variability over Coco Island, an extremely strong ENSO phase may affect the local moisture field.

Atlantic Related Climate Indices
Much like those of the Pacific and Indian Oceans, Atlantic large-scale climate indices also modulate moisture variability over the analyzed stations. In contrast to the ENSO, their influence is more limited to the area of occurrence of SST anomaly and consequently we were able to distinguish several groups of stations that are affected by various phenomena.
Distinct inter-annual signals were found in the Canaries, at the MAS1 (Maspalomas, Gran Canaria island) and LPAL (Roque de los Muchachos, La Palma island) stations. Although Gallego et al. [73] have found the ENSO signal in the subtropical Atlantic, in results obtained by us of PWV anomalies, did not show a direct dependency on this phenomenon. Therefore, we followed conclusions obtained by Sánchez-Benítez et al. [74], who pointed out that tropical Atlantic SST is responsible for the largest variability of precipitation in this region. Figure 7a presents PWV non-linear changes with reference to the Tropical Northern Atlantic index (TNA) [75], that describes anomalies in the monthly SST from 15 • W-57.5 • W and 5.5 • N-23.5 • N basin. Both stations show overall consistency with the TNA and small discrepancies between LPAL and MAS1 may be the effect of altitude difference, which exceeds 2000 m. The values of cross-correlation coefficients amounted to 0.48 and 0.45 with −3 and −4 month lags for LPAL and MAS1, respectively.
Atmosphere 2021, 12, x FOR PEER REVIEW 12 of 24 gersdorp, South Africa). Since in general the climate of south Africa is strongly associated with the South Atlantic SST [76], the HARB and HRAO stations were investigated for occurrence of oscillations resulted from it. For this purpose, the Tropical Southern Atlantic (TSA) index [75], which describes SST anomalies in the 30° W-10° E and 20° S-0° basin, was used. As can be observed (Figure 7b), both PWV non-linear trends tend to be modulated by the TSA, although the absolute values of PWV variations do not exceed 0.95 mm. Both these stations, however, are located inland at an altitude of more than 1388 m, which generally reduced the average amount of water vapor. The strong negative anomaly up to −1.1 mm in 2018 is linked to the severe drought that occurred in South Africa [77]. Similarly, as in case of previously investigated climate indices, again we calculated cross-correlation coefficients between PWV non-linear trend and the TSA index. Despite the overall compliance resulting from Figure 7, the values of cross correlations coefficient were small and amounted to 0.32 and 0.35 (with +3-and +5-months lag), for HARB and HRAO, respectively. These results may arise from the fact that there are visible oscillations in PWV long term trends, most likely driven by local environmental conditions and generally low levels of PWV (similar were observed in LPAL and MAS1). To reduce their impact on the cross-correlation process we adopted, for both the TSA index and PWV trend, the annual moving average. As a result, the value of cross-correlation coefficient increased to 0.52 and 0.59 for HARB and HRAO, respectively (a similar procedure adopted for LPAL and MAS1 cases resulted in an increasing cross-correlation coefficient up to 0.54 and 0.55, respectively). Long-term variations in SST in the Atlantic were also tested to verify their impact on Central American-Caribbean region (CARIB) PWV variability, in which three stations are located, namely CRO1 (Christiansted, Virgin Islands, United States), GUAT (Guatemala City, Guatemala) and MANA (Managua, Nicaragua). Martinez et al. [78] pointed out that Atlantic Multidecadal Oscillations (AMO) [79] should be considered during investigation On the opposite side of Equator, Atlantic SST anomalies also seemed to affect PWV long term variability at the HARB (Pretoria, South Africa) and HRAO stations (Krugersdorp, South Africa). Since in general the climate of south Africa is strongly associated with the South Atlantic SST [76], the HARB and HRAO stations were investigated for occurrence of oscillations resulted from it. For this purpose, the Tropical Southern Atlantic (TSA) index [75], which describes SST anomalies in the 30 • W-10 • E and 20 • S-0 • basin, was used. As can be observed (Figure 7b), both PWV non-linear trends tend to be modulated by the TSA, although the absolute values of PWV variations do not exceed 0.95 mm. Both these stations, however, are located inland at an altitude of more than 1388 m, which generally reduced the average amount of water vapor. The strong negative anomaly up to −1.1 mm in 2018 is linked to the severe drought that occurred in South Africa [77]. Similarly, as in case of previously investigated climate indices, again we calculated cross-correlation coefficients between PWV non-linear trend and the TSA index. Despite the overall compliance resulting from Figure 7, the values of cross correlations coefficient were small and amounted to 0.32 and 0.35 (with +3-and +5-months lag), for HARB and HRAO, respectively. These results may arise from the fact that there are visible oscillations in PWV long term trends, most likely driven by local environmental conditions and generally low levels of PWV (similar were observed in LPAL and MAS1). To reduce their impact on the cross-correlation process we adopted, for both the TSA index and PWV trend, the annual moving average. As a result, the value of cross-correlation coefficient increased to 0.52 and 0.59 for HARB and HRAO, respectively (a similar procedure adopted for LPAL and MAS1 cases resulted in an increasing cross-correlation coefficient up to 0.54 and 0.55, respectively).
Long-term variations in SST in the Atlantic were also tested to verify their impact on Central American-Caribbean region (CARIB) PWV variability, in which three stations are located, namely CRO1 (Christiansted, Virgin Islands, United States), GUAT (Guatemala City, Guatemala) and MANA (Managua, Nicaragua). Martinez et al. [78] pointed out that Atlantic Multidecadal Oscillations (AMO) [79] should be considered during investigation of Caribbean inter-annual variability, particularly in terms of rainfall. To verify whether the air humidity is also under the influence of this oscillation we conducted a cross-correlation process. Although AMO describe long-duration changes in SST throughout the entire North Atlantic basin (0 • -80 • N, 80 • W-0 • E), we obtained significant values of crosscorrelation coefficient that amounted to 0.47 (+1 months lag) and 0.45 (+1 month lag) for CRO1 and MANA, respectively. At the same time the GUAT station was characterized by only a 0.31 cross-correlation coefficient. Nevertheless, a much higher consistency was obtained when instead of AMO, the Caribbean (CAR) index, that describes SST anomalies in the Caribbean Sea basin, was used. Figure 8a presents PWV non-linear trends for CRO1, GUAT and MANA stations compared to the CAR index. As can be noticed, the stations are in good agreement with SST variability in this region. More precisely, cross-correlation coefficient values were calculated at 0.48 (+7 months lag), 0.63 (+4 months lag) and 0.64 (+7 months lag) for CRO1, GUAT and MANA, respectively. The SST anomaly in all CARiB region does not, however, explain all PWV variations, e.g., negative anomaly for CRO1 station in March 2015 (−1.3 mm). Despite this, the overall dependency is clear, much as in the case of the GUAT and MANA stations. In the case of CARIB stations, it is also worth to emphasizing two aspects: (1) a clear upward trend in all the stations analyzed in this region and (2) occurrence of the highest and lowest anomaly at CRO1 (Christiansted, Virgin Islands, U.S.) and MANA (Managua, Nicaragua) stations that tend to appear at a frequency of about 5-6 years.
term of coefficient and lag values). Negative anomalies also occurred at a similar time, although they are shifted by several months. Surprisingly, the highest similarity between the CAR index and PWV variability was found for the NKLG station. In Figure 8b it is clearly visible how well matched these two signals are. The calculated cross-correlation coefficient amounted to 0.73 with only +1 month lag, which indicates strong teleconnections between SST in northern-east Atlantic and PWV long-term variability over central-west Africa.

Non-Linear Trend Variability in Relation to Local SST
Although the impact of climate indices analyzed so far had a global or regional range, not all stations were characterized by statistically significant variations related to them. In addition, the investigated influence of e.g., ENSO or IOD did not explain all PWV non-linear trend changes over the stations. Therefore, we conducted an analysis of the impact of local SST variability on the air moisture. The local SST was determined as an average value from the SST (GHRSST project) in a 2-degree radius around the station, and in this form subjected to the SSA process, with the aim of estimating a non-linear trend. Consequently, we obtained SST decadal variability for 38 out of 43 stations, that is, for stations which were in sufficient proximity to the oceans. Stations excluded due to these criteria were: BOGT (Bogota, Colombia), BRAZ (Brasilia, Brazil), HARB (Pretoria, South Africa), HRAO (Krugersdorp, South Africa), MBAR (Mbarara, Uganda) and UNSA (Salta, Argentina). In contrast to the formerly investigated climate indices, in this case we only calculated correlation at lag 0, since local SST variations do not involve large-scale teleconnections in affecting local PWV. Detailed information about the obtained results is given in Table 2. For better comparison with previously discussed results, the table was SST anomalies in the Caribbean Sea were also found to affect stations located much further to the south and east. We found out that the PWV non-linear trend estimated for KOUR (Kourou, French Guiana) and NKLG (Libreville, Gabon) (Figure 8, right) represents a similar signal as the CAR index. For the KOUR station cross-correlation coefficient amounted to 0.56 (+9 months lag), which is similar to results obtained for CRO1 (both in term of coefficient and lag values). Negative anomalies also occurred at a similar time, although they are shifted by several months. Surprisingly, the highest similarity between the CAR index and PWV variability was found for the NKLG station. In Figure 8b it is clearly visible how well matched these two signals are. The calculated cross-correlation coefficient amounted to 0.73 with only +1 month lag, which indicates strong teleconnections between SST in northern-east Atlantic and PWV long-term variability over central-west Africa.

Non-Linear Trend Variability in Relation to Local SST
Although the impact of climate indices analyzed so far had a global or regional range, not all stations were characterized by statistically significant variations related to them. In addition, the investigated influence of e.g., ENSO or IOD did not explain all PWV nonlinear trend changes over the stations. Therefore, we conducted an analysis of the impact of local SST variability on the air moisture. The local SST was determined as an average value from the SST (GHRSST project) in a 2-degree radius around the station, and in this form subjected to the SSA process, with the aim of estimating a non-linear trend. Consequently, we obtained SST decadal variability for 38 out of 43 stations, that is, for stations which were in sufficient proximity to the oceans. Stations excluded due to these criteria were: BOGT (Bogota, Colombia), BRAZ (Brasilia, Brazil), HARB (Pretoria, South Africa), HRAO (Krugersdorp, South Africa), MBAR (Mbarara, Uganda) and UNSA (Salta, Argentina). In contrast to the formerly investigated climate indices, in this case we only calculated correlation at lag 0, since local SST variations do not involve large-scale teleconnections in affecting local PWV. Detailed information about the obtained results is given in Table 2. For better comparison with previously discussed results, the table was supplemented with information about the climate index with which we investigated the highest correlation. Results show that most of the analyzed stations present a high agreement between PWV and SST non-linear trend. The absolute mean value of correlation coefficient, calculated based on all 38 stations, was 0.57. The highest agreement was found for YAR2 (Yarragadee, Australia) (0.87), while the lowest for the TCMS (0.25) station. Figure 9 (top) presents examples of stations that are in high agreement with local SST, namely YAR2 and the THTI (correlation 0.80) station. In their case it can be clearly noticed how GNSS PWV non-linear trend extracted through the SSA method varied simultaneously with SST: the highest and lowest PWV anomalies occurred at almost the same time as the highest and lowest SST anomalies. This is especially evidenced, again, during the extremely strong La Niña (2010-2012), that was responsible for occurrences of the highest anomaly at the YAR2 station (2.1 mm), and the lowest at the THTI station (−3.9 mm). Although the THTI station is highly correlated with ENSO (cross-correlation 0.70), such results were less expected at the YAR2 station, which is poorly correlated with ENSO (cross-correlation below 0.4). However, since this ENSO phase was record-breaking and caused an increase in humidity throughout almost all of Australia [80,81], we found its distinct impact also at other stations located in west Australia, namely YARR (Yarragadee, Australia) and KARR (Karratha, Australia) (not shown). Such an impact of the strong 2010-2012 La Niña episode on moisture distribution over the Australian region (including the previously discussed TOW2 and DARW stations) was also recognized in GPS-derived PWV by Choy et al. [82]. Table 2. Correlation between PWV non-linear trend and local SST non-linear trend, supplemented with information on the climate index correlated with given station the most (p-positively, n-negatively). Another group of stations were characterized by a slightly lower, but still statistically significant similarity between PWV and SST variability. ASPA (correlation 0.58) and TWTF (Taoyuan, Taiwan) (correlation 0.59) may serve as examples (Figure 9, middle). Like YAR2 and THTI, the ASPA station shows a clear relationship between PWV and SST variability. However, on the basis of this station it can also be noticed how this method is sensitive to a lack of original data. Despite the fact that all GNSS PWV time series were pre-processed, in the way of filling gaps, the negative peak in SST in September 2008 is not reflected in PWV, thus showing how data gaps negatively affect the reliability of non-linear trend estimation using the SSA method. On the other hand, there was also the TWTF station, that represents a distinct trend that can only partly be explained by raising SST (since 2011). Nevertheless, such a trend is consistent with the positive PWV trend between 2006-2012 in Taiwan discussed by Yeh et al. [83]. This station, is also characterized by a distinct oscillation, with a 3-5-year period. Almost identical signal, both in terms of trend and oscillations, were found at the neighboring stations TCMS and TNML (both Hsinchu, Taiwan).

Site
The last group states for the stations reflects the absolute value of the correlation coefficient less than 0.4 (Figure 9, bottom). In this group we can distinguish stations like REUN (Le Tampon, France) (correlation 0.36), whose PWV non-linear trend presents some similarity to the SST anomalies, although due to the altitude difference (1552 m) and local environment conditions they are expressed by a statistically lower correlation to the SST. Reduced similarity between PWV over the station and local SST can also be related to the lack of GNSS data, like in the case of the FALE station (correlation 0.33) (Figure 9, bottom). This station had 72% of the originally delivered GNSS solutions, while the remaining 28% were calculated using the SSA method. The lack of real data affects not only the PWV nonlinear trend during the period for which data were artificially created, but may also affect PWV before and after data gaps. Consequently, longer and more numerous breaks results in less reliable trend estimation, which have been found at the FALE station since 2013. Therefore, any analyses and interpretations should taking into account the presence of filled breaks in the time series. This effect also contributes to reducing the value of the correlation coefficient for the formerly mentioned TCMS and TNML stations, especially compared to the TWTF station. Although they are all in the similar location and have analogous PWV anomalies, TCMS and TNML are characterized by much a lower correlation (0.25 and 0.27) than for TWTF. However, in contrast to TWTF (89% data completeness), both TCMS (74% data completeness) and TNML (70% data completeness) suffer from a lack of GPS observations, especially since 2008. pre-processed, in the way of filling gaps, the negative peak in SST in September 2008 is not reflected in PWV, thus showing how data gaps negatively affect the reliability of non-linear trend estimation using the SSA method. On the other hand, there was also the TWTF station, that represents a distinct trend that can only partly be explained by raising SST (since 2011). Nevertheless, such a trend is consistent with the positive PWV trend between 2006-2012 in Taiwan discussed by Yeh et al. [83]. This station, is also characterized by a distinct oscillation, with a 3-5-year period. Almost identical signal, both in terms of trend and oscillations, were found at the neighboring stations TCMS and TNML (both Hsinchu, Taiwan). The last group states for the stations reflects the absolute value of the correlation coefficient less than 0.4 (Figure 9, bottom). In this group we can distinguish stations like REUN (Le Tampon, France) (correlation 0.36), whose PWV non-linear trend presents

Linear Long-Term Variability
Next to the non-linear, we also investigated linear variability, to evaluate the general trend of changes using least square adjustment. Figure 10 represents PWV linear trend estimated for the January 2001-December 2018 period. The weighted mean value, estimated based on all stations (where the individual uncertainties were treated as weights), was 0.08 mm/year, with an uncertainty (1-sigma) equal to 0.01 mm/year (note that the KOKV station was excluded from this analysis due to insufficient length of the time series). Although comparative analysis of tropospheric linear trend for different time spans is not fully reliable [84], we decided to refer our results to those obtained by Chen and Liu [36], who also studied GPS-based trend of PWV in tropics during a similar period. For the 2000-2014 time span they verified the occurrence of a trend equal to 0.023 mm/year, which is much lower than the trend estimated by us. Such significant difference, however, can be mostly explained by the number of stations involved in the analysis. In our case we used 43 stations, which provide broader coverage than the number of stations used in the aforementioned study (17). In addition, Chen and Liu clearly pointed out that inhomogeneous distribution and the small number of stations reduces the reliability of the estimated GNSS-based linear trends. On the other hand, they have also estimated trends using other techniques and obtained similar results to those presented by us: 0.066 mm/year (ECMWFbased), 0.083 mm/year (NCEP-based) or 0.073 mm/year (radiosonde-based). Although in both these studies the mean linear trend value was positive, it is worth underlining that none of the stations analyzed by us were characterized by a negative trend. The highest trend was found for the BAKO stations and amounted to 0.18 ± 0.02 mm/year, while the smallest one was found for TUVA and was equal to 0.00 ± 0.02 mm/year. According to Figure 10 it can be stated that the largest changes occurred in the area from Australia up to Taiwan. Next to the BAKO, stations with other high rates of PWV linear trends can be found here, including DARW (0.14 ± 0.02 mm/year), KARR (0.14 ± 0.02 mm/year), NTUS (0.15 ± 0.01 mm/year) or TNML (0.17 ± 0.02 mm/year). Distinct trends were also observed in Hawaii and in Central America, while smaller ones were observed in the southern-west Pacific and South America. In general, there was no correlation between rate of PWV change (expressed through a percentage ratio between the linear trend and the mean PWV) and station altitude. This is most likely related to the insufficient number of analyzed stations, which does not allow for a reliable estimation of such dependency, especially considering various climate phenomena affecting local PWV across global tropics. An exception to this is the Hawaii region, in which four stations are located close to each other, therefore under a similar climate phenomenon, but with large altitude differences. In this case the highest changes were observed at the highly seated MAUI (3044 m) and MKEA (3728 m), and amounted to 3.03% and 2.69% per year, respectively. Clearly lower results, equal to 0.44% were obtained for the KOKB station (1150 m), while the lowest, 0.36%, for the HNLC station (6 m). According to Figure 10 it can be stated that the largest changes occurred in the area from Australia up to Taiwan. Next to the BAKO, stations with other high rates of PWV linear trends can be found here, including DARW (0.14 ± 0.02 mm/year), KARR (0.14 ± 0.02 mm/year), NTUS (0.15 ± 0.01 mm/year) or TNML (0.17 ± 0.02 mm/year). Distinct trends were also observed in Hawaii and in Central America, while smaller ones were observed in the southern-west Pacific and South America. In general, there was no correlation between rate of PWV change (expressed through a percentage ratio between the linear trend and the mean PWV) and station altitude. This is most likely related to the insufficient number of analyzed stations, which does not allow for a reliable estimation of such dependency, especially considering various climate phenomena affecting local PWV across global tropics. An exception to this is the Hawaii region, in which four stations are located close to each other, therefore under a similar climate phenomenon, but with large altitude differences. In this case the highest changes were observed at the highly seated MAUI (3044 m) and MKEA (3728 m), and amounted to 3.03% and 2.69% per year, respectively. Clearly lower results, equal to 0.44% were obtained for the KOKB station (1150 m), while the lowest, 0.36%, for the HNLC station (6 m).

Discussion
The main goal of this study is to demonstrate a usefulness of the GNSS technique in monitoring long-term, climate-related changes in tropospheric moisture patterns. Variability in tropical climate patterns triggers various changes in the prevailing weather conditions over the globe, for which their in-depth examination is crucial. One such canonical tropical phenomenon is ENSO. The GNSS technique may provide novel insight into shifts in regional and global distribution of atmospheric moisture during ENSO initiation, evolution and decay. This can be achieved because GNSS techniques are able to represent changes in PWV associated with ENSO. In fact, PWV variability derived from GNSS data compares well with RS data across a range of temporal scales ( Figure 11). Furthermore, GNSS data allow direct analysis of sub-diurnal to intra-annual time scale. This provides an opportunity to investigate multi-scale interactions based on GNSS data, without additional resources, which is a necessity for RS networks. Such data can help to leverage detailed field studies with targeted observations. As an example, the GNSS station in Tahiti, which has been operational during NOAA ENRR field program, was able to monitor evolution of the atmospheric moisture content throughout the 2016 El Niño event, that is beyond the IOP between 26 January and 29 March (Figure 11). In order to fully leverage GNSS observations, one need to break the location constraint-just as all operational RS stations, GNSS observations are land-locked. Most of the global tropics, including areas where phenomena such as ENSO or IOD evolve, is open ocean. However, recently PWV data have been retrieved from a small, regional network of ship-borne GNSS receivers [85]. This proves that GNSS technology can be expanded to oceans, taking advantage not only of research fleet but also densely populated shipping routes across global oceans. This would be neither a simple nor an easy objective, but appropriate frameworks already exist. For example, the Ship of Opportunity Program (SOOP) leverages cargo ships to conduct automated oceanographic observations. Integration of GNSS observations in such frameworks is definitely feasible and would provide a tremendous improvement in observations of atmospheric moisture content over open ocean.

Conclusions
This paper addressed the subject of inter-annual long-term variability of the tropical In order to fully leverage GNSS observations, one need to break the location constraint-just as all operational RS stations, GNSS observations are land-locked. Most of the global tropics, including areas where phenomena such as ENSO or IOD evolve, is open ocean. However, recently PWV data have been retrieved from a small, regional network of ship-borne GNSS receivers [85]. This proves that GNSS technology can be expanded to oceans, taking advantage not only of research fleet but also densely populated shipping routes across global oceans. This would be neither a simple nor an easy objective, but appropriate frameworks already exist. For example, the Ship of Opportunity Program (SOOP) leverages cargo ships to conduct automated oceanographic observations. Integration of GNSS observations in such frameworks is definitely feasible and would provide a tremendous improvement in observations of atmospheric moisture content over open ocean.

Conclusions
This paper addressed the subject of inter-annual long-term variability of the tropical PWV derived from GNSS observations. First, non-linear changes were identified by means of the SSA method, showing how significant PWV anomalies in tropical regions can be dependent on the occurrence of various climate phenomena. At most of the analyzed stations, PWV was found to be affected by the ocean-atmosphere conditions caused by ENSO. A significant cross-correlation (the coefficient exceeded 0.4) between PWV nonlinear trend and MEI index (describing ENSO phases) was found at 21 out of 43 stations (49%), located in the Pacific Ocean basin, Australia, the Indian Ocean basin, Africa, and South America. Twelve of them were positively correlated to the ENSO, thus showing an increase in PWV anomalies during El Niño episodes and a decrease during La Niña, while nine of them were negatively correlated. Although the ENSO signal was observed across all global tropics, it was a dominant mode of PWV long-term variability only for stations located in the vicinity of the west Pacific basin (11 stations). Nevertheless, its extremely strong phases affected PWV situated outside the west Pacific, including stations that statistically are not correlated with the ENSO.
Beside the ENSO, we found a distinct influence of five other phenomena, expressed through regional climate indices, that were identified as a primary driver of PWV nonlinear variability over the analyzed stations. These cases concern among others the Hawaii region, for which stations show significantly higher consistency with the HAW and NPGO (north Pacific Gyro Oscillation) index than with the MEI, but also stations in the Indian Ocean that were reflecting changes related to the different phases of IOD. The DMI index (east and west variant) turned out to be the most consistent with the PWV anomalies at five stations, including the MBAR station located in east-central Africa. PWV non-linear trends over the Atlantic basin were dominated by Atlantic-related climate indices. Here, TNA and TSA were mainly affecting the subtropical stations, located on the Canary Islands and in South Africa. At the same time, variations related to the CAR index were found in the Central American-Caribbean region and northern-east South America. The highest consistency between this index and PWV long-term variability, however, was found for the NKLG station, located on west coast of Africa, thus indicating strong influence of the northern-east Atlantic on the long-term moisture variations in equatorial Africa.
In total, nine stations were found not to be significantly influenced by any of the analyzed climate indices. These stations were located in South America (AREQ (Arequipa, Peru), BOGT (Bogota, Colombia), BRAZ (Brasilia, Brazil)), west Australia (KARR, YAR2, YARR) and Taiwan (TCMS, TNML, TWTF). In the case of the South American stations, the lack of correlation between climate indices and PWV non-linear trend is most likely due to their locations. They are situated inland at high altitudes (2449 m, 2554 m and 1119 m, for AREQ, BOGT and BRAZ, respectively), which limits the impact of oceanic-related phenomena on local PWV. Such a situation, however, did not occur regarding Australian stations, which are situated on the west coast, in the immediate vicinity of the Indian Ocean, or regarding Taiwanese stations that are also located on the coast. PWV non-linear trends from these stations, however, showed a correlation to the long-term variability of local SST. In general, the mean value of the cross-correlation coefficient between PWV and the local SST non-linear trend, calculated based on 38 stations, was 0.60. For more than 70% of the stations the correlation value exceeded the threshold of significance (0.4), reaching a maximum value up to 0.87 (YAR2 station).
There was also an upward linear trend in the mean GNSS PWV in the tropics. While the weighted mean value was estimated at 0.08 ± 0.01 mm/year, there was visible spatial variability in the results obtained for particular stations, e.g., distinct trends were observed in the regions of the Maritime Continent, Australia and the Caribbean, while more subtle changes were observed for South America, Africa and the Pacific. Generally, not a single station was characterized by a negative trend value, and the highest linear trend was found for the BAKO station (Cibinong, Indonesia) and amounted to 0.18 ± 0.02 mm.
The analysis conducted here showed that long-term homogeneous GNSS PWV time series delivered through PPP processing can be successfully used in inter-annual climate variability studies. By analyzing GNSS PWV using the SSA method, we obtained results with a clear PWV response to different climate indices and local SST. In addition, the SSA method allowed us to investigate the total influence of many factors on the inter-annual PWV variability, including those unrelated to the large-scale climate phenomena, such as local droughts. The data completeness, however, is one of the major factors that affect proper analysis results.