Ocean tide loading displacements in western Europe: 1. Validation of kinematic GPS estimates

GPS has been extensively used to estimate tidal ground displacements, but the accuracy of this has not been systematically verified. Using more than 20 sites distributed across western Europe, we show that postprocessed kinematic precise point positioning GPS with appropriately tuned process noise constraints is capable of recovering synthetic tidal displacements inserted into real data, with a typical accuracy of 0.2 mm depending on the time series noise. The kinematic method does not result in erroneous propagation of signals from one coordinate component to another or to the simultaneously estimated tropospheric delay parameters. It is robust to the likely effects of day‐to‐day equipment and reference frame changes, and to outages in the data. A minimum data span of 4 years with at least 70% availability is recommended. Finally, we show that the method of reducing apparent coordinate time series noise by constraining the tropospheric delay to values previously estimated in static batch GPS analysis, in fact, results in the suppression of true tidal signals. Using our kinematic GPS analysis approach, periodic displacements can be reliably observed at the 0.2 mm level, which is suitable for the testing and refinement of ocean tide and solid Earth response models.


Introduction
Periodic ground displacements due to solid Earth body tides (EBT) and ocean tide loading (OTL) must be removed from GPS, satellite laser ranging, very long baseline interferometry (VLBI), and Doppler orbitography and radiopositioning integrated by satellite (DORIS) geodetic measurements in order to maximize their extensive use in geophysical studies such as monitoring changes in mean sea level, measuring glacial isostatic adjustment (GIA), elastic displacements due to present-day ice loss, tectonics, Earth rotation and satellite orbits, and reference frame determination. Similarly, satellite laser altimetry measurements over grounded ice (from the CryoSat-2 mission and previously ICESat) must be corrected for tidal displacement, otherwise observed surface elevation change rates will be biased and will result in erroneous estimates of ice sheet thinning and mass balance [e.g., Yi et al., 2000]. EBT effects are routinely modeled, with an often stated [e.g., Yuan and Chao, 2012] accuracy of around 1%, while modeling (predicting) OTL displacement requires both an ocean tide model and a Green's function that describes the Earth's elastic and anelastic response to the tidal load. Errors due to ocean tide models have for many years been the accuracy limitation in OTL predictions [e.g., Baker and Bos, 2003;Bos and Baker, 2005;Penna et al., 2008], with negligible Green's function errors assumed.
GPS has been used to estimate tidal displacements since around 2000 in order to validate both EBT models [Ito et al., 2009] and OTL displacement predictions (and hence different ocean tide models). Besides correcting space geodetic measurements, identifying appropriate ocean tide models is critical for obtaining accurate estimates of geophysical processes from monthly gravity field estimates from the Gravity Recovery and Climate Experiment satellite mission, including GIA, ocean circulation, and seasonal/climatic hydrological and ice sheet mass change. In regions where ocean tide models are inaccurate, GPS estimates of OTL displacement could be used instead of model predictions. VLBI also offers the ability to estimate OTL displacements [e.g., Petrov and Ma, 2003] but is only available at around 35 sites worldwide, whereas there are now thousands of continuous and quasi-continuous GPS sites from which tidal displacements may be estimated, and the number of sites is continually increasing. the harmonic parameter estimation approach. Second, the postprocessed estimation of amplitudes and phases from the harmonic analysis of coordinate time series with high temporal resolution (typically finer than 2 h), hereafter termed the kinematic approach.
The first GPS results of the harmonic estimation approach were reported by Schenewerk et al. [2001], following the method's introduction in VLBI data analyses by Schuh and Moehlmann [1989] and then Sovers [1994]. Their estimated M 2 OTL height displacement amplitudes were quality controlled by comparison with predictions based on the Schwiderski [1980] ocean tide model, with 90% of the 353 global sites (of which over half were in North America) considered agreeing to within 5 mm. Allinson et al. [2004] used a similar approach and reported 0.5 mm M 2 height-amplitude agreements between the harmonic displacement parameters estimated using 1000 days of GPS data and OTL predictions from the FES99 ocean tide model for the UK Ordnance Survey inland station LEED. King et al. [2005] used such GPS harmonic displacement parameters as reference to validate different ocean tide models around Antarctica, and King [2006] reported a 0.8 mm M 2 height vector difference compared with FES99 for AMUN near the South Pole. Thomas et al. [2007] compared GPS, VLBI, and modeled M 2 OTL height values at 25 global (mainly inland) sites and found 1.5 mm RMS agreement between GPS and VLBI, but sub-mm agreement between GPS and OTL predictions using the NAO.99b and TPXO6.2 ocean tide models. Similarly, Yuan et al. [2009] obtained M 2 height RMS agreements of around 1 mm with predictions from the NAO.99b model for 12 sites around Hong Kong, while Yuan and Chao [2012] demonstrated M 2 height RMS residuals of around 0.5 mm for over 600 sites in western USA with spatial coherence, after modeling in the GPS processing both EBTs and OTL using the FES2004 ocean tide model. Yuan et al. [2013] extended this to the global scale, finding M 2 height RMS agreements (with FES2004) of 0.4 mm for 307 inland sites with continental-scale spatial coherence, and suggested that GPS harmonic displacement parameters may potentially be estimated with sufficient accuracy to provide constraints on the internal structure of the Earth. However, common to all these studies is the lack of quality control of the GPS tidal displacements with modeled values; in no case have truly independent accuracy assessments been performed to ensure geophysical interpretation can be carried out with confidence.
The kinematic GPS tidal displacement estimation concept was first demonstrated for relative displacements by Khan and Tscherning [2001] and then Khan and Scherneck [2003], who processed 49 days of GPS data from a baseline from FAIR to CHI3 in Alaska, estimating coordinates every hour and obtaining relative M 2 OTL height amplitudes and phase lags which agreed with predictions from different ocean tide models to 1.1-3.4 mm and 10-15°, respectively. King [2006] extended this relative approach to determine absolute values of OTL height displacement, using coordinate time series generated every 5 min using kinematic precise point positioning (PPP) GPS for AMUN near the South Pole. However, the vector difference of the M 2 height component compared with that estimated using the harmonic parameter estimation method was 1.8 mm in amplitude. Yun et al. [2007] analyzed 7 weeks of relative hourly GPS positions around Korea and found M 2 height relative amplitude RMS errors with respect to NAO.99b, GOT00.2, and FES99 ranged from 1.0 to 4.6 mm, and 6-10°for the phase lags. Melachroinos et al. [2008] analyzed 105 days of data from eight sites in Brittany as part of a European network, estimating hourly relative positions using the GINS software; the resulting M 2 height complex misfits with respect to OTL model predictions were 1.8-4.0 mm, although differences of 10-20 mm arose at one site. Vergnolle et al. [2008] undertook a similar study, using the same European campaign data set (but analyzing data from 14 sites and using the GAMIT software), with RMS M 2 height misfits of 2-4 mm obtained with respect to FES2004 OTL predictions. Ito et al. [2009] analyzed 4 months of Japan GNSS Earth Observation Network System (GEONET) data using kinematic PPP GPS (GPS Tools version 0.6.3 with no EBT or OTL modeled) and formed coordinate time series at a 30 s interval. At the Tajimi GPS station, the M 2 height discrepancy with the modeled EBT plus OTL (NAO.99Jb) was 1% (~1.0 mm) in amplitude and 0.19°in phase. Ito and Simons [2011] used kinematic GPS to infer elastic parameters and structure of the crust/upper mantle of the western United States having observed residual M 2 height signals of 1-2 mm from a yearlong coordinate time series corrected for EBT and OTL. However, they did not present any quality control on their GPS estimates of tidal displacement, with the only previous evidence of attainable quality using GPSTools being Ito et al. [2009].
The disparity in the reported M 2 height displacement accuracy between the harmonic parameter estimation and kinematic GPS approaches, the lower reported accuracies of the kinematic GPS method, and the lack of truly independent testing of the accuracy attainable motivate this study in order to enable the confident geophysical interpretation of tidal displacement from GPS. Using controlled input periodic displacements, we investigate the accuracy with which kinematic PPP GPS may be used to estimate tidal displacement without the reliance on comparisons with necessarily imperfect geophysical models. We assess the displacement accuracy attainable with kinematic PPP GPS, the length of time series required, and the sensitivity of the estimated tidal displacement to time series discontinuities, gaps, and GPS processing settings such as parameter process noise. We therefore seek to determine whether kinematic PPP GPS may be used to infer the elastic and anelastic internal properties of the Earth based on (sub) millimeter level changes in tidal displacement, and the level at which it may be used to validate (and discriminate between) ocean tide models.

