4D Tropospheric Tomography using GPS Estimated Slant Delays

Tomographic techniques are successfully applied to obtain 4D images of the tropospheric refractivity in a local dense network. In the lower atmosphere both the small height and time scales and the non-dispersive nature of tropospheric delays require a more careful analysis of the data. We show how GPS data is processed to obtain the tropospheric slant delays using the GIPSY-OASIS II software and define the concept of pseudo-wet delays, which will be the observables in the tomographic software. We then discuss the inverse problem in the 3D stochastic tomography, using simulated refractivity fields to test the system and the impact of noise. Finally, we use data from the Kilauea network in Hawaii and a local 4x4x41-voxel grid on a region of 400 Km$^2$ and 15 Km in height to produce 4D refractivity fields. Results are compared with ECMWF forecast.


Introduction
One goal of ongoing research on the application of tomographic techniques to the atmosphere using GPS signals is to obtain 4D images of the refractivity in the troposphere using a ground network. In this paper we present our first results on the subject, showing that a careful combination of GPS data processing together with tomographic Kalman-filtering can reproduce the state of the refractivity in the neutral atmosphere.
The effect of the atmosphere on GPS signals is measured as an extra delay. The ionospheric electron content produces a dispersive delay and can hence be corrected by using two different frequencies and linearly combining them. The neutral atmosphere, however, induces a delay independent of frequency. This delay (∆L) can be related to the refractivity by ∆L = L 10 −6 N dl, (1.1) N ≈ 77.6 P T + 3.73 · 10 5 P w T 2 , (1.2) where P is the total atmospheric pressure (mbar), P w is the water vapor pressure (mbar), and T is the atmospheric temperature (K) ( [Smith et al.. 1953]). The extraction of this slant delay from GPS measurements requires accurate modeling. To take into account the dependence of the tropospheric delay on the satellite elevation, mapping functions are used, and a Zenith Total Delay (ZTD) is estimated at each station. To consider azimuthal variability, a horizontal gradient is also estimated. The time evolutions of these parameters is modeled as random walk stochastic processes. The total atmospheric delay can be partitioned into two components: the hydrostatic delay, due to the dry gases in the troposphere and the nondipole component of water vapor, and the wet delay, due to the dipole component of water vapor. The contribution of the hydrostatic component to the ZTD is larger than that of the wet and it can be estimated from surface pressure measurements. We remove the zenith hydrostatic delay (ZHD) from the total slant delays through an accurate estimation of the pressure according to site measurements, and then form the pseudo-wet delays (PWD): we map the time series zenith wet delay (ZWD) and total horizontal gradients back to each ray direction and add the residuals to account for any data mismodeling. The refractivity associated is termed pseudo-wet refractivity (N ). The PWD's are the observables in the tomographic approach. The reason for removing of the ZHD is two-fold. On the one hand, it reduces the scale height of the problem (typically, the hydrostatic component has a scale height of 10 Km while the wet component has a 2 Km scale height) and thus eases the numerical solution. In addition, theN field obtained is closely related to the water vapor distribution. In the paper we will first discuss the approach taken to estimate the PWD's using the JPL software GIPSY-OASIS II (GOA II) (see [Webb et al. 1997]) and the network used in the study. Then we briefly discuss the inverse problem and analyze the impact of the needed constraints through simulations. Finally, tomographic solutions with GPS data measured in the Kilauea network in Hawaii during 1st February 1997 and their comparison with independent measurement sources will be presented.

Tropospheric Slant Delay Estimation
The GPS signals L 1 and L 2 are first combined to remove the effect of the ionosphere. The resulting observable, L C , is modeled in GOA II using a Point Positioning Strategy in terms of tropospheric effects (ZWD+ZHD+Gradients), geometrical factors and receiver clocks. Satellite clock corrections and orbits are provided by JPL, as well as earth orientation parameters. Tropospheric delays are mapped to the zenith using Niell mapping functions ( [Niell 1996]) and related to the gradients following the expression ( [Davis et al. 1993]): (2.1) The time dependence of the tropospheric component is handled using a kalman filter and a random walk stochastic process, with a drift rate of 0.25 cm/ √ h for the zenith delay and 0.03 cm/ √ h ( [Bar-Sever et al. 1998]) for the gradients. We then form the PWD from the solution as discussed above.

Network and data considered
GPS data was tracked with a sampling of 30 s in the Kilauea network in Hawaii. A map with the stations is shown in Figure 3. It can be seen that they cover an area of about 400 Km 2 . The heights of the stations are distributed in a range of about 2000 m, which seems particularly well-suited for vertical resolution in the tomographic solution.

