An observation operator for ground-based GPS slant delays

An observation operator for the Global Positioning System (GPS) Slant Total Delay (STD) observations is developed. The operator is implemented in the High Resolution Limited Area Model (HIRLAM) for data assimilation of tropospheric humidity information contained in the geodetic estimates of microwave signal phase delays. The observation operator integrates refractivity along the signal path between the satellite and the receiver to obtain STD. The signal path is solved using the geometrical path, as given by the geometrical azimuth and zenith angles of the satellite as viewed at the receiver, as the starting point. The tropospheric effect on the signal path is accounted for. A preliminary validation of the operator has been performed against observational data. Model counterparts are compared with the geodetic estimates of Zenith Total Delay (ZTD) and STD. The bias and standard deviation between modelled and observed STD are −1.0 and 1.7 cm, respectively, at zenith, and −7.2 and 8.6 cm, respectively, at zenith angle 80◦.


Introduction
Meteorological data assimilation combines heterogeneous observations with a priori knowledge on the atmospheric state and produces an a posteriori estimate of the atmospheric state (e.g. Lorenc, 1986). The state vector can be used as an initial condition (analysis) of numerical weather prediction (NWP) models. Due to the deterministic nature of the NWP model forecast and to the chaotic nature of atmospheric flow, high quality analysis is essential for a useful model prediction.
Observation modelling is a necessary step for utilizing observations in meteorological data assimilation. In essence, it consists of expressing the observed quantity in terms of model variables, stored in the model state vector. The observation model is usually called observation operator and it produces the model counterpart of the observed quantity.
Observation modelling of in situ measurements is different from observation modelling of remote sensing measurements. In situ measurements usually represent a physical quantity which is also represented by the model variables. In contrast, remote sensing data and model data usually represent different physical quantities. Therefore, the simplest observation operator for in situ measurements consists of interpolation of the model state to the observation location. An observation operator for remote sensing data requires, in addition to interpolation, also observation modelling (e.g. Andersson et al., 1994).
An observation operator for ground based signal delays of Global Positioning System (GPS) is presented and evaluated in this article. The observing system is presented first.
The nominal constellation of GPS consists of 24 navigation satellites continuously broadcasting microwave signals to ground based receivers. The basic principle of high accuracy geodetic positioning is to measure the phase shift of the carrier wave as the signal propagates from the satellite to the receiver. The atmosphere affects the measurements by two mechanisms. First, pressure and humidity influence the navigation signals by decreasing their speed of propagation. Consequently, the apparent distance between the satellite and the receiver is longer than the true geometric distance: in GPS literature, this phenomenon is called as tropospheric refraction or tropospheric delay (e.g. Hofmann-Wellenhof et al., 2001). Second, free electrons in the higher levels of atmosphere cause the ionospheric refraction, which is not discussed here. For a signal coming from the zenith direction, the tropospheric delay corresponds to an excess path length of about 2.5 m at the mean sea level. The contribution of the humidity is usually less than 10% of the total delay. Several advantages can be attributed to ground-based GPS as a meteorological observing system. Networks of receiver stations can be built with relatively low cost to produce humidity observations with high spatial and temporal resolution in all weather conditions. The major disadvantage is the heavy preprocessing needed for production of meteorological observations. The preprocessing leads to complicated observation errors including station-dependent biases and spatial and temporal error correlations, which are difficult to estimate and properly account for in the data assimilation process (Eresmaa and Järvinen, 2005).
The tropospheric delays can be estimated either for signals coming from the zenith direction (Zenith Total Delay, ZTD; e.g. Stoew et al., 2001) or for the slanted paths between each satellitereceiver pair (Slant Total Delay, STD; Ware et al., 1997;Alber et al., 2000;de Haan et al., 2002). From the NWP point of view, ZTD is a desirable observation since it is a linear function of vertically Integrated Water Vapour (IWV) above the GPS receiver (Bevis et al., 1992). IWV is directly related with the model humidity variable and it is ideally suited to the model geometry.
In a fine-scale NWP model, assimilation of STD instead of ZTD is supported by a number of aspects. First, since there are several satellites visible from each receiver, the amount of available STD observations is roughly 10-fold the amount of available ZTD observations. Second, STD observations contain potentially information on the atmospheric anisotropy of humidity in the vicinity of receiver, which information has been averaged out from ZTD observations in the preprocessing. Third, ZTD can be assimilated as a special case of STD; in fact, ZTD is derived from raw measurements along the slanted signal paths. Finally, STD observations provide potential for water vapour tomography in kilometric scale modelling. On the other hand, the procedure for production of STD observations is still in a stage of evolution (Elosegui and Davis, 2004). It is not clear if the accuracy of near real time STD observations is adequate for the demands set by the operational NWP systems. Also, due to the geometry of the slanted signal paths, assimilation of STD is not as straightforward as the assimilation of ZTD. Moreover, it is currently not known how to properly account for the complicated observation error correlations of STD.
This article focuses on the description of an observation operator for STD developed at the Finnish Meteorological Institute (FMI). The observation operator is implemented in the High Resolution Limited Area Model (HIRLAM; Undén et al., 2002) 3D-Var assimilation system Lindskog et al., 2001). The characteristics of the HIRLAM model are discussed in Section 2. Section 3 deals with the alternative approaches for designing the observation operator. The algorithm is discussed in detail in Section 4, and the validation aspects are covered in Section 5. Conclusions are summarized and discussion is given in Section 6.

