Estimating trends in atmospheric water vapor and temperature time series over Germany

Ground-based GNSS (Global Navigation Satellite System) has efficiently been used since the 1990s as a meteorological observing system. Recently scientists have used GNSS time series of precipitable water vapor (PWV) for climate research. In this work, we compare the temporal trends estimated from GNSS time series with those estimated from European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA-Interim) data and meteorological measurements. We aim to evaluate climate evolution in Germany by monitoring different atmospheric variables such as temperature and PWV. PWV time series were obtained by three methods: (1) estimated from ground-based GNSS observations using the method of precise point positioning, (2) inferred from ERA-Interim reanalysis data, and (3) determined based on daily in situ measurements of temperature and relative humidity. The other relevant atmospheric parameters are available from surface measurements of meteorological stations or derived from ERA-Interim. The trends are estimated using two methods: the first applies least squares to deseasonalized time series and the second uses the Theil–Sen estimator. The trends estimated at 113 GNSS sites, with 10 to 19 years temporal coverage, vary between −1.5 and 2.3 mm decade−1 with standard deviations below 0.25 mm decade−1. These results were validated by estimating the trends from ERA-Interim data over the same time windows, which show similar values. These values of the trend depend on the length and the variations of the time series. Therefore, to give a mean value of the PWV trend over Germany, we estimated the trends using ERA-Interim spanning from 1991 to 2016 (26 years) at 227 synoptic stations over Germany. The ERA-Interim data show positive PWV trends of 0.33± 0.06 mm decade−1 with standard errors below 0.03 mm decade−1. The increment in PWV varies between 4.5 and 6.5 % per degree Celsius rise in temperature, which is comparable to the theoretical rate of the Clausius– Clapeyron equation.

Abstract. Ground-based GNSS (Global Navigation Satellite System) has efficiently been used since the 1990s as a meteorological observing system. Recently scientists have used GNSS time series of precipitable water vapor (PWV) for climate research. In this work, we compare the temporal trends estimated from GNSS time series with those estimated from European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA-Interim) data and meteorological measurements. We aim to evaluate climate evolution in Germany by monitoring different atmospheric variables such as temperature and PWV. PWV time series were obtained by three methods: (1) estimated from ground-based GNSS observations using the method of precise point positioning, (2) inferred from ERA-Interim reanalysis data, and (3) determined based on daily in situ measurements of temperature and relative humidity. The other relevant atmospheric parameters are available from surface measurements of meteorological stations or derived from ERA-Interim. The trends are estimated using two methods: the first applies least squares to deseasonalized time series and the second uses the Theil-Sen estimator. The trends estimated at 113 GNSS sites, with 10 to 19 years temporal coverage, vary between −1.5 and 2.3 mm decade −1 with standard deviations below 0.25 mm decade −1 . These results were validated by estimating the trends from ERA-Interim data over the same time windows, which show similar values. These values of the trend depend on the length and the variations of the time series. Therefore, to give a mean value of the PWV trend over Germany, we estimated the trends using ERA-Interim spanning from 1991 to 2016 (26 years) at 227 synoptic stations over Germany. The ERA-Interim data show positive