Kinematic GPS
where i is the tidal constituent, A k,i and ϕ k,i are the amplitude and Greenwich phase lag, respectively, per component and constituent, ω i is the constituent's angular frequency, and χ i (t 0 ) is the constituent's astronomical argument at reference time t 0 .
Equation (1) represents the harmonic expansion of the tidal potential, which theoretically comprises an infinite number of tidal constituents. Cartwright and Tayler [1971] computed a 505 term harmonic expansion that represents around 99.9% of the total tidal potential [Agnew, 2007] and, of these terms, the lunar semidiurnal constituent with period 12.42 h (M 2 ) usually dominates OTL and EBTs. Hence, real tidal displacement results and examples in this paper will focus on M 2 . We will focus on the height component, whose displacements are typically 3 times larger than those in the horizontal components [Baker, 1984].
If GPS coordinate time series are of sufficient length and temporal resolution (higher sampling than the Nyquist frequency), then harmonic analysis enables the amplitudes and phase lags of the tidal displacement constituents to be estimated. The dominant constituents M 2 , S 2 , K 1 , O 1 , K 2 , and N 2 can, according to the Rayleigh criterion [e.g., Pugh and Woodworth, 2014], be separated from significant constituents that are nearby in frequency with around 1 year or more of data. However, unless at least 18.61 yearlong time series are available and used to estimate the effects of minor constituents simultaneously, nodal corrections (±3% for M 2 amplitude and ±2°for M 2 phase) must be applied to the estimated displacements at lunar frequencies [Agnew, 2007].

Kinematic Precise Point Positioning GPS Key Concepts
The accurate estimation of tidal displacement constituent amplitudes and phase lags via harmonic analysis of high-rate GPS coordinate time series is directly dependent on the quality of the GPS positions. To ensure absolute amplitudes and phase lags are estimated in a global reference frame, the PPP GPS technique, pioneered for postprocessed static positioning by Zumberge et al. [1997] and later extended to kinematic positioning [e.g., Kouba and Héroux, 2001], is preferred. If relative positioning techniques are used, "levering" or reference values are needed, which assumes that the displacement values for at least one site are perfect (e.g., from a tidal model). Kinematic PPP GPS positioning requires dual-frequency carrier phase GPS data and highly accurate satellite orbit and clock data. The receiver coordinates are estimated at high temporal frequency in a batch least squares process or as a random walk time-varying parameter as part of a Kalman or square root information filter. When using a Kalman or square root information filter, the values for the process noise, observational weights, and data sampling interval influence the variation in the state vector that can arise between the previous and current epoch on incorporating the current epoch's measurements. With the ambiguities treated as time invariant and the receiver clock as a white noise parameter, the random walk process noise values assigned to the receiver coordinates and tropospheric delay are critical to the ability of kinematic PPP GPS to accurately measure tidal ground displacements. If a priori EBT and OTL displacement models are applied, few millimeter magnitude residual tidal displacements are expected, necessitating low-process noise values and long data arcs (e.g., 24 h), in order to confidently separate the coordinates, tropospheric delay, receiver clock, and ambiguity terms (which should be resolved to integer values). The coordinates from the back-smoothed long data arcs are then concatenated to form coordinate time series, with the amplitude and phase lag estimated per constituent using harmonic analysis. The back smoothing from session end to beginning ensures the parameter estimates at the start of the session are commensurate with those at the end, i.e., they are not degraded by the PPP forward processing slow convergence times [e.g., Cai and Gao, 2013]. The accuracy with which tidal displacements may be estimated is therefore dependent on the process noise, the time series sampling and length, noise properties, offsets, and data gaps. A further consideration is how well each of the GPS-estimated parameters can be decoupled from each other, in that there is not always a direct propagation of component tidal displacement to the estimated component coordinate. For example, King et al. [2003] showed that height tidal displacements can propagate into the horizontal components when using 1 h batch least squares sampling; Penna et al. [2007] showed that unmodeled horizontal tidal displacement propagates to the height coordinate component when using 24 h GPS sampling; Dragert et al. [2000] showed that unmodeled tidal displacement can be absorbed by the tropospheric delay estimates in 24 h GPS sampling.

