On the inclusion of GPS precipitable water vapour in the nowcasting of rainfall

Introduction Conclusions References


Introduction
Atmospheric water vapour is a very heterogeneous and often rapidly varying meteorological field, being difficult to monitor.Direct observations of that variable by synoptic meteorological stations provide only local surface measurements, while radiosondes provide a continuous vertical profile although at poor temporal resolution due to its expensive operational cost.Remote sensing of water vapour by different space-borne platforms, namely geostationary and polar or-biting sensors, also suffers from severe time and/or spatial sampling problems.In recent years, GNSS (Global Navigation Satellite System) data, mostly that associated with the more popular GPS (Global Positioning System) system, obtained from ground-based stations, started to offer an alternative image of the distribution of the vertically integrated water vapour (precipitable water vapour, PWV), and that information was found to contribute positively to atmospheric analysis for weather forecast (Bevis et al., 1994;Baker et al., 2001;Vedel et al., 2004;Bock et al., 2007;Seco et al., 2012).The rapid increase in the density of such networks in some regions has been motivating the exploration of GNSS data for other meteorological applications.Some examples include 3-D water vapour tomography from short-term field experiments (Champollion et al., 2005), the analysis of countrywide networks (Bender et al., 2011) and even the deployment of high-density GNSS systems specifically designed for the monitoring of atmospheric convection (Adams et al., 2011).It is also expected that with the improvement of the current GNSS systems and the development of future new ones, this technique should provide more refined information about the atmospheric state, allowing also an interoperability in what concerns the water vapour estimation through multi-GNSS processing techniques (Li et al., 2015a) and multi-sensor approaches such as SAR (Benevides et al., 2015a).
Some studies revealed the potential of the GNSS meteorology to analyse the water vapour distribution at near real time and its applicability for weather nowcasting (Haase et al., 2003;Brenot et al., 2006;Yan et al., 2009;Karabatic et al., 2011).Assimilation of GPS-derived water vapour data into weather forecast models has already been studied, resulting in better forecast of severe rain events (Vedel et al., 2004;Cucurull et al., 2004), even at the convective scale (Yan et P. Benevides et al.: On the inclusion of GPS precipitable water vapour al., 2009).The technique allows us to estimate water vapour measurements with high sampling rate with less than 1 h of temporal delay.However a major drawback is the estimation of the near-real-time satellite orbit with high accuracy, which has an impact on the GPS atmospheric measurements (Karabatic et al., 2011).Nevertheless, recent developments have improved the GPS PWV precision through near-real-time orbit, clock and phase delay corrections, using a large network of stations and multi-GNSS processing techniques (Li et al., 2014(Li et al., , 2015a)).Champollion et al. (2004) show the relevance of GNSS data in the analysis of a case study of torrential precipitation in France, where retrieved PWV fields indicate that smallscale variability of water vapour may have a crucial role in the triggering of deep convective processes.In mountainous areas, both shallow and deep-convection processes may be enabled by the redistribution of water vapour in local circulations, which may equally be assessed by GNSS data (Iwasaki and Miki, 2001).Coastal areas, especially with topographic features, are also known to be characterized by important spatial gradients of the water vapour distribution, leading to specific atmospheric circulations, with large horizontal gradients in the cloud and precipitation fields (Champollion et al., 2004;Brenot et al., 2006;Bastin et al., 2007;Bock et al., 2007;Van Baelen et al., 2011).Most of the torrential rainfall events in coastal regions are driven by the advection of highhumidity air masses from the ocean, especially near warm waters, as can be found in the Mediterranean coast (Vedel et al., 2004).If this warm and moist air encounters cold polar air transported southward in the mid and high troposphere, as frequently occurs in that region in the autumn, these are conditions that favour static instability.This will likely lead to convection, especially in places with topographic forcing, initiating deep convective storms (Champollion et al., 2004).In all these, and in many other cases of interest, a better analysis of the distribution of water vapour is a key factor to improve our understanding of the initiation of precipitation events and our ability to forecast such events.
In addition to GPS, the most popular GNSS systems are GLONASS, BeiDou and the future Galileo.The availability of multi-GNSS sensors represents an opportunity to improve the precision and reliability of GNSS measurements, including the atmospheric water vapour sensing.The multi-GNSS data processing increases the number of satellite observations throughout the atmosphere, improving the geometry, the processing time, the redundancy and consequently the robustness of the integrated water vapour measurements (Li et al., 2015a;Benevides et al, 2015b).In particular, the new signals frequency provided by the GPS modernisation and the future Galileo will result in better ionospheric noise estimations reducing the tropospheric delay error (Karabatic et al., 2011).Nevertheless, currently BeiDou and Galileo are not at their full potential, affecting the technique performance in some regions, and most of the GNSS receivers have to be updated in order to collect data from all of the present sys-tems, implying that this technique has room to evolve in the coming years.
The objective of this paper is to analyse the water vapour content given by GPS atmospheric processing, with complementary meteorological data, in order to study precipitation events and provide useful information for its nowcasting.The aim is finding a relationship between the time evolution of the GPS delay and precipitation, for different rainfall events, with emphasis on those associated with torrential rainfall.Three years (2010 to 2012) of hourly accumulated precipitation data were collected from several classical meteorological stations in the vicinity of Lisbon, Portugal, and the nearest GNSS stations were processed in order to obtain a PWV time series.This region has a Mediterranean climate with strong Atlantic influence, with frequent episodes of heavy rainfall associated with large horizontal gradients in the humidity fields.Four case studies are analysed in some detail to justify a simple conceptual model between the PWV signal and precipitation.A full year of continuous hourly observations of both GPS delay and precipitation is then analysed, and a simple forecast model is tested against these data.