Introduction
Water vapor is considered the most active greenhouse gas that permanently affects the Earth's climate. Due to its high temporal and spatial variations, the precipitable water vapor (PWV) content in the atmosphere has to be regularly and accurately determined for meteorological and climatological purposes. PWV is the amount of water (in millimeters) that would result from condensing a column of water vapor that extends from the measurement point to altitudes of about 12 km. Water vapor mainly resides in the lowest 3 km of the atmosphere and its content generally increases with air temperature. While other observation systems such as radiosondes and microwave radiometers provide PWV measurements that are limited in the temporal and (or) spatial resolutions, ground-based GNSSs provide time series of accurate PWV estimates with 15 min (for this work) sampling at dense GNSS networks, without significant additional costs. Since Bevis et al. (1992) presented the Global Positioning System (GPS) as an efficient meteorological tool, GNSS data have been increasingly used for estimating atmospheric parameters, particularly precipitable water vapor (Gendt et al., 2004;Luo et al., 2008;Jade and Vijayan, 2008;Bender et al., 2008;Alshawaf et al., 2015). GNSS-Published by Copernicus Publications on behalf of the European Geosciences Union. 3118 F. Alshawaf et al.: Estimating trends in atmospheric water vapor and temperature time series based estimates of zenith total delay (ZTD) or PWV have been assimilated into numerical weather prediction models to improve the quality of the output (Bock et al., 2005;Bennitt and Jupp, 2012). They have also been used to improve the performance of high-resolution atmospheric models (Pichelli et al., 2010). Besides meteorology, GNSS estimates of PWV have been employed over Scandinavia for climatological research (Elgered and Jarlemark, 1998;Gradinarsky et al., 2002;Nilsson and Elgered, 2008). The authors found that PWV shows an increase of 1.2-2.4 mm decade −1 . Haas et al. (2003) used ground-based GPS, very long baseline interferometry, radiosonde, and microwave radiometer data to assess long-term trends in PWV time series over Sweden. An increase of about 0.17 mm year −1 within the period 1980-2002was observed. Hausmann et al. (2017 analyzed a decadal time series of PWV (2005PWV ( -2015 from mid-infrared Fourier transform infrared measurements above the mountain Zugspitze. For that time period, they did not observe statistically significant trend in PWV time series. The PWV time series from ground-based GNSS and the European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA-Interim) data might show temporal inconsistencies due to, for example, hardware replacement or inconsistent processing methods (Ning et al., 2016). Therefore, homogenization of the atmospheric data is indispensable for climatological research to properly estimate climatic longterm trends. Vey et al. (2009) andNing et al. (2016) analyzed PWV time series estimated at global GNSS sites to detect and correct for inhomogeneities in the data. Atmospheric reanalysis models such as ERA-Interim have also been employed for climate research. The analysis fields are produced based on 4D-Var assimilation of regular and irregular meteorological data, including surface and upper-air atmospheric fields (Dee et al., 2011). Bengtsson et al. (2004) observed an increasing long-term trend with a slope of 0.16 mm decade −1 in the water vapor data set of ERA 40 within the period of 1958-2001. They suggested applying corrections for the changes in the observing system when using the data for PWV analysis to achieve trend values comparable to GNSS. ERA-Interim and MERRA (Modern Era Retrospective-Analysis for Research and Applications) were also used for trend analysis (Simmons et al., 2007;Rienecker et al., 2008).
Typically, climate scientists consider a period of 30 years as an appropriate time necessary to average variations in weather and evaluate climatic effects for a particular site, as described by the World Meteorological Organization (Arguez and Vose, 2011). Data collected and averaged or summed in some way over 30 years are referred to as climate normals. A 30-year period is recommended, as it is sufficiently long to filter out the interannual variations or anomalies but at the same time short enough to show climatic trends. It is then obvious that the GNSS temporal span is still too short for estimating a reasonably proper climatic trends in this sense. The previous studies using GNSS-based PWV time series for assessing the trends show highly variable estimates for different time windows as well as different research regions. In this paper, we present the PWV trends estimated using GNSS sites over Germany and compare them with the trends estimated from other data sets. The current climate normal period should cover the period from 1 January 1991 to 31 December 2020 (Arguez and Vose, 2011). So, we analyzed time series available from 1 January 1991 to June 2016 (26 years) to provide more robust information about the climatic trends. These data sets are the ERA-Interim reanalysis and surface meteorological data from the German Meteorological Service (DWD). The former data set provides global PWV grids while the latter does not. However, different studies have used the dew point that is computed using surface measurements of temperature and relative humidity to approximate the total column PWV (Reitan, 1963;Bolsenga, 1965;Smith, 1966;Tuller, 1977). The formula presented to obtain the PWV from surface measurements is described in Sect. 4. This empirical relation requires only information that can accurately be determined on the ground. The accuracy of dew-point-based PWV approximations depends of course on the atmospheric conditions and the variability of the moisture profiles. It is, however, obvious that PWV estimations based just on atmospheric conditions at the Earth's surface would not always be in complete agreement with, for example, PWV values from balloon soundings integrated through the atmosphere. Since the possibility for obtaining a data set with long time series and high spatial resolution for estimating PWV trends is very limited, we evaluated the potential of this method for climate analysis. We first obtained the PWV based on dew point temperature measurements and evaluated the quality of the time series. Then we used them to estimate the PWV as well as temperature trends. In this work, we apply a preprocessing step to evaluate the quality and homogeneity of the time series ahead of the trend estimation. For checking the homogeneity of the time series, we use the ERA-Interim as a reference. We apply the technique of singular spectrum analysis to detect possible change points followed by a t test to identify the significance thereof. The description of the approach for homogeneity check is beyond the scope of this paper and details are found in Wang (2008) and Ning et al. (2016).
In this paper, we want to show that GNSSs are, and will remain, a promising data source for climate research, particularly in the future when the time series are adequately long. We first extract PWV time series by processing GPS data and evaluate their quality using ERA-Interim as a reference. These PWV time series are then used to estimate the water vapor trends using both least squares and Theil-Sen estimator. The same procedure is applied to the ERA-Interim data over the same time period, for validation. Since the GNSS sites are independently mounted, the length of the PWV time series at the network sites is variable. For a reasonable trend estimation, we used time series with lengths of 10 to 19 years. Since the length of time series is variable and it is still below the normal (30 years), suggested by climatologists, we extended the research using other data sets (e.g., ERA-Interim and synoptic data). First, it is more reasonable to give a mean value of the change per decade (trend) using longer time series; second, using time series of the same length at all sites in the research region enables us to observe specific spatial features of the trend, as presented later in the paper. We also calculate the change in PWV per degree Celsius rise in temperature to check if it is consistent with the Clausius-Clapeyron relation.
This paper is organized as follows. In Sect. 2, we describe the method for PWV determination using GNSS data and the comparison with ERA-Interim. In Sect. 3, we present the methods for estimating the atmospheric trends followed by the results of estimating the decadal rate of change in Sect. 4. Then, the conclusions of this research are presented.