Validation Method and Data Set
To independently validate the accuracy with which tidal displacements may be estimated using kinematic PPP GPS, displacements with known amplitude and phase lag are required. An appropriate method is to introduce into the GPS data a synthetic periodic displacement of controlled, known amplitude and phase and assess how accurately the synthetic signal introduced can be recovered from the coordinate time series generated from kinematic GPS processing. Sites everywhere in the world are subject to EBT and OTL displacement (at the poles, the semidiurnal and diurnal constituent EBTs are zero, but OTL is not); and therefore, unless the models are perfect, some residual signal at the major tidal periods will prevail in the coordinate time series and will contaminate any artificial displacement introduced at such a tidal period. Instead, we introduced a synthetic signal with period 13.9585147 h, in the approximate semidiurnal to diurnal range in order to be commensurate with tidal signals, but sufficiently distinct from the major and minor semidiurnal and diurnal tidal constituents so to be readily distinguishable in harmonic analysis. This chosen period is not entirely arbitrary but in fact corresponds to a minor tide with Doodson number (2,À6,2,4,0,0), written without the +5 bias. The potential amplitude of this wave is less than 10 À6 of the M 2 wave and the true amplitude of OTL and EBT displacement at this period is therefore negligible, but the presence of the constituent in standard tidal potential catalogs facilitates its estimation in readily available tidal analysis software such as ETERNA [Wenzel, 1996]. The amplitude of the synthetic signal should be chosen to be commensurate with those residual tidal signals expected in GPS data, e.g., due to errors arising from ocean tide models, OTL Green's functions, or EBT models. Penna et al. [2007] suggest that depending on the global region, M 2 height errors of 0-5 mm in OTL displacements may be expected, based on RMS intermodel agreements. For the EBTs, the dominant M 2 constituent has a maximum amplitude of about 15 cm at the equator, decreasing to about 7 cm at latitude 45°, to zero at the poles, so an assumed 1% EBT model error leads to expected residual M 2 height displacement errors of 0.0-1.5 mm. Other recent studies report anomalous M 2 height signals of around 1-3 mm from comparison with GPS displacement estimates [Ito and Simons, 2011], and spatially coherent M 2 height anomalies of 0.2 mm in the horizontal and 0.4 mm in height reported [Yuan and Chao, 2012].
To be deemed valid and reliable for geophysically interpreting residual tidal displacements of around 0-5 mm, kinematic GPS solutions should satisfy the following.
1. The synthetic signal should be recovered with quantifiable and sufficient accuracy to detect and distinguish between the aforementioned residual tidal signals, i.e., for height better than 0.4 mm. 2. The measurement accuracy should be similar across the 0-5 mm expected residual displacement range.
For geophysical interpretation, it is important to use GPS to estimate residual signals of a few millimeters rather than the entire signal (which can be up to 50 mm amplitude for the M 2 OTL height component), in order to reduce unmodeled errors due to nodal corrections which may amount to a few tenths of a millimeter for the total signal but are negligible for the residual case. 3. Unmodeled ground displacements should propagate only into the kinematic GPS coordinate parameters, the propagation should be 100%, and there should be no "leakage" per orthogonal component (i.e., East displacement should propagate solely into East, North into North, and height into height).

Journal of Geophysical Research: Solid Earth
10.1002/2015JB011882 4. Tropospheric zenith wet delays (ZWDs) should be estimated as reliably as possible and absorb none of the unmodeled ground displacement. These can be quality controlled by comparison with ZWDs from colocated radiosondes and also from daily static GPS processing, which have been extensively validated against radiosondes, water vapor radiometers, numerical weather models, and used for the generation of the International GNSS Service (IGS) tropospheric product [e.g., Byun and Bar-Sever, 2009]. 5. The RMS of the GPS carrier phase residuals should be minimized.
These criteria are collectively considered in the rest of this paper, together with the appropriate process noise, sensitivity to GPS coordinate time series offsets, time series noise, data gaps, and length.

Data Set
In order to develop an optimum analysis strategy, data were first used from four sites. Three GPS sites from Britain, namely Camborne (CAMB), Lerwick (LERW), and Herstmonceux (HERT) were chosen, as they experience almost the complete global range of M 2 OTL displacement over a small region, i.e., height amplitudes ranging from 42 mm at CAMB, to 7 mm at LERW, and 5 mm at HERT. Furthermore, CAMB is located in Cornwall, southwest England, where anomalous M 2 height signals of around 3 mm arise consistently for multiple sites in kinematic PPP GPS coordinate time series, having corrected for OTL and EBTs [Bos et al., 2015]. CAMB, LERW, and HERT were also chosen as they are colocated with radiosondes (two launches per day), enabling GPS-estimated tropospheric delays to be validated, and have continuous data sets spanning at least 4 years, commensurate with the typical length of time series available for the majority of GPS sites across western Europe considered by Bos et al. [2015], thereby providing a representative indication of the tidal displacement accuracy before geophysical interpretation is carried out. The Swiss GPS site Zimmerwald (ZIMM) was also selected since it is an inland site at which OTL effects are small (M 2 height amplitude 7 mm), and any errors in OTL displacement predictions due to ocean tides are expected to be very small; and hence, any residual tidal signals at ZIMM are more likely to be due to EBT errors. No radiosonde data were available for ZIMM. To provide further accuracy assessment of the kinematic method and confirm the derived optimum analysis strategy, an additional 17 GPS sites across western Europe were selected, as shown in Figure 1, with the 21 sites in total forming a sample of the 259 sites considered by Bos et al. [2015]. The additional 17 sites were chosen to provide a sample that included a geographical spread (hence encompassing natural variations in OTL, EBTs, and tropospheric delay) with minimal data gaps over a 4 year time span (so to directly compare and assess the effects of different time series noise). The sample also incorporated sites which had long time series (up to 12 years from 1999.0 to 2011.0) with no known offsets, in order to assess the effect of time series length on the harmonic displacement measurement accuracy. Four years of data were collated from 2007.0 to 2011.0 for all 21 sites (in addition to the 2003.0-2007.0 data span for CAMB, HERT, LERW, and ZIMM), in order to test any sensitivity to data spans covering different years.

