Long GPS coordinate time series: Multipath and geometry effects

Within analyses of Global Positioning System (GPS) observations, unmodeled subdaily signals propagate into long‐period signals via a number of different mechanisms. In this paper, we investigate the effects of time‐variable satellite geometry and the propagation of a time‐constant unmodeled multipath signal. Multipath reflectors at H = 0.1 m, 0.2 m, and 1.5 m below the antenna are modeled, and their effects on GPS coordinate time series are examined. Simulated time series at 20 global IGS sites for 2000.0–2008.0 were derived using the satellite geometry as defined by daily broadcast orbits. We observe the introduction of time‐variable biases in the time series of up to several millimeters. The frequency and magnitude of the signal is dependent on site location and multipath source. When adopting realistic GPS observation geometries obtained from real data (e.g., including the influence of local obstructions and hardware specific tracking), we observe generally larger levels of coordinate variation. In these cases, we observe spurious signals across the frequency domain, including very high frequency abrupt changes (offsets) in addition to secular trends. Velocity biases of more than 0.5 mm/yr are evident at some sites. The propagated signal has noise characteristics that fall between flicker and random walk and shows spectral peaks at harmonics of the draconitic year for a GPS satellite (∼351 days). When a perfectly repeating synthetic constellation is used, the simulations show near‐negligible time correlated noise highlighting that subtle variations in the GPS constellation can propagate multipath signals differently over time, producing significant temporal variations in time series.


Introduction
[2] Geophysical interpretation of GPS coordinate time series most commonly involves the determination of rates of crustal motion and amplitudes and phases of periodic motions caused by seasonal mass loads (for example, atmospheric, hydrological, or oceanic) [Dong et al., 2002;van Dam and Wahr, 1998].Each geophysical parameter of interest derived from long GPS time series is biased at some level by residual systematic error, particular those that manifest as long-period spurious trends or (quasi) periodic signals.In the case of GPS, these systematic errors and their propagation into GPS time series are not yet all well understood.Recent literature has highlighted that longperiod systematic errors may occur in coordinate time series due to one of two primary mechanisms.
[3] First, spurious signals may occur directly due to unmodeled long-period signals, such as satellite antenna modeling errors that propagate differently as the satellite constellation changes [Ge et al., 2005].The second originates in the presumption that all subdaily signals are modeled at the observation level, either completely within the functional model or through partial mitigation from the stochastic model (e.g., elevation-dependent weighting).However, this is not yet the case and residual subdaily (systematic) errors remain which have been shown to propagate into time series [e.g., Penna et al., 2007] due to a combination of different mechanisms [e.g., Stewart et al., 2005].
[4] One well-studied example of how subdaily signals propagate into longer-period signals is the case of an unmodeled subdaily tidal signal.Unmodeled semidiurnal and diurnal signals at certain tidal periods have been shown to propagate into fortnightly, semiannual and annual periods [Penna and Stewart, 2003;Stewart et al., 2005] with admittances exceeding 100% in some cases [Penna et al., 2007].The propagated signal was shown to be highly sensitive to input signal frequency, coordinate component of the unmodeled signal and site location.King et al. [2008] showed that, even after modeling for solid earth tides and ocean tide loading displacements, substantial signal at subdaily periods remained in GPS coordinate time series and these propagate into annual and semiannual periods in 24 h solutions with median amplitudes of ∼0.5mm, but reaching several millimeters at several sites.As indicated, studies of seasonal geophysical loading phenomena [e.g., Blewitt et al., 2001;Wu et al., 2003] are therefore adversely affected, as are estimates of linear velocity, by the presence of such systematic errors.
[5] Systematic errors in least squares solutions, such as used in GPS data analysis, propagate according to the least squares design matrix.Therefore, the propagation mechanism is controlled by the observation geometry (as defined by the receiver location(s), satellite constellation and local obstructions) plus the chosen parameterization of the solution.Parameters typically include (but are not limited to) site coordinates, adjustments to tropospheric zenith delay, atmospheric gradients, phase ambiguity terms (when not fixed to integers) and receiver clocks (depending on the data differencing approach adopted).Temporal variations in the observation geometry or the number or type of parameters estimated will therefore likely produce temporal variations in the propagated signal.
[6] It is notable, therefore, that even in the trivial yet unrealistic case of the GPS antenna location and number of estimated parameters remaining constant with time, the GPS satellite constellation is constantly evolving as satellites are commissioned and decommissioned, or removed from the solution due to satellite eclipse or maneuver.Furthermore, site specific obstructions such as vegetation or man-made structures may change with time, producing a further change in the observation geometry.The consequence is that, even if an unmodeled signal remains completely constant in time, the way in which it will propagate is likely to change with time.If these temporal variations are suitably large, then the ensuing systematic error will likely bias the GPS time series significantly, resulting in erroneous interpretation of geophysical signals such as tectonic velocity, glacial isostatic adjustment, vertical motion of tide gauges or seasonal geophysical loading signals.Furthermore, such errors would degrade the GPS contribution to the International Terrestrial Reference Frame [Altamimi et al., 2007].
[7] In this paper we test this hypothesis, using one source of presently unmodeled subdaily signal: carrier phase multipath, taking "multipath" to mean signal reflections from planar surfaces not part of the antenna itself, as adopted by Georgiadou and Kleusberg [1988].Errors of this type are similar in nature to antenna phase center variation mismodeling and antenna imaging (changes in the antenna phase pattern induced by conducting material in the vicinity of the antenna) [Georgiadou and Kleusberg, 1988], in that they exhibit an elevation dependency.We do not consider these additional sources of systematic error here, yet simply note that the propagation mechanism is likely to be somewhat similar.Early studies investigating multipath found very small effects on geodetic time series, leading to the thought that multipath effects are mitigated by averaging when sufficient observational time spans are used [e.g., Davis et al., 1989;Lau and Cross, 2007a].However, years of experience with continuous GPS time series, combined with analysis advances leading to improved signal-to-noise ratio of long time series, have shown that GPS time series are highly sensitive to GPS hardware changes (including receiver, receiver firmware and antenna), suggesting that multipath may be playing an important role; yet the mechanism for this is not well understood.As a contribution to understanding the effect of carrier phase multipath on multiyear GPS coordinate time series, we perform several trials using simulated and real data, using various simulated carrier phase multipath signals.
[8] We commence in section 2 by introducing the adopted model of carrier phase multipath that we use throughout this paper to perturb a multiyear, simulated GPS coordinate time series.The simulation approach is introduced (section 2.2) before detailing the different satellite constellation configurations adopted to assess the different characteristics of the multipath propagation (section 2.3).First, in order to show the influence of time variability on the propagation, we start with a theoretical constellation that has a fixed orbit repeat time, a constant number of satellites and no obstructions above the elevation mask at each site (section 2.3.1).Second, we adopt the same clear horizon but with a realistic time variable satellite constellation taken from the broadcast orbits, (section 2.3.2).Finally, we detail the most realistic constellation that includes site specific time variable changes to the observation tracking as observed in real data (for example hardware changes and physical obstructions on the horizon, section 2.3.3).
[9] In section 3 we compare and discuss the simulated time series generated using the three constellation configurations, in addition to investigating the effect of changes to the adopted functional and stochastic models within the evolving constellation scenario.In section 4, we provide a comparison against time series computed using real data in a PPP approach with GIPSY 5 software [Webb and Zumberge, 1995], and in section 5 we present two possible mitigation strategies that involve novel weighting strategies of the input observations in order to minimize the time-and geometry-dependent propagation of the multipath signal.