Determination of atmospheric PWV from GNSS data
We used GPS data collected in central Europe, mainly in Germany as shown in Fig. 1. The research region is well covered by 351 permanent GNSS sites with an average separation distance of 30 km. Homogeneous time series with length from 10 to 19 years are available from 119 sites. The second data set we used is the ERA-Interim reanalysis with a spatial resolution of 79 km in longitude and latitude, 60 vertical levels with the model top at 0.1 hPa (about 64 km), and 6 h temporal resolution. Additionally, there are 326 meteorological stations operated by the DWD with data profiles spanning more than 60 years at a temporal rate of 1 h. The climate data center created by the DWD provides long homogeneous time series for climate studies (http://www.dwd. de/EN/climate_environment/cdc/cdc_node.html). They provide surface measurements of temperature, pressure, water vapor pressure, precipitation, snow cover, and other meteorological parameters for climate research. In this section, we briefly describe the methods for PWV determination using GNSS phase observations and a comparison between the different data sets.
Based on the method of precise point positioning (Zumberge et al., 1997), GNSS observations are processed to produce site-specific atmospheric ZTD. The ZTD is an estimate of the total propagation delay caused by the dry gases and water vapor of the atmosphere. Employing meteorological data measured directly at the GNSS site or interpolated from the adjacent meteorological station, the zenith hydrostatic delay (ZHD) is calculated. For each GNSS site, the nearest meteorological station triangle is used to interpolate the measurements at that site (Gendt et al., 2004). The ZHD, in meters, at the GNSS site is then calculated according to the model of Saastamoinen (1973) reported in Davis (1986, pp. 51): where H is the orthometric height in kilometers and φ is the latitude of station. P is the corresponding air pressure at the station in hPa. The air pressure P at the GNSS site in Eq. (1) is obtained by vertically interpolating the surface pressure P s using the barometric formula: where T s is the surface air temperature at the meteorological station in kelvin, z and z s are, respectively, the altitude in kilometers of the GNSS and meteorological station above mean sea level (a.m.s.l.), L is the temperature lapse rate in K km −1 , R is the universal gas constant (8.31447 J mol −1 K −1 ), M is the molar mass of Earth's air (0.0289644 kg mol −1 ), and g is the average Earth's gravitational acceleration (9.80665 m s −2 ). In the lower atmosphere, the temperature is related to the elevation change using the following linear relation: By analyzing ERA-Interim temperature profiles over Germany, we found that the lapse rate changes between summer and winter and in space. The value of L varies between 3 and 7 K km −1 for this research region. These values result in 2 mm change in the ZHD at altitude difference of 1 km. Similarly, the change in PWV is below 0.2 mm, which can be neglected. Once the ZHD is calculated, the zenith wet delay (ZWD) is obtained by and it is converted into PWV using the empirical factor (Bevis et al., 1994): For more details on the GNSS data processing, the reader is referred to Gendt et al. (2004) and Bender et al. (2011). We compared the PWV obtained from GNSS with ERA-Interim data. Figure 2 shows the results for three sites at different altitudes as well as the mean and standard deviations of the time series difference. ERA-Interim grid provides values of PWV at grid points separated by about 79 km in longitude and latitude. The ECMWF provides software to horizontally interpolate the current ERA-Interim grid at different locations of the GNSS stations as described in Heise et al. (2009). We did not account for altitude difference, which has a significant impact in mountainous areas. For the sites located in flat terrain, the two data sets show strong correlation with the bias values below 1 mm and standard deviation of less than 2 mm (Fig. 2). The bias between the data sets increases for sites in mountainous regions. The time series of the site 0285 (Garmisch, Germany; 1779 m a.m.s.l.), for example, show a larger bias between GNSS and ERA-Interim data, which is explained as follows: we average PWV of four distant grid points around the GNSS site. With the rough spatial resolution, the variability of surface topography is not well captured in the reanalysis data, which significantly increases the height difference between GNSS and the model, and hence the PWV difference. Additionally, the daily mean in ERA-Interim is obtained by averaging four PWV values per day, while using GNSS there are 96 PWV estimates per day. This positive bias was also observed by Morland et al. (2006) and Camy-Peyret (2002, 2003) when comparing ERA 40 and GPS in the Alps for the site Jungfraujoch at 3584 m a.m.s.l.
For accurate determination of the PWV from GNSS measurements, it is required to have measurements of mainly air pressure and temperature at the GNSS sites or within a short spatial range. In the absence of meteorological measurements, would the interpolation of pressure and temperature from reanalysis data be a good replacement? To answer this question, we compared the PWV time series extracted from the ZTD by using both measurements at the meteorological stations and ERA-Interim data. To calculate the ZHD, the in situ measured pressure and temperature are horizontally interpolated to the GNSS site and then vertically interpolated to the GNSS antenna phase center. For GNSS sites below the lowest ERA-Interim level, the pressure and temperature are extrapolated at the site altitude as described in Heise et al. (2009). The ZWD is then extracted and converted into PWV. Figure 3 shows the scatterplots of PWV obtained using surface measurements and ERA-Interim data. We found that in regions of smooth topography, the ERA-Interim data and the measurements provide almost the same values of PWV and pressure. In regions of steep topographic gradients, however, the ERA-Interim data show slightly different results, which is mainly related to the pressure data as observed from Fig. 3. The deviations between the measured pressure and the ERA-Interim pressure increase in mountainous regions, which affects the calculation of the ZHD and hence the obtained PWV.
Besides station pressure, an important factor for an accurate determination of PWV is the conversion factor , which should be calculated using measurements of surface temperature. Askne and Nordius (1987) determined the conversion factor as follows: where ρ w is the density of water and R w is the specific gas constant of water vapor (461.5 J kg −1 K −1 ). In our research, we used the values of the physical constants k 3 and k 2 given by Bevis et al. (1994), T m was given by Davis et al. (1985) as where T is the air temperature and P wv water vapor pressure at vertical levels. Davis et al. (1985) suggested the use of water vapor pressure and temperature profiles from radiosondes; however, it is easier to get these profiles from numerical atmospheric models. In this work, we obtained T m as described in Heise et al. (2009) using the ERA-Interim model that covers 60 vertical levels extending from the Earth's surface up to 0.1 hPa. T m can be well approximated based on air surface temperature by the following formula (Bevis et al., 1992): where T s is the surface temperature in kelvin. For our research region, we compared T m obtained from both methods (7) and (8) as shown by the scatterplot of Fig. 4. The surface temperature and vertical profiles of water vapor pressure and temperature in Eq. (7) from ERA-Interim were employed. The difference between the T m calculated from both methods at the GNSS site 0522 (Pirmasens, Germany; 399 m a.m.s.l.) has a mean value of 0.97 K and a standard deviation (SD) of 2 K. Repeating the calculations for the site 0285 (Garmisch, Germany; 1779 m a.m.s.l.), the bias increases to 3.02 K and the SD is 1.83 K. Not only surface pressure grids are inaccurate in mountainous regions (Fig. 3d) but also pressure profiles, which might be related to the coarse grid of ERA-Interim. Also, the temperature profiles have inaccuracies, although less than those for the pressure. By using the integration in Eq. (7), the accumulated error in the calculated T m will be higher, and the bias between this T m and that calculated using only the surface temperature will increase, as observed from the right plot in Fig. 4. However, 3 Decadal variability in time series of atmospheric variables

Estimating the trend using least squares regression
Econometricians have developed reasonably simple models that are capable of interpreting, testing hypotheses, and fore-  ). The bias is 0.97 K for the first site and 3.02 for the second, and the SD is 2 K for the first and for the second 1.83 K.
casting economic data. The method was to decompose the time series into a trend, a seasonal, a cyclic, and an irregular component (Enders, 1995). The trend component represents the long-term behavior of the time series, while the seasonal and the cyclic components represent the regular and periodic movements. The time series also contain a stochastic irreg-ular component. Time series of PWV and temperature, for example, have different temporal variations that can be reasonably modeled using this approach. Here holds an additive model, such that the time series y t can be extended as where T t is a deterministic trend component with slow temporal variations, S t represents the seasonal component with known periodicity (e.g., 12 months for PWV and temperature), and I t represents the irregular (stationary) stochastic component with short temporal variations. We did not observe a regular signal that lasts longer than 1 year, so we excluded the cyclic component for the model. The presence of seasonality might mask the small changes in the linear trend. Therefore, for proper trend analysis, the seasonal component has to be estimated and removed from the time series, which is known by seasonal adjustment (Enders, 1995). The deseasonalized data are useful for extracting the long-term trend and exploring the irregular component of a time series. The seasonal adjustment is applied as an iterative procedure as follows. To best estimate the seasonal component, the linear trend has first to be estimated and removed from the time series. There are different methods to estimate the trend such as using moving average or parametric trend estimation. Here, we used the method of moving average with a window length of 1 year that is able to smooth out seasonal and irregular signals. We employ time series of PWV and temperature with daily values (the GNSS-based estimates of PWV have a temporal resolution of 15 min, but we average them to get mean daily values for climatological studies). The trend is estimated as follows: T t = y t−q + y t−q+1 + · · · + y t+q−1 + y t+q d .
Since the time series are daily and the seasonal signal is annual, the value of d is 365 and q = (d − 1)/2. For d = 366, q = d/2 and the trend is estimated from T t = 0.5y t−q + y t−q+1 + · · · + y t+q−1 + 0.5y t+q d .
The estimated trend component is subtracted from the original time series and the detrended signal is averaged to estimate the seasonal componentŜ t as follows. We first obtain with n the number of data samples. Then w t is centered, i.e., we derive a seasonal signal with a zero mean.
For an additive model,Ŝ t should fluctuate around zero to avoid any influence from the trend. The estimated seasonal component is subtracted from the original time series to obtain a seasonally adjusted time series dy t , i.e., Figure 5 shows an example of the trend, seasonal, and irregular components of PWV time series at site 0896 in Berlin.
To estimate the slope of the trend, we fit a straight linê T =b +m t to the trend component produced by the moving average step. The standard deviation of the estimated slope (called standard error) is calculated as (Wigley et al., 2006) where n−2 is the degree of freedom for n data points. The approximate 95 % confidence interval is expressed asm ± 2 sm. Weatherhead et al. (1998) presented another way to calculate the standard deviation of the estimated slope.
where σ I denotes the standard deviation of the irregular component and n y denotes the number of years of the data. φ I represents the one-lag autocorrelation of the irregular component.