GPS Processing and Introduction of Synthetic Displacement Signals
The continuous GPS (CGPS) data were processed using GNSS-Inferred Positioning System (GIPSY) v6.1.2 software in kinematic PPP mode, using a square root information filter as outlined in section 2. Jet Propulsion Laboratory (JPL) "repro1" fiducial final precise satellite orbits (in the IGS08 reference frame), clocks, and Earth orientation parameters were held fixed, estimating receiver site coordinates, ZWDs, east-west and north-south tropospheric gradients (using the gradients definition of Bar Sever et al. [1998]), and receiver clock offsets every 5 min, with ambiguities fixed to integers according to Bertiger et al. [2010]. The VMF1 gridded tropospheric mapping function was used, together with the associated European Centre for Medium-Range Weather Forecasts-based a priori zenith hydrostatic and wet delays [Boehm et al., 2006], with an elevation cutoff angle of 10°. EBTs were modeled according to the International Earth Rotation Service 2010 Conventions [Petit and Luzum, 2010]. OTL displacement amplitudes and phase lags were computed for 11 principal constituents Journal of Geophysical Research: Solid Earth 10.1002/2015JB011882 (M 2 , S 2 , O 1 , K 1 , N 2 , K 2 , Q 1 , P 1 , M f , M m , and S sa ) using the FES2004 [Lyard et al., 2006] ocean tide model, the SPOTL software v3.2 [Agnew, 1997], and Green's functions based on the Gutenberg-Bullen Earth model and then converted from the center of mass of the solid Earth (CE) to center of mass of the Earth system (CM) frame using the values available from http://holt.oso.chalmers.se/loading/ CMC/FES2004.cmc. The OTL displacement was modeled in the GIPSY processing using D. Agnew's hardisp routine (obtainable through Petit and Luzum [2010]), which accounts for companion tides using spline interpolation of the tidal admittances of the 11 principal constituents to form 342 constituents. A synthetic periodic displacement with period~13.96 h was applied when modeling the GPS observations at each site (phase referenced to zero defined at GPS time frame epoch J2000) for each of the East, North, and height components, with amplitudes that varied from 0 to 6 mm depending on the test considered. This was done as for the OTL and EBT displacement modeling, i.e., by correcting the approximate coordinates per measurement epoch. The phase lag of the 13.96 h signal was applied using a convention of lags negative, although except for one set of tests, the phase lag was kept constant as zero. To minimize daybreak edge effects, the data were processed in 30 h data arcs centered daily on 12:00 (GPS time), with the 5 min parameter estimates corresponding to each unique GPS day extracted and concatenated to form time series. The discrete 5 min parameter estimates were then averaged to one value every 30 min to reduce both the time series noise and harmonic analysis computational time, while remaining well below the Nyquist period of the 13.96 h and semidiurnal/diurnal tidal signals.