GPS atmospheric processing
The propagation of the GPS signal is affected by atmospheric effects both at high levels, in the ionosphere, and at lower levels, in the troposphere.The ionosphere is a dispersive medium causing a delay that can be practically cancelled through a linear combination of the GPS signal frequencies (L1, 1575.42 MHz andL2, 1227.60 MHz).The remaining atmospheric delay is mainly caused by the water vapour content present in the tropospheric layer.Its pathway departs from the satellites to the station receivers at the speed of light, but alterations are induced in its course because of density differences in the medium causing refraction.Clouds, moisture, temperature and pressure gradients, mainly caused by the water vapour, provoke delay and bending in the signal path travel.Such delay is estimated by remote sensing space techniques, including those directly provided by GPS data at ground stations, in the form of numerical tropospheric parameters in the zenith direction, often complemented by an estimate of their local horizontal gradients.In order to improve the precision in GPS positioning, one needs accurate and precise atmospheric parameters.Conversely, very accurate geodetic positions of the fixed stations and a guarantee of small errors from all sources in the inversion process are prerequisites for a successful estimation of the tropospheric parameters from GPS data processing (Tregoning et al., 1998).
The most influential GPS errors not directly originated by the atmosphere are the antenna phase centre variation, multipath, inaccurate satellite orbits and clock delays (Herring et al., 2010).The tropospheric parameters computed in GPS processing are determined through all sets of slant path observations of delay from each station to all observed satel-lites in its horizon at a single instant, defining the zenith total delay (ZTD) and representing conditions of propagation in a volume of the troposphere resembling an inverted cone (Benevides et al., 2013).The extent of that volume depends on the cut-off elevation angle of the station observations and corresponds roughly to a 25-50 km radius at an altitude of 5 km (Bock et al., 2007;Koulali et al., 2012) around its zenith direction, despite the low fraction of water vapour expected above that height (Brenot et al., 2006).The computation of ZTD is accomplished by the application of mapping functions that project the original slant path delay observations into the zenith position at the station (Niell, 2001;Boehm et al., 2006) and can be separated in the dry (or hydrostatic) and humid (or non-hydrostatic) components.The former is defined in GPS processing as the zenith hydrostatic delay (ZHD) and is due to the effect of dry air, representing a percentage not below 90 % of the total tropospheric delay.The latter, also known in the literature as the wet component, and defined as the zenith wet delay (ZWD), represents less than 10 % of the signal, being related with the tropospheric water vapour distribution.The hydrostatic component varies very slowly in time and can be very accurately estimated, while the wet component is highly variable and hard to model, given the variability of the water vapour and the current deficiencies in its observation.Estimation of the ZHD can be done accurately using surface pressure measurements at the GPS station, following Saastamoinen (1972): where P 0 is the surface pressure (hPa), φ is the geodetic latitude and H ref represents the height (km) above the geoid.Therefore, ZHD is proportional to P 0 .A priori zenith hydrostatic delay error estimates are about 0.2 mm for a pressure measurement accuracy of 1 hPa (Tregoning and Herring, 2006).After the determination of ZHD, ZWD is computed by subtracting ZHD from ZTD. ZWD computed from the GPS delay, can be related to a more common representation of the integrated humidity profile, the previously referred PWV, representing the total mass of water vapour in an atmospheric column with unit area.This quantity is measured in kg m −2 , although it is generally expressed in the equivalent practical unit of millimetres (the height of an equivalent column of liquid water).The empirical relation between PWV and ZWD was proposed by Bevis et al. (1992): where k 3 and k 2 are empirical constants, R v is the specific gas constant for water vapour, ρ is the liquid water density and T m is a mean temperature of the atmospheric column.
In practice T m is often computed from the observed surface temperature through an empirical relation constrained by a large set of radiosonde or reanalysis data (Bevis et al., 1994).
In the present study, the estimation of T m was based on a yearly set of radiosondes in Lisbon (Mateus et al., 2014).
Several studies have shown that PWV estimates obtained from GPS are typically achieved with an accuracy of the order of 1 to 2 kg m −2 (mm) when compared with classical measurements like radiosonde, radiometer, lidar or radar (Rocken et al., 1995;Baker et al., 2001;Niell, 2001;Haase et al., 2003;Brenot et al., 2006;Snajdrova et al., 2006;Bastin et al., 2007;Jin et al., 2007;Byun and Bar-Sever, 2009;Champollion et al., 2009;Li et al., 2014), although some studies reach a bias of 3 kg m −2 or more (Tregoning et al., 1998;Bock et al., 2007;Boniface et al., 2012).The sensitivity of PWV estimation relatively to the uncertainty in κ has been estimated by Bevis et al. (1994) to be of the order of 1 or 2 %.Brenot et al. (2006) also evaluated different conversion factors and estimated that they lead to discrepancies in PWV smaller than 0.3 kg m −2 , which is less than the expected error from direct meteorological PWV measurements.As for the estimation of the PWV in near real time for nowcasting proposes, a compromise in the GPS processing between the adoption of ultra-rapid orbits and its coarser precision has to be taken into account (Byun and Bar-Sever, 2009).Furthermore, the application of precise point position GPS processing instead of a network solution speeds up the processing time but at the cost of degrading the PWV accuracy up to 2 cm (Karabatic et al., 2011).