Signal Multipath
[10] The space surrounding an antenna may be subdivided into three regions including the reactive near field in the region nearest the antenna, the radiating near field and the far field out to infinity [e.g., Balanis, 2005].The boundaries of these regions are not sharply defined although criteria have been developed in order to delineate them in practice.Given an antenna with maximum dimension D, and signal with wavelength l, and for D > l [Balanis, 2005] defines the first and second boundaries occurring at distances from the antenna surface, respectively.For a GPS choke ring antenna with D = 0.38 m, l L1 = 0.19 m and l L2 = 0.24 m, then R 1 In this paper we consider multipath sources solely within the radiating near field, or distances in the range ∼0.03 m to ∼1.2-1.5 m.For an examination of the effect of reactive near field sources on GPS time series, see Dilssner et al. [2008], and for the effect of phase center modeling errors on site velocities, see Steigenberger et al. [2009].
[11] To date there is no widely accepted model for multipath, partly due to the complexity of real world GPS antenna environments.One relatively simple model shown to approximate observed multipath has been described by Elosegui et al. [1995], based on the earlier work of Georgiadou and Kleusberg [1988] and Young et al. [1985].This model is based on the assumption that multipath is caused solely by a horizontal reflector at some height, H, below the GPS antenna phase center causing an attenuation, a, of the signal voltage amplitude.So, for a satellite with elevation angle, ", phase bias d L (in units of meters) of the L1 or L2 phase due to multipath may be modeled as [Elosegui et al., 1995] where l is the carrier phase wavelength for the L1 or L2 carrier phase signal.Since most geodetic GPS positioning is performed using the ionosphere-free linear combination (LC) of the raw carrier phase observables (L1 and L2), the LC phase delay can be computed as [12] We show in Figure 1 (left) d L LC for a range of heights above the reflector (H = 0.1, 0.2, and 1.5 m) and for two different attenuation values (a = 0.05 and 0.1).As noted by Elosegui et al. [1995], the effect of changing the height above the reflector is to change the rate at which d L LC varies with elevation (increased rate of variation with increasing height).Within a reasonable range of attenuation values (0.01 to 0.1), a change in attenuation has the effect of approximately uniformly scaling the bias as seen by comparing Figure 1 (left).
[13] Despite this model providing a useful approximation, it is based on geometric ray optics which is not appropriate in the radiating near field.In addition, values of d L LC will be modified by the antenna gain pattern.To improve on this, we adopt a model as developed by T. A. Herring (personal communication, 2009, hereafter denoted HMM) that extends the Elosegui et al. [1995] model through the use of Fresnel equations relating to electric field amplitudes, and attempting to take into account antenna gain properties: with the antenna gain pattern, consisting of the direct (g d ) and reflected (g r ) gain, modeled using a simple modified dipole model expressed as a function of rate of change (G) of the antenna gain with signal zenith angle, such that also, a = Sg r R a is the amplitude of the reflected signal for a given surface roughness (S), with being the Fresnel equation for an electric field perpendicular to the plane of incidence.[14] This depends on refractive indices n 1 and n 2 , with n 1 = 1 for air and n 2 appropriate for the reflective medium.In practice this can be derived from the square root of tabulated dielectric constants, and for concrete this is typically taken as 4 (n 2 = 2).
[15] Here d L LC may be formed analogously to equation ( 2).
[16] For the purposes of this paper we adopt values which approximate the amplitude of the signal in GPS phase residuals (T. A. Herring, personal communication, 2009), namely S = 0.3, G = 1.1 and n 2 = 2, although these are not definitive, and indeed some sites may require a different value of S, as we shall demonstrate.We show in Figure 1 (right) d L LC for S = 0.3 and 0.5 based on HMM.The amplitude of the modeled signals scales approximately linearly with variations in S. Increasing either G to 1.2 or n 2 to 100 result in increased model signal amplitude of about 50%.Only small changes in frequency or phase occur.Two differences to the model of Elosegui et al. [1995] can be identified.First, the signal in HMM is substantially reduced at high elevations, which is in general agreement with the pattern seen in GPS carrier phase residuals.Second, both amplitude and frequency of the HMM is proportional to H, whereas the Elosegui et al. [1995] model does not show the same sensitivity of amplitude to variations in H (compare Figure 1, left versus right).

Simulation Approach
[17] Our GPS simulator is based on the undifferenced GPS observable, and is conceptually similar to that of Santerre [1991].Observation-level biases, such as d L LC , are entered via the "observed minus computed" term (b) of the batch least squares adjustment: [18] The design matrix, A, is defined by the partial derivatives of the functional model with respect to the parameters.The effect of the unmodeled signal on these parameters is reflected in the estimated values x.W is the inverse of the observation variance-covariance matrix.
[19] Initial station coordinates are taken from the ITRF2005 coordinate set and satellite positions are taken from one of three orbit configurations (section 2.3).Satellite orbits and clocks are assumed known and fixed.Unless otherwise specified below, estimated parameters are the corrections to initial station coordinates, tropospheric zenith delays, receiver clocks (normally every epoch) and real valued ambiguity terms (normally one per satellite pass).In these simulations, nonzero values of the estimated parameters represent parameter bias.Correct integer ambiguity fixing may be simulated by simply removing these terms from the design matrix [Santerre, 1991].
[20] This simulator has previously been used by King et al. [2003] to study propagation of tidal signals in subdaily coordinate estimates; the observed biases in output time series were accurately reproduced in the simulator.We have also verified the simulator by reproducing the propagated periodic biases observed in the GIPSY solutions of Penna et al. [2007] to within very small errors.We therefore assume the simulator is capable of reproducing the effects of systematic errors on real GPS solutions (that use an undifferenced observation strategy), but with the advantage of controlling all systematic error sources.
[21] For the tests described here, we used a typical 24 h "observation" session and estimated adjustments to topocentric station coordinates (north, east and up) once per day (i.e., once per session), tropospheric zenith delay parameters once per hour, and used an elevation cutoff angle of 7°.By default we applied uniform weighting to all observations and include a single real value ambiguity parameter per satellite pass (we assess the impact of elevation-dependent weighting and ambiguity fixing throughout section 3).We used one measurement epoch every 300 s (to reduce the computational burden) and considered data over the years 2000-2007 inclusive for a sample of 20 sites (Figure 2) in the International GNSS Service (IGS) [Dow et al., 2005].
[22] By fixing the satellite orbits, clocks and earth orientation parameters to predetermined values we are potentially oversimplifying the problem.However, since multipath is generally highly localized (e.g., monument types, heights above reflectors and reflector conditions are all site specific) they are unlikely to systematically bias parameters with large spatial scales, and hence we argue that a site-by-site analysis is reasonable.This remains to be verified in further work.
[23] Our simulations represent a substantial advance on the work of Elosegui et al. [1995], which was limited in part by the available computational power at the time.First, they performed a simplified adjustment, considering mainly a one-dimensional (vertical) coordinate with, at most, one tropospheric zenith delay term per day and they did not examine the effects of ambiguity fixing.Second, they report on simulations for only a small number of (noncontinuous) days.Third, they reported on the propagation at only a single midlatitude site.Fourth, they used a multipath model not entirely appropriate to the near field.Each of these differences affects the least squares design matrix (or b in the last case), and hence could alter the propagation into long-term GPS time series.2.3.1. Clear Horizon: Constant Constellation [24] In order to provide a non time variable reference orbit, we adopted an artificial configuration that used a fixed 24 satellite constellation, each with a fixed orbit repeat period of 24 h minus 246 s, close to the average GPS satellite orbital repeat time [Agnew and Larson, 2007].To achieve this constellation, we adapted the "perfect GPS orbit" scenario as developed in a simulator implemented by Penna and Stewart [2003].The horizon at each site was assumed clear for this configuration; that is, satellites are observed without obstruction down to the elevation cutoff angle of 7°.