Process Noise Optimization
To determine the optimum coordinate and tropospheric process noise values for the kinematic PPP GPS estimation of tidal ground displacements of a few millimeters in amplitude, data from CAMB were first considered, as it is located in Cornwall, southwest England, where anomalous M 2 signals in GPS height time series are seen throughout the region [Bos et al., 2015], and it is crucial to ascertain their authenticity before attempting geophysical interpretation. Coordinate and ZWD time series were generated at a 30 min temporal resolution from 2003.0 to 2007.0 for a series of varying coordinate and tropospheric process noise values applied in the GIPSY PPP solutions. Apart from the synthetic displacement inserted, all other models and options were kept fixed, including the use of a carrier phase observation standard deviation of 10 mm. First, the coordinate process noise was fixed at the GIPSY-recommended value for slow-moving platforms, i.e., 0.57 mm/√s, and time series generated, applying in turn a different ZWD process noise value which ranged from 10 À5 through to 10 2 mm/√s, either side of the GIPSY-suggested value of 0.05 mm/√s. Each test was repeated for four different input synthetic 13.96 h periodic displacements, namely 0, 2, 4, and 6 mm amplitudes, applied simultaneously in each of the East, North, and height components.
To quantify the effect of the varying ZWD process noise, the kinematic ZWDs were compared with ZWDs estimated from both static GPS processing and radiosondes. The standard deviation of the differences between the kinematic GPS-estimated ZWDs and those from a "standard" GIPSY 30 h static PPP (estimated every 5 min and then averaged to one value every 30 min) is shown in Figure 2, for the 0 and 6 mm input amplitude displacement cases. The static PPP processing followed that outlined in Williams and Penna [2011], except GIPSY v6.1.2 software and JPL repro1 fiducial orbits and clocks were used, the coordinates were in the IGS08 reference frame, and 30 h not 24 h data arcs were used, while a ZWD process noise of 0.05 mm/√s was used. Also shown in Figure 2 are the standard deviations of the differences between the kinematic GPS-estimated ZWDs and those from twice daily colocated radiosonde launches. It can be seen that the ZWD differences with respect to GIPSY static and radiosonde values are minimized for a tuned tropospheric process noise of 0.10 mm/√s, not for the GIPSY-suggested value of 0.05 mm/√s (which was also that applied in generating the GIPSY static ZWDs). Also shown in Figure 2 is the magnitude of the vector difference [e.g., Penna et al., 2008, equation (4)] between the height component amplitude and phase of the output 13.96 h signal and its controlled, known input amplitude and phase; the amplitude of the residual M 2 height signal; the standard deviation of the coordinate time series height component and the median of the daily RMS values of the carrier phase residuals.
While the resulting time series cannot be considered for geophysical interpretation without also tuning the coordinate process noise, it can be seen from Figure 2 that the tuned ZWD tropospheric process noise also gives a near minimum error (vector difference) in the estimated 13.96 h synthetic output signal of about 0.1-0.2 mm, compared with its controlled known input value. Similar trends are found for the East and North components, but with smaller magnitudes. In other synthetic tests (not shown) for frequencies within the semidiurnal band, we observe frequency-dependent variations in the magnitude of the vector differences across the ZWD process noise range but obtain similar minimum errors. We also observe slight frequency dependent variations, of a factor of 3 or less, in the optimum ZWD process noise but with no appreciable effect on signal recovery (minimum error) compared with the chosen value of 0.10 mm/√s.   Using a coordinate process noise of 0.57 mm/√s, there appears to be negligible absorption of the unmodeled synthetic/M 2 displacement by the ZWDs, for all amplitudes, since the agreement with respect to static GPS and radiosondes is consistent for both the 0 and 6 mm input synthetic signal amplitude cases (and also for the 2 and 4 mm cases which we tested but have not shown). It can be seen that the M 2 height displacement is not minimized using the optimum ZWD process noise, but this is not desired, as the process noise settings are tuned so to accurately detect real signals, rather than to minimize or maximize estimated ground displacement. Similarly, we do not expect to necessarily minimize the RMS of the residuals without also tuning the coordinate process noise, as unmodeled ground displacement will propagate into the residuals.
Having determined an initial optimum ZWD process noise value of 0.10 mm/√s, the tests were repeated in order to obtain the optimum coordinate process noise. The ZWD process noise was now held fixed at 0.10 mm/√s but the coordinate process noise (same value applied to each of the X, Y, and Z components as set in GIPSY v6.1.2) varied from 0.3 × 10 À2 to 0.3 × 10 3 mm/√s, including the GIPSY-suggested slow-moving kinematic value of 0.57 mm/√s. The variations in the same six quantities considered for the ZWD process noise tuning are shown in Figure 3, and it can be seen that the RMS of the carrier phase residuals and the synthetic signal error are minimized using the same coordinate process noise (3.2 mm/√s and looser, although ZWD errors with respect to GIPSY static values increase very slightly for the looser values). This suggests an optimum (actually lower bound) tuned coordinate process noise of 3.2 mm/√s. It is also apparent that, for coordinate process noise values of 3.2 mm/√s and looser, the unmodeled ground displacement propagates to only the coordinates (i.e., not ZWDs also), as the "Ht synth err" values are near zero. Meanwhile, the coordinate standard deviation is not minimized when the synthetic displacement error is minimized, which is to be expected as real displacements (which will increase the scatter) must be detectable. When using the tuned coordinate process noise, an M 2 height amplitude of 2.9 mm arises (East and North both 0.6 mm), but this signal is suppressed and the height standard deviation reduced if tighter process noise values are used, with the undesirable consequence that the unmodeled signal propagates into larger carrier phase residuals. It can also be seen from Figure 3 that the value for the tuned process noise is the same for both the 0 and 6 mm amplitudes of the input synthetic periodic displacement and, as expected, we also found this for the 2 and 4 mm input amplitudes, i.e., the tuned process noise applies across the 0-6 mm amplitude range considered. While an improved agreement with ZWD static GIPSY estimates is obtained on tightening the coordinate process noise from 3.2 mm/√s, this is to be expected and is not a concern (or necessarily desired), as the kinematic solution is simply tending toward the static case as the kinematic coordinate process noise is reduced (although it does not exactly reach the static value because the ZWD process noise used here is different to the GIPSY static default). Likewise, the lack of a clear minimum in discrepancy with respect to radiosonde-derived ZWD is not a concern, because the discrepancy is dominated by random noise. Coordinate process noise (mm/sqrt(s)) CAMB 0 mm amp: 0.10 mm/sqrt(s) ZWD process noise  To check if the initial ZWD optimum process noise of 0.10 mm/√s was still valid, having tuned the coordinate process noise from an initial value of 0.57 to 3.2 mm/√s, the ZWD process noise tests were repeated, holding the coordinate process noise fixed at 3.2 mm/√s. The results obtained are shown in Figure 4, and it can be seen that the synthetic signal and ZWD errors are still minimized when the previously tuned ZWD and coordinate process noise values of 0.10 mm/√s and 3.2 mm/√s, respectively, are used. The coordinate precision (standard deviation) is also minimized, suggesting the optimum tropospheric and coordinate process noise values have been found. It is apparent that the RMS residual errors may be minimized if looser tropospheric process noise values are used, but both the coordinate and tropospheric delays degrade substantially if so. There is also a small increase in the magnitude of the ZWD errors compared with using a coordinate process noise of 0.57 mm/√s, possibly since a looser constraint results in increased coordinate and ZWD correlation. Crucially, however, M 2 displacement amplitudes with the optimum tuned settings remain at East 0.55 mm, North 0.55 mm, and height 2.88 mm, regardless of the amplitude of the input synthetic displacement (the maximum error, i.e., vector difference, obtained is 0.16 mm), strongly suggesting that these tuned coordinate and ZWD process noise values enable accurate 0-6 mm amplitude semidiurnal and diurnal periodic tidal ground displacements to be detected to better than 0.2 mm accuracy, and that thẽ 3 mm amplitude M 2 height signals in southwest England, considered by Bos et al. [2015] are genuine, are not caused by GPS processing error, and may be interpreted geophysically. The tests illustrated for CAMB were repeated for LERW, HERT, and ZIMM, with the same trends and same optimum process noise values found.
All the results presented (and in the rest of this paper) for amplitude, phase, and vector difference values shown were computed estimating 13.96 h and M 2 signals at the exact periods using least squares harmonic analysis, with the 6 mm input displacement scenario amplitudes agreeing to those with ETERNA to 0.05-0.16 mm across the four sites (within the ETERNA error estimates of 0.20 mm). As a further check, the 2007.0-2011.0 data of the 21 sites shown in Figure 1 were processed using the optimum settings, and the amplitudes and Greenwich phase lags of the M 2 residuals were computed and compared with ETERNA estimates (which had median, minimum, and maximum amplitudes of 0.55, 0.23, and 2.97 mm, respectively). Median amplitude agreements of 0.02 mm (minimum 0.00 mm, maximum 0.26 mm) and median phase agreements of 2.9°(minimum 0.1°, maximum 21.0°) were obtained. For the ETERNA computations, its standard high-pass filter for half hour data (n30m30m2.nlf) was used, which removes the long-term drift so no drift parameters were estimated. The Tamura tidal potential [Tamura, 1987] was used, and 18 wave groups were estimated: 10 daily (Q 1 , O 1 , M 1 , P 1 , S 1 , K 1 , Ψ 1 , Φ 1 , J 1 , and OO 1 ); six semidiurnal (2N 2 , N 2 , M 2 , L 2 , S 2 , and K 2 ); and the wave groups containing M 3 and M 4 .