Estimating the trend using Theil-Sen estimator
The Theil-Sen estimator presented by Theil (1950) and Sen (1968) aims to robustly find the linear fit of a data set despite containing outliers. If (t 1 , y 1 ), . . ., (t n , y n ) represent the data points, then the Theil-Sen estimator determines the slope of the line that connects each data pair. The median among the slopes of all pairs in the slope of the fit, i.e., The standard error of the estimated slope is calculated as in Eq. (15). We compared the two methods of trend estimation using PWV time series at the site Lindenberg (14 • 6 E, 52 • 12 N), where GNSS, ERA-Interim, synoptic, and radiosonde data are available, as shown in Fig. 6. The bias between PWV from synoptic data and GNSS is 0.04 mm, while that to the ERA-Interim is −0.21 mm. The bias of both GNSS and synoptic PWV to the radiosonde PWV is 0.95 mm, which is because the former are daily values while the radiosonde provides an instant measurement (at 12:00 UTC). This, however, does marginally affect the estimation of the trend. For the least squares method, we estimate and remove the seasonal component and filter out the irregular component to provide the trend shown in Fig. 6b. Table 1 shows the slope of the linear trend estimated at the site Lindenberg using the PWV time series in Fig. 6 computed using the least squares and Theil-Sen methods. Applying both methods to three different data sets shows a positive trend of about 0.5 mm decade −1 with standard deviation of 0.04 mm decade −1 .   We also analyze the change in PWV in relation to the change in temperature. As temperature rises, the air capacity to hold moisture increases at the Clausius-Clapeyron rate. The water vapor pressure e is related to temperature T as follows: where H v is enthalpy of vaporization and R is the universal gas constant. This relationship indicates that 1°C rise in the temperature increases the vapor pressure by 7 %. Based on this formula, the change in the PWV can theoretically be related to the change in the temperature. The PWV is linearly related to the vapor pressure as presented in Tuller (1977), i.e., PWV = 2.3e. By substituting this into Eq. (18), the increment in PWV should, in theory, be the same as the increment in the vapor pressure (approximately 7 %) per degree Celsius rise in temperature. This was also observed by analyzing the temperature, water vapor pressure, and PWV data sets. We obtained the change in PWV and vapor pressure per 1 • rise in temperature as shown in Fig. 12. The increase in the water vapor pressure at 227 stations is in the range of 4.5 and 6.5 %, which is close to the Clausius-Clapeyron rate. We observed a similar rate of change for the PWV with the temperature, or more precisely, PWV 2 PWV 1 = 1.003 e 2 e 1 .