Clear Horizon: Evolving Constellation
[25] To investigate the propagation impacts of a more realistic time variable orbit in our simulations, this "evolving constellation" configuration adopted orbits from the daily broadcast ephemerides, with satellites set as unhealthy eliminated from the simulations.Again, this configuration assumed no obstructions above the elevation mask at each site.
[26] The time variability of this configuration can be clearly seen in Figure 3 (blue) where we show constellation statistics as a function of elevation over time.The minimum, maximum and standard deviation of epoch by epoch satellite elevations (computed in daily bins) are shown for a representative set of sites.Changes in the blue lines reflect the changing constellation geometry over time.The minimum elevation angle is governed by the 7°elevation cutoff used in the simulations.The maximum elevation angle and standard deviation of all elevation angles are governed by the actual satellite constellation as observed from each site location.
[27] The temporal variability of the satellite geometry in Figure 3 is striking.MAW1 exhibits a strong secular term in maximum satellite elevation over time, although this is a feature of all high-latitude sites to some extent.More equatorial sites exhibit high-frequency variation in maximum elevation and low-frequency fluctuations in elevation standard deviation.A clear trend in elevation standard deviation between 2000 and 2002 is also a common feature to many sites.

Obstructed Horizon: Evolving Constellation
[28] The final orbit configuration assessed within the simulations takes into consideration time variable changes that occur due to the influence of local tracking issues, hardware changes and site obstructions above the elevation cutoff angle that are also potentially time variable and highly site specific.We derived this constellation by using only those satellite observations that are used in analyses of real world GPS data in a conventional GIPSY PPP analysis [Webb and Zumberge, 1995;Zumberge et al., 1997].Constellation statistics from this constellation (shown in brown, Figure 3) show considerable variation from the clear horizon, constant constellation configuration discussed previously.Notable differences are observed at sites such as MCM4 and DUBO which exhibit both secular and clear quasi-periodic characteristics.

