Spatial and temporal correlations of the GPS estimation errors Borys Stoew and Gunnar Elgered Onsala Space Observatory, Chalmers University of Technology, 439 92 Onsala, Sweden Onsala March 16, 2005 This report summarizes the work carried out within Work Package 3000 in the TOUGH Project. It contains the final versions of Delivery D19 on spatial correlation of errors and Delivery D22 on temporal correlation of errors. TOUGH (Targeting Optimal Use of GPS Humidity Measurements in Meteorology) is a project supported by the 5th Framework Programme of the European Commission. Summary We have studied the spatial correlation structure of the estimation errors when the GPS technique is used to derive the atmospheric propagation delay in the zenith direction. These error correlations are parameterized using an analytical function. For spatial scales smaller than 200 km the GPS error correlations become significant; for larger spatial scales, the error correlations are small and slowly decreasing with distance. These results are based on a study of the differences between the propagation delays obtained from a numerical weather model, and the corresponding GPS estimates. The temporal correlations of the propagation delay estimation errors are also examined and the decorrelation times are of the order of 1–2 days. The temporal correlation results are derived for the Swedish and Finnish GPS sites, using the relation between the estimation errors for the vertical coordinate and the zenith propagation delay. Contents 1 Introduction 1 2 Measuring and estimating the ZTD 2.1 Radiosonde data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 GPS data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Microwave radiometer data . . . . . . . . . . . . . . . . . . . . . . . . . . 1 3 3 4 3 Extracting ZTD time series from HIRLAM 3.1 The non-linear observation operator for the ZWD . . . . . . . . . . . . . . 3.2 The tangent-linear observation operator . . . . . . . . . . . . . . . . . . . 3.3 ZTD Innovation vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4 5 5 4 Spatial correlations using analytic models 5 5 Spatial correlations based on NWP and GPS data 7 6 Spatial correlations using vertical coordinate residuals 12 7 Temporal correlations of slant delay residuals using ray tracing 7.1 Hydrostatic and wet mapping functions . . . . . . . . . . . . . . . . . . . 7.2 Ray tracing using radiosonde ascents . . . . . . . . . . . . . . . . . . . . . 13 13 14 8 Temporal correlations of ZTD from GPS position residuals 16 9 Discussion 18 References 19 ii 1 Introduction In short term weather forecasting, various sources of data describing the atmospheric systems and processes are used. The measurement techniques providing information about the atmospheric parameters vary not only in sampling rate and spatial distribution, but also in accuracy. The potential of using atmospheric delays estimated from GPS data in near real time and the high temporal/spatial resolution of the measurements make the GPS data assimilation worthwhile (Gutman et al., 2004). Traditionally, the errors in the radiosonde measurements of humidity are assumed to have no spatial correlation in the numerical weather prediction (NWP) model setup. Spatially correlated errors lead to a non-diagonal structure of the error covariance matrices, whose inverting is computationally demanding. The use of a proper error-covariance model ensures the positive definiteness of the respective covariance matrices, so that they have no eigenvalues near zero and are non-singular. The errors in the estimates of atmospheric propagation delay, both from GPS processing and from NWP models, exhibit spatial correlation structures which deserve special attention (Berre, 1997; Stoew et al., 2001). A statistical description of the temporal correlations of errors in the atmospheric delay estimated from GPS data is needed for the proper operation of an NWP model — in particular, for the weighting of past zenith total delay (ZTD) estimates in an observation bias reduction scheme (Lindskog, 2001; Daley, 1992). The spatial and temporal correlations of the estimation errors are assessed in the following sections of this report. Section 2 makes a brief overview of the relations between relevant atmospheric parameters and the ZTD. We introduce the radiosonde data which we later use to estimate mapping function errors through the method of ray tracing. We also review the estimation of the ZTD using GPS and microwave radiometer data. Section 3 describes how the ZTD time series are obtained from NWP models. Thereafter we focus on the spatial (horizontal) correlation of ZTD errors. In Section 4, we discuss the analytic functions used to model the horizontal correlations of the GPS estimation errors in the ZTD; we provide the resulting parameters in Section 5. Using the residuals of the estimated vertical coordinates, we assess the spatial error correlations in Section 6. The temporal correlation of ZTD estimation errors is studied in Section 7 where we discuss two different atmospheric delay mapping functions (MFs). The residual estimates of the propagation delay, obtained by ray tracing of radiosonde data, are used to estimate the temporal correlations of the errors contributed by the MFs used in the GPS processing. The mapping function is only one component introducing a temporal correlation of the ZTD estimation errors. For example, satellite orbits and clocks also contribute to the effect. An attempt to estimate the total effect is made in Section 8 which presents temporal correlations of the ZTD errors based on the assumption that these estimation errors are strongly correlated with those of the estimated vertical positions of the GPS sites. Section 9 ends the report with a discussion concerning the obtained results. 2 Measuring and estimating the ZTD The refractivity in an air parcel is often written N = k1 pd e e + k2 + k3 2 , T T T (1) where the parameters k1 , k2 and k3 are empirically determined from laboratory experiments, T is the absolute temperature, e represents the water vapour pressure, and p d is the sum of the partial pressures of the “dry” gases. The refractivity is usually split into 1 a hydrostatic and a non-hydrostatic (wet) part: N = N h + Nw . (2) Accordingly, the ZTD can be presented as a sum of a zenith hydrostatic delay (ZHD) and a zenith wet delay (ZWD): ℓ = ℓ h + ℓw . (3) The ZHD is due to the refractivity associated with the gasses in hydrostatic equilibrium, obeying the hydrostatic equation for an ideal atmosphere: dp = −ρg . (4) dz Above, p is the partial pressure of the gas at an altitude z, ρ is the gas density, and g is the gravitational acceleration. The ZWD in (3) is mostly due to water vapour (e.g., Bevis et al., 1992). Davis et al. (1985) proposed an approximate formula for estimating the ZHD at a given site with latitude θ in degrees and height ho in km above the ellipsoid from pressure measurements: P , (5) ℓh = (2.2768 ± 0.0024) f (θ, ho ) where ℓh is in mm, the ground pressure P is in hPa, and f (θ, ho ) = 1 − 0.00266 cos(2θ) − 0.00028 ho . (6) The amount of water vapour expressed as the height of the column equivalent liquid water is termed integrated precipitable water vapour (IPWV) and is defined as Z 1 ∞ I= ρv dz , (7) ρ l ho where ρl and ρv are the densities of liquid water and water vapour. The IPWV is typically measured in mm. The quotient ℓw , (8) Q= I where both ℓw and I are in mm, has been described by Askne and Nordius (1987): µ ¶ k3 ′ −8 Q = 10 ρl Rv k2 + . (9) Tm Above, Rv = 461.5 J kg−1 K−1 is the specific gas constant for water vapour; we assume k2′ ≈ 22 K/hPa and k3 ≈ 3.7·105 K2 /hPa after Emardson and Derks (2000); the parameter Tm is defined as the mean temperature of the water vapour: Z ∞ ρv dz ho . (10) Tm = Z ∞ ρv dz ho T Based on radiosonde profiles acquired over the US, Bevis et al. (1992) proposed the linear regression Tm ≈ 70.2 + 0.72 To , where To is the surface temperature in Kelvin. Other relations for Europe are reported by e.g. Emardson and Derks (2000). The ZTD time series are dominated by the temporal variation of the ZWD, even though the absolute values of the ZHD are about an order of magnitude larger those of the ZWD. From the viewpoint of numerical weather modelling, the ZTD can easily be calculated using integration over the different levels of the atmospheric model. Conversely, the ZTD estimates based on the GPS measurements carry information about atmospheric pressure and humidity. As the current amount of humidity measurements assimilated into NWP models is insufficient for many applications, the GPS estimates of the ZTD are likely to become valuable for the weather forecasting community. 2 2.1 Radiosonde data Balloons carrying measurement sensors can be used to acquire vertical profiles of temperature, pressure, humidity, and wind. The regular launching of balloon-borne radiosondes (RS) is costly; this limitation has an impact on the amount of the acquired data and on their quality. Recent studies over large regions of the northern hemisphere have shown that a low RS launch frequency (e.g. only twice a day) may introduce bias of up to 3% to the estimates of average amount of precipitable water (Dai et al., 2002). The data obtained from the atmospheric soundings are often reduced to profiles containing the values at a set of standard pressure levels (surface, 1000, 850, 700, 500, 400, 300 hPa), along with the levels representing significant changes in the observed atmospheric parameters. The radiosondes used now in the Nordic countries are of the Vaisala RS-80 type. Recent investigations on the performance of the Vaisala RS-80 relative humidity (RH) sensors suggest an age-related contamination by the packaging of the sondes (Niell et al., 2001; Wang et al., 2002). As the sensors typically underestimate the actual RH values by as much as 5% in very humid conditions, correction algorithms should be applied (Lesht, 1999; Wade and Schwartz, 1993). The vertical profiles of RH, and the air temperature and pressure can be used to derive the partial pressure of water vapour and subsequently the wet refractivity. The RH expressed in percent is defined as: U% = e 100 , esw (11) where e is the partial pressure of water vapour; esw is the saturation vapour pressure over water surface, and characterizes the conditions for a given volume and temperature when water can be found in equilibrium in both its liquid and gaseous phases. Empirical formulas for calculating esw (T ) have been derived by fitting a model to tabulated measurements (e.g. Crane, 1976; List, 1958). Crane (1976) introduced the model ¶ µ ¶¸ · µ T T − 273 − 5.31 ln (12) esw (T ) = 6.105 exp 25.22 T 273 using tabulated values of esw in hPa and T in Kelvin degrees. The model has an uncertainty within ±0.4% over the temperature range from 243 K to 303 K (−30 ◦ C to +30◦ C). Note that the partial pressure of water vapour equals the saturation pressure at the dew point temperature Td , i.e. e = esw (Td ) . (13) 2.2 GPS data GPS is a navigation system using several satellite-borne radio transmitters. At present, the constellation consists nominally of 24 active satellites in 6 orbital planes, allowing for a simultaneous visibility of at least 4 satellites for any receiver on the surface of the earth (Hofmann-Wellenhoff et al., 1997). For the purposes of meteorology (water vapour in the atmosphere), high measurement accuracy is achieved by using state of the art GPS receivers capable of measuring the carrier phase of the received signal (Blewitt, 1993). Information about the moisture of the atmosphere can be extracted from the ZTD fields, using a network of GPS receivers. The data processing can be carried out with the GIPSY-OASIS software (Webb and Zumberge, 1993) using the Precise Point Positioning solution (Zumberge et al., 1997). This processing is now routinely performed at the analysis centre at the Onsala Space Observatory for a receiver network (Fig. 1) (e.g., Emardson et al., 1998). 3 KEVO KEVO KIRU KIRU SODA SODA OULU OULU ROMU ROMU SKEL SKEL VILH VILH KUUS KUUS OVER OVER ARJE ARJE UMEA UMEA OSTE OSTE JOEN JOEN VAAS VAAS SVEG SVEG SUND SUND OLKI OLKI KARL KARL VANE VANE GOTE GOTE RR O UO TTU LEKS MART LEKS MART VIRO VIRO METS METS VAST VAST LOVO LOVO NORR NORR A RA OR BBO SAA S NN OO JONK JONK VISB VISB OSKA OSKA 10˚E HASS HASS MALM MALM Fig. 1. Locations of the GPS sites in Sweden and Finland. 2.3 Microwave radiometer data A microwave radiometer is used at the Onsala Space Observatory for independent assessment of atmospheric water vapour variability in general and the ZWD in particular. In the following, it is referred to as a Water Vapour Radiometer (WVR). It operates at 21.0 and 31.4 GHz measuring the sky brightness temperatures with a long-term repeatability (order of years) of approximately 1 K. This corresponds to an absolute uncertainty in the ZWD of 7 mm; the repeatability of the measured brightness temperatures becomes better for shorter time scales. Typically, RMS differences obtained when WVR data are compared to other independent measurements, such as radiosondes, are of the same order. More details on the performance of the WVR have been documented elsewhere (Elgered and Jarlemark, 1998). 3 3.1 Extracting ZTD time series from HIRLAM The non-linear observation operator for the ZWD The non-linear observation operator H(·) is used to project the background state onto the observation space. The background IPWV estimate in the current HIRLAM version is derived from (7): N X 1 I≈ qi (Pi+1/2 − Pi−1/2 ) , (14) ρl gs (θ) i=1 where gs (θ) is the gravity acceleration as a function of latitude of the site; qi is the background specific humidity at model level i; N denotes the number of model levels. The 4 background pressure Pi±1/2 is evaluated at the model layer boundaries, and approximates the partial pressure of the dry gases in the hydrostatic equation; the other parameters are evaluated at the center of each grid cell. The ZWD part of the operator H(·) is based on (9), and is defined by the numerically calculated IPWV for the model levels: µ ¶ N 10−8 Rv X k3 ′ ℓw ≈ qi k2 + (Pi+1/2 − Pi−1/2 ) , gs (θ) Ti (15) i=1 where Ti is the background temperature at model level i. 3.2 The tangent-linear observation operator The tangent-linear form H is the derivative of H(·) with respect to the state x of the numerical weather model, evaluated at the background xb : H = ∇x H(x) |x=xb . (16) The linear form of the observation operator for the IPWV over a single site, based on (14), is: " N # N X X 1 qi (P̄i+1/2 − P̄i−1/2 ) , (17) q̄i (Pi+1/2 − Pi−1/2 ) + I¯ ≈ ρl gs (θ) i=1 i=1 where the tangent-linear pressure and specific humidity are denoted P̄i±1/2 and q̄i , respectively, and are the derivatives of the corresponding parameters at the background state. 3.3 ZTD Innovation vector ZTD estimates for one year were acquired, starting with 1 May 2000, at a sampling rate of 5 minutes, acquired using 35 receivers on the territory of Sweden and Finland. These time series were subsequently decimated at 6-hour sampling intervals in order to match the HIRLAM forecast runs. These forecast runs were performed off-line at the Swedish Meteorological and Hydrological Institute, and covered the same period. For the calculation of the high resolution short term forecasts the HIRLAM 3DVAR was used (Källén, 1996; Gustafsson et al., 2001; Lindskog et al., 2001). The differences between the short term ZTD forecast fields and the estimates of ZTD from GPS (the innovations) are shown in Fig. 2. Note that the HIRLAM–GPS innovation time series (ℓH − ℓG) appear to have higher variability in the summer. For all the GPS sites (Fig. 1) there is a persistent over estimation of ZTD by the GPS technique of about 1 cm. When a rectangular “sliding window” was applied to the differenced time series for each of the GPS sites, there was no apparent seasonal dependence in the resulting slowly varying bias between the HIRLAM and the GPS ZTDs. We do not correct the innovation data for this time-varying component, and we only remove the mean difference when the horizontal correlations of the innovations were calculated for each pair of GPS sites. 4 Spatial correlations using analytic models When modelling atmospheric parameters, it is a common practice to assume horizontally homogeneous, isotropic random field, locally time-stationary over a given period. Our approach to spatial correlation modelling is based on the statistical analysis of turbulent atmosphere (Tatarskii, 1971). We are interested in the spatial structure of the correlation 5 0.06 HIRLAM−GPS innovations in ZTD, m 0.04 0.02 0 −0.02 −0.04 −0.06 0 50 100 150 200 250 Days since May 1, 2000 300 350 400 Fig. 2. Time series of the differences between the HIRLAM and GPS zenith total delay estimates (i.e. the innovations) for the site at Onsala. The corresponding time series of the innovations were produced for all GPS sites, and then used to calculate the correlations of the GPS estimation errors. coefficient ρ which is the normalized two-dimensional auto-covariance function of a random field. We use a sampled version (both in time and in space) of this field to calculate ρ as a function of the site separation r (e.g., Stoew et al., 2001). Different authors have proposed various analytical expressions which describe the spatial correlation of a random field as a function of distance (Daley, 1991; Tatarskii, 1971). These expressions satisfy the requirement that any proper 2D correlation model should yield a non-negative power spectral density (PSD) function. We consider the Gaussian and the von Kármán analytical functions in order to model spatial correlations in polar coordinates: µ 2¶ r (18) ̺1 (r) = a · exp − 2 r1 µ ¶ν µ ¶ b r r ̺2 (r) = ν−1 Kν (19) 2 Γ(ν) r2 r2 where Kν (x) is the modified Bessel function of the second kind with parameter ν > 0; the parameters a and b are scaling factors; r1 , r2 are characteristic lengths for the respective function fitted to the data. Any linear combination of (18) and (19) forms a proper 2D correlation function. To describe 2D spatial correlations, a second-order autoregressive function is sometimes used in the literature (e.g., Daley, 1991); it is, however, a special case of (19). For a 2D isotropic field, the Fourier transform is equivalent to a zero-order Hankel transform (Erdélyi et al., 1954). Thus we can write the normalized power spectral densities 6 corresponding to (18) and (19): µ ¶ g1 (κ) κ 2 r12 = exp − g1 (0) 4 1 g2 (κ) =¡ ¢ν+1 g2 (0) 1 + κ 2 r22 (20) (21) where κ corresponds to a set of 2D spatial wave number pairs (k, l), for which κ 2 = k 2 +l2 . This approach has its limitations which stem from the ignoring of the curvature of the earth surface. Hollingsworth and Lönnberg (1986) described the isotropic component of model-prediction error covariances using finite sums of Bessel functions, taking into account the spherical harmonic expansion of the field. We refrain from using their representation of spatial correlations, as it yields large variations of the correlation model for baselines where no data are available (Stoew, 2001). 5 Spatial correlations based on NWP and GPS data It is possible to separate the contributions from two estimation error fields, assuming that the spatial correlation structure is known for one of them. We denote the ZTD true values as ℓi and ℓj for the locations i and j; the ZTD estimated using the HIRLAM and GPS are G H H G G ℓH i and ℓi , respectively. The differences ei = ℓi − ℓi and ei = ℓi − ℓi represent the time series of the estimation errors for the two methods. The covariance of the differenced time series (the innovations) carries information about the spatial structure of the e H and eG fields: i h G H G − ℓ cov ℓH − ℓ , ℓ i i j j h i G H G = cov ℓH − ℓ + ℓ − ℓ , ℓ − ℓ + ℓ − ℓ i i j j i i j j h i G H G = cov eH i − ei , ej − ej h i h i H G G = cov eH , e + cov e , e (22) i j i j , where in the last line we use the fact that there is no correlation between the HIRLAM and the GPS estimation errors in ZTD (no GPS data are assimilated into HIRLAM). H An analytical model of the term cov[eH i , ej ], derived using the NMC method (Parrish and Derber, 1992), is used to fit a curve to the remaining term in order to obtain the horizontal correlations of the GPS errors. The NMC method was applied at SMHI for three different set-ups of HIRLAM, namely for the 44-, 33-, and 22-km horizontal resolution. The resulting models of the HIRLAM prediction error spatial correlations are presented in Fig.3. The spatial correlations depend on the chosen HIRLAM resolution, and as the model resolution increases, the correlation spatial scales become smaller (i.e., the correlation curves become steeper). In the same plot we present also the spatial correlations derived using the ensemble assimilation method for a 44-km set-up of HIRLAM. The parameter values for the corresponding models, ̺H(r), are summarized in Table 1. The spatial correlations for the innovations ℓH − ℓG are shown in Fig. 4. These correlations can be expressed by the normalized pairwise covariances: h i G H G cov ℓH i − ℓi , ℓj − ℓj ̺i,j(r) = , (23) σi σj 7 Analytic function fit (NMC/44km) Analytic function fit (NMC/33km) Analytic function fit (NMC/22km) Analytic function fit (ENS/44km) 0.9 0.8 Correlation coefficients 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 200 400 600 800 1000 Baseline length, km 1200 1400 1600 1800 Fig. 3. Analytic functions modelling the horizontal correlations of forecast errors for ZTD calculated for HIRLAM. The solid line presents the model of the forecast error correlations, estimated using the NMC method for a 44-km resolution NWP model. For the 33-km and the 22-km resolution model, the NMC method reveal even smaller correlation scales (# and ▽, respectively). The results from the use of the ensemble assimilation method applied to the 44-km resolution version of HIRLAM are also shown (×). where r is the site separation; σi , σj are the standard deviations of the innovation time series at sites i and j, respectively. They are of the order of 1.1–1.3 cm, and decrease slightly with site latitude and height. Note that because spatially uncorrelated estimation/forecast errors contribute to the standard deviation of the innovations, ̺ i,j(0) ≤ 1 . A composite model, based on a linear combination of three von Kármán curves (19), Table 1. Parameters of the composite von Kármán correlation models (see text). The values of the parameters ν and r2 for the ̺H(r) component were estimated using the NMC method at a horizontal resolution of 44, 33, and 22 km. The ensemble assimilation method was used to derive the parameters in the fourth group in the table. Model ai Parameter ν r2 ̺H(r) ̺G 1 (r) ̺G 2 (r) 0.27 0.53 0.07 1.05 3.33 2.98 192.80 37.02 685.33 NMC, 44-km resolution ̺H(r) ̺G 1 (r) ̺G 2 (r) 0.37 0.40 0.14 1.12 2.65 1.01 135.76 38.01 714.69 NMC, 33-km resolution ̺H(r) ̺G 1 (r) ̺G 2 (r) 0.36 0.42 0.13 0.90 2.21 1.01 127.10 42.31 726.27 NMC, 22-km resolution ̺H(r) ̺G 1 (r) ̺G 2 (r) 0.33 0.41 0.16 1.82 2.68 0.80 84.93 36.23 717.73 Ensemble assimilation, 44-km resolution 8 was fitted to the correlation values: ̺(r) = a1 ̺G (r) + a2 ̺G (r) + a3 ̺H(r) , 1 2 (24) where a1 , a2 , a3 are scaling factors; ̺H(r) represents the signature of the HIRLAM ZTD forecast error correlations, and is plotted using the dash-dotted line in Fig. 4. The main impact of the HIRLAM forecast errors is for baselines r < 1000 km. The other two von Kármán curves, ̺G (r) and ̺G (r), model the contribution of GPS errors to the inno1 2 vation correlations (the solid line in Fig. 4). The composite model (24) is presented by the dashed line, and has 7 degrees of freedom: three for each ̺G (r) and ̺G (r), and one 1 2 H for ̺ (r). The parameter values for the model are shown in Table 1. The slowly declining spatial correlations for baselines r > 400 km in the same graph would have to be attributed to large scale spatial correlations in GPS estimation errors. This behavior is in agreement with earlier simulation results (Jarlemark et al., 2001). However, for small separations, the results depend strongly on the correctness of the model ̺H(r), and on the availability of ZTD data from nearby GPS sites. For the short baseline separations, the ZTD error correlations derived using the NMC method decline at a slower rate compared to the HIRLAM–GPS difference data. This result indicates that, in reality, the forecast error correlations may have shorter length scales than those estimated by Berre (1997). It is further supported by the notion that these correlations are calculated for the +6 hour forecasts, utilizing differences between forecasts valid at the same time, but generated at different ranges (+24 and +48 hours respectively). We speculate that the spatial scales of the forecast errors should decrease if the forecast range is decreased. These scales also depend on the chosen horizontal resolution of the NWP model, as demonstrated in Figures 4, 5, and 6. A new method, used by the European Centre for Medium-Range Weather Forecasts, based on ensemble assimilation indicates smaller length scales than those determined by the NMC technique. The newer ensemble assimilation approach uses sets of randomly perturbed observations and parameters of the numerical model to derive the error correlation structures in the forecasts (Houtekamer et al., 1996). The results from the ensemble assimilation approach are shown in Fig. 7. 9 Horizontal correlations of the HIRLAM−GPS ZTD innovations Innovation correlations Model fit (innovations) Model fit (Hirlam forecasts) Model fit (GPS estimates) 0.9 0.8 Correlation coefficients 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 200 400 600 800 1000 Site separation, km 1200 1400 1600 1800 Fig. 4. Horizontal correlations of the HIRLAM–GPS zenith total delay innovations for a one year period (∗). The dash-dotted line represents the contribution of the forecast errors correlations for ZTD, estimated using the NMC method for HIRLAM at a resolution of 44 km. The dashed line shows the best fit to the innovation data of a composite model of the correlations. The solid line represents the remaining term, corresponding to the correlations of GPS errors. Horizontal correlations of the HIRLAM−GPS ZTD innovations Innovation correlations Model fit (innovations) Model fit (GPS estimates) Model fit (Hirlam forecasts) 0.9 0.8 Correlation coefficients 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 200 400 600 800 1000 Site separation, km 1200 1400 1600 1800 Fig. 5. Horizontal correlations of the ZTD innovations for a one year period (∗). The dash-dotted line represents the contribution of the forecast errors correlations for ZTD, estimated using the NMC method for HIRLAM at a resolution of 33 km. The dashed line shows the best-fit composite model for the innovation correlations. The solid line represents the remaining term, corresponding to the GPS error correlations. 10 Horizontal correlations of the HIRLAM−GPS ZTD innovations Innovation correlations Model fit (innovations) Model fit (GPS estimates) Model fit (Hirlam forecasts) 0.9 0.8 Correlation coefficients 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 200 400 600 800 1000 Site separation, km 1200 1400 1600 1800 Fig. 6. Horizontal correlations of the ZTD innovations for a one year period (∗). The dash-dotted line represents the contribution of the forecast errors correlations for ZTD, estimated using the NMC method for HIRLAM at a resolution of 22 km. The dashed line shows the best-fit composite model for the innovation correlations. The solid line represents the remaining term, corresponding to the GPS error correlations. Horizontal correlations of the HIRLAM−GPS ZTD innovations Innovation correlations Model fit (innovations) Model fit (GPS estimates) Model fit (Hirlam forecasts) 0.9 0.8 Correlation coefficients 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 200 400 600 800 1000 Site separation, km 1200 1400 1600 1800 Fig. 7. Horizontal correlations of the ZTD innovations for a one year period (∗). The dash-dotted line represents the contribution of the forecast errors correlations for ZTD, estimated using the ensemble assimilation method. The dashed line shows the best fit to the innovation data of a composite model of the correlations. The solid line represents the remaining term, corresponding to the correlations of GPS errors. 11 6 Spatial correlations using vertical coordinate residuals It is well known that there is a strong correlation between the estimation errors in the ZTD and the vertical coordinate. In order to assess the results in the previous section it is worthwhile to study the spatial error correlations of the vertical coordinate. Vertical coordinate residuals have been calculated using GPS data from 1996 to the beginning of 2001 from daily solutions (Johansson et al., 2002; Scherneck et al., 2003). The data have been collected from the SWEPOS and FinnRef networks (see Fig. 1). We calculated the cross-correlations of the time series for each pair of sites, and the results are shown in Fig. 8. The errors in the estimated station height have common sources with the errors in the ZTD estimated using GPS. The plot in Fig. 8 supports the simulation results by Jarlemark et al. (2001), suggesting a slow decrease of the horizontal correlations of the GPS estimation errors. For the short baselines (∼200 km), the correlation values remain below 0.9. This may indicate a random, uncorrelated component in the data, which has a scaling effect on all correlation coefficient values. Although the spatial correlations of these residuals are clearly present, the scatter of the data points prevents us from providing a reliable fit of an analytical model of the type discussed in Section 4. 1 0.9 0.8 Correlation coefficient 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 200 400 600 800 1000 Site separation, km 1200 1400 1600 1800 Fig. 8. Horizontal correlations of the vertical coordinate residuals for the period 1996–2001, using data from SWEPOS and FinnRef. 12 7 7.1 Temporal correlations of slant delay residuals using ray tracing Hydrostatic and wet mapping functions The slant propagation delay is usually converted to an equivalent zenith delay using mapping functions (Davis et al., 1985). As an example, the wet-delay MF (WMF) m w (ε) relates the equivalent zenith wet delay to the slant delay ℓw (ε): ℓzw (ε) = ℓw (ε) , mw (ε) (25) where ε is the elevation angle. In GPS processing, the ZWD estimate is a weighted average of the equivalent zenith delays ℓzw (εi ), estimated for the directions towards the satellites visible at a given epoch. It is somewhat different from the actual value of the delay in the zenith direction, ℓ w (90◦ ), due to (a) inhomogeneities of the atmosphere, and (b) imperfections in the mapping function. The hydrostatic mapping functions (HMFs) and delays are related in a similar way. The ray-tracing software described below assumes the atmosphere to be horizontally isotropic above the RS launch site. Most often, the total slant delays are expressed by ℓ(ε) = ℓh mh (ε) + ℓw mw (ε) , (26) where ℓh and ℓw denote the ZHD and the ZWD, respectively. The GPS estimates utilize mapping functions which at present introduce an rms error in total delay of less than 6 mm at satellite elevation angles larger than 10◦ (Bisnath et al., 1997; Niell, 1996). The Ifadis (1986) MFs are used for ε > 2◦ , and are in the form: 1 mI(ε) = a sin ε + sin ε + , (27) b sin ε + c The Ifadis HMF/WMF utilize surface measurements of the temperature, water vapour-, and total pressure. Time series of the Ifadis HMF values for ε = 10◦ at the Landvetter site are shown in Fig. 9. The Niell (1996) HMF and WMF are based on statistics covering RS ascents from 26 locations and have the general form: a 1+ 1+ mN(ε) = sin ε + b 1+c a sin ε + , (28) b sin ε + c where the parameters a, b, and c are obtained using only information about the site location. These MFs are useful for ε > 3◦ . For the Niell WMF, these parameters depend only on site latitude, while the behavior of the HMF is related also to the day of the year and the station height (see Fig. 9). The differences between the parameterizations of the Niell and the Ifadis mapping functions are summarized in Table 2. 13 5.564 Niell hydrostatic Ifadis hydrostatic 5.562 Hydrostatic MF value at ε = 10° 5.56 5.558 5.556 5.554 5.552 5.55 5.548 5.546 0 50 100 150 200 250 Days since Jan 1, 2000 300 350 400 Fig. 9. Time series of the Niell and Ifadis hydrostatic mapping functions, calculated for ε = 10 ◦ at Landvetter. The systematic differences between these two MFs can reach 0.1% in the winter for this site. The seasonal variations in the Niell MF are only related to site latitude and the day of the year, while the Ifadis HMF utilizes surface measurements. 7.2 Ray tracing using radiosonde ascents The ray-tracing software uses vertical profiles of RH, air temperature and pressure to calculate the wet and hydrostatic propagation delays, and the geometric delays (bending), for a given elevation angle. The results from the ray tracing procedure are then used to produce the normalized wet/hydrostatic atmosphere thickness (slant delays divided by the corresponding zenith values) in order to obtain the residuals for the Niell and Ifadis mapping functions, for ε = (5◦, 10◦, 15◦ ). Ideally, the time correlations of these residuals should represent the temporal behavior of the GPS slant delay estimation errors. The results for the Niell HMF/WMF are shown in Fig. 10, spanning the year 2000. Note that the temporal correlations of the Niell WMF residuals appear uncorrelated on scales larger than one day, while the HMF residuals are correlated on scales of several days. The results for the Ifadis HMF/WMF are shown in Fig. 11. Note the slowly-varying bias in the HMF residuals (the HMF underestimates the ray-trace results in the winter); this bias causes the correlations of the Ifadis HMF residual to be rather high on monthly scales. Furthermore, there is an apparent diurnal signature in these correlations, which should be attributed to the differences in the day-time and night-time variations of the total pressure and temperature in the vertical direction. The values of the temporal correlations for these structure function residuals change insignificantly, only by about 1% for the different elevation angles. Table 2. Summary of the properties of the Niell and the Ifadis mapping functions. Hydrostatic MF Wet MF Niell Depends on height, latitude and day of the year. Depends only on station latitude. Ifadis Parameterized using surface temperature, pressure and humidity. 14 Slant delay residuals, m 0.02 Niell wet Niell hydrostatic a) 0.01 0 −0.01 −0.02 0 50 100 150 200 250 Days since Jan 1, 2000 300 350 400 1 Correlation coefficient b) 0.8 0.6 0.4 0.2 0 −0.2 0 5 10 15 Time lag τ, days 20 25 30 Fig. 10. Residuals of the Niell HMF, calculated for ε = 10◦ at Landvetter (upper graph). The rms variations of the HMF residuals are about 5 times larger than those of the Niell WMF. The temporal correlations of the Niell mapping function residuals are shown in the lower graph. The WMF residuals (#) are practically uncorrelated on scales larger than one day, while the Niell HMF residuals (×) are correlated on scales of several days. Slant delay residuals, m 0.02 Ifadis wet Ifadis hydrostatic a) 0.01 0 −0.01 −0.02 0 50 100 150 200 250 Days since Jan 1, 2000 300 350 400 1 Correlation coefficient b) 0.8 0.6 0.4 0.2 0 −0.2 0 5 10 15 Time lag τ, days 20 25 30 Fig. 11. Residuals of the Ifadis HMF, calculated for ε = 10◦ at Landvetter (upper graph). The rms variations of the HMF residuals are about 6 times larger than those of the Ifadis WMF; there is also an apparent systematic bias in the HMF for this site. The temporal correlations of the Ifadis MF residuals are shown in the lower graph. The Ifadis WMF residuals (#) seem to be uncorrelated for time scales larger than one day. The Ifadis HMF residuals (×) are correlated on scales of several days and exhibit a diurnal pattern. 15 8 Temporal correlations of ZTD from GPS position residuals We calculated an estimate of the normalized auto-covariance function for a 1-year long time series of differences between ZWDs, estimated using a WVR and GPS data at Onsala. The results are shown in Fig. 12. These correlation values decrease to e −1 at a time lag of about 1 day. The temporal structure of the WVR errors is unknown; therefore, the GPS error temporal correlations cannot be deduced from the GPS–WVR differences without reservations. Correlation coefficient 1 0.8 0.6 0.4 −1 e 0.2 0 0 5 10 15 Time lag τ, days 20 25 30 Fig. 12. The temporal correlation values for the differences between GPS and WVR data at Onsala, for the year 2001. The time correlation of these differences decreases to about 1/e at a time lag τ ≈ 1 day. Another way to analyze the temporal correlations of the errors in GPS-estimated ZWD is to study the differences between time series of integrated water vapour derived radiosonde- and GPS measurements. Although the errors in the RS data can be assumed temporally uncorrelated, the variance of the RS errors actually makes the timecorrelations of the differences appear smaller (see Fig. 13). Furthermore, there are few colocated RS and GPS sites, and using similar differenced time series from nearby sites would introduce a component due to the atmosphere itself. Corr. coefficient 1 0.8 0.6 0.4 −1 e 0.2 0 0 5 10 15 Time lag τ, days 20 25 30 Fig. 13. The temporal correlation values for the differences between GPS and RS data at Sundsvall, for the year 1997. The time correlation of these differences decreases to about 1/e at a time lag τ ≈ 1/2 day, partly due to the variance of the random measurement errors in the RS data. The temporal behavior of the errors in GPS-estimated ZWD/ZTD can be described indirectly, using the residual of the vertical coordinate of a permanent GPS site, as these two parameters have common error sources. Their relation further depends on the data processing strategy, the used zenith delay MFs, the antenna phase center modelling, and the modelling of the ocean/atmosphere loading effects. Time series of the vertical coordinate residuals, discussed earlier, are shown in Fig. 14, together with the temporal correlations for the Onsala GPS site. The decrease in the correlations of these residuals is consistently slow for all studied SWEPOS GPS sites; the decorrelation time is about 1–2 days, and is summarized in Table 3. The variances of the individual time series were found to be less than 100 mm2 . Apparently, the northern16 Vertical coordinate residuals, m 0.04 a) 0.02 0 −0.02 −0.04 −0.06 1996 1997 1998 1999 Calendar year 2000 2001 2002 Correlation coefficient 1 b) 0.8 0.6 0.4 −1 e 0.2 0 0 5 10 15 Time lag τ, days 20 25 30 Fig. 14. Time series of the vertical coordinate residuals for Onsala (upper graph), and the respective temporal autocorrelation values (lower graph). The time correlation typically decreases to 1/e at a time lag τ ≈ 2 days. most sites have slightly longer decorrelation times and slightly higher variances. This can further be attributed to snow and ice, accumulating on the antenna radomes during the winter (Jaldehag et al., 1996). Lastly, the data are generated at different rates (daily coordinate estimates vs. 5-minute intervals for the ZTDs). These considerations imply that the presented vertical residual correlations should be used as an approximation of the ZWD/ZTD temporal error correlations cautiously, as the former over-estimate the latter. Correlation coefficient 1 0.8 0.6 0.4 −1 e−1 0.2 0 0 5 10 15 Time lag τ, days 20 25 30 Fig. 15. Temporal correlation values for the vertical coordinate residuals from 2001 at Onsala (#). The time correlations of the GPS–WVR differences for 2001 are plotted for comparison (dashed line), and show consistent correlation values. Although the linear trends, annual and sub-annual oscillations (up to 3 cycles per year), and air-pressure loading have been removed from these time series (Scherneck et al., 2003), there are unmodelled components on time scales larger than one month, causing the correlation estimates to remain rather high for time lags τ > 10 days for some GPS sites, including the one at Onsala. In particular, the years 1996 and 2000 appear to have short-term trends (see Fig. 14). Figure 15 presents the temporal correlations obtained 17 Table 3. Decorrelation of the residuals of the vertical coordinates for the SWEPOS sites. The data sets span the years 1996–2001. The temporal correlation function is shown for time lag τ = 2 days. The GPS sites are sorted by increasing latitude values. Site HASS OSKA ONSA VISB BORA JONK NORR VANE LOVO KARL MART LEKS SVEG SUND OSTE UMEA VILH SKEL OVER ARJE KIRU Latitude [◦ N] Correlation ̺(τ ) at τ = 2 days Variance [mm2 ] 56.09 57.06 57.40 57.65 57.72 57.74 58.59 58.69 59.34 59.44 60.59 60.72 62.02 62.23 63.44 63.58 64.70 64.88 66.31 66.32 67.88 0.28 0.30 0.33 0.34 0.32 0.28 0.30 0.30 0.33 0.32 0.36 0.35 0.35 0.32 0.39 0.34 0.38 0.39 0.33 0.36 0.40 71 72 70 72 60 63 64 62 67 64 69 76 68 69 80 69 75 80 72 85 100 from a subset without an apparent trend in the vertical residual data for Onsala (the year 2001), together with a plot of the correlations of GPS–WVR differences from Fig. 12. We also performed a study of smaller, seasonal subsets of these data and found no seasonal dependence of the time correlation coefficients. 9 Discussion The estimation of the atmospheric delay using GPS is associated with errors whose behavior is assumed to be random. The spatial structure of the error correlations in the ZTD is derived most reliably using innovation time series. The spatial scales of the horizontal correlations are of the order of several hundred kilometres, and the proposed analytical models are described in Table 1 of Section 5. An implementation of the von Kármán function for the purpose of a 3DVAR data assimilation should be simple. Some contribution to the ZTD estimation errors can be systematic, due to modelling (imperfect mapping functions). The temporal correlations of the corresponding errors have time-scales of days, and this is comparable to the time-scales of the atmospheric processes. Furthermore, temporally correlated errors can be introduced into the ZTD estimates by the GPS data processing strategy itself. In particular, a Kalman filter based data processing software assumes the ZTD to be a random walk process with chosen parameters, which reflect the average behavior of the atmosphere. However, a rapid change in the atmospheric state cannot be modelled well enough, resulting in short-term systematic differences. Such systematics have time scales of less than an hour, and are shorter than the typical range of the numerical forecasts of today. These systematics will become important in the future 4D data assimilation, when GPS data will be used at the precise acquisition epochs. A typical solution to the problem of time-correlated estimation errors is “thinning” of the Kalman filter estimates, i.e. using results that are sufficiently separated in time (Brown and Hwang, 1997). The comparisons to WVR data in Section 8 suggest that the contribution of the atmosphere related systematic errors is small. The assumption that the daily vertical coordi18 nate residuals and the ZTD estimation errors have a similar time-correlation structure is supported by the WVR-GPS comparison for the site at Onsala. This indicates that longterm, unmodelled trends in the satellite orbits and tidal motions of the sites propagate into the ZTD estimation errors, with a variance of less than 100 mm2 . Acknowledgements. The work on the spatial correlations of the zenith delay estimation errors was carried out in close collaboration with Nils Gustafsson. We are grateful to Alan Rogers, Jim Davis, Tom Herring, and Arthur Niell for providing us with the ray-tracing software. References Askne, J. and H. Nordius (1987). Estimation of tropospheric delay for microwaves from surface weather data. Radio Sci., 22:379–386. Berre, L. (1997). Non-separable structure functions for the HIRLAM 3DVAR. Tech. Rep. No. 30, Swed. Meteorol. and Hydrol. Inst. Master Thesis. Bevis, M., S. Businger, T.A. Herring, C. Rocken, R.A. Anthes and R.H. Ware (1992). GPS meteorology: Remote sensing of atmospheric water vapor using the Global Positioning System. J. Geophys. Res., 97(D14):15,787– 15,801. Bisnath, S.B., V.B. Mendes and R.B. Langley (1997). Effects of tropospheric mapping functions on space geodetic data. In Proceedings of the IGS Analysis Center Workshop, Jet Propul. Lab., Pasadena, Calif., 12–14 March. Blewitt, G. (1993). Advances in global positioning system technology for geodynamics investigations: 1978–1992. In Smith, D.E. and D.L. Turcotte, eds., Contributions of Space Geodesy to Geodynamics: Technology, vol. 25, pp. 195–213. AGU Geodynamics Series, Wiley & Sons, Inc. Brown, R.G. and P. Hwang (1997). Introduction to Random Signals and Applied Kalman Filtering. John Wiley & Sons. Crane, R.K. (1976). Refraction effects in the neutral atmosphere. In Meeks, M.L., ed., Methods of Experimental Physics, vol. 12B, pp. 186–200. Academic Press, N.Y. Dai, A., J. Wang, R.H. Ware and T. van Hove (2002). Diurnal variation in water vapor over North America and its implications for sampling errors in radiosonde humidity. J. Geophys. Res., 107(D10), 4090, 10.1029/ 2001JD000642. Daley, R. (1991). Atmospheric Data Analysis. Cambridge atmospheric and space science. Cambridge University Press, Cambridge, U.K. Daley, R. (1992). The effect of serially correlated observation and model error on atmospheric data assimilation. Mon. Wea. Rev., 120(1):164–177. Davis, J.L., T.A. Herring, I.I. Shapiro, A.E.E. Rogers and G. Elgered (1985). Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. Radio Sci., 20(6):1593–1607. Elgered, G. and P.O.J. Jarlemark (1998). Ground-based microwave radiometry and long term observations of atmospheric water vapor. Radio Sci., 33(3):707–717, 98RS00488. Emardson, T.R. and H.J.P. Derks (2000). On the relation between the wet delay and the integrated precipitable water vapour in the European atmosphere. Meteorol. Appl., 7(1):61–68. Emardson, T.R., G. Elgered and J.M. Johansson (1998). Three months of continuous monitoring of atmospheric water vapor with a network of Global Positioning System receivers. J. Geophys. Res., 103(D2):1807–1820, 10.1029/97JD03015. Erdélyi, A., W. Magnus, F. Oberhettinger and F.G. Tricomi (1954). Tables of Integral Transforms, vol. 2. McGrawHill Book Co. Gustafsson, N., L. Berre, S. Hörnquist, X.Y. Huang, M. Lindskog, B. Navascués, K.S. Mogensen and S. Thorsteinsson (2001). Three-dimensional variational data assimilation for a limited area model. Part I: General formulation and the background error constraint. Tellus, 53A(4):425–446. Gutman, S.I., S.R. Sahm, S.G. Benjamin, B.E. Schwartz, K.L. Holub, J.Q. Stewart and T.L. Smith (2004). Rapid retrieval and assimilation of ground based GPS precipitable water observations at the NOAA Forecast Systems Laboratory: impact on weather forecasts. J. Met. Soc. Japan, 82(1B):351–360. Hofmann-Wellenhoff, B., H. Lichtenegger and J. Collins (1997). Global Positioning System, Theory and Practice. Springer, N.Y., 4th revised ed. Hollingsworth, A. and P. Lönnberg (1986). The statistical structure of short-range forecast errors as determined from radiosonde data. Part I: The wind field. Tellus, 38A(2):111–136. Houtekamer, P.L., L. Lefaivre, J. Derome, H. Ritchie and H.L. Mitchell (1996). A system simulation approach to ensemble prediction. Mon. Wea. Rev., 124(6):1225–1242. Ifadis, I. (1986). The atmospheric delay of radio waves: Modeling the elevation dependence on a global scale. Tech. Rep. No. 38L, Chalmers Univ. Tech., Göteborg, Sweden. ISBN: 99-0605353-4. Jaldehag, R.T.K., J.M. Johansson, J.L. Davis and P. Elósegui (1996). Geodesy using the Swedish permanent GPS network: Effects of snow accumulation on estimates of site positions. Geophys. Res. Lett., 23(13):1601–1604, 96GL00970. Jarlemark, P., J. Johansson, B. Stoew, L. Gradinarsky and G. Elgered (2001). Spatial error correlation of GPS atmospheres as determined from simulations. Phys. Chem. Earth, 26(6–8):451–456. Johansson, J.M., J.L. Davis, H.G. Scherneck, G.A. Milne, M. Vermeer, J.X. Mitrovica, R.A. Bennett, B. Jonsson, G. Elgered, P. Elósegui, H. Koivula, M. Poutanen, B.O. Rönnäng and I.I. Shapiro (2002). Continuous GPS measurements of postglacial adjustment in Fennoscandia. 1. Geodetic results. J. Geophys. Res., 107(B8), 2157, doi:10.1029/2001JB000400. Källén, E. (1996). HIRLAM Documentation Manual. System 2.5. Tech. rep., SMHI, Norrköping, Sweden. Lesht, B.M. (1999). Reanalysis of radiosonde data from the 1996 and 1997 water vapor intensive observation periods: Application of the Vaisala RS-80H contamination correction algorithm to dual-sonde soundings. In 19 Proc. of the 9th Annual Atmospheric Radiation Measurements (ARM) Science Team Meeting, San Antonio, Texas. [ Available online http://www.arm.gov/pub/docs/documents/technical/confp9903/lesht-99.pdf ]. Lindskog, M. (2001). Impact of observation processing on numerical weather prediction. Phil. Lic. Thesis, Dept. Meteorology, Stockholm Univ. Lindskog, M., N. Gustafsson, B. Navascués, K.S. Mogensen, X.Y. Huang, X. Yang, U. Andræ, L. Berre, S. Thorsteinsson and J. Rantakokko (2001). Three-dimensional variational data assimilation for a limited area model. Part II: Observation handling and assimilation experiments. Tellus, 53A:447–468. List, R.J. (1958). Smithsonian Meteorological Tables. Smithsonian Inst., Washington, D.C. Niell, A.E. (1996). Global mapping functions for the atmosphere delay at radio wavelengths. J. Geophys. Res., 101(B2):3227–3246, 95JB03048. Niell, A.E., A.J. Coster, F.S. Solheim, V.B. Mendes, P.C. Toor, R.B. Langley and C.A. Upham (2001). Comparison of measurements of atmospheric wet delay by Radiosonde, Water Vapor Radiometer, GPS, and VLBI. J. Atmos. Oceanic Technol., 18(6):830–850. Parrish, D.F. and J.C. Derber (1992). The National Meteorological Centre’s spectral statistical interpolation analysis system. Mon. Wea. Rev., 120(8):1747–1763. Scherneck, H.G., J.M. Johansson, H. Koivula, T. van Dam and J.L. Davis (2003). Vertical crustal motion observed in the BIFROST project. J. Geodynamics, 35:425–441. Stoew, B. (2001). On the Temporal and Spatial Structure of the Estimation Errors in GPS Meteorology. Tech. rep. No. 411L, Dept. of Radio and Space Science, Centre for Astrophysics and Space Science, Chalmers Univ. of Technol., Göteborg, Sweden. Stoew, B., J.M. Johansson and G. Elgered (2001). An assessment of estimates of IWV from ground-based GPS data. Meteorol. Atmos. Phys., 77(1–4):99–107. Tatarskii, V.I. (1971). The Effects of the Turbulent Atmosphere on Wave Propagation. Israel Program for Scientific Translations, Jerusalem. Wade, C.G. and B. Schwartz (1993). Radiosonde humidity observations near saturation. 8th Symposium on meteorological observations and instrumentation, Anaheim, California, Jan 17–22. Wang, J., H.L. Cole, D.J. Carlson, E.R. Miller, K. Beierle, A. Paukkunen and T.K. Laine (2002). Corrections of humidity measurement errors from the Vaisala RS-80 radiosonde – Application to TOGA-COARE data. J. Atmos. Oceanic Technol., 19:981–1002. Webb, F.H. and J.F. Zumberge (1993). An introduction to the GIPSY/OASIS-II. JPL Publ. D-11088, Jet Propul. Lab., Pasadena, California. Zumberge, J.F., M.B. Heflin, D.C. Jefferson, M.M. Watkins and F.H. Webb (1997). Precise point positioning for the efficient and robust analysis of GPS data from large networks. J. Geophys. Res., 102(B3):5005–5017. 20