Estimating the trends using GNSS-based PWV
In this section, we show the estimated trends using three data sets, GNSS, ERA-Interim, and synoptic data of PWV and temperature. First, we estimated the trends of PWV at 351 GNSS sites with time series of 4 to 19 years long and the corresponding standard deviations of the estimated slope as shown in Fig. 7. The size of the marker is proportional to the length of the time series (small squares indicate short time series). As observed from the figure, there are high trend values, particularly at sites with short time series. Therefore, in Fig. 8a, we eliminated all sites with time series shorter than 10 years. At the remaining 119 sites the PWV trend varies between −1.5 and 2.3 mm decade −1 (except for six sites) with precision of the estimated trends below 0.25 mm decade −1 .
To validate these estimates, we analyzed ERA-Interim data over the same times where GNSS data are available (Fig. 8c).
The results from concurrent ERA-Interim time series show high similarity in the trend values and the variations of the trend in space.

Estimating the trends using longer time series
Since the trend is estimated from GNSS time series of different length, it is reasonable to provide a mean value for the whole region or observe spatial features of the trends. Therefore, and in order to get more insight and more reasonable conclusions about the long-term temporal variations of PWV, it is necessary to analyze time series spanning one predefined period for all stations. Since the last climate normal extends from 1991 to 2020, we analyzed time series of 26 years (January 1991-June 2016) from ERA-Interim and synoptic data. We investigated time series at 227 meteorological stations where the ERA-Interim is horizontally interpolated at the synoptic station using bilinear interpolation. Figure 9 shows the estimated trends using ERA-Interim PWV time series by, first, applying the least squares to the seasonally adjusted data and, second, using the Theil-Sen method. Both methods show similar values of the trend, positive with values of 0.34 ± 0.06 mm decade −1 . As observed from Fig. 9, the trend tends to increase in the direction to northeastern Germany. In order to validate these results, it is necessary to have a long data set, which is not available for this research. However, DWD provides surface measurements of atmospheric parameters that are accurate and homogenous so that they are proper for climate studies. It is not possible to accurately determine the total column water vapor using surface meteorological observations alone. However, it was shown in the 1960s that it is possible to approximate the atmospheric PWV based on dew point temperature measurements, which is considered an indicator of the amount of moisture in the air (Reitan, 1963). The dew point temperature in turn is determined based on the air temperature and relative humidity. Reitan (1963) presented a basic relationship between the mean monthly PWV and mean monthly surface dew point temperature by the following regression form: where PWV is in centimeters and T d is the dew point temperature in degrees Fahrenheit. a and b are estimated to have the values of −0.981 and 0.0341 (Reitan, 1963). The standard error in the PWV estimate was 0.18 cm. Following the same procedure, Bolsenga (1965) obtained slightly different estimates for a and b using hourly and mean daily observations. Smith (1966) obtained a similar regression equation with the coefficient a not being constant. It depends on the vertical distribution of the atmospheric moisture, i.e.,   with the value of λ dependent on the site latitude and the season of year (Smith, 1966).
In this work, we estimated the coefficients a and b at each meteorological station by fitting the curves in Eq. (19) to the ERA-Interim PWV data. The median values for a and b using daily PWV are −1.346 and 0.039, which are close to the values −1.249 and 0.0427 presented by Bolsenga (1965). For monthly PWV, the median values are −1.224 and 0.037 for a and b, respectively.
We used measurements of surface dew point temperature to obtain the daily PWV and time series for the whole network are evaluated using the ERA-Interim data. The PWV value at the meteorological station is computed by applying bilinear interpolation to the ERA-Interim PWV at four grid points around that station. The altitude difference was not accounted for. Figure 10a shows the bias and standard deviation values of daily PWV for 227 stations as well as the bias against the altitude difference of the two data sets (ERA-Interim height−station height). The bias is centered around 0.15 mm and the standard deviation around 2.5 mm. From Fig. 10b we observe that the higher the altitude difference, the larger is the mean PWV difference.
Next, we estimated the trends using the time series of dew-point-based PWV after removing local environment effects, which are presented in Fig. 11. The trend values vary in the range of 0.48 ± 0.13 mm decade −1 over the research region. From the figure, we observe the increase in the estimated trend when moving towards northeastern Germany. The color gradient in this figures is similar to that shown by ERA-Interim in Fig. 9. However, the values of the slopes estimated from ERA-Interim and synoptic data are different, which is not surprising -first because of the coarse resolution of the ERA-Interim data and second due to altitude difference, which might result in different trends. In order to justify these results, a data set with a higher spatial resolution than that of ERA-Interim is required.  The same procedure is applied to estimate the trends from temperature and dew point temperature time series. The estimated temperature trends from surface measurements at 227 stations shown in Fig. 13a fluctuate in the range of 0.39 ± 0.1 K decade −1 . In Fig. 13c the trend esti-mated for dew point temperature time series, in the range of 0.48 ± 0.11, are shown. We calculated the change in PWV per 1 • Celsius rise in temperature, and the results are shown in Fig. 12 Figure 12. The increase in atmospheric water vapor pressure and PWV per 1 • Celsius rise in temperature using data at 227 stations. Figure 13. Estimated temperature trends using surface measurements of (a) temperature, (c) dew point temperature, and the corresponding standard error of the estimated slope (b, d). and 6.5 %, which is comparable to the theoretical rate of the Clausius-Clapeyron equation. Also, the estimated trend in PWV is correlated with that from the dew point temperature, which is exhibited in Fig. 13c. Using ERA-Interim temperature and dew point temperature leads to the same observation; however, the trend values are slightly different. We also observed that the trends of dew point temperature are almost in the same range as those for PWV, which makes time series of the dew point temperature proper to provide reasonably adequate information about PWV trends.