Location and data
The Greater Lisbon area is characterized by large interannual variability and spatial precipitation heterogeneity (Soares et al., 2012).Annual mean precipitation ranges from 600 to 800 mm and is mostly concentrated in the colder months (November to February) and essentially absent in summer.Located in the western sector of the Iberian Peninsula, it has a Mediterranean climate, strongly modulated in wintertime by the North Atlantic Oscillation (Trigo et al., 2005).In this experiment, the GPS data were collected from a GNSS mesoscale network composed by 15 permanent stations regionally distributed around Lisbon, covering an area of about 100 × 150 km 2 .Such stations are part of the permanent national networks of the Portuguese Geographic Institute (IGP) and of the Army Geographic Institute of Portugal (IGeoE).Meteorological data were provided by the Portuguese Institute of Sea and Atmosphere (IPMA).Figure 1 shows the geographical distribution of the stations belonging to the GNSS network with a closer zoom around the Lisbon centre, jointly with a digital terrain model with darker to lighter tons corresponding to lower to higher elevations, and it also shows the location of the meteorological stations as well as the radiosonde site used for the experiment.The altitudes of the GNSS stations vary from 22 m in PACO to 356 m at ARRA station.Regional relief is not steep but exhibits some complexity, being relatively flat only in the areas around the Tagus basin, which are located eastward of ALCO station.The coastal configuration is rather complex and is marked by the large Tagus (south of FCUL, IGP0, IGEO and west of ALCO) and Sado estuaries (south-east of ARRA and PAML).The geographical setup favours the advection of westerly marine air masses into the coastal areas often characterized by complex topographic relief, occasionally leading to torrential rainfall (Champollion et al., 2004;Bastin et al., 2007).

Methods
The GPS data were processed using the GAMIT/GLOBK software (v10.5;Herring et al., 2010) in daily doubledifference processing sessions on a weekly basis at least, tightly constraining the stations to the ITRF08 reference frame, adding several International GNSS Service (IGS) stations and using IGS precise final orbits to constrain the system.IGS stations around the world are selected to provide a better coordinate estimate and to reduce the correlation of tropospheric parameters between nearby stations within the network (Rocken et al., 1995).A second step is performed only running GAMIT, fixing the precise coordinates estimated in the first step and configuring some software options toward the enhancement of the atmospheric estimates.To reduce an effect of the ZTD discrepancy between daily estimates, a time overlap window strategy (Haase et al., 2003) was adopted in the second step of the GPS processing.The first and last values of the daily window of the GPS data processing have an edge effect caused by the stochastic parameter determination, which can result on several millimetres gap discrepancy (Jin et al., 2007).The time overlap win-dow strategy consists of performing four processing sessions with a 12 h duration (e.g.21:00-09:00, 03:00-15:00, 09:00-21:00, 15:00-03:00 UTC) on each day and using only the central 6 h on each time window, therefore resulting in a complete 24 h set of a day.To avoid the orbit adjustment starting from the final hours of the previous day, orbits adjusted in the first daily processing step are used.ZTD is parameterized on GAMIT as a stochastic function variation of the ZHD using piecewise linear interpolation between temporal nodes, constrained to a Gauss-Markov process that is parameterized with a power density function of 1 cm h 1/2 for the atmospheric zenith parameters and 2 cm h 1/2 for the gradient parameters.Therefore, ZTD is generated within a time interval of 15 min and the respective horizontal variation gradients within 30 min, thus creating equilibrium between the temporal variation coverage of the atmospheric delays and the software processing time consumption (Herring et al., 2010).A 7 • elevation angle cut-off was fixed together with the VMF1 global mapping functions for mapping the slant path directions of the GPS from each satellite elevation angle to the station zenith position, taking into account different factors such as the Earth curvature at different latitudes and seasonal changes (Boehm et al., 2006).Measurements of the atmospheric pressure on the GPS station provide superior precision on the delay determination and minimisation of height station errors throughout the process (Tregoning and Herring, 2006), but unfortunately most of those in this data set do not possess coupled meteorological sensors (except CASC).As a consequence, ZHD was modelled throughout a global grid containing surface pressure data with 6 h frequency values, calculated from the ECMWF meteorological reanalysis.The temperature at the surface is obtained from the Global Pressure and Temperature model, which has a coarser temporal resolution than the pressure values obtained from the global grid data (Boehm et al., 2007).An ocean loading model derived from tides is also used (FES2004), being recommended for the precise calculation of station heights and tropospheric delays, together with an atmospheric pressure loading model from NCEP (Tregoning and van Dam, 2005).