Recovery of Synthetic Periodic Input Displacements
The process noise tests focussed on the ZWD and the coordinate height component, although the respective 0, 2, 4, and 6 mm amplitude 13.96 h synthetic periodic ground displacements considered were applied simultaneously to all three of the East, North, and height coordinate components. The very small 13.96 h signal vector difference errors (shown in Figure 4 as Ht synth err) strongly suggest that the input height displacement propagates to the estimated height coordinate only, with close to 100% admittance (ratio of output amplitude to input amplitude). To test this further, and that unmodeled ground displacement per East, North, and height component propagates only into the respective coordinate components, the optimum ZWD and coordinate process noise runs were repeated for each of CAMB, HERT, LERW, and ZIMM, but with input 13.96 h synthetic period displacement amplitudes of 2 mm in East, 4 mm in North, and 6 mm in height introduced in the same processing run.
The recovered amplitudes and phase lags of the output 13.96 h signal obtained from the harmonic analysis of the coordinate time series are listed in Table 1, together with the errors (vector differences) with respect to the controlled input displacement (denoted in Table 1 as the "simultaneous tropospheric solution"). Also shown in Table 1 are amplitudes, phase lags, and vector differences for the 0 mm 13.96 h synthetic ground displacement case, to provide a quantification of the lower bound of the errors, while the quantities considered in section 4, namely ZWDstatic, ZWDsonde, and the median RMS residuals, are also listed. It can be seen from the vector differences listed in Table 1 that a measurement accuracy of better than 0.20 mm is obtained at all four sites for all coordinate components (East and North errors are typically a factor of two smaller than the height component), and that unmodeled ground displacement per component propagates directly into the matching coordinate component. This confirms that the spurious signals found for static GPS by King et al. [2003] and Penna et al. [2007] do not arise and provides further confirmation that geophysical interpretation of ground displacements with approximate semidiurnal and diurnal periods and amplitudes as small as 0.2 mm may be undertaken with confidence using the kinematic PPP GPS approach. The RMS of the residuals varies between 3.5 and 5.5 mm depending on the site and is the same per site for the case when the RMS residuals 5.12 5.51 3.54 4.68 a Shown for sites CAMB, LERW, HERT, and ZIMM, for 13.96 h synthetic periodic displacement input amplitudes of 0 mm, and then the introduction of 2, 4, and 6 mm in the East, North, and height components, respectively. Also shown are values from the a priori fixed tropospheric solution, in which gradients and ZWDs were not estimated but instead fixed as a priori values from a static GPS solution, in which the displacements were also introduced. All values are in mm, except for phases, which are given in degrees.

10.1002/2015JB011882
synthetic displacement is inserted and when none is inserted, as are the quantities ZWDstatic and ZWDsonde. This provides further confirmation that the unmodeled ground displacement propagates into the coordinates directly and is not absorbed by the estimated tropospheric delay. This can also be seen from Figure 5 ("kinematic" lines), which shows for CAMB the one-sided power spectra for the east-west and north-south tropospheric gradients and ZWD (for the 2 mm in East, 4 mm in North, and 6 mm in height 13.96 h synthetic periodic displacement case), calculated for this figure and in the rest of the paper according to Welch [1967] with seven segments (each comprising 17,520, 30 min data points, i.e., 365 days) that are 50% overlapping. In addition, a Hann window was applied to all of each segment. There are no spikes in the power spectral density (PSD) above the noise level at the 1.719 cycles per day (cpd) frequency of the synthetic periodic input displacement, for any of the gradients or ZWD.

Tropospheric Delay Mitigation
The tests described thus far have comprised simultaneously estimating tropospheric delay and coordinate parameters at each measurement epoch (5 min) and then averaging them to single values every 30 min. However, Geng et al. [2009] found that the proportion of hourly PPP solutions that could be improved using ambiguity fixing (rather than leaving as float) increased if tropospheric delay values estimated from a previous daily static GPS solution were introduced and fixed as a priori values in the 30 min solutions. Using a similar strategy and, in addition, introducing previously estimated east-west and north-south tropospheric gradients as fixed, Reuveni et al. [2012] improved the standard deviation of GIPSY-estimated 30 s kinematic PPP height estimates by nearly 50%, suggesting that improved estimates of subdaily strain will arise. To test whether such a strategy also leads to improved estimates of tidal displacement from kinematic GPS, the 2 mm East, 4 mm North, and 6 mm height input displacement  Figure 5. Power spectra of east-west and north-south tropospheric gradients and ZWDs for CAMB (2003.0-2007.0, 30 min data point averages) for kinematic (simultaneously estimating coordinates, ZWDs, gradients, and receiver clocks every 5 min, shown in blue) and static (estimating one coordinate over 30 h, and ZWDs, gradients, and receiver clocks every 5 min, shown in red) GPS solutions. In each solution, 13.96 h synthetic displacements of 2, 4, and 6 mm have been introduced into the East, North, and height components, respectively.
Journal of Geophysical Research: Solid Earth 10.1002/2015JB011882 runs were reprocessed. In the reprocessing, the data were first processed as daily (30 h, i.e., 3 h overlap at each of the start and end of the day) static PPP ambiguity-fixed solutions, estimating one set of coordinates per day but estimating ZWDs and tropospheric gradients every 5 min. The same 30 h data arc was then processed kinematically as an integer-fixed solution, introducing the previously estimated ZWDs and gradients as fixed a priori values, and estimating only coordinates, together with the (white noise) receiver clock offset, both estimated every 5 min and then averaged to one value every 30 min. This solution is termed the "a priori fixed tropospheric solution," whereas the original solution in which the ZWD, gradients, and coordinates are estimated simultaneously is termed the "simultaneous tropospheric solution." The coordinate power spectra for CAMB are shown in Figure 6 for both the simultaneous and a priori fixed tropospheric solutions. It can be seen that with the a priori fixed tropospheric approach, there are substantial reductions in noise in all three coordinate components for frequencies in the range of about 0.3 to 8 cpd, including at the introduced synthetic displacement frequency of 1.719 cpd. However, it can be seen from Table 1 that for the a priori fixed tropospheric solution, the resulting amplitudes of the known signal for the East, North, and height components are only 1.10, 2.13, and 1.58 mm, respectively, compared with the 2, 4, and 6 mm truth values. Thus, substantial increases in vector difference compared with the truth values arise with the a priori fixed rather than simultaneous tropospheric approach, i.e., 4.4 mm vector differences for the height component for the 6 mm input amplitude case (similar changes arise at LERW, HERT, and ZIMM, as can be seen from Table 1). The use of the a priori fixed tropospheric delays has erroneously oversmoothed the solution, with the consequence that erroneous geophysical interpretation can arise. For example, using the simultaneous tropospheric approach, at CAMB for "simultaneous" (blue lines) and "a priori fixed" (red lines) tropospheric solutions, with 13.96 h synthetic displacements of 2, 4, and 6 mm introduced into the East, North, and height components, respectively. In the simultaneous tropospheric solution, coordinates, ZWDs, gradients, and receiver clocks are all estimated every 5 min, whereas in the a priori fixed tropospheric solution, only coordinates and receiver clocks are estimated.
Journal of Geophysical Research: Solid Earth 10.1002/2015JB011882 a residual 12.42 h period M 2 height signal of 2.9 mm amplitude is obtained, which is explained geophysically by Bos et al. [2015] to arise from deficiencies in the Gutenberg-Bullen Green's function used here, whereas it has been erroneously reduced to only 0.5 mm amplitude in the a priori fixed tropospheric solution. It can therefore be concluded that while coordinate time series noise reductions may be obtained using the a priori fixed tropospheric approach of Reuveni et al. [2012], it results in the failure to detect actual subdaily periodic ground displacement signals of geophysical interest which are present in the time series and is not appropriate.
The failure of the a priori fixed tropospheric approach to detect the 13.96 h periodic height displacement can be explained by first considering the power spectra of the east-west and north-south tropospheric gradients and ZWDs from the daily static GPS solution, which are shown for CAMB in Figure 5, in addition to those from the simultaneous (kinematic) tropospheric approach. The spectra for the daily static GPS solution (in which the coordinate does not vary over the 30 h batch) show how some of the 13.96 h displacement propagates into the tropospheric parameters. Conversely, the spectra for the (kinematic) simultaneous tropospheric solution show no discernible signal above the noise level at the 13.96 h period, i.e., as discussed above the unmodeled height displacement propagates entirely into the estimated heights and not the ZWDs. For the a priori fixed tropospheric solution, the introduction in the kinematic run of the ZWDs and gradients from the 30 h static solution (which has absorbed the synthetic displacements introduced to the data), as fixed values, results in the canceling out of the synthetic periodic displacement. Thus, it is effectively modeled (together with residual M 2 ground displacement), with the erroneous outcome that only very small signals prevail in the kinematic height time series, as can be seen from Table 1