Conclusions
In this paper, we aimed at estimating climatic trends from GNSS-based precipitable water vapor time series and surface measurements of air temperature in Germany. First, we compared PWV time series obtained from GNSS and ERA-Interim, which show strong correlation with an average bias of −1.7 mm (ERA-Interim-GNSS) and standard deviations of 2.63 mm.
By comparing the GNSS-based PWV with those from ERA-Interim, the results show small bias values in flat terrain, while the bias increases in mountainous regions. This is mostly caused by the coarse spatial resolution of the ERA-Interim data and hence the inability to properly represent the topography.
To evaluate the temporal evolution of PWV and temperature, we modeled the time series with an additive model that contains trend, seasonal, and stochastic irregular components. The time series are seasonally adjusted to remove the periodic signal, and the trend component is then analyzed after filtering out the irregular component caused mainly by weather variations. The comparison of this method with the Theil-Sen estimator shows insignificant differences in the estimated trends. The GNSS-based estimated PWV trends change between −1.5 and 2.3 mm decade −1 for time series that are 10 to 19 years long. Since the PWV time series at different GNSS sites are not concurrent, we could not draw specific conclusions about the mean trend or spatial features of the trend over the whole research region. Therefore, we extended the research to analyze 26year (1991-2016) time series from ERA-Interim and synoptic stations. Using dew point temperature, we could produce PWV time series at 227 stations with a bias below 1.2 mm to the ERA-Interim data. By analyzing time series of 26 years from ERA-Interim and synoptic data, the PWV trends are observed to be positive and in the range of 0.34 ± 0.06 and 0.48 ± 0.13 mm decade −1 , respectively. The ERA-Interim PWV shows lower trend values of the trend. We found that the trends estimated, using 26 years of data for each station, tend to show a positive gradient when moving from southwestern to northeastern Germany. This was observed in ERA-Interim and synoptic data for both PWV and temperature time series.
The increment in PWV varies between 4.5 and 6.5 % per degree Celsius rise in temperature, which is comparable to the theoretical rate of the Clausius-Clapeyron equation. The magnitude of the PWV trend slightly differs from that of the dew point temperature. Hence, we can consider the trends estimated from the dew point temperature as a measure for the PWV trends in case of lack of observations. It would be illuminating to validate the results of this research using a data set that has a higher spatial resolution than the ERA-Interim.
Data availability. In this paper, we used (1) ERA-Interim data, which can be downloaded from the ECMWF servers (http://apps. ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/; Dee et al., 2011); and (2) synoptic data, which can be downloaded from the DWD data center ftp://ftp-cdc.dwd.de/pub/CDC/. The GNSS-based data can be acquired by sending a request to the GFZ.
Competing interests. The authors declare that they have no conflict of interest.