Clear Horizon: Constant Versus Evolving Constellation
[29] Figure 4 shows the simulation output highlighting the propagation of the HMM d L LC into the height component of three representative sites using H = 0.1 m, H = 0.2 m, and H = 1.5 m, for both the ambiguity free and fixed scenarios.The effect of a time-varying constellation may be seen by comparing the constant and evolving constellation scenarios (Figures 4 (left) and 4 (right), respectively).The three sites are NTUS, an equatorial site; POTS, a midlatitude site; and MCM4, a high-latitude site.Summary statistics including the offset, slope and RMS (white noise) are provided for all sites in Table 1.
[30] On first inspection, the most noticeable difference between these simulations is a clear shift from a stationary time series with high-frequency periodic variability for the constant constellation, toward clear low-frequency variability present within the evolving constellation output.Also present in the later output are occasional transient events.For example, the H = 1.5 m NTUS ambiguity fixed example (Figure 4) shows what could be interpreted as a significant offset in the coordinate time series during late 2000.Closer examination of the time series shows that the offset is ∼1 mm, and occurs over no more than 2-3 days, and possibly as quickly as from one day to the next.This offset is purely an artifact of the propagation of the systematic error introduced into the simulations as a function of the time variable satellite constellation.
[31] In relation to the north and east coordinate components, we observe a similar pattern of variability yet approximately an order of magnitude smaller in magnitude.This finding agrees with intuition given the simulated multipath signal is azimuthally symmetrical and hence one would expect the effects on north and east components to be relatively small.In the main, we show below only Figures 4-11 pertaining to the height component, including the results for north and east components in the supplementary material where warranted.
[32] Returning to the vertical component, and now considering the height bias (offset in Table 1), we observe biases reaching just over 7 mm for all tested values of H, and note that the offset is approximately equal for the constant and evolving constellation cases.Biases for H = 0.2 m tend to be smaller, but the mean bias will change depending on elevation cutoff angle (see below).Table 1 confirms the output for the constant constellation has a zero trend computed over the entire period (2000.0-2008.0).The evolving constellation case shows trends reaching ± 0.16 mm yr −1 in some instances; the effect over shorter periods could be substantially larger, as is suggested by Figure 4, and it will also scale with input multipath amplitude.
[33] Considering next the variability of the simulation output, (RMS in Table 1), we observe a progressive increase in scatter from H = 0.1 m to H = 1.5 m, and from constant to evolving constellations (mean values of 0.01,0.07and 0.15 mm for the constant case, and 0.04,0.21and 0.49 mm for the evolving case, H = 0.1, 0.2, and 1.5, respectively).The RMS shows a small yet interesting correlation with latitude for the evolving constellation case, with the maximum RMS present near the equator, reducing to a minimum at approximately ±50°, and subsequently increasing toward the poles.This pattern is expressed across all values of H, increasing significantly from H = 0.1 to H = 1.5.
[34] The difference between the ambiguity free and fixed cases is interesting for both the constant and evolving constellation scenarios.In the time domain (Figure 4), there appears to be a clear reduction in the energy of the periodic components within the time series when moving from the ambiguity free to fixed cases.To assist in the analysis of these periodic components and potentially time correlated noise structures, we compute global stacks of power spectra using the Lomb-Scargle periodogram [Scargle, 1982] as implemented according to Press et al. [1992] with an oversampling factor of 4. We repeated the analysis with an over sampling factor of 1 and found no significant differences in the expression of various harmonic peaks present in the spectra.In order to aid comparison between different simulations, prior to stacking in the frequency domain, the power spectra for each time series are normalized by an arbitrary variance of 1 mm 2 for each coordinate component.We generate our stacked spectra by computing the median normalized power estimate across each frequency.
[35] The stacked spectra for the constant constellation (Figure 5, left) confirm that these time series are composed of a stationary series with a complex harmonic comb of energy superimposed.The energy within this harmonic comb is maximum for the H = 0.2 m case, and a minimum for the H = 0.1 m scenario.It is interesting that dominant peaks occur at harmonics of a period close to one year.On closer inspection, these harmonics are in fact multiples of 351.4 days, equivalent to the so-called GPS satellite draconitic year (abbreviated here as dy) which defines the average repeat period of the GPS satellites in inertial space relative to the sun.In a study of stacked spectra taken from the weekly solutions of the IGS, Ray et al. [2008] noticed    energy centered at the same fundamental frequency.Ray et al. [2008] were unable to find similar energy present in time series from other techniques, nor geophysical data, and concluded that the signals were likely to be spurious GPS errors occurring as a result of aliasing and/or orbital errors and/or propagation of site-dependent effects such as multipath.Our results shown in Figure 5 confirm that a simple multipath model is capable of propagating in such a way as to excite a very broad harmonic comb whose fundamental frequency is a function of the repeat orbit period of the satellites.However, comparing spectra for the ambiguity fixed and free solutions (Figure 5), it is evident that ambiguity fixing reduces the propagation of multipath signal to these frequencies, as also found by Tregoning and Watson [2009] when dealing with unmodeled tidal signals.We show here that multipath is therefore likely to be contributor to these harmonics, at some level.We return to the effects of ambiguity fixing shortly.
[36] The relative energy of the periodic terms present in Figure 5 is reduced when moving to the evolving constellation scenario.While clearly still present in the spectra (see in particular H = 0.2m in Figure 5 and also the north and east component in the supplementary material), the periodic energy is partially masked by an increased time correlated noise structure not present in the constant constellation scenario.At frequencies lower than 20 cpdy in the evolving orbit case, we observe a spectrum that is between flicker noise and random walk (with a typical spectral index of approximately -1.5 for the evolving constellation case in Figure 5).In a separate study, we have instead propagated white noise as the multipath signal through the simulator using the evolving constellation configuration.The output consists of simply white noise which therefore suggests that temporal correlation in real GPS time series is not introduced by a time-variable constellation alone, but by timevariable satellite constellation combined with unmodeled systematic errors, at least of the kind considered here.This finding is of interest given the observed presence of flicker noise with unknown origin in real GPS time series [e.g., Langbein, 2008;Williams et al., 2004].We note that short baseline studies have revealed noise with a spectral index higher than pure flicker noise [Hill et al., 2009;King and Williams, 2009], similar to what we show here.Multipath of the kind modeled here is therefore likely to be a secondary contributor toward the dominant flicker noise observed in global GPS analyses.Its importance will clearly increase with networks of smaller spatial scale and as models relevant to larger-scale analyses improve.
[37] At frequencies higher than 20 cpdy, the evolving constellation case shows more systematic periodic features with relatively clearly defined peaks surrounded by quite broad cusps of energy, centered on periods of (in order of decreasing power in the H = 0.2 m evolving constellation scenario) ∼14.6 days, ∼7.3 days, and ∼5.4 days (∼24, ∼48, and ∼65 cpdy, respectively).The peak at ∼5.4 days has a particularly broad cusp, similar to a much less energetic cusp surrounding a broad peak at ∼2.4 days.These features are harmonics of the orbit repeat period and are evident in both the evolving constellation and the constant constellation.They persist (with some variation) in kinematic solutions and so cannot be due to propagation of the signals due to the 24 h observation session.Their frequencies are altered if alternative satellite repeat periods are used in the constant constellation simulations.These broad cusps of energy appear to be related to aliasing of the multipath signal as it is sampled from satellites on the same orbital plane.For instance, switching from 600 s sampling to 300 s removed an energetic peak at 2.4 days and reduced signal in the other bands.Noise in these frequency bands in real GPS time series is currently substantially higher than these signals, so their detection is unlikely to be simple, although they may be important contributors to the high-frequency component of their spectra.
[38] The difference in the stacked spectra between the ambiguity free and fixed case is dramatic for both the constant and evolving constellation scenarios (Figure 6).Fixing ambiguities in these simulations has the effect of reducing noise and minimizing the relative energy of the periodic components that dominate the ambiguity free simulations.We do not observe as much reduction in the higher-frequency periodic terms, although the reduction is still substantial.As is typical, ambiguity fixing has a dramatic effect on the east coordinate component, and an important effect on the north coordinate component [Blewitt, 1989] as shown in Figures S1 and S2 of the auxiliary material.1 Further, our results suggest that ambiguity fixing does not only suppress random errors, but fixing has an important role to play in mitigating systematic errors, at least of the type modeled here [see also King et al., 2003].
[39] In relation to the analysis presented by Ray et al. [2008], we note the cumulative IGS solution combination investigated in their study included solutions computed without ambiguity fixing or with only partial fixing.Our results indicate this may have exacerbated the expression of the draconitic periodic components seen within the IGS data set, although it is not clear to the level of contribution from multipath related errors; indeed tidal signals also appear to contribute [Tregoning and Watson, 2009].Homogenously reprocessed time series with high levels of ambiguity fixing are required before the exact influence of multipath and/or tidal effects may be quantified.Tregoning and Watson [2009] have presented such a solution and shown less prominent harmonics in some cases, although the harmonics still clearly exist.Therefore, it is unlikely that multipath of the type modeled here is the dominant source of harmonics reported by Ray et al. [2008] and Tregoning and Watson [2009].Regardless, multipath does exist in GPS time series, and this signal propagates into coordinate time series noise even after ambiguity fixing, as we have shown here.We focus next on the sensitivity of our findings to variations in the functional and stochastic models.