HIRLAM model
HIRLAM is a hydrostatic limited area NWP system developed as a Nordic co-operation together with the meteorological institutes of the Netherlands, Ireland and Spain. Operational HIRLAM model configurations typically cover Europe, the Northern Atlantic and parts of the North America and Arctic Ocean in a transformed latitude-longitude grid at 31 or 40 vertical levels. The HIRLAM model variables are temperature (T), specific humidity (q), horizontal wind components (u, v) and logarithm of surface pressure (ln p s ). Lateral boundary fields are provided by the European Centre for Medium-Range Weather Forecasts (ECMWF).
The HIRLAM model makes use of hybrid vertical coordinate which is a mixture of geometric height and pressure (Simmons and Burridge, 1981

General approach for the signal path determination
The STD observation operator calculates the tropospheric delay corresponding to a given model state vector, a given GPS satellite and a given GPS receiver. A geodetic software is used for estimation of positions of the satellite and the receiver ( s and r ; notation follows Hofmann-Wellenhof et al., 2001) simultaneously with the tropospheric and ionospheric delays. The signal path between the satellite and the receiver depends on the atmospheric state and cannot be accurately deduced from the raw GPS measurements alone. However, accurate signal path is essential for modelling STD using the model state variables. A necessary task of the observation operator is therefore to determine the signal path across the model grid with sufficient accuracy. In terms of signal path determination, the modelling of STD resembles modelling of radio occultation measurements, which also provide meteorological information and are widely studied (e.g. Zou et al., 1999;Healy, 2001;Poli and Joiner, 2004).
There are at least two alternative approaches for determination of the signal path. One option is to utilize the Fermat's principle of wave propagation. This leads to an observation operator based on the Least Travel Time (LTT) principle, including a ray-tracing algorithm. By an iterative approach using NWP model fields of pressure p, T and q, such an observation operator implicitly accounts for the tropospheric effect on the signal path. Another option uses the Geometrical Path (GP) approximation as the starting point: the assumption is that the tropospheric effect on the signal path can be accounted for by explicitly modifying the Tellus 58A (2006), 1 straight GP connecting the satellite and the receiver. For given s and r , GP follows from geometry, and it is not dependent on the atmospheric state. The geometrical azimuth (α g,r ) and zenith angles (ζ g,r ) at the receiver location are used for defining the GP.
As far as the accuracy of the signal path determination is concerned, LTT principle is expected to be better than the GP approximation. On the other hand, calculation of the signal path is fast with an observation operator based on GP approximation, such as the observation operator described in this article.
There are three factors that deserve special attention in the signal path determination methodology. These factors are explained in the following.
(1) Due to the curvature of the Earth, α g,r and ζ g,r , as calculated for the receiver coordinates, are not the same as the geometrical azimuth and zenith angles as viewed in the upper parts of the signal path. More precisely, as the viewing point ν is moved along the signal path towards the broadcasting satellite, the geometrical zenith angle ζ g,ν decreases.
(2) The geometrical zenith angle ζ g,ν differs from the apparent zenith angle ζ a,ν . This is because refractivity increases from near-zero values outside the atmosphere to gradually larger values deeper in the atmosphere, as the GPS signal enters the atmosphere.
(3) Signal path determination is further complicated by the horizontal model level height gradients. Due to these gradients, the model level heights at the receiver location are not sufficient for the accurate determination of the signal path.
At the cost of modelling accuracy, it is possible to gain some computational performance by neglecting these three factors thus simplifying the signal path determination algorithm. The effect of the simplifications on observation modelling will be studied later in this article. It is expected that the modelling errors caused by applying these simplifications generally increase with increasing ζ g,r .

Algorithm description
The STD observation operator consists of separate algorithms for determination of the signal path across the model grid, pro-

Determination of the signal path across the model grid
The signal path is uniquely defined by a set of three-dimensional coordinates indicating the intersections of the signal path with the NWP model levels. Let κ k denote the intersection with the model level k. Intersection κ k is defined by latitude ϕ k , longitude λ k and height above the mean sea level z k .
The number of NWP model levels is denoted by K. For computational reasons, also model orography (K + 1) and receiver height (K + 2) are counted here as model levels. Therefore, K + 2 intersections κ k are needed for defining a signal path, i.e. k ∈ {1, 2, . . . , K + 2}. The receiver is at κ K +2 , which may in some cases lay between the NWP model levels. Similarly, κ K +1 is defined as the intersection of the signal path with the model orography.
Starting from the receiver at κ K +2 and proceeding upwards, the intersections are calculated applying GP by repeating the following algorithm kernel.
(2) Calculate the height difference z k between height of the model level k at (ϕ k , λ k ) and height of κ k+1 . The height of the model level k is calculated with the hydrostatic equation.
(3) Using z k , α g,r and ζ g,r , calculate the increments of latitude ϕ k and longitude λ k between κ k and κ k+1 .
(4) Add ϕ k , λ k and z k to the coordinates of κ k+1 to obtain κ k .
Applying the geometry as illustrated in Fig. 1, meridional and zonal displacements ( n and e) are obtained by s ≈ z k tan ζ g,r , n = s cos α g,r ≈ z k tan ζ g,r cos α g,r , e = s sin α g,r ≈ z k tan ζ g,r sin α g,r , where s is the horizontal displacement. Note that these equations are not valid in the vicinity of the poles. The meridional and zonal displacements are further transformed into latitude and longitude increments ( ϕ k and λ k ) by where a k+ 1 2 is the local radius of the curvature of the Earth, and subscript k + 1 2 refers to mean of the values at k and k + 1. Note that a k+ 1 2 varies with ϕ because the Earth is not exactly spherical.
4.1.2. Application of the algorithm kernel. Due to the three specific factors described in Section 3, blind application of the algorithm kernel is not sufficient for accurate signal path determination. The algorithm kernel corresponds to a GP observation operator applying all three simplifications. The deviations of the signal path determination from the algorithm kernel are discussed in the following.
To account for the effect of the curvature of the Earth, law of sines is used. In (1), z ν is the local geometrical height of ν above the mean sea level. It is assumed that the geometrical azimuth angle is the same as viewed either from ν or from the receiver, i.e. α g,ν = α g,r .
The refractive bending is accounted for by using apparent zenith angle ζ a,ν instead of ζ g,ν . The apparent zenith angle is calculated by the Snell's law sin ζ a,ν = sin ζ g,ν 10 −6 N + 1 , where refractivity N is a function of p, T and q. As will be seen later, these model variables are interpolated to the intersections for the calculation and integration of refractivity, and therefore accounting for the bending effect is computationally easy. It is assumed that the refractive bending has no effect on the azimuth angle, i.e. the apparent azimuth angle α a,ν = α g,ν = α g,r .
The model level height gradients are taken into account by iterating the steps 2-4 of the algorithm kernel. At each iteration step, the height difference z k is recalculated using height of the model level k at κ k obtained on the previous iteration step. Iteration convergence limit of one metre of z k is reached after 2, occasionally after 3, iteration steps.

Projection of the model variables on the signal path
Calculation of tropospheric delay requires values of p, T and q along the signal path at the intersections. These values are obtained through horizontal bilinear interpolation from the four model grid points surrounding each intersection. For the orography (κ K +1 ), vertical extrapolation is performed using constant lapse rate for T (6.5 • C km −1 ) and constant q (lowest model level value). In case the receiver is between two NWP model levels, the values at the receiver (κ K +2 ) are interpolated linearly with respect to height; otherwise vertical extrapolation with constant lapse rates is used. Receiver level pressure is calculated using the hydrostatic equation.
In case the receiver station height deviates significantly from the model orography height, the observation operator is unaware of the actual balanced structure of the planetary boundary layer. One could think of shifting the boundary layer structure to the observation height in a similar manner as the vertical profiles of the lateral boundary data are shifted for a limited area model (Majewski, 1985). However, this is not included in the observation operator described here.

Integration of refractivity along the signal path
where n = 10 −6 N + 1. In the microwave frequency domain, N is obtained by where p d and e are partial pressures of dry air and water vapour, and Z d and Z w are the compressibility factors of dry air and water vapour, respectively (Bevis et al., 1992). The empirical coefficients are k 1 = 77.60 K hPa −1 , k 2 = 70.4 K hPa −1 and k 3 = 3.739 × 10 −5 K −2 hPa −1 (Bevis et al., 1994). Deviations of Z d and Z w from unity are a few parts per thousand and will be neglected here. By substituting p d = p − e and e = qp 0.622+0.378q and arranging the terms, (4) can be rewritten as Refractivities N k at the intersections κ k are calculated using (5).
The observation operator calculates STD as the sum where STD mod is the contribution from the signal path between the receiver at k = K + 2 and the model top at k = 1 and STD top is the contribution from the signal path above the model top.
Integration of STD mod is done by where STD k represents the contribution of STD mod from the layer between the model levels k and k + 1. The numerical integration scheme adopted here is as follows. The behaviour of N with respect to height is quasi-exponential and is therefore approximated by N = exp(a + bz).
Tellus 58A (2006), 1 Since N and z at the intersections are known, the parameters a and b can be solved separately for each layer by using N k , N k+1 , z k and z k+1 . The contribution STD k can then be integrated analytically by where notation k + 1 2 refers to mean of the values at k and k + 1, and ζ a,k+ 1 2 is the apparent zenith angle. The Saastamoinen model (Saastamoinen, 1972;Bevis et al., 1992) is applied for calculation of STD top by where p 1 , ζ a,1 , ϕ 1 and z 1 are pressure (hPa), zenith angle, latitude and height (m) at κ 1 i.e. at the model top, respectively, and the denominator accounts for variations of gravity with ϕ 1 and z 1 . Since the top of the HIRLAM model grid is at 10 hPa, the magnitude of STD top is about 1% of STD.

Validation of the STD observation operator
The STD observation operator is validated by two methods. First, an experiment with standard atmospheric data is performed to ensure that modelled ZTD and STD are not sensitive to the NWP model grid configuration. Second, quantitative estimates of the accuracy of modelled ZTD and STD are made by comparison with observational ZTD and STD data.

The experiment with standard atmospheric data
Tropospheric delays modelled by the observation operator should not depend on the discretization details, i.e. the NWP model configuration, of the continuously stratified atmosphere. An experiment using standard atmospheric data is performed in order to check that modelled ZTD and STD are sufficiently independent on the number of the model levels K. Subarctic summer standard vertical profiles (McClatchey et al., 1971) of p, T and q are used as the reference atmosphere. The experiment is started with K = 31 using vertical hybrid coordinates of the operational HIRLAM at FMI. Model level pressures are calculated, and the standard profiles of T and q are interpolated to the model levels linearly in logarithm of p. ZTD and STD corresponding to the standard atmosphere are calculated by the observation operator. For STD, geometrical zenith angle 80 • is used. K is then increased to 2K − 1 by introducing new model levels between the original model levels. Standard profiles are interpolated to the new levels and ZTD and STD are recalculated. This procedure is repeated until K = 1921.  Table 1 shows the resulting behaviour of ZTD and STD as a function of K. As K is increased from 31 to 61, ZTD decreases by 0.03 cm, and further increase of K does not affect on ZTD. Behaviour of STD is qualitatively very similar. Assuming the delays modelled with K = 1921 to be the closest to the truth, a positive bias of 0.03 cm at zenith and 0.09 cm at 80 • zenith angle will be attributed to inaccuracies of the numerical integration. These biases are very small compared with the expected observation and modelling errors. It is concluded that the method described in Section 4.3 is numerically adequate for refractivity integration.

Comparison with observational data
Modelled STD is compared with post-processed ZTD and STD observations derived from the geodetic measurements. The operational HIRLAM NWP model, with 31 levels and horizontal resolution of 0.2 • , is used for the comparison. Model counterparts are calculated using 6-hr forecasts as NWP model fields. No attempt is made to obtain the model background at appropriate observation time.
5.2.1. Comparison with observed ZTD. The observational ZTD data set consists of ZTD estimates for April 2001 from 35 permanent GPS receiver stations in Finland and Sweden as postprocessed by the Chalmers University of Technology. Nearly instantaneous (average over 5-min interval) values closest to the valid time of the 6-hr forecast are used. The geodetic processing strategy used for ZTD production has been described by Stoew and Elgered (2004). The receiver stations are denoted by open squares in Fig. 2. Figure 3 shows a scatterplot consisting of 3662 ZTD observations and their model counterparts. Modelled ZTD is close to observations, but there is a negative bias indicating that the modelled ZTD is systematically lower than the observed ZTD. Systematic difference (bias), standard deviation of the difference (sdev) and root-mean-squared difference (rms) between modelled and observed ZTD, as calculated over the whole data set, are −0.84, 1.04 and 1.34 cm, respectively. For individual receiver stations, bias varies between −1.6 and −0.28 cm, sdev varies between 0.72 and 1.5 cm, and rms varies between 0.93 and 2.0 cm. These results are consistent with the comparisons reported by a number of authors (e.g. Gustafsson, 2002;Haase et al., 2003;Pacione and Vespe, 2003). The NWP forecast field and the observation operator as the origin of the bias cannot be excluded. It is, however, likely that a part of the bias is related to the systematic observing system errors, such as receiver multipath, antenna phase center variations, mapping function errors, unmodelled vertical movements of the Earth's crust and systematic assumptions of the geodetic network solution.  These receiver stations are circled in Fig. 2. Modelled ZTD follows closely observed ZTD variations at both receiver stations, although the time series for Virolahti is notably biased. Fig. 2) in the Western Europe for the period of 1-24 May 2003. STD is estimated for all satellite-receiver pairs with 30-s time interval using 80 • zenith angle cut-off. The algorithm for STD production is as follows. First, Zenith Hydrostatic and Wet Delays (ZHD and ZWD) are estimated for time slots of 20 min. ZHD and ZWD are then projected to geometrical zenith angles by Niell hydrostatic and wet mapping functions (Niell, 1996) to yield Slant Hydrostatic and Wet Delays (SHD and SWD). The isotropic part of STD is obtained as the sum of SHD and SWD. The least squares residuals of the geodetic network solution, corresponding to each satellite-receiver pair at each observation epoch, are corrected for multipath and used for inclusion of the anisotropic component to STD observations. The processing method is discussed in more detail by de Haan et al. (2002).

Comparison with observed STD. The observational STD data set is provided by the Technical University of Delft and it contains observations from 16 receiver stations (black dots in
The methodologies for production of meteorologically useful STD observations are evolving. The STD observations produced by the standard method (Alber et al., 2000) may lack part of the anisotropic information due to the least squares estimation procedure (Elosegui and Davis, 2004). Also, the observations may include station-dependent biases varying with zenith angle. These errors are not estimated for this study. Therefore the validation of the observation operator using geodetic STD estimates can only be regarded as a preliminary validation step. A more exhaustive validation could be based on e.g. comparison against ground based water vapour radiometer measurements (Braun et al., 2001).
In order to assess the implications of each of the three simplifications described in Section 3, STD is modelled using five Tellus 58A (2006), 1 Included Included different versions of the observation operator ( Table 2). The versions GP and ITER use only the geometrical zenith angle at the receiver (ζ g,r ). However, the effect of the model level height gradients is taken into account in ITER, but neglected in GP. CURV uses (1) to account for the effect of the curvature of the Earth. BEND uses (2) to account for the effect of the refractive bending of the signal path. ADVV makes no simplifications and is thus the most advanced observation operator version presented in this article. Figure 5 shows statistics of the difference between modelled and observed STD as a function of ζ g,r , calculated for all receiver stations in a subset of data consisting of 360 857 observations at 10-min time interval. Each curve corresponds to one observation operator version, as summarized in Table 2. The amount of observations increases uniformly with ζ g,r . The observation binning interval for statistics calculation is 1 • of ζ g,r . Figure 5a shows the systematic difference (bias) between modelled and observed STD. At small zenith angles (up to 40-50 • ) the bias is about −1 cm and it is of similar magnitude as is the bias between modelled and observed ZTD. This bias is considered in this article as "the near zenith bias". The two overlaying curves (GP, ITER; thin solid, dotted) show large positive bias values at large zenith angles. This bias is considered in this article as "the modelling bias". The dashed curve (BEND) also reaches large positive bias values but, very importantly, the larger is the zenith angle, the more BEND lags behind GP and ITER; that is, the modelling bias due to neglecting bending effect increases with zenith angle. These three curves are produced with observation operator versions neglecting the effect of the curvature of the Earth. Note that at about the zenith angle 62 • , the modelling bias of GP, ITER and BEND exactly counterbalances the near zenith bias (point of cross-over with zero bias line). The two other curves (CURV, ADVV; heavy solid, dashdotted) remain closer to the near zenith bias. The modelling bias of CURV starts to increase at zenith angles larger than about 75 • . Note that CURV neglects modelling of the effect of the refractive bending, which is exactly the difference between the curves for CURV and ADVV (and between GP or ITER and BEND). The remaining modelling deficiencies increase the negative bias of ADVV below the value of the near zenith bias. (2) the larger is the zenith angle, the more important it is to model the effect of the refractive bending on the signal path, and (3) it is not necessary to model the effect of horizontal model level height gradients, at least as long as the NWP model resolution is a few tens of kilometers. Figure 5b shows the random difference (standard deviation) between modelled and observed STD. Standard deviation increases with zenith angle and there is practically no difference between different observation operator versions. For the observation operator versions accounting for the curvature of the Earth (CURV, ADVV), the random difference between modelled and observed STD is comparable in magnitude with the systematic difference. It is concluded that the observation modelling alternatives presented in this article contribute equally to the random difference between modelled and observed STD, i.e. there is no particular skill in one observation operator version over another version. An explanation to this counter-intuitive result is as follows. All observation operator versions presented in this article are able to determine the signal path with nearly equal accuracy and thus there are only insignificant NWP model state differences in the calculation of refractivity. Once the signal path is determined, there are however systematic differences in different observation operator versions (i.e. systematically different ζ a,k+ 1 2 in (8)) on how to calculate STD. Therefore, the bias of different observation operator versions varies significantly (Fig. 5a), whereas the random difference does not (Fig. 5b). A corollary is that further improvements in the accuracy of the signal path determination are expected to only modestly decrease the random difference between modelled and observed STD. Figure 5c shows the total perceived difference (rms) between modelled and observed STD. Since the random differences are the same for all observation operator versions, differences between the rms curves in Fig. 5c are solely due to the systematic differences between different observation operator versions. Rms increases with zenith angle. At the zenith angle 62 • , the near zenith bias and the modelling bias of GP, ITER and BEND (all missing the effect of the Earth's curvature) are counterbalanced. This results in rms of GP, ITER and BEND at zenith angles between 50 • and 65 • being smaller than rms of CURV and ADVV. Similarly, at the zenith angles larger than 70 • , the near zenith bias and the modelling bias of CURV (missing the effect of the refractive bending) partly cancel each other. As a result, rms of CURV is smaller than rms of ADVV, despite the fact that ADVV is physically the most justified observation operator version presented in this article.
STD bias, sdev and rms at zenith direction, i.e. at ζ g,r = 0, can be compared with bias, sdev and rms of ZTD. The estimates of these are obtained by least squares fitting a straight line into the curves in zenith angle range 0-20 • . The resulting values are −1.00, 1.69 and 1.97 cm for bias, sdev and rms, respectively. This indicates that ZTD observation accuracy is better than STD observation accuracy at zenith direction, which is understood here as a consequence of the ZTD data processing which averages out the random errors of STD observations. Temporal averaging over 5-min time interval for ZTD production is another explanation for better statistics. Figure 6 shows time series of ζ g,r (solid curve) and the difference between modelled and observed STD (dots) for a signal path corresponding to receiver station Wettzell (49.14 • N, 12.88 • E; circled in Fig. 2) and GPS satellite SVN 15 for 12 May 2003 05-09 UTC. STD is modelled by observation operator version ADVV. During the overpass, observed STD decreases from 6.84 to 2.30 m as the satellite ascents to zenith. Since the NWP field is constant during this 4-hr period, the short time scale variations of Fig. 6 are interpreted as variations in observed STD. During the first 40 min of the time series, ζ g,r is over 50 • and random variations of STD are relatively large (of the order of 2 cm). At 6-9 UTC the random variations are smaller than 1 cm. The discontinuities in the time series result from the 20-min time slotting of the raw data in the geodetic processing algorithm. One interpretation of Fig. 6 is as follows: the accuracy of STD observations is not better than a few millimeters even at small zenith angles, and the accuracy is worse than a centimeter at large zenith angles. This interpretation is in contrast with the STD accuracy reported by Braun et al. (2001).

Summary and discussion
This article describes the observation operator for assimilation of STD observations by HIRLAM 3D-Var system. The observation operator consists of separate algorithms for determination of the GPS signal path from satellite to receiver across the model grid, projection of model variables to the signal path, and integration of neutral refractivity along the signal path to obtain STD. The signal path determination algorithm is based on the geometrical Tellus 58A (2006), 1 path (GP) approximation rather than on the theoretically more justified Least Travel Time principle.
As revealed by an experiment using standard atmospheric data, ZTD and STD modelled by the observation operator are fairly insensitive to the number of NWP model levels. The observation operator is able to closely reproduce ZTD observed at geodetic ground-based receiver stations. It is concluded that the refractivity integration is correctly performed.
Three possible simplifications of the signal path determination algorithm are studied in this article. Adopting all of the simplifications leads to neglection of the effects of the curvature of the Earth, the refractive bending and the horizontal model level height gradients (GP-version of the observation operator; see the main text). This is relevant for zenith angles smaller than about 40 • -50 • . For larger zenith angles the effect of the curvature of the Earth needs to be taken into account (CURVand ADVV-versions). The significance of accounting for the refractive bending of the signal path (BEND-and ADVV-versions) increases with zenith angle. However, the horizontal gradients of the model level heights (ITER-version) can safely be neglected when zenith angles are smaller than 80 • .
The STD modelled by the physically most justified observation operator version is systematically lower than the observed one. As measured by standard deviation of the difference between modelled and observed STD, the accuracy of the NWPbased STD estimates is about 1.7 cm at zenith direction. Standard deviation increases with increasing zenith angle and reaches 8.6 cm at zenith angle 80 • . While the systematic difference between modelled and observed STD is strongly dependent on the STD modelling details, the random difference is not.

Acknowledgments
This work is carried out within the EU FP5 project "TOUGH". TOUGH is a shared-cost project (contract EVG1-CT-2002-00080) co-funded by the Research DG of the European Commission within the RTD activities of the Environment and Sustainable Development sub-programme (5th Framework Programme). The observation operator presented in this article is a public deliverable D38 of TOUGH and it can be downloaded from http://web.dmi.dk/pub/tough/. Funding from TEKES, the National Technology Agency of Finland, is also thankfully acknowledged. The authors are grateful for Nils Gustafsson (Swedish Meteorological and Hydrological Institute), Jan Johansson (Onsala Space Observatory, Chalmers University of Technology, Sweden), Hans van der Marel (Technical University of Delft, the Netherlands) and John de Vries (Royal Netherlands Meteorological Institute) for fruitful discussions and for provision of the ZTD and STD observations. Three anonymous reviewers are kindly acknowledged for a number of suggestions for improving the manuscript and the numerical scheme of the observation operator.