Results
In order to verify the possible contribution of GPS PWV estimates to nowcast severe rainfall episodes, GPS and meteorological data from 2010 to 2012 are here analysed.The latter year was analysed in a continuous mode to justify a statistical analysis of a full annual cycle.For the remaining 2010 and 2011 years the analysis was performed on 12 series of data, each containing one or more rainfall episodes in consecutive days, including all days with accumulated precipitation above 25 mm, in at least one of the observed meteorological stations.
The IDL station, a very well-maintained meteorological observatory that includes a classical rain gauge and a modern digital system, is the closest to the IGP0 GPS receiver, at a distance of about 1 km, and constitutes the best a priori match for the rainfall analysis.However, the verification of the PWV data can only be done at the radiosonde station (Fig. 1) located in the Lisbon airport, which is distanced at about 6 km from the receiver.

The annual cycle of PWV
A characterisation of the annual PWV variability is presented in Fig. 2, showing the 2012 yearly series of the hourly GPS atmospheric data for the IGP0 station.PWV is characterized by a clear annual cycle, the amplitude and behaviour of which are a function of the local climate (Haase et al., 2003;Jin et al., 2007;Byun and Bar-Sever;2009).Hourly PWV values range from 5 mm in winter to almost 45 mm in summer and early autumn.Monthly mean PWV ranges between 7 mm in February and around 20 mm from June to November.The record signal also shows evidence of strong synoptic variability, responding to the passage of water vapour carrying meteorological systems.
In the Lisbon area, most of the seasonal cycle of ZTD is associated with the wet component, ZWD, as ZHD shows little variability and a weak seasonal cycle (Fernandes et al., 2013).However, the level of variability shown in Fig. 2 is consistent with that found in coastal Mediterranean regions (Jin et al., 2007), where due to its low altitude and consequently high mean humidity and humidity variability, large values of standard deviation of the GPS measurements are generally found (Haase et al., 2003).In the region of Lisbon, higher values of PWV are expected in the presence of southwesterly flows carrying maritime tropical air masses over a warmer summer-to-autumn ocean.

Events of severe rain
To characterise the behaviour of PWV retrievals in the presence of severe rainfall, all episodes with observed rain above 25 mm day −1 in the period 2010-2012 were selected, corresponding to a total of 12 time segments, each containing one or more rain episodes and the adjacent days.Here the analysis will be restricted to IGP0 GPS data compared with all the meteorological stations following a longitudinal axis around Lisbon, except for Setubal, which is located further south (Fig. 1).Four cases are analysed in more detail; three with the highest recorded rain rate and a fourth with also a high rain rate but where an important feature of some PWV signals can be explained.
Figure 4a shows data from 28 to 30 October 2010, which includes the strongest event in that period with rain rate up to 38 mm h −1 .In this case, there is a very clear and simple relation between the evolution of the GPS PWV and rainfall.Two episodes of intense rain occur, each after a build-up of PWV, from about 15 to 35 mm, and both initiate a sharp reduction in PWV during the following hours.In each episode, the build-up phase takes much longer than the subsequent PWV reduction, and the peak of rain lags the peak in PWV.Two cold fronts passages were observed during this period, with embedded convection provoking severe flooding in the centre of Lisbon on 29 October.
The second strongest event occurred on 23 September 2012, with maximum observed rain rate around 22 mm h −1 , as shown in Fig. 4b.Here the analysis is extended to 9 days, because there is significant structure in the time series of PWV along the full period.The severe rain event occurred just after the largest PWV peak (at 00:00 UTC on 23 September), with a maximum PWV around 43 mm.Rain in Lisbon (IDL) occurred 2 h after the peak, whereas at Cabo Raso (at the coast, to the west of Lisbon, see Fig. 1) it coincides with the IGP0 peak, an indication that this was an eastward moving storm.A second PWV peak of about 37 mm, is observed late on 25 September, but with much smaller amounts of rain followed by three faster oscillations with little or no observed rain.
The third strongest rainstorm led to 20 mm h −1 in Lisbon, on 21 April 2011, as shown in Fig. 4c.The rain event also follows a local PWV peak with a lag of about 1 h, but that peak was not exceptional (about 22 mm).During the 7-day period (17 to 23 April 2011) the PWV signal and that of rain are rather complex, and there is no clear link between the evolution of PWV and the several observed rain episodes.
The last selected case, presented in Fig. 4d, corresponds to a cold front originated by a polar air mass transported from the North Atlantic, observed in the Lisbon area from 30 to 31 August 2011.It was associated with two very large PWV peaks approaching 40 mm, and observed rain spikes approaching 10 mm h −1 .Three PWV curves are presented referring to stations PACO, IGP0 and ALCO, located from west to east respectively.In addition to the small bias due to different station heights, it is noticed that the signal patterns are similar and seem to occur first in the western locations and later in the eastern ones, possibly indicating an eastward movement of the water vapour.For the IGP0 station, the first PWV peak (38 mm) registered at 14 h of day 30, after a growth of 63 % from 24 mm and lasting 14 h and was followed by a fast decrease down to 22 mm in about 6 h but with no significant observed rain.The second peak (40 mm), observed at 21 h on day 31, after an increase of 18 mm in about 18 h, was followed by strong precipitation peaking at 8 to 10 mm h −1 .Furthermore, some stations registered sev-eral precipitation spikes of similar amplitude in the following day, mostly during negative trends in PWV, but without further significant PWV build-up.