Effect of Changing the Functional and Stochastic Models
[40] We consider from here on only ambiguity free height time series.We focus on the height component since the modeled multipath signals have the most important impact on this component.Furthermore, as NTUS, POTS and MCM4 are representative of the larger set of sites, we primarily show results for these sites only.We show ambiguity free solutions for ease of comparison with GIPSY PPP solutions presented below.For the sake of clarity, we refer to the already presented evolving constellation simulations as our "reference simulations." [41] Figure 7a shows the effect of increasing the reflectivity from S = 0.3 to S = 0.5, as illustrated in Figure 1.Increasing the magnitude of the input reflectivity nearly doubles the input systematic error which in turn results in an apparently uniform near doubling of the coordinate bias across all frequency bands, including the mean bias.This linear relationship begins to break down for small H, but is suitably linear for at least 0.1 < S < 1.0, which enables simple scaling of the S = 0.3 results in our subsequent tests below.In a worst case scenario (without considering changes to other parameters within the multipath model), the subsequent results could therefore be scaled up by a factor of at least 3.3.
[42] Adopting S = 0.3 we now investigate how the biases change with alternative stochastic or functional models.We show the effect of (1) elevation-dependent weighting (Figure 7b), (2) 24 h batch coordinate solutions (Figure 7c), (3) elevation cutoff angle (Figure 8) and ( 4) the removal of clock terms (as occurs in a double difference solution) (Figure 7d) by comparing them in turn with the reference simulation time series.
[43] First, to assess the impact of elevation-dependent weighting against the reference simulation that used a uniform weighting strategy, we apply a commonly used stochastic model in geodetic positioning [e.g., Jin et al., 2005]: where w i is the inverse of the variance of each observation, i.
Figure 7b provides a comparison between simulations using uniform weighting (more saturated colors) and elevationdependent weighting (less saturated colors), considering only the height component.The general pattern is a reduction in mean bias at all three sites for H = 0.1 m and H = 1.5 m when adopting elevation-dependent weighting.The exception is an increase in bias magnitude of 5-10 mm for H = 0.2 m.The variation in the detrended time series is reduced by up to 70% for H = 1.5 m, 0-25% for H = 0.2 m.The change in variability for H = 0.1 m is small, but does increase or decrease depending on the site.
[44] Second, we consider the effect of 24 h batch coordinate solutions.Penna and Stewart [2003] and Stewart et al. [2005] found that in the case of unmodeled periodic signals, 24 h batch coordinate solutions were a dominant propagation mechanism.To examine this effect in the case of unmodeled carrier phase multipath, we repeated the simulations using uniform observation weighting, but estimating site coordinate parameters kinematically, that is every measurement epoch (300 s). Figure 7c compares the kinematic solutions, after applying a 24 h running mean and sampling to one measurement per 24 h to align with the conventional 24 h solution time base.Applying a running mean in this way produces an aliased signal of about 1% of the input signal and is hence negligible.In contrast to the case of propagation of periodic signals, the 24 h batch solutions appear to somewhat dampen the effects of multipath on coordinate time series.This is most obvious for the H = 0.2 m solutions, with smaller scatter but greater bias in the 24 h solutions.Figure 6.Effect of ambiguity fixing on the height time series in the form of stacked (median) power spectra of difference time series (ambiguity free minus ambiguity fixed).See Figure S2 for the north and east components.
[45] Third, we consider the effect of elevation cutoff angle.Elosegui et al. [1995] found substantial variation in coordinate bias depending on elevation cutoff angle, as would be expected when Figure 1.To examine this further, we repeated the simulations using elevation cutoff angle of 15°instead of 7°, presenting the difference in Figure 8.In addition to the mean bias effect discussed by Elosegui et al. [1995], the temporal variability of the time series depends, interestingly, on the elevation cutoff angle.The solutions all show substantial variations in the mean bias (although these would be reduced using elevationdependent weighting).
[46] Fourth, we consider the case where the receiver clock terms are removed from the least squares solutions.This simulation is somewhat analogous to a double-difference solution with multipath at only one of two closely located stations.Figure 7d shows that effect of removing the clock terms is small, with the biggest differences seen at NTUS.The time series continue to exhibit substantial mean and time-variable bias.As a result of the agreement with the reference simulations, it may be that our simulations are also representative of double-difference solutions, but this remains to be tested.
[47] Perhaps the most striking aspect of the simulations shown up to this point is the large magnitude, broadband, temporal variation in the time series of all three coordinate components.The only time-variable component of the simu-lation is the time-variable orbits contained in the broadcast ephemerides; the simulated multipath signal is time constant for each simulation and only a function of signal elevation.This strongly suggests, therefore, that the time-variable satellite constellation is resulting in time-variable systematic error propagation, in contrast to one widely held belief that such signals tend to "average out" [e.g., Davis et al., 1989] and simply produce a time-constant bias [see also Dilssner et al., 2008].It is the temporal variability that we focus on for the remainder of this paper.

Clear Horizon: Evolving Constellation Versus Obstructed "Real" Horizon
[48] To verify that the simulated signals are indeed similar to those observed using real data with the same introduced systematic error source, we analyzed raw RINEX data from the same sites in GIPSY v5 [Webb and Zumberge, 1995] using PPP with JPL legacy fiducial orbits and clocks, and applied the Vienna Mapping Function (VMF1) [Boehm et al., 2006] and nominal hydrostatic zenith delays [Tregoning and Herring, 2006] based on ECMWF as given in the VMF1 grids.We model ocean tide loading displacements with TPXO6.2 [Egbert and Erofeeva, 2002] and solid Earth tides according to IERS2003 [McCarthy and Petit, 2004].As with the reference simulations, and as is typical in GIPSY PPP solutions, a 7°elevation cutoff angle and uniform observation weighting were adopted.For the sake of this comparison, ambiguities were not fixed to integers.We first computed three sets of coordinate time series using the same observation-level biases as presented in section 2.1.In practice, this was achieved by modifying the L1 and L2 phase center variation file, sampling the output of equation (3) in 1°intervals.We then computed a fourth solution per site without introducing any bias.This reference set of time series are then used to difference the perturbed time series against to test the influence of introducing the bias.Since common geophysical and other common mode signals cancel, the remaining signal can be considered to be the effect of the systematic propagation of the unmodeled input signals.
[49] The results for all considered sites are shown in Figure 9. Examining Figure 9 (left) where the GIPSY results are plotted, it is immediately obvious that the effect on site coordinates is generally substantially much larger than in the reference simulations (Figure 9, right).In general, the GIPSY results show larger variation across a wide frequency range, with the notable exceptions of ONSA and POTS, where the GIPSY and simulated time series show similar levels of temporal variability throughout their time series.We attribute these differences to the GIPSY solutions adopting a different satellite geometry caused by site specific obstructions, tracking characteristics and data editing.As discussed in section 2, our "clear horizon, evolving constellation" simulations presume every satellite above the 7°elevation cutoff mask was observed.This rarely occurs in practice due to site obstructions, receiver tracking problems, data editing in the software or the hardware elevation cutoff mask being set to a value higher than 7°( as shown in Figure 3).
[50] To confirm that the difference between the clear and obstructed "real" horizon geometries introduces this change in propagation, we reran the reference simulation using only the satellites present on an epoch-by-epoch basis, in the GIPSY residual files.The output is shown in Figure 9 (right).Comparing Figure 9 (left and right), the simulated and GIPSY coordinate time series now match remarkably well.
[51] The most striking example of the effects of satellite tracking is MCM4 (Figure 9, noting the different scale to the other sites).For much of the time series, annual signals with amplitude >0.01 m are evident in the H = 0.2 m case.These signals end in mid-2006, corresponding to a change in receiver make (from an AOA SNR-12 ACT to an ASHTECH Z-XII3).This matches the timing of the end of annual signals evident in the satellite elevation standard deviation at MCM4 in Figure 3.We examined the effect of perturbing our simulator by white noise (instead of the HMM) and found that while the scatter decreased over time, the output MCM4 time series was stationary without any systematic structure.DUBO, GOUG, GRAS, and NYAL also show annual signals for sections of their time series.Offsets and transient events are also evident in the time series.We recall that these effects are purely due to changing observation geometry in the presence of unmodeled systematic errors.Offsets reach several cm (e.g., KERG, MALI, MCM4, NYAL).All of the spikes seen in the time series are of the same origin suggesting that at least some of the "outliers" present in GPS time series are the result of systematic error propagation.
[52] Of the sites with relatively smoothly varying time series, one site, NKLG, shows a distinct drift in the coordinate time series for H = 0.2 m, of about 1.5 mm over 2000.0-2007.0,equivalent to a bias of 0.2 mm yr −1 .LHAS also shows a drift, equivalent to 0.75 mm yr −1 over 2002.0-2006.0.We observe that most other sites have velocity biases >0.1 mm yr −1 , even after offsets are removed.
[53] Global stacked power spectra of the simulations computed using the obstructed "real" evolving constellation show a marked increase in the noise floor when compared to the clear horizon equivalent (Figure 10, "uniform" weighting).Of significant interest here in the spectra computed using the real constellation is a substantial reduction in variance for the H = 0.1 m scenario versus H = 0.2 m and H = 1.5 m (i.e., comparing Figure 10, right), particularly at the low frequency end of the spectrum.Here, we observe normalized power estimates almost 2 orders of magnitude lower than the other cases, given otherwise identical multipath conditions (notably the reflectivity coefficient, S).
[54] Together, these results suggest that changes to satellite tracking (caused by changes to hardware, firmware, or local obstructions) can be highly detrimental to GPS coordinate time series when unmodeled systematic errors are present in the GPS solutions.Comparing Figures 3 and 9 suggests that it is more than just the minimum elevation cutoff angle that affects the propagation into coordinate time series; the actual distribution of observations across the sky must be temporally stable for the highest precision results.