. Time Series Noise and Offsets
While the synthetic signals at the four sites considered thus far for the process noise tuning and error propagation tests have been recovered to better than 0.2 mm, the error in an estimated sinusoid is dependent on the signal-to-noise ratio (SNR), which is dependent on the time series PSD. Therefore, to obtain an indication of the typical and maximum tidal displacement errors for the western Europe sites considered by Bos et al. [2015] and thus ensure confident geophysical interpretation could be undertaken, a sample 17 of their 259 sites was considered, i.e., the sites denoted by a plus sign in Figure 1. These were chosen as they encompassed the full range of time series noise experienced by the 259 sites (determined from the integrated PSD over the spectral band 12.92-15.00 h, i.e., close and either side of 13.96 h but without incorporating any tidal signals); they had 4 years of data during the window 2007.0-2011.0 with no gaps and a geographical spread to ensure any orbital and region-wide phenomena affected all sites similarly. For each site, the GPS data were processed with a 4 mm amplitude 13.96 h synthetic displacement introduced into each of the East, North, and height components, and the phase varied from 0 to 330°, with the vector difference error computed per run and then the RMS computed. Figure 7 shows the RMS error against integrated PSD for the 17 sites, with a high correlation of 0.78, confirming the dependence of the quality of the estimated tidal displacement on time series noise close to the tidal period. Importantly, it can be seen from Figure 8 that the maximum vector difference is only 0.43 mm, while the median vector difference RMS error is 0.14 mm, strongly suggesting that height component geophysical displacements of greater than 0.2-0.4 mm can be measured using 4 years of kinematic PPP GPS data.
We next considered the contribution to the time series noise caused by offsets or daybreak effects. These arise as the GPS data have been processed kinematically in 30 h sessions (with coordinates from the middle 24 h extracted), i.e., with a 6 h overlap per day to minimize daybreak and offset edge effects. Hence, daybreak/offset effects can arise from small differences in the reference frame of the fixed fiducial satellite orbits from day to day if different sites are held fixed or due to equipment changes (antenna or receiver hardware, or receiver firmware). Longer data arcs are theoretically possible when processing GIPSY PPPs which could help to mitigate such daybreak effects, but this is not straightforward. To assess whether any daybreak or offset effect prevails, we simulated a synthetic daybreak effect. Each 24 h constant offset was selected independently from a Gaussian distribution with fixed standard deviation throughout the time series (2007.0-2011.0), again with a 4 mm amplitude 13.96 h signal applied simultaneously to all three coordinate components. We found that for independent daily offsets of 10 mm RMS (i.e., day-to-day offsets up to 14 mm RMS) and 5 mm RMS, the noise in the 13.96 h band of the offset time series was respectively, factors of approximately 5 and 20 less than the noise of the GPS time series themselves and had negligible effect on the recovered synthetic signal. 6.1.1. Tidal Cusps As explained in section 2, the 13.96 h signal was chosen to assess the tidal displacement accuracy as it is commensurate with the semidiurnal to diurnal tidal periods but does not coincide with any dominant tidal constituent of geophysical interest. The accuracy of better than 0.2 mm thus far demonstrated for the 13.96 h signal can be deemed representative of the accuracy of the dominant M 2 tidal constituent provided the time series noise at the two frequencies are similar. An increase in noise at M 2 compared with 13.96 h potentially arises due to tidal cusps [Munk et al., 1965], which are gradual increases in power around the tidal peaks, due to the nonlinear interaction of tides and weather, and are visible in some sea level records. An increased noise level near M 2 would mean a lower SNR and therefore a larger error of the estimated M 2 amplitude and phase lag. Figure 8 shows an enlargement of the height component simultaneous tropospheric solution PSD of CAMB that was shown in Figure 6, for the frequencies near M 2 , and shows a very slight gradual increase in power at M 2 . However, note that near the synthetic 13.96 h (1.719 cpd) signal, a very slight gradual increase in power can also be seen. This strongly suggests that these arise only due to spectral leakage in the Welch method that was used to compute the PSD, because in reality this synthetic signal cannot interact with any other physical signal to produce a tidal cusp. We can conclude that the noise levels around the M 2 and 13.96 h periods are very similar; and that therefore, our predicted accuracy by which the synthetic signal at the 13.96 h period can be estimated is representative of that obtained at the M 2 period.