A linear broken line analysis
The case studies previously analysed indicated that one often observes a pattern of heavy rain occurring after a peak in PWV, leading to a subsequent sharp decrease of the latter variable.Despite the observation of this simple characteristic, the relationship between variations of PWV and precipitation may be rather complex, as one would expect due to the heterogeneous distribution of precipitation in both space and time, and especially during severe weather events.The cause effect of this relation is not reversible, meaning that large variations of PWV often occur without a rainfall registered event.However these large variations could indicate cases where the horizontal transport of PWV is responsible for its local reduction, or where rain occurred in non-observed locations.To assess the relative importance of these different situations, one will now look at 1 continuous year of GPS data and their relation with the rainfall meteorological data at a group of paired stations.The pairs GNSS meteorological stations are IGP0-IDL (distanced 1 km), CASC-Cabo Raso (distanced 7.5 km), PAML-Setubal (distanced 3 km) and FCUL-airport (distanced 3.5 km; see Fig. 1 for geographical distribution within the region of study).The analysis uses the 2012 data.
In this method, the evolution of the hourly PWV retrieved from the GPS signal is analysed by linear fitting of the PWV signal to a broken line.PWV temporal evolution at a station is continuously verified, grouping the hourly increments or decrements in line segments.Linear least-squares fitting of the previous 6 h of PWV is performed to evaluate the trend signal.If the signal is reversed, the line is terminated, resulting in a broken line segment with an increasing or decreasing trend.Thereby the PWV hourly variations are approximated by segments of broken lines (or ramps) monotonically changing in the significant peaks, where a ramp signal inversion is observed.A 6 h time interval for the minimum broken line length was set empirically in order to enable the linear fitting algorithm to perform a better discretization of the PWV signal characteristics without being affected by the noisy features that are often verified between two consecutive hourly measurements.This also provides a reasonable time interval to be meaningful in the framework of the rain nowcasting.Observing the example of the pair IGP0-IDL, a total of 8636 of hourly GPS observations at IGP0 were included in the analysis with some hourly gaps between the start and ending of the year due to GPS processing issues discussed above.The aggregated time series contains 1186 ascending and descending ramps.
The case studies presented in the previous section suggest that rain events are likely to occur in descending PWV ramps, a few hours after a significant PWV peak.Such behaviour is consistent with a build-up of PWV in atmospheric convective systems, in either shallow or deep-convection processes, leading to the initiation of rain.However, one must note that such a process is not directly observed in the ZWD data at one single station, since the technique is not capable of reproducing the 3-D nature of the water vapour field.Keeping that in mind, one may test the time correlation between the ZWD signal and rain.Figure 5 shows the fraction of rain events which were observed at the chosen station in the 6 h period after a PWV maximum, as a function of the maximum attained PWV (Fig. 5a), the maximum increase in PWV ( PWV; Fig. 5b) and the maximum rate of change of PWV (∂PWV/∂t; Fig. 5c), all computed during the previously ascending ramp.Rain is here defined at the resolution of standard meteorological observations, e.g.0.1 mm h −1 .Data in Fig. 5a show that rain only occurred when PWV peaks exceeded 12 mm in an ascending ramp and that in general the probability of rain increased with the PWV peak value.Figure 5b shows that ramps with PWV changes above 15 mm are generally associated with rain, with a rain probability above 60 %.An apparently stronger message is given by Fig. 5c, showing a quasi-linear relation between the probability of rain and the ∂PWV/∂t in the preceding ascending ramps.All figures indicate a clear positive correlation between the probability of rain and the corresponding measure of the intensity of the previously ascending ramp.
A more complete view of the 2012 data is shown in Table 1, where observed rain is aggregated in four different PWV ramp classes.These correspond to rain occurring (A) during ascending PWV ramps which last for more than 6 h; (B) in descending ramps up to 6 h after those ascending ramps; (C) during shorter-lived ascending ramps, up to 6 h; (D) during descending ramps not included in (B).Class B, corresponding to 10 % of the total time, is responsible for 50 % of the rain, leading to the highest mean rain rate at 0.34 mm h −1 , more than 5 times the mean annual rain rate.In contrast, class A shows values corresponding to 37 % of the total time and is responsible for only 15 % of the rain, with a rain rate 2 times smaller than the annual average.Between class A and class B there is a change in mean rain rate by a factor of 11.Classes C and D show intermediate results, with class C, corresponding to 11 % of total time, presenting a mean rain rate 50 % larger than the annual average and class D, the longest class with 41 % of the time, at less than half of the annual mean rain rate (the same as class A).
Table 2 uses the same 2012 data to look at the relation between the ramp type and rain intensity, already suggested by Fig. 5.The first main result, not unexpected, is that most of long ascending ramps (657 in class A) did not give rise to observed rain; on average these ramps followed by no rain were, however, weak, leading to 74 % of false positive cases when measured by the tendency of PWV.The remaining 26 % of the long ascending ramps (227 in 884) were followed by rain; the rain exceeded 5 mm in 6 h only in about 5 % (41 in 884).Nevertheless, almost all the rain events larger than 10 mm  (82 % or 22 in 27) occur after long ascending ramps.Considering all events of severe rain in the 3-year period 2010-2012 (last row in Table 2) for the IGP0-IDL stations case, 92 % (24 in 26) of those events of torrential rainfall occurred after long ascending ramps with large prior tendencies of PWV.
Both Table 1 and 2 indicate that there is important information in the local PWV signal directly retrieved from GPS data, at the high temporal resolution that it permits, even if such information is incomplete as it overlooks the 3-D structure of the water vapour field.