The inverse problem
In tomography, one wants to obtain the solution fields (refractivity in the tropospheric case) from the integrated value along a ray path. The amount and distribution of data, however, is often insufficient to determine a unique solution. Thus, some additional information has to be added to the system. To this end, we rewrite the pseudo-wet refractivity asN ( is any set of basis functions and ǫ( r, t) is the quantization error. In our tomographic approach, three dimensional voxels are used. Then, our observations are modeled by This is a set of linear equations of the form Ax = y. As mentioned, the system may not have a solution and thus, we seek to minimize the functional χ 2 (x) = (y − Ax) T · (y − Ax). Even so, the unknowns may not be fixed by the data, although the number of equations is usually much greater than the number of unknowns. This is so because we are trying to obtain a solution with more degrees of freedom (more voxels) than the resolution carried by the data. In previous work, (see [Ruffini et al. 1998] and [Rius et al. 1997 ]) we have discussed the use of a correlation functional to confine the spatial spectrum of the solution to the low portion of the frequency space. The same concept can be expressed by adding new equations (constraints) that impose that the density in a voxel be a weighted average of its neighbours. Let us now take a closer look to the constraints. Horizontal smoothing is needed to account for the nonuniform distribution of the rays. Water vapor is mostly concentrated in the first 2 Km in height. Now, if we aim for a tomographic solution with a resolution of 350 m (a trade-off between data resolution and computer load), some smoothing constraints have to be added in the vertical direction. Finally, we impose a zero value constraint in the highest layer. The constraints are represented by the equation B · x = 0. To smooth out the time variability, a Kalman filtering was used, modeling the troposphere as a random walk stochastic process with a drift rate of δ = 0.14 cm/(Km √ h). The tomographic process is integrated in our home-made software package LOTTOS (LOcal Tropospheric TOmography Software).

Results with simulated data
In order to test the software for the network's geometry and tune the different parameters, simulations play an essential role. We have simulated a 3D pseudo-wet refractivity field following the expression (following [Davis et al. 1993 where h w = 2 Km, h d = 10 Km, g d and g w are the horizontal gradients up to a multiplicative constant in the hydrostatic and wet component, respectively and ρ is the horizontal coordinate. We have set N W 0 = 150 mm/Km and N D 0r = 2 mm/Km, the latter to account for any residual hydrostatic component, g d = 0 and two different non-zero g w values, applied depending on the latitude of the station. The geometry of the rays for the 1st February 1997 in the Kilauea network has been used to generate the simulated rays. The PWD for each ray has been formed according to L sim = voxN i · l i , whereN i is the pseudo-wet refractivity associated to voxel i and l i is the length of the ray across this voxel. The weight of the constraints was tuned to give the best fit with simulations, successively adding the constraints in the horizontal plane, then the vertical smoothing constraints, and finally adding a zero value constraint for the highest layer. In Figure 1 we show the error in reconstructing the simulated field. Vertical reconstruction has also been verified adding a perturbation at different heights to the simulated fields of Equation 3.1; good vertical resolution is achieved thanks to the distribution in height of the stations. We also note (see Figure 1) that vertical reconstruction is gradually lost as the perturbation is placed at higher altitudes (above 2 Km, resolution starts degrading). However, when the boundary constraint kicks in, the solution agrees again with the simulated field, slowly decaying to zero.
The impact of noise in the rays has also been considered. Following the data processing described above, we have seen that post-fit residuals are dependent on elevation (roughly as 1/sin(e), where e is the elevation); therefore, we can assume a noise value for the ZWD, then map it to the slant directions, and add it to the total simulated delay. We show in Figure 2 the results of the noise analysis. As it can be seen, a ZWD noise of 0.5 cm represents a noise in the solution of less than 3 mm/Km rms (about 2% of the maximum value).

Results from the Kilauea network
Pseudo-wet delays from the Kilauea network in Hawaii have been input to the system. We have processed data from 18 stations for February 1st, 1997. We have used LOTTOS considering 4x4 voxels of 6 ′ in latitude and 7 ′ in longitude, 41 layers of 350 m height and 30-minutes batches for the Kalman filter. The tomographic fields have then been used to generate the reconstructed PWD which, in turn, have been processed to calculate the horizontal gradients and ZWD for each station. These have been compared to those obtained with GOA II. In Figure 3 we show the magnitude and direction of the 24-h mean value of the gradient as calculated using GOA II (yellow) and using the tomographic solution (green). In Figure 4 we show the time series of these magnitudes for a particular station (PGF3) and in Table 1 we show the correlation values over time for each station. This shows that the reconstructed 4D field is consistent with the data. Vertical distribution is shown in Figure 5 for a given latitude and time. We have computed the values ofN ( r, t) profiles with data from the European Center for Medium-Range Weather Forecasts (ECMWF) analysis and compared with our tomographic profiles, averaging them to meet the ECMWF maps resolution of 0.5 degrees. Results are presented for three different hours (06h, 12h and 18h) in Figure 6 showing good agreement.

Conclusions
We have successfully obtained tomographic 4D images of the pseudo-wet refractivity for a local dense network. The concept of slant pseudo-wet delays has been introduced. Simulated fields have been used to validate and tune our LOTTOS software package. Finally, data from the Kilauea network have been processed and the reconstructed 4D fields of the refractivity have been compared with on-site estimated horizontal gradients and ECMWF vertical profiles, yielding good correlation in both cases. Although more work is needed in this area, our results provide the proof-of-concept of the tomographic inversion of tropospheric GPS data.