Reprocessed height time series for GPS stations

Precise weekly positions of 403 Global Positioning System (GPS) stations located worldwide are obtained by reprocessing GPS data of these stations for the time span from 4 January 1998 until 29 December 2007. The processing algorithms and models used as well as the solution and results obtained are presented. Vertical velocities of 266 GPS stations having a tracking history longer than 2.5 yr are computed; 107 of them are GPS stations located at tide gauges (TIGA observing stations). The vertical velocities calculated in this study are compared with the estimates from the co-located tide gauges and other GPS solutions. The formal errors of the estimated vertical velocities are 0.01– 0.80 mm yr−1. The vertical velocities of our solution agree within 1 mm yr−1 with those of the recent solutions (ULR5 and ULR3) of the Universit́ e de La Rochelle for about 67– 75 per cent of the common stations. Examples of typical behaviour of station height changes are given and interpreted. The derived height time series and vertical motions of continuous GPS at tide gauges stations can be used for correcting the vertical land motion in tide gauge records of sea level changes.

Abstract. Precise weekly positions of 403 Global Positioning System (GPS) stations located worldwide are obtained by reprocessing GPS data of these stations for the time span from 4 January 1998 until 29 December 2007. The processing algorithms and models used as well as the solution and results obtained are presented. Vertical velocities of 266 GPS stations having a tracking history longer than 2.5 yr are computed; 107 of them are GPS stations located at tide gauges (TIGA observing stations). The vertical velocities calculated in this study are compared with the estimates from the co-located tide gauges and other GPS solutions. The formal errors of the estimated vertical velocities are 0.01-0.80 mm yr −1 . The vertical velocities of our solution agree within 1 mm yr −1 with those of the recent solutions (ULR5 and ULR3) of the Université de La Rochelle for about 67-75 per cent of the common stations. Examples of typical behaviour of station height changes are given and interpreted. The derived height time series and vertical motions of continuous GPS at tide gauges stations can be used for correcting the vertical land motion in tide gauge records of sea level changes.