A simple forecast algorithm
The results shown in Fig. 5 and in the tables suggest the possibility of using the time evolution of GPS PWV data as an element in the nowcasting of rain.Rain occurs in events of different intensity and duration, and a forecaster is mostly interested in those that either contribute to a large fraction of the total annual rain or exceed relevant thresholds in rain rate.In the case of the Lisbon centre (IGP0-IDL station case) 2012 observations, a total of 614 mm of rain was observed (climate normal is near 750 mm), of which 594 mm was accumulated in 108 events with rain rates above 0.2 mm h −1 , while only 178 mm was accumulated in 12 events with rain rate above 5 mm h −1 .
One will look at the success of very short time forecasting system of rain in two complementary ways: first each rain event (of the two families previously described) is analysed and considered a forecast success if the rate of change of PWV in the previous 6 h exceeded a given threshold (cf.Fig. 5c); then, for the same threshold, the full time series of ∂PWV/∂t is analysed to locate all hours when that threshold is exceeded and the forecast is considered a false alarm, if there was no rain observed in the 6 following hours.This particular analysis was performed with the 2012 data over the single pair located on Lisbon centre, IGP0-IDL, to evaluate the rain forecast potential seen from a single GNSS station.
Figure 6 shows the results of the first forecast test, for three values of the ∂PWV/∂t threshold (0.6, 1.5 and 2.5 mm h −1 ).
For the lowest threshold the system only fails to forecast 36 mm (in 614), recovering more than 90 % of the total annual rain, or 16 mm in 178 of severe rain (also more than 90 %), and most rain events are forecasted up to 2 h or less before.When one increases the threshold, the fraction of non-forecasted rain increases (227 mm for 1.5 mm h −1 , 343 mm for 2.5 mm h −1 ), although for these values of ∂PWV/∂t, the fraction of forecasted severe rain decreases much slowly (84, 75 %).However, as one would expect, the high rate of success of this simple method implies a significant rate of false alarms.
Figure 7 shows the full behaviour of the method as a function of the ∂PWV/∂t threshold, for a wider range of this parameter (0.3 to 3.5 mm h −1 ), summing for all times in the forecast window (1 to 6 h).For the lowest values of the threshold, most rain events are predicted including all severe events, but the rate of false alarms is very high, approaching 85 % when the method forecasts slightly more than 85 % of the total rain.Both the fraction of correctly forecasted rain and the fraction of false alarms decrease monotonically with the chosen threshold for ∂PWV/∂t.The success rate for severe rain events decreases rather slowly until ∂PWV/∂t = 2.5 mm h −1 but drops abruptly beyond that value.When 1.5 mm h −1 is the rain warning threshold, the system detects about 75 % of the total rain and more than 90 % of the severe rain, with a rate of false alarms of 60 to 70 %.At 2.5 mm h −1 , the system still detects 75 % of the severe rain, with a false alarm ratio of about 40 %, but only forecasts 41 % of the total rain.Note that the false alarm rate was computed in two slightly different ways, as explained in the caption of Fig. 7. Figure 7 also shows the mean accumulated rain in the well-forecasted events (black line with triangles), indicating that values of ∂PWV/∂t between 2 and 2.5 mm h −1 lead to the highest values of mean forecasted rain.Considering that the proposed algorithm only relies on the GPS data in a single station, this seems a promising result, meriting further studies with more data in more locations.