Length and Gaps
Many sites have intermittent data gaps caused by equipment failure or difficulties in data processing. We tested the sensitivity to this using data from all 21 sites (Figure 1) with the 4 year time series (2007.0-2011.0) as considered above with 4 mm amplitude 13.96 h synthetic periodic displacements input (in each of the three coordinate components). We artificially removed data and assessed the change in vector difference of the output versus input 13.96 h height signal. For each site, we reduced the data availability to ratios of 0.95, 0.85, 0.75, and 0.65 of the full time series, by imposing data gaps of various lengths at random points in the time series until the aforementioned availability ratios were achieved. To simulate longer equipment failures, we used data gaps with random, exponentially distributed lengths with means of 2, 7, 30, and 90 days, whereas to simulate short equipment or data processing failures we used fixed-length data gaps of 30 min and 12 h. Figure 9 shows that, taken over all availability ratios, we observe little impact of the data outage pattern on the recovery error, which remains below 0.2 mm (at 1σ). However, as the availability ratio is reduced tõ 0.65, the recovery error approaches 0.2 mm. We conclude that at availability ratios of~0.7 and higher, our results remain robust for the typical minimum time series length of 4 years. We also tested input amplitudes of 2 and 6 mm and obtained very similar results.
Finally, we tested the effect of data span using the 17 sites with long time series (up to 12 years from 1999.0 to 2011.0) denoted by the plus signs in Figure 1, with input displacements as for the data gap tests. We subsampled the full time series from these sites to various lengths (1.5, 2, 2.5, 3, 3.5, 4, 5, and 6 years) and observed the change in recovered 13.96 h signal (quantified using the vector difference). For each data span tested, 25 start times were chosen randomly, ensuring that the availability ratio over the subsequent chosen interval remained at least 0.75, notwithstanding any naturally occurring outages within the time series. Figure 10 shows that data spans of 4 years or more allow estimation of the 13.96 h signal to better than about 0.2 mm (at 1σ) of the 12 year values, which increase to 0.2-0.4 mm if less than 2.5 years of data are used.

Conclusions
By assessing the ability of GIPSYestimated kinematic PPP GPS 30 min averaged positions to recover controlled input displacements with periods of 13.96 h (similar to the major semidiurnal and diurnal OTL constituent periods) Figure 9. For 21 sites with up to 4 year long height time series (2007.0-2011.0), differences in recovered 13.96 h vector difference between time series with and without data artificially removed by a given proportion ((top) availability ratio) and distribution ((bottom) gap length < 1 day, fixedlength gaps; gap length > 1 day, exponentially distributed gap lengths with mean as indicated). The inverted triangles show results for individual sites, 25 realizations of each combination of availability ratio and distribution, for 13.96 h input signals of amplitude 4 mm. The grey error bars are one sigma. The bold symbol to the left of each group shows the overall mean bias and standard deviation for all realizations at a given availability ratio or distribution. Data span (years) Figure 10. Differences in recovered 13.96 h vector difference between 12 year height time series and those reduced to a given span (1.5, 2, 2.5, 3, 3.5, 4, 5, and 6 years), for 17 sites. The inverted triangles show results for individual sites, 25 realizations of each data span, and for the 13.96 h input signals of amplitude 4 mm. The grey error bars are one sigma. The bold symbol to the left of each group shows the mean and standard deviation for all realizations at a given data span.

10.1002/2015JB011882
and amplitudes of 0-6 mm (commensurate with typical errors in OTL predictions and EBT models), we have shown a periodic vector displacement height and horizontal accuracy of better than 0.2 mm. This was obtained by first extensively testing the accuracies attainable at four western Europe sites using 4 years of CGPS data and tuning the process noise applied in the square root information filter. Then at 17 sites distributed across western Europe with different time series noise, inserted 13.96 h periodic displacements of 4 mm amplitude and phase lags varying from 0 to 330°were recovered with a median RMS vector difference height error of 0.14 mm. These were obtained using an optimum zenith tropospheric delay process noise of 0.10 mm/√s and a coordinate process noise of 3.2 mm/√s, evaluated by collectively minimizing the RMS of the carrier phase residuals, the coordinate time series standard deviation, the zenith tropospheric delay errors (from comparison with radiosonde estimates and those from static GPS processing), and the controlled periodic input displacement error. We confirmed that the errors remained commensurate across the 0-6 mm amplitudes considered. Therefore, kinematic PPP GPS may be used to validate and discriminate between OTL and EBT models at these accuracy levels in geophysical studies. Previous studies [e.g., Yuan and Chao, 2012;Yuan et al., 2013] have reported an RMS~0.4-0.5 mm M 2 height agreement between parameterized GPS harmonic displacements and geophysical models, whereas this study demonstrates both an improved accuracy and also the use of truly independent validation data, i.e., without using any geophysical models whose quality is later desired to be tested. We also obtain a periodic displacement accuracy better than the previously reported M 2 height displacement accuracies obtained using kinematic GPS, which were typically 1-5 mm in amplitude and 0-10°in phase, based on comparisons with GPS harmonic displacements or geophysical models [e.g., Khan and Tscherning, 2001;Khan and Scherneck, 2003;King, 2006;Yun et al., 2007;Melachroinos et al., 2008;Vergnolle et al., 2008;Ito et al., 2009]. With our kinematic PPP GPS method, we found that unmodeled periodic displacement propagates directly into the matching coordinate component with around 100% admittance and without any arising spurious signals as reported for static GPS by King et al. [2003] and Penna et al. [2007].
Besides tuning the process noise values appropriately, we have shown that the harmonic displacement accuracy is sensitive to the time series noise, with an upper bound vector difference RMS height error of 0.4 mm obtained for sites in western Europe. Offsets at the daybreaks of up to 14 mm, arising due to reference frame changes or equipment changes, have negligible effect on the height estimation accuracy, and we have demonstrated that time series of 4 years or longer are advisable for attaining the aforementioned measurement accuracies. If shorter time series are used, for example 2.5 years, amplitude errors increase slightly to around 0.2-0.4 mm. As long as data are available for at least 70% of the time series span considered, there is negligible accuracy degradation on the estimated displacement. We have also demonstrated that it is imperative that the zenith tropospheric delay is simultaneously estimated as a parameter together with the receiver coordinates. If the method of Reuveni et al. [2012] is used, whereby previously estimated ZWDs from a static GPS position solution are inserted to the kinematic PPP solution and held fixed, real periodic displacements will be erroneously removed from the kinematic PPP coordinate time series and hence not detected.
This paper has described how tidal displacements may be estimated using kinematic PPP GPS using the GIPSY v6.1.2 software. If other geodetic GPS softwares are used to apply the method, similar process noise tuning tests should first be carried out, as the method is sensitive to these values and to carrier phase weighting and data sampling. With this appropriate preparatory work, we expect that scientific GPS software should be capable of recovering periodic displacements with 0.2 mm accuracy, sufficient for the validation and refinement of ocean tide and solid Earth response models used to compute ocean tide loading and Earth body tides.