Introduction
Satellite radar altimetry and tide gauge measurements are the primary techniques for sea level change investigations. Satellite radar altimetry measures absolute sea level using data obtained during the last 35 yr from altimetry satellite missions GEOS3, SEASAT, GEOSAT, ERS-1, ERS-2, GFO, TOPEX/Poseidon, Jason-1, Envisat, Jason-2, CryoSat-2 and recently HY-2A. Over 1750 tide gauge stations located worldwide measure relative sea level providing long time se-ries (totaling longer than 120 yr). Therefore, analysis of tide gauge measurements with the purpose of long-term sea level change research requires a well-defined reference frame. Such reference frame can be realized through precise positions and velocities of Global Positioning System (GPS) stations located at or near tide gauges. One of the purposes of the GPS Tide Gauge Benchmark Monitoring (TIGA) Working Group, former Pilot Project (http://adsc.gfz-potsdam.de/ tiga/index TIGA.html) (Schöne et al., 2009) of the International Global Navigation Satellite System (GNSS) Service (IGS, (Dow et al., 2009)), is to create such a reference frame. The accuracies required by the oceanographic community for sea level studies are about 5-10 mm for station positions and less than 1 mm yr −1 for vertical motions (Schöne et al., 2009). Solutions of positions of GPS stations located near to tide gauges were derived by different research groups in (e.g. Zhang et al., 2008 using relative models for antenna phase centre variations (PCV) and obsolete processing models. However, the switch within the IGS from using a relative to an absolute PCV model mainly affecting the station height, use of new processing software, models and strategies, inclusion of new TIGA GPS stations in the solutions required and made possible a reprocessing within the TIGA project. Within this reprocessing three TIGA analysis centres (a consortium of University of Canberra, University of Tasmania and Australian National University (CTA); German Research Centre for Geosciences (GFZ) and Université de La Rochelle (ULR)) computed global station network solutions. The GFZ and ULR TIGA analysis centres also contributed to the first IGS data reprocessing campaign (IGS repro1, http://acc.igs.org/reprocess.html) and ITRF2008 .
Published by Copernicus Publications on behalf of the European Geosciences Union.

S. Rudenko et al.: Reprocessed height time series for GPS stations
Various authors used GPS measurements to investigate crustal motions at tide gauges and sea level changes regionally (e.g. Buble et al., 2010;Sanchez and Bosch, 2009) and globally (e.g. Wöppelmann et al., 2007Wöppelmann et al., , 2009). Recently, several studies have been made to investigate the influence of various types of systematic errors on station position and vertical velocity estimates. Thus, Munekane and Boehm (2010) investigated troposphere-induced errors in GPS-derived geodetic time series. King and Watson (2010) showed the influence of multipath and geometry effects on long GPS coordinate time series. Effects of azimuthal multipath asymmetry on long GPS coordinate time series were studied by Goebell and King (2011). King et al. (2012) examined the simulated effect of the electromagnetic coupling of a GPS antenna-monument on GPS coordinate time series longer than 2.5 yr. Fu et al. (2012) have found 0.3 mm coordinate differences between solutions using ocean tidal loading (OTL) computed using two different centres of mass: the centre of mass of the Earth system and the centre of mass of the solid Earth. The influence of non-tidal ocean loading effects on geodetic GPS heights has been studied by Williams and Penna (2011). The quality assessment of the recently reprocessed GPS realization of the terrestrial reference frame  shows that the GPS-derived origin is at the centimeter level consistent with the Satellite Laser Ranging (SLR) one with a drift lower than 1 mm yr −1 . From the recently reprocessed ULR solution, Bouin and Wöppelmann (2010) found agreement within 2 mm yr −1 of the tide gauge measurements and vertical velocities for 84 % of the GPS stations co-located with the tide gauges by analyzing 10 yr of continuous GPS data of more than 200 permanent GPS stations distributed worldwide. Global sea-level rise estimates calculated from tide gauge records corrected using GPS data depend on the terrestrial reference frame used. Thus, errors in the reference frame scale rate and origin rate influence the estimated global sea level rise (Collilieux and Wöppelmann, 2011). That is why generation of a stable reference frame containing precise positions of GPS stations located near to tide gauges is very important.
In this paper, we describe the procedure, models used and the results obtained from the analysis of continuous GPS data from a global network of 403 GPS stations for about a 10-yr time span (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007). Vertical velocities of 266 GPS stations with time series longer than 2.5 yr are computed. We compare our solution with the estimates from co-located tide gauges as well as with the GPS-derived vertical velocities from GFZ previous and external solutions.
The remainder of the paper is organized in the following way. The GPS data processing algorithm and the input data, reference frames and models used are described in Sect. 2. The main results of GPS data reprocessing are presented in Sect. 3. The methodology of station vertical velocity computation and some typical and interesting examples of station height changes are given in Sect. 4. Finally, the results ob-2 S. Rudenko et al.: Reprocessed GPS station height time series tions with time series longer than 2.5 yr are computed. We 70 compare our solution with the estimates from colocated tide gauges as well as with the GPS-derived vertical velocities from GFZ previous and an external solutions. The remainder of the paper is organized in the following way. The GPS data processing algorithm and the input data, 75 reference frames and models used are described in Section 2. The main results of GPS data reprocessing are given in Section 3. The methodology of station vertical velocity computation and some examples of station height changes are given in Section 4. Finally, the obtained results are discussed (Sec-80 tion 5), conclusions are drawn and outlook is given.

GPS data processing algorithm
GPS data of the global network of 403 stations covering time span 4 January 1998 -29 December 2007 (GPS weeks 939 to 1459) were analysed using EPOS-Potsdam software (ver-85 sion 7) (Gendt et al., 1994) recently elaborated. The network includes 187 continuous GPS at tide gauges (TIGA) stations and 216 IGS stations. The network was split in two subnetworks of up to 216 globally stations each. The subnetworks are combined to form daily solutions using up 90 to 30 distributed worldwide IGS05 reference stations in the ITRF2005 datum (cluster connectors). A global map of the station network used for reprocessing is given in Fig. 1. The subnetwork of IGS stations was used to estimate GPS satellite orbits, clocks and Earth rotation parameters that were 95 needed for processing GPS data of the TIGA station subnetwork. Daily solutions are combined into three-day solutions by applying orbit continuity constraints (Beutler et al., 1996).
The observation data, reference frames, measurement 100 models and orbit models used are described in Table 1. Terrestrail reference frame was defined in the following way. Initial coordinates of stations present in ITRF2005  were taken from ITRF2005 and estimated for remaining stations. Initial values of station veloc-105 ities were used from ITRF2005, if available, and computed using NNR-NUVEL1A model (McCarthy et al., 2003) for the rest stations. GPS data from a few hundred stations were processed using data processing strategy for huge GNSS global networks 110 (Ge et al., 2006). The following parameters are estimated in the least-square adjustment. The X, Y, Z station coordinates of all stations are estimated weekly using free network strategy with constraints 1-1000 m to a priori values, no station is fixed. Receiver and satellite clocks are solved for 115 at each epoch assuming white noise process. One receiver clock is fixed and used as a time reference. Satellite initial position and velocity, solar radiation pressure scale, y-bias, sine/cosine terms and stochastic impulses (at noon) are estimated for all satellites once per arc. Yaw rate is adjusted 120 for BLOCK II/IIA satellites during shadow crossing. sphere zenith delay coefficients are solved for each station at 1 hour intervals. Troposphere gradients in elevation and azimuth are estimated every 12 hours for each station. Ambiguities are fixed according to (Ge et al., 2005). The x and 125 y Earth pole coordinates and their rates, Length-of-day and GPS satellite phase center offsets are adjusted daily. In a weekly solution, UT1 is fixed for the first day and estimated for the remaining days. week in the range from 0939 till 1459. The number of stations included in the weekly solutions reaches within the time interval up to 300 as it can be seen from Fig. 2, which shows also the rapidly growing observation network in the first years. tained are discussed (Sect. 5), conclusions drawn and outlook provided.

GPS data processing algorithm
GPS data of a global network of 403 stations covering time span 4 January 1998 to 29 December 2007 (GPS weeks 939 to 1459) were analysed using EPOS-Potsdam software (Gendt et al. (1994), version 7) recently elaborated. The global network of GPS stations ( Fig. 1) used for reprocessing was split in two subnetworks. The first subnetwork includes 216 IGS stations; the second subnetwork includes 187 continuous GPS at tide gauges (TIGA) stations: 107 TIGA observing stations and some other new stations. The subnetworks are combined to form daily solutions using up to 30 distributed worldwide IGS reference stations (cluster connectors) using the procedure described in Zhang et al. (2007) and taking into account the global distribution of the reference stations. All available reference stations were used that were available over the time span. The IGS station subnetwork was used to estimate GPS satellite orbits and clocks that were introduced as fixed values when processing GPS data of the TIGA station subnetwork. To stabilize the GPS satellite orbits, daily solutions are combined into three-day solutions by applying orbit continuity constraints (Beutler et al., 1996). Three-day solutions are combined into weekly solutions using the algorithm outlined in Zhang et al. (2007).
The observation data, reference frames, measurement and orbit models used are described in Table 1. The terrestrial reference frame was defined in the following way. The IGS05 (being an IGS realization of ITRF2005  for GPS stations) was used as a priori terrestrial reference frame. The initial coordinates of stations present in the IGS05 were taken from the IGS05 and estimated for remaining stations. Coordinates of the reference stations were estimated with tight constraints to their initial values, and loose constraints were used for all other stations. The initial values of station velocities were used from the IGS05, if available,

Measurement models
Satellite antenna to centre of mass offsets spacecraft-specific z-offsets and block-specific x-and y-offsets from file igs05.atx (Schmid et al., 2007) Phase centre variations (PCV) absolute model PCV for receiver and satellite file igs05.atx (Schmid et al., 2007) Antenna radome calibrations applied, if given in file igs05.atx (Schmid et al., 2007); otherwise the radome effect is neglected and instead standard antenna model (radome => NONE) is used GPS satellite attitude model GPS satellite yaw attitude model (Bar-Sever, 1996) based on nominal yaw rates RHC phase rotation correction phase wind-up applied (Wu et al., 1993) Marker to antenna eccentricity dN, dE, dU eccentricities from site logs applied Troposphere modelling empirical Global Pressure and Temperature (GPT) model (Boehm et al., 2007), Saastamoinen "dry" and "wet" model for zenith delay, Global Mapping Function (GMF) (Boehm et al., 2006 and computed using the NNR-NUVEL1A model (McCarthy et al., 2003) for the remaining stations. GPS data from up to 300 stations per week were processed using the data processing strategy for huge GNSS global networks (Ge et al., 2006). The following parameters are esti-mated in the least-square adjustment. The Cartesian station coordinates are estimated weekly using free network strategy with constraints 1-1000 m to a priori values; no station is fixed. Receiver and satellite clocks are solved for at each epoch assuming white noise process. One receiver  sphere zenith delay coefficients are solved for each station at 1 hour intervals. Troposphere gradients in elevation and azimuth are estimated every 12 hours for each station. Ambiguities are fixed according to (Ge et al., 2005). The x and 125 y Earth pole coordinates and their rates, Length-of-day and GPS satellite phase center offsets are adjusted daily. In a weekly solution, UT1 is fixed for the first day and estimated for the remaining days. week in the range from 0939 till 1459. The number of stations included in the weekly solutions reaches within the time interval up to 300 as it can be seen from Fig. 2, which shows also the rapidly growing observation network in the first years. clock is fixed and used as a time reference. Satellite initial position and velocity, solar radiation pressure scale, y-bias, sine/cosine terms and stochastic impulses (at noon) are estimated for all satellites once per arc. Yaw rate is adjusted for BLOCK II/IIA satellites during shadow crossing. Troposphere zenith delay coefficients are solved for each station at 1 h intervals. Troposphere gradients in elevation and azimuth are estimated every 12 h for each station. Ambiguities are fixed according to Ge et al. (2005). The Earth rotation parameters (ERPs), namely, x and y Earth pole coordinates and their rates and length of day (LOD), are adjusted daily. In a resulting weekly solution, Universal Time UT1 is fixed for the first day of the week and the daily estimates for UT1 and LOD are constrained to obtain a continuous Earth rotation throughout the week.

Results of GPS data reprocessing
The GFZ TIGA GT1 solution contains weekly coordinates of GPS stations, daily values of x and y Earth pole coordinates and their rates and LOD. The solution is available in the Solution (Software/technique) Independent Exchange (SINEX) format via anonymous FTP at TIGA archive (ftp://ftp.gfz-potsdam.de/pub/transfer/kg igs/igstiga/ solutions/) as files /wWWWW/gftWWWW7.snx.Z and at the Crustal Dynamics Data Information System (CDDIS) (ftp://cddis.gsfc.nasa.gov/gps/products/WWWW/repro1/) as files gt1WWWW7.snx.Z, where WWWW stands for GPS week in the range from 0939 till 1459.
We have up to 300 stations included in the weekly solutions ( Fig. 2) for the time interval we considered in this study. This figure also illustrates the rapidly growing observation network in the first years.
As an indicator of relative solution stability, the coordinate repeatabilities of daily solutions with respect to the weekly solution were calculated. The averaged values per week over all stations are depicted in Fig. 3. It is clearly visible that the north and east coordinate components improve over the entire time span and reach even 1 mm level. The up component looks relatively stable from 2001 onwards and reaches      Vertical trends were determined from the GT1 SINEX files 175 by extracting the vertical coordinates and converting them to longitude, latitude and height using the WGS84 geoid model. The time series were then fitted to extract a linear trend using the standard deviation values as reciprocal weights to account for measurement errors. Trend changes were deter-180 mined using the BFAST package (Verbesselt et al., 2010). The number of breakpoints (minimum sequence length) was adjusted so as to obtain reasonable estimates of the trend changes, mirroring the assumed underlying processes. Antenna changes and other events influencing the verti-185 cal trend component were taken from the GPS log files for the respective stations. Some of the GPS stations examined here are located at tide gauges. Here, the trends from the tide gauges were compared with the GPS trends and sea level radar altimetry data from the TOPEX mission to separate the 190 origin of the trend signal where possible. TOPEX/Poseidon sea level anomaly data was provided by Saskia Esselborn, GFZ; a seasonal component was extracted using the Loess algorithm from the STL package (Cleveland et al., 1990).

Plate tectonics: Neah Bay NEAH, Canada
Neah Bay, Washington lies on the Juan de Fuca Strait in the Cascadia Subduction zone, which is subject to crustal uplift as the North American Plate is shifted over the Juan 215 de Fuca Plate. The GPS time series shows a secular trend of 3.96 ± 0.02 mm/yr (4.0 ± 0.0 mm/yr with a seasonal signal removed) (Fig. 6). This is in accordance with the results obtained by Verdonck (Verdonck, 2006), who arrives at a 4.0 ± 0.1 mm/yr trend through the analysis of 220 tide gauge data. Sea level from altimetry yields a trend of standard deviation of solution residuals weighted average is shown in Fig. 4. The accuracy of North and East components is for the investigated time span of about 1 mm and the Up 160 component reaches about 3 mm at the end. Based on these results it can be stated that the required position accuracies are fullfiled needed for determinantion of accurate station height time series and that a precise reference frame was realized.
The daily adjusted values of x and y Earth pole co-165 ordinates and Lenght-of-day (LOD) are compared to the combined solution of the first IGS Reprocessing campaign (IG1). The good agreement between the two solutions in x and y pole coordinates with mean and standard deviation of -0.012 ± 0.056 and -0.026 ± 0.049 mas and LOD with origin of the trend signal where possible. TOPEX/Poseidon sea level anomaly data was provided by Saskia Esselborn, GFZ; a seasonal component was extracted using the Loess algorithm from the STL package (Cleveland et al., 1990). de Fuca Plate. The GPS time series shows a secular trend of 3.96 ± 0.02 mm/yr (4.0 ± 0.0 mm/yr with a seasonal signal removed) (Fig. 6). This is in accordance with the results obtained by Verdonck (Verdonck, 2006), who arrives at a 4.0 ± 0.1 mm/yr trend through the analysis of 220 tide gauge data. Sea level from altimetry yields a trend of Based on these results it can be stated that the required position accuracy needed for the determination of accurate station height time series is fulfilled and that a precise reference frame was obtained.
The daily adjusted values of x and y Earth pole coordinates and length of day of the GT1 solution are compared to those of the IG1 combined solution. The good agreement between the two solutions in x and y pole coordinates with mean values computed using original daily values and standard deviation (scatter in time series of the ERP differences) of −0.012 ± 0.056 and −0.026 ± 0.049 mas and LOD with 0.000 ± 0.023 ms d −1 can be seen in Fig. 6.

Station vertical velocity computation methodology
Vertical trends were determined from the GT1 SINEX files by extracting Cartesian station coordinates and converting them to longitude, latitude and height using the WGS84 geoid model. The time series were then fitted to extract a linear trend using the standard deviation values as reciprocal weights to account for measurement errors. Trend changes were determined using the Breaks For Additive Seasonal and Trend (BFAST) package (Verbesselt et al., 2010). The BFAST algorithm uses a four-step iterative procedure to detect breakpoints in time series. First, the ordinary least squares (OLS) residuals-based moving sum (MOSUM) test is used to detect whether breakpoints do occur in the time series. If the test indicates a significant change (at the significance level p < 0.05), the breakpoints are estimated from the seasonally adjusted data. In the second step, the trend is estimated using robust regression. The OLS-MOSUM test is then applied again in the third step to test for breakpoints in the seasonal component of the time series. In the final step, the seasonal component is estimated from the detrended data. The above steps are iterated until the number and position of breakpoints are unchanged. Details regarding the procedure are described in Verbesselt et al. (2010). The number of breakpoints was adjusted so as to obtain reasonable estimates of the trend changes mirroring the assumed underlying processes.
Antenna changes and other events influencing the vertical trend component were taken from the GPS log files for the respective stations. Some of the GPS stations examined here are located near tide gauges. Here, the trends from the tide gauges were compared with the GPS station height trends and sea level radar altimetry data from the TOPEX/Poseidon mission to separate the origin of the land movement trend signal where possible. All tide gauge data were retrieved from the Permanent Service for Mean Sea Level (PSMSL) (Woodworth and Player, 2003). Linear trends at tide gauges were determined using a standard linear model with ordinary least square adjustment. TOPEX/Poseidon sea level anomaly data were provided by Saskia Esselborn, GFZ. A seasonal component was extracted using the Loess algorithm from the STL package (Cleveland et al., 1990).
No atmospheric loading corrections were applied to the GPS data, except for BRAZ station, where the hydrological seasonal cycle was the main point of interest, and VAAS and MAR6 stations, for which also hydrological issues were reviewed. That is why the data of BRAZ, VAAS and MAR6 stations were corrected for the atmospheric loading using Green's functions as described in detail in Sect. 4.3.4.
While most stations show consistent, steady height trends, a minority demonstrate deviations caused by antenna or receiver changes. Following the recommendations by Blewitt and Lavallée (2002), stations with the total length of time series shorter than 2.5 yr were not considered for the trend estimation. The vertical velocities of GPS stations located at tide gauges and some IGS stations computed by us are presented in Table 2. The seasonal component was removed from the time series before estimating the trend leading to small error estimates given in Tables 2-3 for our solution.
The time series of height changes of all GPS stations of GT1 solution are available at ftp://ftp.gfz-potsdam.de/pub/ home/ig/nana/GPS station heights/. In the following, a few interesting examples are treated in detail.

Plate tectonics: Neah Bay (NEAH), Canada
Neah Bay, Washington, lies on the Juan de Fuca Strait in the Cascadia subduction zone, which is subject to crustal uplift as the North American plate is shifted over the Juan de Fuca plate. The GPS time series provides a secular trend of 3.96 ± 0.02 mm yr −1 (4.0 ± 0.0 mm yr −1 with a seasonal signal removed), and the tide gauge measurements give a secular trend of −1.8 ± 0.1 mm yr −1 (Fig. 7). This agrees well with the results of Verdonck (2006), who computed an estimated rate of land movement of 4.0 ± 0.1 mm yr −1 . In his estimate, he subtracts a mean eustatic sea level rise 2.0 mm yr −1 from the trend calculated from the tide gauge data (−2.0 mm yr −1 ). These three stations, all located around the Gulf of Bothnia, are subject to land movement processes due to glacial isostatic adjustment (GIA). All stations show large secular up-lift rates with Skellefteå leading at 10.80±0.06 mm yr −1 , followed by Vaasa at 8.66 ± 0.04 mm yr −1 and Mårtsbo (7.62 ± 0.01 mm yr −1 ). The GPS station Skellefteå lies in the vicinity of Furuogrund tide gauge (approximately 11 km distance). The tide gauge trend computed by us using full data and for the case, when the annual signal is removed, is −8.1 ± 0.2 mm yr −1 . So, the residual sea level rise from the sum of the tide gauge and vertical land movement trends is 2.7 ± 0.3 mm yr −1 .

Glacial isostatic adjustment
Nedre Gavle tide gauge situated at approximately 10 km distance from the GPS station MAR6 has a trend of −6.0 ± 0.2 mm yr −1 (1896-1986) computed by us (Fig. 8, upper  panel). Since the gauge stopped operating in 1986, there is no common period with the GPS station. The comparison with the land movement trend at Mårtsbo (7.62 ± 0.01 mm yr −1 ) indicates that a large part of the tide gauge signal can be explained by the land movement through GIA. Assuming that the GIA trend was constant over the last 100 yr, the residual sea level trend from the sum of the tide gauge relative sea level trend and land movement signals is 1.62±0.2 mm yr −1 .
Vaasa tide gauge trend estimated by us is −7.6 ± 0.2 mm yr −1 . The vertical land motion at the co-located GPS station VAAS with atmospheric loading correction applied is 8.6 ± 0.2 mm yr −1 . So, the sum of the tide gauge trend and vertical land motion leads to the absolute sea level rise at Vaasa equal to 1.0 ± 0.2 mm yr −1 .
Despite rather large distance (>340 km) between GPS stations Vaasa and Mårtsbo, the height time series for these stations are strongly correlated at R 2 = 0.93. When reduced for atmospheric loading, the correlation drops only marginally to R 2 = 0.92 (Fig. 8, lower panel). A large part of this is due to the trend. However, even the detrended time series score a R 2 = 0.45 correlation (Fig. 9, lower panel). The deviations between the two series consist mainly in an annual multiweek drift at Vaasa, which is caused by snow cover on the antenna. This figure demonstrates that the drift signals coincide with the periods of increased snow coverage in the area. The snow cover data were taken from Robert Dill's hydrological land surface discharge model (LSDM, Dill (2008)). Its influence on the correlation is checked by decomposing the detrended time series and determining the correlations for the respective components. When the linear trend is removed, the time series correlate at R 2 = 0.73 for the residual nonlinear trend, and R 2 = 0.59 for the full series with the nonlinear trend removed. Only the snow-dominated seasonal cycle shows no correlation at R 2 = −0.09. It is worth to remark that the snow disturbance is correctly removed by the STL algorithm when removing the seasonal cycle. This example illustrates that hydrological errors occurring at a seasonal frequency can be easily and automatically removed when the seasonal cycle is calculated by taking the multi-year monthly mean, instead of fitting a sine. The tide gauge time series exhibit the same strong correlation, with    Esquivel et al. (2006) .   both cases for the common period   (Fig. 9, upper panel).

Glacial isostatic adjustment: Canada (Churchill, CHUR, and Kuujjuarapik, KUUJ)
Both Churchill and Kuujjuarapik lie in the Hudson Bay, which is also subject to the glacial isostatic adjustment. Both GPS stations show strong secular uplift trends: 9.67 ± 0.03 mm yr −1 (CHUR) and 11.52 ± 0.06 mm yr −1 (KUUJ) (Fig. 10, middle panel). Despite their distance (1067 km), the two GPS stations display significant correlation (R = 0.92) of their height time series, also of the detrended ones (R = 0.47, Fig. 10, lower panel). The secular trend in the GPS data is nicely visible also in the Churchill tide gauge data. For 1940, the Churchill tide gauge yields a trend of −9.7 ± 0.1 mm yr −1 (−9.54 ± 0.2 mm yr −1 with a seasonal component removed). It is clear even from the visual inspection of the Churchill tide gauge (Fig. 10, upper panel) that the overall negative regional sea level trend has flattened considerably since the beginning of the 1990s. For the 1998-2007, the trend at Churchill yields −2.8 ± 3.5 mm yr −1 . Removing the seasonal signal reduces the standard error and results in a substantially larger trend estimate of −5.18±0.97 mm yr −1 . This means that, for 1998-2007, the full data trend is around two-thirds smaller than for the whole period of coverage, and about half when the seasonal signal is removed. An explanation for this behaviour can be found in Déry et al. (2011), who mention that there is a notable positive trend in river discharge into Hudson Bay, starting around the 1990s and continuing until 2008, the end of their time series. At the same time, there is a shift in the seasonality of river discharge into the bay. With a positive trend in winter and negative trend in summer streamflow, the time series variance is expected to increase. This explains the large impact of the seasonal signal on the trend for the 1998-2007 time series. The tide gauge signal, which has long been dominated by GIA, is apparently now influenced much more strongly by the changing hydrological processes in Hudson Bay.

Bogotá (BOGT), Colombia
This station is a peculiar example. It is located near the building of the Instituto Geográfico Agustín Codazzi (IGAC) in Bogotá and has been continuously subsiding during the past years. The station was treated in Kaniuth et al. (2001) together with another station, BOGA, which is positioned on the top of the IGAC building. In their treatment, the authors suggest that construction work may have worsened the subsidence processes caused by sedimentation in the area.
The vertical trend at BOGT has almost doubled since it was first estimated in their 2001 paper (−25.2 ± 1.4 mm yr −1 ), now reaching the alarming rate of −44.21 ± 0.19 mm yr −1 (Fig. 11). Obviously, the station cannot be used for any purposes within the TIGA framework.

Trend changes in GPS station data: Arequipa (AREQ), Peru
Arequipa lies in southwestern Peru, in the subduction zone where the oceanic Nazca plate is subducted under the South American plate along the Peru-Chile trench. It is an interesting example of a trend change, as the station recorded the 2001 Arequipa-Camana-Tacna area (M:8.4W) earthquake (Utsu, 1990) and its effects on land movement. The trend change is visible in the time series (Fig. 12)   . The correlation between the two is very high at R 2 > 0.9. The lower panel shows the detrended GPS time series for the co-located stations. The correlation is lowered significantly ( R 2 = 0.45) by the errors caused through snow cover on the antenna at VAAS station. The green shading depicts the modeled mean snow cover (in meters) at the four neighboring grid points surrounding the station in the LSDM hydrology model. with exponential behaviour that can be expected for a postseismic event.

Trend change in tide gauge data: Kushimoto (P208), Japan
Kushimoto (GPS station P208) is located at the southern tip of the Kii Peninsula, Wakayama prefecture. It is situated within the Nankai trough subduction zone and was affected, among others, by the 1946 Nankai earthquake (M = 8) (Utsu, 1990). Sato et al. (1998) mention that, prior to the 1944 earthquake, the southern part of the peninsula is supposed to have subsided with a constant rate of 5-6 mm yr −1 . The Kushimoto tide gauge record indicates a clear trend change coinciding with the minor (M = 5.6) Wakayama prefecture earthquake of 9 May 1987 (Utsu, 1990). The relative sea level trend from the tide gauge before the earthquake is a moderate +0.6 ± 0.2 mm yr −1 , increasing by almost ten times to +6.5 ± 0.2 mm yr −1 after the earthquake. The vertical land movement trend determined from the GPS station P208 over the 2005-2010 period is −6.3 ± 0.1 mm yr −1 . A comparison with the TOPEX/Poseidon (1994-2001) sea level anomaly trend (−1.8 ± 1.7 mm yr −1 ) strengthens the assumption that the major part of the tide gauge trend is caused by the land movement (Fig. 13). Apparently, after the 1987 earthquake, the southern part of the Kii peninsula resumed its original subsiding motion. Still, (Isoda et al., 2004) remarked that the Kushimoto tide gauge is located on the (uplifting) Eurasia plate. Emery and Aubrey (1991) mentioned submergence caused by groundwater extraction as another possible cause for subsidence.

Trend change in tide gauge and GPS data:
Aburatsu (P211), Japan Another example of a trend change is Aburatsu tide gauge located at Nichinan, Miyazaki prefecture, Japan. One trend S.   (Utsu, 1990). An analysis of the height time series for the co-located GPS station P211 shows (Fig. 14

Variations caused by hydrological processes: Brasilia (BRAZ), Brazil
The time series of this GPS station (Fig. 15)  To separate the atmospheric pressure loading effects from the hydrological variations, a pressure correction was applied following van Dam et al. (1994). NCEP 6-hourly surface pressure data were convolved with Farrell's Green's functions using the programme provided on Tonie van Dam's website (van Dam, 2010) and the technique outlined in van Dam et al. (1994).
From 2006 on, a strong negative trend is visible in the GPS time series. This decrease can be also found in GRACE data (Kurtenbach et al., 2009). For the comparison, the GRACE data were resampled to monthly values and filtered, then converted to equivalent water height. The gridded product was provided by Henryk Dobslaw, GFZ. For a qualitative comparison, the GRACE water column data are scaled by factor 150 in Fig. 15. The comparison with water mass data from the LSDM model (Dill, 2008) also shows a strong decrease in water column height. The cause of the subsidence is apparently the drought that hit this region in late summer and fall of 2006 (EM-DAT, 2011).

Stations with zero height trend
Tuktoyaktuk (TUKT) station in northwestern Canada is an example for a site with a height trend very close to zero (−0.01 ± 0.02 mm yr −1 ). Located close to the Alaskan border, the station is apparently affected neither by the positive Fig. 13. The apparent sea level rise (+6.5 ± 0.2 mm yr −1 ) at Kushimoto tide gauge after the 1987 Wakayama prefecture earthquake and its comparison with the height trend (−6.3 ± 0.1 mm yr −1 ) at the co-located GPS station P208.
GIA signal centered on Eastern Canada, nor by the land movement signals caused by the recent ice loss (Larsen et al., 2004) commonly found in south Alaska. The moderate length of the time series (4.3 yr), however, makes this assessment a preliminary one. The seasonal cycle produces a spurious trend (0.5 ± 0.3 mm yr −1 ), if not removed (Fig. 16).
Some more GPS stations with longer time series demonstrate the height trend being very close to zero. These are, for example, GPS station CKIS (Rarotonga, Cook Islands, in free association with New Zealand) with the height trend of 0.04 ± 0.03 mm yr −1 obtained over the 6.3 yr time series; PLO3 (Point Loma, CA, USA) for which the vertical velocity of −0.01±0.09 mm yr −1 was computed at the 8.5-yr time interval; RWSN (Rawson, Argentina) with the height trend of −0.01 ± 0.02 mm yr −1 calculated at the 7.9-yr time span; TUVA (Funafuti, Tuvalu) for which the vertical velocity of −0.01 ± 0.05 mm yr −1 was estimated using the 6.1-yr time series and some other stations.
We have compared the vertical velocities of our solution with the previous GFZ solution (Zhang et al., 2008) obtained by processing GPS data from a global network of 370 GPS stations from January 1994 till December 2006 and containing the vertical velocities of 335 stations aligned to ITRF2000 (Altamimi et al., 2002) and two solutions of Université de La Rochelle: ULR3 solution (Bouin and Wöppelmann, 2010)

Conclusions
We have reprocessed GPS data from a global network of 403 GPS stations over the 10-yr time span (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)   and ULR5) solutions of Université de La Rochelle indicates that the accuracy of the vertical velocities is below 1 mm yr −1 for the most GPS stations of our solution. The height time series for the GPS stations co-located at tide gauges and vertical velocities at these stations can be used to correct the estimates of regional and global sea level changes based on tide gauge data. It is planned to perform a new reprocessing of continuous GPS data at longer time span (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) for an increased number of GPS stations at tide gauges by using newer models and ITRF2008 as a priori terrestrial reference frame within the second IGS data reprocessing campaign. This will allow computation of a longer time series of station coordinates for a larger number of GPS stations with an increased accuracy of vertical velocity.