Discussion and conclusions
Previously published case studies by different authors have reported consistent and similar results, supporting the existence of a positive correlation between PWV and rain, with atypical positive water vapour content before and during precipitation, but always pointing out that not all PWV peaks lead to precipitation (Champollion et al., 2004;Bastin et al., 2007;Yan et al., 2009;Brenot et al., 2014).Some of the apparent mismatch between PWV data and precipitation may be due to inadequate sampling of a very heterogeneous rain field by a relatively sparse rain gauge network, e.g. to cases when precipitation eventually did occur in some locations but was not observed by the network.However, it is very likely that most of those discrepancies are due to the 3-D nature  of the PWV field, which is not measurable in a single GPS station.
Some hints of a possible way of addressing that issue are displayed in Fig. 8, showing the hourly evolution of 2-D fields of PWV retrieved from interpolated GPS observations over a wider area, based on the larger set of GNSS stations in Fig. 1, during the two main PWV peaks observed in Fig. 4d.In the top row, corresponding to the first non-precipitating peak (without observed precipitation) there is evidence of eastward advection of the vapour field from Lisbon, associated with the south-eastward displacement of the PWV maximum between 15 (peak) and 18 h, fading away in the following hours.However, in the bottom row the maximum of PWV remains rather stationary in the Lisbon area and shows evidence of progressive weakening, consistent with the observed regional rain.
Results shown in Fig. 8 are preliminary, being constrained by the small number of GPS receivers, implying a poor coverage especially in the north-eastern and south-western regions.Nevertheless, they suggest that there is scope for significantly improving the use of PWV maps produced by the GNSS systems in meteorological nowcasting.For example, the use of accurate high-resolution tropospheric gradients from multi-GNSS processing could allow us to identify more clearly strong humidity gradients in severe weather situations (Li et al, 2015b).Work on PWV tomography from GPS data (Champollion et al., 2009;Adams et al., 2011;Van Baelen et al., 2011) suggests that one may also look at the time evolution of 3-D water vapour fields, directly addressing the impact of deep convection in those fields, or the relevance of surface fluxes in the moistening of the atmospheric boundary layer in specific storms (Iwasaki and Miki, 2001;Champollion et al., 2004;Brenot et al., 2014).In the case of events which are dominated by horizontal advection as in inflows of moist warm air in coastal regions (e.g.Cucurull et al., 2004), a simpler 2-D approach plus time as shown in Fig. 8, but with a good horizontal resolution and automatic post-processing possibly taking into account terrain effects (Van Baelen et al., 2011), would be relevant to provide support for weather nowcasting.
Some limitations of GPS PWV data, leading to noisy signals in the retrieved PWV, have been identified.These include antenna scattering effects raising the GPS delay and its residuals and creating a stronger multi-path effect that hinders the separation of the atmospheric components (Emardson et al., 1998;Champollion et al., 2004).Integrated water vapour measurements by GPS are also affected by hydrometeors, such as hail, fog or dew (Van Baelen et al., 2011).Brenot et al. (2006) evaluated the ZHD error for large ice hydrometers formed in hail storms and quantified an overestimation of ZHD up to 18 mm, especially inside strong convective cells, therefore resulting in an underestimation of ZWD and PWV on these locations with an error that can reach 15 %.
In spite of all these limitations the present study showed that, at least in the Lisbon climate, the GPS signal has relevant information concerning the onset and the intensity of rainstorms.This result is consistent with Inoue and Inoue (2007) study, where data from "summer thunderstorm active days" in a 5-year period were analysed, indicating that the PWV increments registered prior to thunderstorms were significantly higher than those observed in other weather conditions, suggesting a connection between large PWV increments and convective precipitation.A study evaluating 1 month of data in Japan (Shoji, 2013) verified, as in this work, that the precipitation frequency increased rapidly when the PWV reached a certain threshold as a function of the surface temperature, which was not considered in this study but is worth exploring in future work.A much larger amount of continuous data of 9 years in north-eastern Spain was evaluated, being developed a neural network system to forecast rain from GPS PWV data and observed surface pressure in a single point, concluding for the existence of some predictability in a range beyond 2 days (Seco et al., 2012).Their method is not directly comparable with the one proposed here, but it clearly supports the idea of using GPS data to improve rain forecasts.
The present paper introduces a simple algorithm to classify the time evolution of the PWV signal on a single GPS station and to identify conditions favourable for the triggering of precipitation.In spite of the simplicity of the approach and of the limitations of a 1-D (time) view from a 4-D (time and space) PWV field, the results indicate that there is relevant information in the signal at least concerning the probability of precipitation occurrence of most rain events, as almost all of the severe rainfall cases were well predicted.However, it is also clear that such 1-D analysis will not provide by itself a reliable rain prognostic, not only because of a large fraction of false alarm forecasts but mostly because it gives no indication on the rain intensity.The use of GPS PWV data with other meteorological measurements will surely improve the current nowcasting system.Concerning the GNSS data, at least a high-resolution 3-D analysis (2-D space plus time) is required to produce a consistent mass balance algorithm to support the use of retrieved PWV fields for weather nowcasting.A full 4-D tomography analysis may be desirable and should be tested.Longer time series of GPS and rain data for different regional locations need also to be analysed in order to understand the dependence of parameters, such the best threshold for the tendency of PWV on a specific local climate.