Evidence of Multipath Propagation in Real GPS Height Time Series
[55] When comparing our simulated time series with actual height time series computed using GIPSY, we expected little similarity due to the presence of geophysical signal and other noise in the real time series, together with an assumption that the employed model of multipath would be overly simplistic for site-specific similarities to emerge; in fact we made no attempt to replicate the multipath at each site.However, surprisingly, we observe some interesting similarities at some sites.In Figure 9, the actual height time series are shown on the left alongside the simulated time series using the realistic observation geometry.The effect of the multipath signal on GIPSY solutions together with the actual GPS height time series (turquoise) for the station computed using GIPSY without any additional perturbation from the multipath signal (note this will include real geophysical signals in addition to spurious components that may result from multipath or any other unmodeled systematic effect).(right) The reference simulation (all observations) as in Figure 4 together with the simulated time series using only the observations used in the GIPSY solutions (i.e., the obstructed "realistic observations" scenario).Note the different scale for MCM4.
[56] The most striking similarity is seen in the MCM4 time series.The MCM4 height time series shows a high level of agreement with the simulated H = 0.2 m (R = −0.75)time series.The LC phase center of the MCM4 antenna is approximately 0.162 m above a flat concrete slab and hence the agreement between observation and simulation is likely not coincidental.Repeating the simulations with S = 1.0 and H = 0.16 m (see dotted line in Figure 1) effectively scaled the signal to give close agreement with the actual height time series (compare turquoise and purple lines in Figure 11), giving further evidence that the MCM4 time series is dominated by combination of multipath and a time-variable observation geometry.NYAL also exhibits some similarities to the height time series, especially during the pre-2002 period.Due to the presence of real geophysical signal and other noise, it is not possible to exclude similarities at the other sites, and we make no further attempt to investigate this here.

Mitigation
[57] As shown in sections 3 and 4, the propagation of unmodeled multipath signals is sensitive not only to the slowly evolving GPS constellation, but potentially more dramatically to changes in satellite tracking at each station.In the absence of a reliable and generally applicable model for carrier phase multipath, methods for mitigating this effect must therefore be considered in order to maximize the value of the existing global GPS data set.
[58] From the tests reported thus far it is clear that unmodeled carrier phase multipath together with a timeevolving satellite geometry can introduce concerning levels of time correlated noise in the time series and possibly bias estimates of geophysical loading and linear velocity.Even in the best case scenario (clear horizon, constant constellation) the noise is small but not negligible.There are two broad categories of methodology to approach mitigating the temporal variability of these errors.Both attempt to do so by improving the temporal stability of the design matrix.The first of these is to modify the functional model by removing specific observations from the least squares solutions.The second, and perhaps more convenient in practice, is to modify the stochastic model.Either way, the desired effect is to provide daily solutions which use observation geometry largely constant in time.Two relatively simple approaches are examined here, both implemented by modifying the stochastic model.
[59] In an initial step, we processed all data for a given site to determine on a daily basis, the number of observations in each of several elevation/azimuth bins.Here we used 2.5°× 10°bins in elevation and azimuth, respectively.In the first approach, we determine the single day D which has the minimum number of observations across all bins N Az,El , referring to this as "weak day" (WD).In the second approach, we chose the single day D with the least observations below 30°elevation, referring to this as "low elevation" (LE).The subsequent steps are identical for both approaches.
[60] First, for each day we mask out observations in 2.5°× 10°bins where none were found on day D. This is done by applying negligible weight to these observations.This ensures a temporal stability to the observation geometry.Second, we reweight all observations in each bin on each day so that N Az,El is equivalent to N Az,El on day D. This accounts for variation in numbers of observations in various elevations/ azimuths with time (e.g., due to variations in the number of satellites in operation and temporal change in horizon Figure 9. (continued) obstruction).We do this step using bins of 5°× 90°.Given equation ( 3) is azimuthally symmetric we could have used a single 360°bin in azimuth, but we chose to test a mitigation approach which may apply also to azimuthal variation in multipath.We then repeated our reference simulations for each of the two reweighting approaches.
[61] In order to compare results from these mitigation strategies, we return to global stacked power spectra for both the clear and obstructed "real" evolving constellation scenarios (Figure 10).Time series metrics are summarized in Tables 2-4.The most striking result here occurs for the H = 0.2 m scenario, for the obstructed "real" constellation configuration (Figure 10, middle right).Here, we observe a dramatic reduction in the time correlated noise behavior when either mitigation strategy is adopted.Similar, although less significant reductions are observed for H = 0.1 m, and to a lesser extent for H = 1.5 m.When assessing the H = 0.2 m solutions using a purely white noise model, we see a mean RMS reduction from 1.63 mm for the uniform case to 0.71 mm and 0.83 mm for the LE and WD reweighting strategies, respectively (Table 3).We also note a clear reduction in the magnitude and scatter of linear trends present in the reweighted simulation output, particularly for the H = 0.2 m case (Table 3).Here, we observe a reduction from a mean trend of 0.47 mm yr −1 (standard deviation 0.78 mm yr −1 ) to 0.05 mm yr −1 (standard deviation 0.18 mm yr −1 ), for the LE weighting strategy.The reduction in noise resulting from the mitigation approaches is also substantial for the other H, particularly at low frequencies.To further asses these reductions we computed the spectral index using CATS software [Williams, 2008].For the ambiguity fixed uniform weighting case, we see mean spectral indices of −1.3, −1.4 and −1.5 for the H = 0.1, 0.2, and 1.5 m cases, respectively.These drop to between −1.1 to −1.0 for both the LE and Figure 10.Effects of reweighting the solutions as a means of mitigation: Global (median) stacked power spectra for the up component, ambiguity-fixed simulations.(left) Clear horizon, (right) obstructed "real" constellation, (top) H = 0.1 m, (middle) H = 0.2 m, and (bottom) H = 1.5 m.Blue is the standard uniform weighted solution, green is the "LE" strategy, and magenta is the WD strategy.See Figures S3 and S4 for the equivalent north and east components, respectively.WD strategies at H = 0.1 and 0.2, yet show marginal changes for H = 1.5 m (a small increase to −1.67 was observed for the LE case).
[62] Interestingly, the time correlated noise reduction as seen in Figure 10 is not as evident for the clear horizon scenarios.For the H = 0.2 m scenario, the effect of the LE or WD strategies is to reduce low-frequency noise, whereas for H = 1.5 m it increases the low-frequency noise.
[63] To illustrate the effect of these relatively simple approaches, we show in Figure 11 the MCM4 simulations using the LE strategy (for MCM4, the WD approach gives similar performance).Comparing the blue and black lines, we see that the previous spurious signal is clearly and dramatically reduced.