Figure 1 .
Figure1.Location of GNSS, radiosonde and meteorological stations in the regional Lisbon area with a zoom over the city area.Shading according to the terrain elevation.

Figure 2 .
Figure 2. PWV continuous 2012 series at station IGP0 (top).Blue line represents absolute value while the red line is a 30-day average mean.Daily accumulated precipitation in IDL is represented in red (bottom).

Figure 3 .
Figure 3. Comparing precipitable water vapour (PWV) retrieved from GPS data with that computed from radiosonde data in Lisbon (2012 data mostly at 12:00 UTC).

Figure 3
Figure 3 compares the GPS PWV retrievals with direct PWV measurements made by radiosondes from Lisbon, mostly at 12:00 UTC.Results are overall excellent, with a 0.99 correlation between the two data sets, and a very small negative bias in the GPS signal (about 0.1 mm), corroborating other published comparisons (namely, Seco et al. (2012) and Galisteo et al. (2014), both concerning data in Spain).

Figure 4 .
Figure 4. PWV and hourly accumulated rain from (a) 28 to 30 October 2010, (b) 21 to 29 September 2012 and (c) 17 to 23 April 2011; GPS PWV from IGP0 station is in black, hourly precipitation from IDL reference station is in red; other lines represent rain in nearby stations.(d) 29 August to 2 September 2011, with GPS PWV from two additional stations; ALCO in green and PACO in red.Date is presented in hour day/month.

Figure 5 .
Figure 5. Probability of rain in the 2012 data set as a function of PWV maximum (a), PWV maximum increase; PWV (b) and PWV rate of change; ∂PWV/∂t (c).The x axis is represented by the top upper interval limit of the class, where the interval is 3.0 mm, 3.0 mm and 0.35 mm h −1 respectively.

Figure 6 .
Figure 6.Rain forecasted at different time lags as a function of the chosen threshold for the rate of change of PWV.Black lines: total annual rain (614 mm); red lines: accumulated rain in hours with rain rate above 5 mm h −1 (179 mm).Forecast considered successful when ∂PWV/∂t ≥ 0.6 mm h −1 , 1 mm h −1 , 1.5 mm h −1 , at 1 to 6 h prior to the beginning of the rain event; unsuccessful forecasts are computed at forecast ahead time = 0. Data for Lisbon (IGP0-IDL station pair) in 2012.

Figure 7 .
Figure 7. Evaluation of the success of hourly forecasts of rain, when forecasts are issued by the exceedance of a threshold in the rate of change of PWV.The dashed line (solid squares) represents the ratio of hourly false alarms (a forecast of rain leads to no rain in the following 6 h); the dotted line (hollow squares) represents the ratio of false alarm events (consecutive forecasts are aggregated).The blue line (hollow circles) indicates the well-predicted fraction of rain.The red line (solid circles) indicates the well-forecasted fraction of severe rain (> 5 mm h −1 ).The solid black line (triangles) shows the mean accumulated rain in the well-forecasted events.

Figure 8 .
Figure 8. Hourly 2-D fields of GPS PWV around the two main PWV peaks on Fig. 4d.Black triangles represent the stations from the GNSS network.

Table 1 .
Rain accumulated in different PWV change classes, according to 2012 Lisbon data using all GNSS meteorological station pairs.

Table 2 .
Relation observed between rain intensity and PWV ramp data.Data from all Lisbon station pairs in 2012.Last row ( * ) includes all events in the 2010-2012 period for the pair IGP0-IDL.