Discussion
[64] The multipath model we have adopted is undoubtedly simplistic, yet it is known to at least partially replicate observation residuals (T. A. Herring, personal communication, 2009) and here we have shown close similarity between simulated versus real time series at a few sites (e.g., MCM4,  where we needed to significantly increase the surface reflectivity to reach close agreement in amplitude with a correlation R = 0.73).We therefore note that the model is simplistic yet not overly unrealistic.We have also tested the model of Elosegui et al. [1995] and found that, despite differences in some detail, near-identical time-correlated noise and harmonics occur when using this model.Given the significant difference between these models, our major findings relating to the introduction of time correlated noise are therefore likely to be somewhat invariant to model errors.This invariance also suggests that other elevationdependent error sources may propagate similarly, but this requires further investigation.While we have made no attempt to model multipath at individual site locations, simulations of the sort shown here allow us to explore the potential effects of multipath on time series.For instance, examining the time series generated using modeled multipath (S = 0.3), suggests that multipath could be contributing variability of a few mm to coordinate time series for many GPS sites.This metric scales almost linearly with values of reflectivity (S), hence it is likely that this effect reaches several mm at some sites.
[65] The presence of mean biases is concerning for several applications of GPS.First, in the realization of terrestrial reference frames, colocated space geodetic techniques (SLR/VLBI) cannot be expected to match GPS positions if multipath exists with anything other than very small S, regardless of the quality of the local tie.Indeed, any agreement between techniques would be time-variable based on the specific period of GPS data used.Sites with poor intertechnique tie agreements may benefit from GPS relocation, such as was done at FORT by relocating the antenna on a new (higher) monument [Ray et al., 2007].Second, the biases in height introduced by multipath will likely result in biased estimates of tropospheric zenith delay, and hence biased precipitable water vapor, as also noted elsewhere [Meertens et al., 1997].The magnitude of the zenith delay biases will be approximately half to one third of the height bias [e.g., Vey et al., 2002].Third, the absolute calibration of ocean satellite altimeters (e.g., TOPEX/Poseidon and Jason-1 and -2) is completely dependent on determining an in situ estimate of the sea surface height that is free of significant bias in order to compare with altimeter estimates.This is typically undertaken using GPS buoys processed relative to local GPS base station(s) [e.g., Watson et al., 2004].Any bias at a base station introduced by multipath will therefore directly bias the altimeter calibration estimates.Multiple base stations may help mitigate this.This is setup-dependent, however, and agreement between independent solutions to date suggests the effect is relatively small [Bonnefond et al., 2009].Where absolute height accuracy is required, it is not clear based on these results what form of setup is optimal.Adding to the uncertainty which stems from the various simulation results is the fact that we simulated the effect of an infinitely large, homogenous reflecting surface below the antenna, which will not be the exact case in many environments.
[66] The temporal variability we observe is a concern for all studies making use of GPS coordinate time series, particularly since it is not yet straightforward to independently (or retrospectively) determine multipath levels at any given site.We observe velocity biases with rates of up to 2.6 mm yr −1 over the seven year period considered and often larger over shorter periods, sufficient to bias estimates of tectonic motion, three-dimensional glacial isostatic adjustment and vertical rates at tide gauges.Likewise, the long-period spurious signals present in our output simulations may bias geophysical estimates of hydrological, atmospheric or oceanographic loading signals derived using GPS [e.g., Tregoning et al., 2009].Short baseline studies will provide further constraints on the effect of multipath in real GPS time series [King and Williams, 2009] since local site effects are dominant in very short baselines.
[67] The presence of harmonics of the GPS draconitic year [Ray et al., 2008] in our time series demonstrates that an unmodeled subdaily signal, together with a satellite constellation with a repeat time of approximately 10s less than a sidereal day [Agnew and Larson, 2007;Choi et al., 2004] may contribute to these harmonics.Indeed, this is in agreement of a direct evaluation of the analytical expressions of Stewart et al. [2005] computed using this orbital period and an unmodeled signal of the same period as is the case for multipath (Stewart et al. [2005] adopt an orbital period of one sidereal day as defined within the World Geodetic System 1984).We note that the expression of these harmonic terms is partially mitigated through ambiguity resolution which is an important finding in light of the analysis presented by Ray et al. [2008].
[68] In terms of mitigating these signals, two simple approaches adopted here go some way toward the desired outcome of minimizing the propagation effect by standardizing the observational geometry and reducing the level of time correlated noise behavior in the resultant time series.In the absence of precise site-specific multipath modeling techniques, more elegant approaches may be developed that ensure ever greater temporal stability without loss of precision.It is important, however, to note that while these approaches will minimize temporal variability, they will not remove the mean bias that is present; removal of that requires observation level models, of which conventional elevation-dependent weighting is the most common approach adopted in geodetic software.However, we expect that dayby-day elevation-dependent weighting, even those based on daily postfit residuals (e.g., R. W. King and Y. Bock, Documentation for the GAMIT GPS analysis software, version 10.3, Massachusetts Institute of Technology, Cambridge, 2006) will not be sufficient to reduce the temporal variability in the time series.A suitable method of combining day-specific reweighting while preserving a level of temporal stability is required.
[69] We have not yet explored the potential effects of multipath on other parameters estimates, most notably tropospheric zenith delays, although earlier work [Meertens et al., 1997] has reported biases, as would be expected given the similarity of the elevation dependence of low H (Figure 1) and a tropospheric mapping function.Importantly, we note that the signals we simulate here are not the only subdaily signals present in GPS time series [King et al., 2008].Signals of tidal period (both ocean and atmospheric origins) [e.g., Tregoning and Watson, 2009], those related to accumulated snow and thermal expansion [Jaldehag et al., 1996;Meertens et al., 1997;Penna et al., 2007], plus other signal propagation effects (e.g., mapping function, higherorder ionosphere) and orbit modeling errors, will also be propagated into coordinate time series in a time-varying way by the evolving GPS constellation.The importance of the temporal variability within the propagation mechanism in each of these cases remains to be examined.

Conclusions
[70] We have examined the effects of modeled carrier phase multipath on GPS coordinate time series for 20 global sites spanning 2000 to 2007 inclusive.Using a simple multipath model, and a combination of simulation and GIPSY PPP, we have shown that the time-variable GPS satellite constellation produces time-variable biases in coordinate time series.These are site specific and depend on the solution parameterization and the input multipath signal, with the magnitude of the effect highly scalable according to the surface reflectivity, dielectric properties of the reflector and the antenna gain pattern.Harmonics of the GPS satellite draconitic year were identified in time series generated using a synthetic time-constant satellite constellation as well as that defined by the broadcast orbit that evolves over time.However, whereas the time-constant time series contained otherwise white noise, the evolving constellation time series contained time correlated noise broadly compatible with a model between flicker and random walk.Ambiguity fixing was found to reduce the expression of the harmonic signals throughout each coordinate component time series, and the spectral index moved closer to that of flicker noise.
[71] Considering only the actual observations present in GIPSY PPP solutions at each of the sites resulted in dramatic degradation in the quality of many of the simulated time series.Offsets, quasi-periodic (annual) signals, transient events and trends were all evident in a number of sites, purely due to the presence of an unmodeled multipath signal and a time-variable observation geometry.On the other hand, some sites exhibit little change at all, presumably since they observe a near-complete constellation.Somewhat surprisingly given the simplicity of the multipath model, the height time series for one site in particular, MCM4, was found to show close agreement with a multipath signal with origin ∼0.16 m below the antenna.White noise observation errors, propagated in the same way as the multipath signal, do not yield time correlated systematic error structures such as those observed in this study.We conclude that the elevation dependence of the input systematic error is required in order to excite the time correlated structure within the output time series.
[72] In addition to the reflective near field signals studied in this paper, we note that reactive near field and antenna phase center model errors will also affect station coordinate time series [Dilssner et al., 2008;Steigenberger et al., 2009].Visual inspection of the coordinate time series shown by Dilssner et al. [2008] is suggestive of the presence of time-correlated noise, but this requires further investigation, as does the case of antenna phase center model errors and far field multipath [Larson et al., 2007].
[73] Mean height biases were also evident and varied with the antenna reflector distance and elevation cutoff angle.While the level of bias varied with time and stochastic model, for long-period studies, very low setups (with respect to the reflecting surface) appear preferable since these consistently give the smallest bias with the lowest levels of time-correlated noise.This conclusion only applies however to the model considered here and does not consider potential reactive near field effects.We note that some previous work has also found little evidence of multipath where the antenna base was located flush with the reflective surface [Borsa et al., 2007;Meertens et al., 1997].Ultimately, a combination of calibration for reactive near field effects [Dilssner et al., 2008] and improved modeling of the radiating near field using detailed local knowledge [Lau and Cross, 2007b] may prove to be the optimal strategy for dealing with future measurements.However, given the long history of GPS observations in less than ideal settings (e.g., on large diameter pillars), some ability to mitigate the effect is required, either through modification of the functional or stochastic model.We have shown that, for the values of H we tested, it is possible to reduce time-correlated noise using simple modifications to the stochastic model.
[74] In reality, the model of time-constant multipath applied here will be insufficient at many sites.Multipath sources will not typically be azimuthally symmetric, or constant in time.Given the high scalable effect multipath may have on time series, the development of a more generic and historically applicable model or mitigation approach is highly desirable.
[75] Acknowledgments.M.A.K. was funded by a NERC research fellowship and the research was enabled by a Newcastle University visiting fellowship to C.S.W. C.S.W. was otherwise supported under the Australian Research Council's Discovery Projects funding scheme (DP0877381).We thank the IGS community for their data and products and JPL for making the GIPSY software and satellite products available.Some of the figures were generated using GMT [Wessel and Smith, 1998].We thank Tom Herring for kindly suggesting and providing his multipath model.We thank reviewers Jim Ray and Simon Williams and Associate Editor Tim Dixon for constructive comments which helped improve the manuscript.

Figure 1 .
Figure 1.Carrier phase bias due to reflectors using two different multipath models.(left) The Elosegui et al. [1995] model with two different attenuation values, a, at H = 0.1 m (gray), H = 0.2 m (blue), and H = 1.5 m (magenta), sampled every 1°.(right) The HMM for the same values of H and with two different surface reflectivity terms, S, sampled every 1°.Np is the refractive index of the reflective medium.The dotted black line (bottom right) is for S = 1.0 for H = 0.16 m.

Figure 2 .
Figure 2. Site locations from the IGS network used in this study.

Figure 3 .
Figure 3. Daily elevation angle statistics generated using the broadcast constellation (blue) and those observations actually used in a GIPSY PPP analysis (brown).Units are degrees.Note the different scales for each plot.

Figure 4 .
Figure 4. Height bias time series for three sites when propagating d L LC with H = 0.1 m (gray/black), H = 0.2 m (blue) and H = 1.5 m (magenta) for (left) constant constellation and (right) evolving constellation.Colors match those in Figure 1.Ambiguity-free (lighter shades) and ambiguity-fixed (darker shades) solutions are shown.

Figure 5 .
Figure 5. Stacked (median) power spectra using all 20 sites for the up component.(left) Constant constellation, (right) evolving constellation, (top) H = 0.1 m, (middle) H = 0.2 m, and (bottom) H = 1.5 m.Ambiguity-fixed and ambiguity-free solutions are shown.Frequency is cycles per satellite draconitic year (cpdy, taken here as 351.4 days) as described in the text.See text for the normalization strategy.Colors match those of Figure 4. See Figure S1 for the north and east components.

Figure 7 .
Figure 7. Effect of alternative simulations on heights due to (a) alternative surface roughness of S = 0.5, (b) elevation-dependent weighting, (c) using kinematic coordinate solutions (shown after smoothing to 24 h), and (d) removal of clock terms.In each case, the less saturated colors are the reference time series in Figure 4, and the more saturated colors are the result of the respective simulation variation.

Figure 8 .
Figure 8. Difference in simulated height bias using different elevation cutoff angles (7°minus 15°).Colors are the same as for Figure 1.

Figure 9 .
Figure9.Comparison of height time series using clear and obstructed "real" evolving constellation geometries.(left) The effect of the multipath signal on GIPSY solutions together with the actual GPS height time series (turquoise) for the station computed using GIPSY without any additional perturbation from the multipath signal (note this will include real geophysical signals in addition to spurious components that may result from multipath or any other unmodeled systematic effect).(right) The reference simulation (all observations) as in Figure4together with the simulated time series using only the observations used in the GIPSY solutions (i.e., the obstructed "realistic observations" scenario).Note the different scale for MCM4.

Figure 11 .
Figure 11.Comparison of real MCM4 GIPSY height time series (turquoise) with the reference simulation (H = 0.2, S = 0.3, blue), the modified simulation to best fit MCM4 (S = 1.0,H = 0.16 m, purple), and the simulation using the LE mitigation strategy (H = 0.2, S = 0.3, black).The lines have been arbitrarily offset for the sake of clarity.

Table 1 .
Regression Coefficients, Intercept, Trend, and Residual RMS for Site Time Series Simulated Using Broadcast Orbits With Ambiguities Fixed a

Table 2 .
Regression Coefficients for Three Different Weighting Strategies That Each Use Obstructed "Real" Horizon, Evolving Constellation, Ambiguity Fixed Simulations for Multipath Scenario H = 0.1 m

Table 3 .
As per Table 2 but Using Multipath Scenario H = 0.2 m

Table 4 .
As per Table 2 but Using Multipath Scenario H = 1.5 m