- Email: [email protected]
Tomographic determination of the spatial distribution of water vapor using GPS observations
Advances in Space Research 37 (2006) 2211–2217 www.elsevier.com/locate/asr
Tomographic determination of the spatial distribution of water vapor using GPS observations M. Troller a
a,*
, A. Geiger a, E. Brockmann b, J.-M. Bettems c, B. Bu¨rki a, H.-G. Kahle
a
Geodesy and Geodynamics Laboratory, Institute of Geodesy and Photogrammetry, Swiss Federal Institute of Technology, ETH Hoenggerberg, HPV G56, Schafmattstrasse 34, 8093 Zuerich, Switzerland b Swiss Federal Office of Topography, Geodetic Bases and Permanent Networks, Wabern, Switzerland c Federal Office of Meteorology and Climatology, Zuerich, Switzerland Received 10 September 2004; received in revised form 5 July 2005; accepted 6 July 2005
Abstract With the advent of the GPS navigation system, a promising ground based technique has been introduced which makes it possible to estimate the amount of water vapor in the troposphere from operational GPS networks at relatively low additional costs. While the estimation of the integrated amount is currently well established, the determination of the spatial water vapor distribution and its temporal variation are still a major challenge. To account for the vertical resolution, several tomographic approaches were pursued. We developed the software package AWATOS (atmospheric water vapor tomography software) which is based on the assimilation of double differenced GPS observations. Applying a least-squares inversion, the inhomogeneous spatial distribution of water vapor is determined. An extensive investigation has been carried out in Switzerland. GPS measurements are performed by the dense permanent Swiss national GPS network AGNES of the Swiss Federal Office of Topography (swisstopo). A total of 40 equally distributed water vapor profiles have been estimated on an hourly basis. For the purpose of validation, 22 radiosonde profiles were used at the GPS and meteorological station Payerne. Furthermore, data of the numerical weather model aLMo (alpine model in Switzerland, MeteoSwiss) were compared with the tomographic results. An overall good agreement of the three methods with an rms of better than 1.6 g/m3 absolute humidity was achieved. The results show that AGNES can be used as a dedicated network for the purpose of GPS-tomography, using a horizontal resolution of approximately 50 km and height layers of 300–500 m thickness in the lower troposphere. 2005 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: GPS; GPS meteorology; Water vapor; Refractivity; Tomography; Satellite navigation; Numerical weather prediction; Quantitative precipitation forecast
1. Introduction The next generation of numerical weather prediction models currently being developed, with a mesh size of 1–3 km, will require a precise analysis of the threedimensional water vapor distribution. Specifically, quantitative precipitation forecast with explicit simula*
Corresponding author. Tel.: +41 44 633 21 23; fax: +41 44 633 10
66. E-mail address: [email protected] (M. Troller).
tion of convection at the kilometric scale are very sensitive to temperature and humidity structures (Ducrocq et al., 2002; Park, 1999). Moreover, these structures show a high temporal and spatial variability, which must be well captured by the assimilation procedures of the numerical weather prediction model. With the planned use of rapid update cycles, short forecasts frequently repeated, the water vapor analysis should be available many times a day. Currently, the main source of humidity information over land in operational numerical weather prediction
0273-1177/$30 2005 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2005.07.002
2212
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
(NWP) model are radio soundings, with a poor horizontal and temporal resolution, and humidity meters at synoptic stations. GPS meteorology is a promising technique to estimate the total amount of water vapor in the troposphere on a continuous basis using satellite navigation systems (e.g., Bevis et al., 1992). With GPS tomography, it is aimed at determining also the vertical distribution of water vapor in the troposphere with a high temporal resolution (e.g., Flores et al., 2000; Hirahara, 2000; Kruse, 2001; Troller, 2004). In this study, hourly resolution has been investigated, which represents a major improvement compared to time resolution obtained from operational ballon soundings. The development of this approach could be an important source of data for the next generation of numerical weather models, at least over regions with a high density GPS network.
2. Principle of GPS tomography GPS signals are significantly influenced by the atmosphere, especially the ionosphere and troposphere, along their path from the satellite to the GPS antenna. The satellites are transmitting at two different carrier frequencies in the radio L-band. As the ionosphere is a dispersive medium in the radio frequency range, its influence can be eliminated by composing a linear combination of the two GPS carrier frequencies. The remaining effect is caused by the delay of the signal due to the refractivity of the troposphere. Furthermore, this effect can be subdivided into a so-called dry and wet part, the latter being proportional to the integrated precipitable water vapor (e.g., Hofmann-Wellenhof et al., 2001). The main advantages of GPS are that (1) measurements can be carried out on a continuous basis at relatively low expense; (2) GPS is an all-weather observing system, almost insensitive to clouds, and the effect of liquid water can usually be neglected (Elgered, 1993). The wet propagation delay DPD wet of a radio signal from the satellite to the receiver antenna is defined as Z satellite Z satellite 6 DPD ¼ ðn 1Þ ds ¼ 10 N wet ds; ð1Þ wet wet antenna
antenna
where nwet represents the refractive index due to water vapor and Nwet the wet refractivity whereby latter is defined as Nwet = 106 · (nwet 1); ds is the length of a ray path element. The ray bending effect can be neglected, as the cutoff angle of 10 (Mendes, 1999) has been applied throughout the presented evaluation. However, the wet path delay reflects the integrated amount of water vapor, only. To achieve the spatial distribution, a tomographic approach was investigated.
In the tomographic approach, a discretization of the atmosphere with a voxel model is used. The number of layers and the horizontal and vertical size of voxels are correspondingly defined. Ellipsoidal boundary surfaces of the layers account for the earthÕs curvature. At the horizontal boundaries in each layer open voxels are added, i.e., these voxels reach ad infinitum. Consequently, no rays are neither completely nor partially outside the model. Within each voxel, the wet refractivity Nwet is introduced as an unknown constant. Discretizing Eq. (1), the wet slant path delay from satellite r to receiver antenna p as observed is related to the unknowns (refractivity Nwet) as follows: 6 DPD;r wet;p ¼ 10
k X
N wet;i Dsri;p ;
ð2Þ
i¼1
where Dsri;p represents the length of the path through voxel i, and k the number of voxels. Wet slant path delays can be derived from several remote sensing techniques such as GPS, water vapor radiometry and solar spectrometry. In the current approach, we focus on GPS double difference path delays. By following the well-established concept of double differencing (e.g., Beutler et al., 2001), we obtain from (Eq. (2)) the double difference wet slant path delays D2;PD;rs wet;pq of the satellites r and s and the stations p and q PD;r PD;r PD;s PD;s D2;PD;rs ð3Þ wet;pq ¼ Dwet;q Dwet;p Dwet;q Dwet;p . The following approach has been applied to retrieve the double difference path delays: GPS data have been processed using the Bernese GPS Software (Beutler et al., 2001) yielding GPS zenith path delays DPD and double difference phase residuals D2U as processing output. 2;PD;rs The double difference path delays Dpq of the satellites r, s and the stations p, q are reconstructed from the GPSPD derived zenith path delays DPD p ; Dq , which have to be mapped back to the corresponding elevations using the corresponding mapping functions (airmass factors) mðelrp Þ; mðelrq Þ; mðelsp Þ; mðelsq Þ. Furthermore, the double difference phase residual D2 Urs pq is added 2;PD;rs D2;PD;rs ¼ Dpq þ D2 Urs pq ; pq
ð4Þ
where r r PD D2;PD;rs ¼ DPD pq q mðelq Þ Dp mðelp Þ s s PD DPD q mðelq Þ Dp mðelp Þ .
ð5Þ
Using ground meteorological measurements, the zenith dry path delay can be determined by applying the formula of Saastamoinen (Saastamoinen, 1972; Troller, 2004). Subsequently, the double difference dry slant path delay has been constructed using the same approach as in Eq. (5). Finally, the wet part of the double difference
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
slant path delay D2;PD;rs wet;pq is extracted by subtracting the dry part from the total amount. The number of traversing rays per voxel depends on the geometry defined by the number of visible satellites, the distribution of the ground stations and the size of the voxels. GPS-tomography usually generates a partly illposed problem in the sense that only a portion of the unknowns (refractivity per voxel) can be determined whereas the other part is under-determined. Additional information is necessary to solve the equation system. In this approach, the uppermost layer is bounded by 8000 and 15,000 m lower and upper level, respectively. This allows for the assumption of a priori values of Nwet = 0 ppm in that level. In special cases, a priori values can be introduced for other voxels. This procedure is described in Section 4. In addition, the voxels are mutually coupled by using a realistic covariance function to limit the variation of the difference in refractivity of neighboring voxels. However, only direct neighboring voxels are correlated to allow for steep refractivity gradients. Extensive simulations using various weather conditions show the feasibility of this approach and indicate the need of double difference measurements with noise not exceeding 5 mm (Troller et al., 2002). By properly treating the GPS data a double difference noise limit
2213
of 4–5 mm can be achieved in most cases, thus allowing for a tomographic solution with sufficient accuracy.
3. Description of the experiment The observations have been carried out in Switzerland. The data used are collected in the Swiss permanent GPS reference network (AGNES) of swisstopo. It covers the entire Swiss territory with a dense horizontal and vertical resolution (Fig. 1). Large height differences between the GPS stations are suitable for an accurate tomographic solution. An automated near real-time processing provides hourly means of zenith total delays, which are then used to derive instantaneous values of GPS path delays. Data of the meteorological network ANETZ of MeteoSwiss (SMA, 1985) is used to decompose the total delays in its wet and dry part. ANETZ contains a total of 72 well-distributed ground stations covering Switzerland and allows an accurate splitting of the zenith total delay. A period of one week in November 2002 was chosen for the investigations. Rapid weather changes, heavy rainfall, clear conditions and sunny periods occurred during that week.
Fig. 1. GPS stations of the Swiss Permanent GPS Reference Network AGNES (swisstopo). The network contains 30 AGNES stations, with a height distribution from 400–3600 m. For the automated processing, 20 EUREF stations and 23 stations from other networks are included. The figure shows the stations of the Swiss territory only.
2214
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
Fig. 2. 3D view of the evaluation perimeter. The tomographic voxel model consists of 16 layers up to 15,000 m height. The borders of the layers are set at 0, 200, 600, 900, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2700, 3200, 4000, 5000, 8000, and 15,000 m (the figure shows layers up to 5000 m height only). Each layer contains six voxels in longitude and three voxels in latitude (spacing 0.5) plus 22 outer voxels (not shown on the figure). The GPS stations are shown as column according to their station height. Radiosondes are launched from the radiosonde station Payerne (MeteoSwiss), shown as cuboid.
The voxel model above the Swiss territory is chosen with 6 · 3 voxels in horizontal (voxel side lengths 50 km) and 16 layers up to 15,000 m (Fig. 2). Water vapor profiles are retrieved on an hourly basis. The tomographic results are compared with data of the operational high resolution NWP model of MeteoSwiss (aLMo). aLMo is the Swiss implementation of the non-hydrostatic local weather model of COSMO (Doms and Scha¨ttler, 2003). The model is operational at MeteoSwiss since April 2001. It covers most of Western Europe, has a mesh size of about 7 km in the horizontal and uses 45 vertical layers, with a thickness of about 100 m in the lower troposphere. According to the space and time resolution of the tomographic evaluation, a total of 7680 aLMo profiles and radio soundings at the station Payerne (twice a day) are made available by MeteoSwiss for comparison purposes.
is assigned to the uppermost layer (8000–15,000 m) with a regularization factor of 1/900. On the other hand, AWATOS ANETZ contains the same constraints as AWATOS Correlation and in addition one a priori wet refractivity value for each voxel lower than 2000 m.
4. Profile determination and evaluation Double difference residuals and hourly means of zenith path delays are determined using the Bernese software package (Beutler et al., 2001). Reconstructed double-difference wet slant delays are then introduced into the tomographic software package AWATOS. Two types of tomographic profiles are determined, containing different constraints. On one hand, AWATOS Correlation includes inter-voxel constraints between all neighboring voxels. Compared to a double-difference observation, the constraints are downweighted by a factor of 502 (regularization factor = 1/2500). In addition, an a priori wet refractivity of zero
Fig. 3. Cross section of the square root of cofactors at latitude 46.75. The two top layers (5000–15,000 m height) are not shown in the figures. Panel (a) represents the situation of the solution AWATOS Correlation and (b) AWATOS ANETZ. It is seen, that if a priori refractivity is introduced, the precision increases significantly.
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
These refractivity values are based on a collocation and interpolation procedure (COMEDIE software package; Troller et al., 2000), using the operationally available meteorological ground station data (ANETZ, MeteoSwiss). The least-squares adjustment of the tomographic equation system yields the matrices of the cofactors of the unknown parameters from the inverted normal matrices (e.g., Mikhail, 1976). The comparison of the two cofactor matrices of the two solutions reveals the impact on quality of the two different types of solution. The a priori calculation takes no measurements into account. It assesses the geometry and the weighting of the measurements, only. The solution AWATOS Correlation shows a slight increase of the precision with height. Regarding the solution AWATOS ANETZ, a significant improvement of the square root of cofactors is visible in the first 2000 m height but also in the higher voxels (Fig. 3). In Fig. 4, two examples of tomographic profiles, aLMo data and radiosonde profiles are displayed. Furthermore, the profile obtained with COMEDIE is plotted. The lowermost 2000 m of the latter profile are used to constrain the solution AWATOS ANETZ. The
2215
comparison of the profiles confirms the conclusions made of the matrices of cofactors: The degree of agreement of the individual profiles is varying. AWATOS ANETZ fits always more accurately to aLMo than AWATOS Correlation. The 22 radiosonde launches during the investigation week at station Payerne have been used to perform a comparison to the tomographic profiles and the aLMo model (Table 1). The conclusion of Fig. 4 can be verified, the mean rms of AWATOS ANETZ is half as large as of AWATOS Correlation. Fig. 5 shows the mean rms of the tomographic profiles compard to the radiosondes as function of the height. The impact of the a priori refractivity in the first 2000 m height is clearly visible. The comparison of aLMo to the radiosonde profiles is also documented in Table 1. The agreement is within 2.6 ppm (refractivity units) at station Payerne, where radiosondes are available to compare the aLMo model. As radiosondes are available every six or twelf hours and at the station Payerne only, the aLMo data with an hourly resolution have been used to compare the tomographic solutions over entire Switzerland. Thus, a total of 3024 profiles are compared to aLMo for the one-week measurements (Fig. 6). In Fig. 4, two profiles with
5000
5000
GPS Tomography
GPS Tomography
06.11.02 17h
4500
AWATOS Correlation AWATOS ANETZ aLMo solution COMEDIE profile Radiosonde
4000
3500
3500
height [m]
height [m]
3000
2500
2500
2000
2000
1500
1500
1000
1000
500
500
a 0
AWATOS Correlation AWATOS ANETZ aLMo solution COMEDIE profile Radiosonde
4000
3000
0
06.11.02 23h
4500
b
© M. Troller, GGL-ETHZ
10
20
30
40
wet refractivity [ppm]
50
60
0
0
© M. Troller, GGL-ETHZ
10
20
30
40
50
60
wet refractivity [ppm]
Fig. 4. Wet refractivity profiles of the two tomographic solutions, the radiosonde ascending at the station Payerne (MeteoSwiss), data of the numerical weather model aLMo (MeteoSwiss) and the COMEDIE profile. Panel (a) shows a profile with a steep refractivity gradient at 1800 m height. All solutions match accurately together on the order of 10 ppm. Panel (b) shows a situation with an inhomogeneous decrease of refractivity with height. AWATOS Correlation underestimates the wet refractivity in the lowest 1000 m height. However, also the other tomographic solution as well as the radiosonde show differences to the aLMo solution.
2216
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
Table 1 Statistical analysis of the tomographic profiles and aLMo compared to the radiosondes
Mean offset Mean rms Mean r
AWATOS Correlation
AWATOS ANETZ
aLMoa
6.3 12.7 11.0
1.4 5.1 4.9
0.5 2.6 2.6
The analysis contains 22 profiles. The accuracy of the two tomographic solutions varies significantly. The mean rms of AWATOS Correlation is approximately double as the corresponding value of AWATOS ANETZ. Furthermore, the latter has no significant mean offset. a aLMo data for the lowermost voxel (200–600 m height) are not available.
Fig. 5. Mean rms of the two tomographic solutions compared to the radiosondes as a function of the height. The plot shows the comparison of the 22 radiosonde launches during the investigation week at the station Payerne. Generally, the rms is decreasing with increasing height. The impact of the a priori profiles on the first 2000 m height from COMEDIE is clearly visible in the solution AWATOS ANETZ.
18
AWATOS Correlation AWATOS ANETZ
16
mean rms [ppm]
14 12 10 8 6 4
9, coincide with atmospheric situations rapidly changing from heavy rain to dry periods and contrariwise. During these time periods, the differences between GPS-tomography and aLMo are larger. This may mainly be caused by the temporal averaging during each data assimilation period, thus neglecting gradients in time for the assimilation interval.
5. Conclusions and outlook The tomographic approach was successfully performed to process GPS data for the determination of 4 dimensional water vapor data. The investigations showed that the Swiss national reference network AGNES is a dedicated network for GPS tomography if using a grid model with 6 · 3 voxels horizontally. AGNES has a station density of about one GPS receiver per voxel column. Regularization methods are necessary to achieve a stable tomographic solution. Inter-voxel constraints and a priori refractivity information allow usually to achieve a high accuracy. An overall rms of about 5–8 ppm (refractivity units) or 0.8–1.3 g/m3 absolute humidity compared to aLMo has been reached, depending on the tomographic solution. Using a temporal resolution of 1 h, the accuracy remains approximately stable during the whole week, at least for the solution AWATOS ANETZ. Only during rapidly changing atmospheric situations, the rms is slightly decreasing. Further investigations showed (Troller, 2004) that the accuracy achieved at station Payerne, can be expected in whole Switzerland. Nevertheless, in a further development, an improvement of the regularization methods should be aimed at. The use of humidity profiles over land determined by satellites as well as a priori data from NWP should be studied. Further investigations should focus on the augmentation of the tomographic spatial and temporal resolution. The introduction of data from GPS-independent sensors such as radiometers or solar spectrometers should be considered. The impact of the derived information on the quality of NWP forecast should then be thoroughly evaluated.
2 0 3
4
5
6
7
8
9
10
Day of November 2002
Fig. 6. Time series of the mean of all profiles compared to aLMo. Overall, AWATOS ANETZ has a smaller rms than AWATOS Correlation. Three jumps are visible on November 3, 6 and 9.
different levels of agreement, mainly depending on the atmosphereÕs actual state are shown. This behavior can be confirmed by inspecting the time series analysis of Fig. 6. Periods with large rms on November 3, 6 and
References Beutler, G., Bock, H., Brockmann, E., Dach, R., Fridez, P., Gurtner, W., Hugentobler, U., Ineichen, D., Johnson, J., Meindl, M., Mervart, L., Rothacher, M., Schaer, S., Springer, T., Weber, R. in: Hugentobler, U., Schaer, S., Fridez, P. (Eds.), Bernese GPS Software Version 4.2. Astronomical Institute, University of Berne, 2001. Bevis, M., Businger, S., Herring, T., Rocken, C., Anthes, R., Ware, R. GPS meteorology: remote sensing of atmospheric water vapor using the global positioning system. J. Geophys. Res. 97 (14), 15787–15801, 1992.
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217 Doms, G., Scha¨ttler, U. (Eds.). Newsletter No. 3, February 2003. Technical Report, COSMO (Consortium for Small-scale Modelling), 2003. Ducrocq, V., Ricard, D., Lafore, J.-P., Orain, F. Storm-scale numerical rainfall prediction for five precipitating events over France: on the importance of the initial humidity field. Weather Forecast. 17, 1236–1256, 2002. Elgered, G. Atmospheric Remote Sensing by Microwave Radiometry, Chapter Tropospheric Radio-Path Delay from Ground-based Microwave Radiometry, Wiley, New York, pp. 215–258, 1993. Flores, A., Ruffini, G., Rius, A. 4d tropospheric tomography using GPS slant wet delays. Ann. Geophys. 18, 223–234, 2000. Hirahara, K. Local GPS tropospheric tomography. Earth Planets Space 52, 935–939, 2000. Hofmann-Wellenhof, B., Lichtenegger, H., Collins, J. GPS Theory and Practice. Springer, New York, 2001. Kruse, L. Spatial and temporal distribution of atmospheric water vapor using space geodetic techniques. Geoda¨tisch-geophysikalische Arbeiten in der Schweiz, vol, 61, Schweizerische Geoda¨tische Kommision, 2001. Mendes, V.B. Modeling the neutral-atmosphere propagation delay in radiometric space techniques. Ph.D. dissertation, University of New Brunswick, 1999.
2217
Mikhail, E.M. Observations and Least Squares. University Press of America, 1976. Park, S.K. Nonlinearity and predictability of convective rainfall associated with water vapour perturbations in a numerically simulated storm. J. Geophys. Res. (Atmospheres) 104, 31,575– 31,587, 1999. Saastamoinen, J. Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. In: Henriksen, S.W., Mancini, A., Chovitz, B.H. (Eds.). The Use of Artificial Satellites for geodesy, Geophys. Monogr. Ser. AGU, 15, pp. 247–251, 1972. SMA. Charakteristiken der ANETZ-Daten. In: Beitra¨ge zum ANETZDaten Kolloquium vom 17. April 1985 in Zu¨rich, p. 56. Schweizerische Meteorologische Anstalt Zu¨rich, 1985. Troller, M. GPS based determination of the integrated and spatially distributed water vapor in the troposphere. Diss ETH No. 15513, 2004. Troller, M., Bu¨rki, B., Cocard, M., Geiger, A., Kahle, H.-G. 3d refractivity field from GPS double difference tomography. Geophys. Res. Lett. 29 (24), 2149–2152, 2002. Troller, M., Cocard, M., Geiger, A. Modellierung 4 dimensionaler Refraktionsfelder zur Berechnung von Wegla¨ngen-Korrekturen bei Satellitenmessungen. In: Simulation raumbezogener Prozesse: Methoden und Anwendungen, 9, IfGI prints, 2000.
Tomographic determination of the spatial distribution of water vapor using GPS observations M. Troller a
a,*
, A. Geiger a, E. Brockmann b, J.-M. Bettems c, B. Bu¨rki a, H.-G. Kahle
a
Geodesy and Geodynamics Laboratory, Institute of Geodesy and Photogrammetry, Swiss Federal Institute of Technology, ETH Hoenggerberg, HPV G56, Schafmattstrasse 34, 8093 Zuerich, Switzerland b Swiss Federal Office of Topography, Geodetic Bases and Permanent Networks, Wabern, Switzerland c Federal Office of Meteorology and Climatology, Zuerich, Switzerland Received 10 September 2004; received in revised form 5 July 2005; accepted 6 July 2005
Abstract With the advent of the GPS navigation system, a promising ground based technique has been introduced which makes it possible to estimate the amount of water vapor in the troposphere from operational GPS networks at relatively low additional costs. While the estimation of the integrated amount is currently well established, the determination of the spatial water vapor distribution and its temporal variation are still a major challenge. To account for the vertical resolution, several tomographic approaches were pursued. We developed the software package AWATOS (atmospheric water vapor tomography software) which is based on the assimilation of double differenced GPS observations. Applying a least-squares inversion, the inhomogeneous spatial distribution of water vapor is determined. An extensive investigation has been carried out in Switzerland. GPS measurements are performed by the dense permanent Swiss national GPS network AGNES of the Swiss Federal Office of Topography (swisstopo). A total of 40 equally distributed water vapor profiles have been estimated on an hourly basis. For the purpose of validation, 22 radiosonde profiles were used at the GPS and meteorological station Payerne. Furthermore, data of the numerical weather model aLMo (alpine model in Switzerland, MeteoSwiss) were compared with the tomographic results. An overall good agreement of the three methods with an rms of better than 1.6 g/m3 absolute humidity was achieved. The results show that AGNES can be used as a dedicated network for the purpose of GPS-tomography, using a horizontal resolution of approximately 50 km and height layers of 300–500 m thickness in the lower troposphere. 2005 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: GPS; GPS meteorology; Water vapor; Refractivity; Tomography; Satellite navigation; Numerical weather prediction; Quantitative precipitation forecast
1. Introduction The next generation of numerical weather prediction models currently being developed, with a mesh size of 1–3 km, will require a precise analysis of the threedimensional water vapor distribution. Specifically, quantitative precipitation forecast with explicit simula*
Corresponding author. Tel.: +41 44 633 21 23; fax: +41 44 633 10
66. E-mail address: [email protected] (M. Troller).
tion of convection at the kilometric scale are very sensitive to temperature and humidity structures (Ducrocq et al., 2002; Park, 1999). Moreover, these structures show a high temporal and spatial variability, which must be well captured by the assimilation procedures of the numerical weather prediction model. With the planned use of rapid update cycles, short forecasts frequently repeated, the water vapor analysis should be available many times a day. Currently, the main source of humidity information over land in operational numerical weather prediction
0273-1177/$30 2005 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2005.07.002
2212
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
(NWP) model are radio soundings, with a poor horizontal and temporal resolution, and humidity meters at synoptic stations. GPS meteorology is a promising technique to estimate the total amount of water vapor in the troposphere on a continuous basis using satellite navigation systems (e.g., Bevis et al., 1992). With GPS tomography, it is aimed at determining also the vertical distribution of water vapor in the troposphere with a high temporal resolution (e.g., Flores et al., 2000; Hirahara, 2000; Kruse, 2001; Troller, 2004). In this study, hourly resolution has been investigated, which represents a major improvement compared to time resolution obtained from operational ballon soundings. The development of this approach could be an important source of data for the next generation of numerical weather models, at least over regions with a high density GPS network.
2. Principle of GPS tomography GPS signals are significantly influenced by the atmosphere, especially the ionosphere and troposphere, along their path from the satellite to the GPS antenna. The satellites are transmitting at two different carrier frequencies in the radio L-band. As the ionosphere is a dispersive medium in the radio frequency range, its influence can be eliminated by composing a linear combination of the two GPS carrier frequencies. The remaining effect is caused by the delay of the signal due to the refractivity of the troposphere. Furthermore, this effect can be subdivided into a so-called dry and wet part, the latter being proportional to the integrated precipitable water vapor (e.g., Hofmann-Wellenhof et al., 2001). The main advantages of GPS are that (1) measurements can be carried out on a continuous basis at relatively low expense; (2) GPS is an all-weather observing system, almost insensitive to clouds, and the effect of liquid water can usually be neglected (Elgered, 1993). The wet propagation delay DPD wet of a radio signal from the satellite to the receiver antenna is defined as Z satellite Z satellite 6 DPD ¼ ðn 1Þ ds ¼ 10 N wet ds; ð1Þ wet wet antenna
antenna
where nwet represents the refractive index due to water vapor and Nwet the wet refractivity whereby latter is defined as Nwet = 106 · (nwet 1); ds is the length of a ray path element. The ray bending effect can be neglected, as the cutoff angle of 10 (Mendes, 1999) has been applied throughout the presented evaluation. However, the wet path delay reflects the integrated amount of water vapor, only. To achieve the spatial distribution, a tomographic approach was investigated.
In the tomographic approach, a discretization of the atmosphere with a voxel model is used. The number of layers and the horizontal and vertical size of voxels are correspondingly defined. Ellipsoidal boundary surfaces of the layers account for the earthÕs curvature. At the horizontal boundaries in each layer open voxels are added, i.e., these voxels reach ad infinitum. Consequently, no rays are neither completely nor partially outside the model. Within each voxel, the wet refractivity Nwet is introduced as an unknown constant. Discretizing Eq. (1), the wet slant path delay from satellite r to receiver antenna p as observed is related to the unknowns (refractivity Nwet) as follows: 6 DPD;r wet;p ¼ 10
k X
N wet;i Dsri;p ;
ð2Þ
i¼1
where Dsri;p represents the length of the path through voxel i, and k the number of voxels. Wet slant path delays can be derived from several remote sensing techniques such as GPS, water vapor radiometry and solar spectrometry. In the current approach, we focus on GPS double difference path delays. By following the well-established concept of double differencing (e.g., Beutler et al., 2001), we obtain from (Eq. (2)) the double difference wet slant path delays D2;PD;rs wet;pq of the satellites r and s and the stations p and q PD;r PD;r PD;s PD;s D2;PD;rs ð3Þ wet;pq ¼ Dwet;q Dwet;p Dwet;q Dwet;p . The following approach has been applied to retrieve the double difference path delays: GPS data have been processed using the Bernese GPS Software (Beutler et al., 2001) yielding GPS zenith path delays DPD and double difference phase residuals D2U as processing output. 2;PD;rs The double difference path delays Dpq of the satellites r, s and the stations p, q are reconstructed from the GPSPD derived zenith path delays DPD p ; Dq , which have to be mapped back to the corresponding elevations using the corresponding mapping functions (airmass factors) mðelrp Þ; mðelrq Þ; mðelsp Þ; mðelsq Þ. Furthermore, the double difference phase residual D2 Urs pq is added 2;PD;rs D2;PD;rs ¼ Dpq þ D2 Urs pq ; pq
ð4Þ
where r r PD D2;PD;rs ¼ DPD pq q mðelq Þ Dp mðelp Þ s s PD DPD q mðelq Þ Dp mðelp Þ .
ð5Þ
Using ground meteorological measurements, the zenith dry path delay can be determined by applying the formula of Saastamoinen (Saastamoinen, 1972; Troller, 2004). Subsequently, the double difference dry slant path delay has been constructed using the same approach as in Eq. (5). Finally, the wet part of the double difference
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
slant path delay D2;PD;rs wet;pq is extracted by subtracting the dry part from the total amount. The number of traversing rays per voxel depends on the geometry defined by the number of visible satellites, the distribution of the ground stations and the size of the voxels. GPS-tomography usually generates a partly illposed problem in the sense that only a portion of the unknowns (refractivity per voxel) can be determined whereas the other part is under-determined. Additional information is necessary to solve the equation system. In this approach, the uppermost layer is bounded by 8000 and 15,000 m lower and upper level, respectively. This allows for the assumption of a priori values of Nwet = 0 ppm in that level. In special cases, a priori values can be introduced for other voxels. This procedure is described in Section 4. In addition, the voxels are mutually coupled by using a realistic covariance function to limit the variation of the difference in refractivity of neighboring voxels. However, only direct neighboring voxels are correlated to allow for steep refractivity gradients. Extensive simulations using various weather conditions show the feasibility of this approach and indicate the need of double difference measurements with noise not exceeding 5 mm (Troller et al., 2002). By properly treating the GPS data a double difference noise limit
2213
of 4–5 mm can be achieved in most cases, thus allowing for a tomographic solution with sufficient accuracy.
3. Description of the experiment The observations have been carried out in Switzerland. The data used are collected in the Swiss permanent GPS reference network (AGNES) of swisstopo. It covers the entire Swiss territory with a dense horizontal and vertical resolution (Fig. 1). Large height differences between the GPS stations are suitable for an accurate tomographic solution. An automated near real-time processing provides hourly means of zenith total delays, which are then used to derive instantaneous values of GPS path delays. Data of the meteorological network ANETZ of MeteoSwiss (SMA, 1985) is used to decompose the total delays in its wet and dry part. ANETZ contains a total of 72 well-distributed ground stations covering Switzerland and allows an accurate splitting of the zenith total delay. A period of one week in November 2002 was chosen for the investigations. Rapid weather changes, heavy rainfall, clear conditions and sunny periods occurred during that week.
Fig. 1. GPS stations of the Swiss Permanent GPS Reference Network AGNES (swisstopo). The network contains 30 AGNES stations, with a height distribution from 400–3600 m. For the automated processing, 20 EUREF stations and 23 stations from other networks are included. The figure shows the stations of the Swiss territory only.
2214
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
Fig. 2. 3D view of the evaluation perimeter. The tomographic voxel model consists of 16 layers up to 15,000 m height. The borders of the layers are set at 0, 200, 600, 900, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2700, 3200, 4000, 5000, 8000, and 15,000 m (the figure shows layers up to 5000 m height only). Each layer contains six voxels in longitude and three voxels in latitude (spacing 0.5) plus 22 outer voxels (not shown on the figure). The GPS stations are shown as column according to their station height. Radiosondes are launched from the radiosonde station Payerne (MeteoSwiss), shown as cuboid.
The voxel model above the Swiss territory is chosen with 6 · 3 voxels in horizontal (voxel side lengths 50 km) and 16 layers up to 15,000 m (Fig. 2). Water vapor profiles are retrieved on an hourly basis. The tomographic results are compared with data of the operational high resolution NWP model of MeteoSwiss (aLMo). aLMo is the Swiss implementation of the non-hydrostatic local weather model of COSMO (Doms and Scha¨ttler, 2003). The model is operational at MeteoSwiss since April 2001. It covers most of Western Europe, has a mesh size of about 7 km in the horizontal and uses 45 vertical layers, with a thickness of about 100 m in the lower troposphere. According to the space and time resolution of the tomographic evaluation, a total of 7680 aLMo profiles and radio soundings at the station Payerne (twice a day) are made available by MeteoSwiss for comparison purposes.
is assigned to the uppermost layer (8000–15,000 m) with a regularization factor of 1/900. On the other hand, AWATOS ANETZ contains the same constraints as AWATOS Correlation and in addition one a priori wet refractivity value for each voxel lower than 2000 m.
4. Profile determination and evaluation Double difference residuals and hourly means of zenith path delays are determined using the Bernese software package (Beutler et al., 2001). Reconstructed double-difference wet slant delays are then introduced into the tomographic software package AWATOS. Two types of tomographic profiles are determined, containing different constraints. On one hand, AWATOS Correlation includes inter-voxel constraints between all neighboring voxels. Compared to a double-difference observation, the constraints are downweighted by a factor of 502 (regularization factor = 1/2500). In addition, an a priori wet refractivity of zero
Fig. 3. Cross section of the square root of cofactors at latitude 46.75. The two top layers (5000–15,000 m height) are not shown in the figures. Panel (a) represents the situation of the solution AWATOS Correlation and (b) AWATOS ANETZ. It is seen, that if a priori refractivity is introduced, the precision increases significantly.
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
These refractivity values are based on a collocation and interpolation procedure (COMEDIE software package; Troller et al., 2000), using the operationally available meteorological ground station data (ANETZ, MeteoSwiss). The least-squares adjustment of the tomographic equation system yields the matrices of the cofactors of the unknown parameters from the inverted normal matrices (e.g., Mikhail, 1976). The comparison of the two cofactor matrices of the two solutions reveals the impact on quality of the two different types of solution. The a priori calculation takes no measurements into account. It assesses the geometry and the weighting of the measurements, only. The solution AWATOS Correlation shows a slight increase of the precision with height. Regarding the solution AWATOS ANETZ, a significant improvement of the square root of cofactors is visible in the first 2000 m height but also in the higher voxels (Fig. 3). In Fig. 4, two examples of tomographic profiles, aLMo data and radiosonde profiles are displayed. Furthermore, the profile obtained with COMEDIE is plotted. The lowermost 2000 m of the latter profile are used to constrain the solution AWATOS ANETZ. The
2215
comparison of the profiles confirms the conclusions made of the matrices of cofactors: The degree of agreement of the individual profiles is varying. AWATOS ANETZ fits always more accurately to aLMo than AWATOS Correlation. The 22 radiosonde launches during the investigation week at station Payerne have been used to perform a comparison to the tomographic profiles and the aLMo model (Table 1). The conclusion of Fig. 4 can be verified, the mean rms of AWATOS ANETZ is half as large as of AWATOS Correlation. Fig. 5 shows the mean rms of the tomographic profiles compard to the radiosondes as function of the height. The impact of the a priori refractivity in the first 2000 m height is clearly visible. The comparison of aLMo to the radiosonde profiles is also documented in Table 1. The agreement is within 2.6 ppm (refractivity units) at station Payerne, where radiosondes are available to compare the aLMo model. As radiosondes are available every six or twelf hours and at the station Payerne only, the aLMo data with an hourly resolution have been used to compare the tomographic solutions over entire Switzerland. Thus, a total of 3024 profiles are compared to aLMo for the one-week measurements (Fig. 6). In Fig. 4, two profiles with
5000
5000
GPS Tomography
GPS Tomography
06.11.02 17h
4500
AWATOS Correlation AWATOS ANETZ aLMo solution COMEDIE profile Radiosonde
4000
3500
3500
height [m]
height [m]
3000
2500
2500
2000
2000
1500
1500
1000
1000
500
500
a 0
AWATOS Correlation AWATOS ANETZ aLMo solution COMEDIE profile Radiosonde
4000
3000
0
06.11.02 23h
4500
b
© M. Troller, GGL-ETHZ
10
20
30
40
wet refractivity [ppm]
50
60
0
0
© M. Troller, GGL-ETHZ
10
20
30
40
50
60
wet refractivity [ppm]
Fig. 4. Wet refractivity profiles of the two tomographic solutions, the radiosonde ascending at the station Payerne (MeteoSwiss), data of the numerical weather model aLMo (MeteoSwiss) and the COMEDIE profile. Panel (a) shows a profile with a steep refractivity gradient at 1800 m height. All solutions match accurately together on the order of 10 ppm. Panel (b) shows a situation with an inhomogeneous decrease of refractivity with height. AWATOS Correlation underestimates the wet refractivity in the lowest 1000 m height. However, also the other tomographic solution as well as the radiosonde show differences to the aLMo solution.
2216
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217
Table 1 Statistical analysis of the tomographic profiles and aLMo compared to the radiosondes
Mean offset Mean rms Mean r
AWATOS Correlation
AWATOS ANETZ
aLMoa
6.3 12.7 11.0
1.4 5.1 4.9
0.5 2.6 2.6
The analysis contains 22 profiles. The accuracy of the two tomographic solutions varies significantly. The mean rms of AWATOS Correlation is approximately double as the corresponding value of AWATOS ANETZ. Furthermore, the latter has no significant mean offset. a aLMo data for the lowermost voxel (200–600 m height) are not available.
Fig. 5. Mean rms of the two tomographic solutions compared to the radiosondes as a function of the height. The plot shows the comparison of the 22 radiosonde launches during the investigation week at the station Payerne. Generally, the rms is decreasing with increasing height. The impact of the a priori profiles on the first 2000 m height from COMEDIE is clearly visible in the solution AWATOS ANETZ.
18
AWATOS Correlation AWATOS ANETZ
16
mean rms [ppm]
14 12 10 8 6 4
9, coincide with atmospheric situations rapidly changing from heavy rain to dry periods and contrariwise. During these time periods, the differences between GPS-tomography and aLMo are larger. This may mainly be caused by the temporal averaging during each data assimilation period, thus neglecting gradients in time for the assimilation interval.
5. Conclusions and outlook The tomographic approach was successfully performed to process GPS data for the determination of 4 dimensional water vapor data. The investigations showed that the Swiss national reference network AGNES is a dedicated network for GPS tomography if using a grid model with 6 · 3 voxels horizontally. AGNES has a station density of about one GPS receiver per voxel column. Regularization methods are necessary to achieve a stable tomographic solution. Inter-voxel constraints and a priori refractivity information allow usually to achieve a high accuracy. An overall rms of about 5–8 ppm (refractivity units) or 0.8–1.3 g/m3 absolute humidity compared to aLMo has been reached, depending on the tomographic solution. Using a temporal resolution of 1 h, the accuracy remains approximately stable during the whole week, at least for the solution AWATOS ANETZ. Only during rapidly changing atmospheric situations, the rms is slightly decreasing. Further investigations showed (Troller, 2004) that the accuracy achieved at station Payerne, can be expected in whole Switzerland. Nevertheless, in a further development, an improvement of the regularization methods should be aimed at. The use of humidity profiles over land determined by satellites as well as a priori data from NWP should be studied. Further investigations should focus on the augmentation of the tomographic spatial and temporal resolution. The introduction of data from GPS-independent sensors such as radiometers or solar spectrometers should be considered. The impact of the derived information on the quality of NWP forecast should then be thoroughly evaluated.
2 0 3
4
5
6
7
8
9
10
Day of November 2002
Fig. 6. Time series of the mean of all profiles compared to aLMo. Overall, AWATOS ANETZ has a smaller rms than AWATOS Correlation. Three jumps are visible on November 3, 6 and 9.
different levels of agreement, mainly depending on the atmosphereÕs actual state are shown. This behavior can be confirmed by inspecting the time series analysis of Fig. 6. Periods with large rms on November 3, 6 and
References Beutler, G., Bock, H., Brockmann, E., Dach, R., Fridez, P., Gurtner, W., Hugentobler, U., Ineichen, D., Johnson, J., Meindl, M., Mervart, L., Rothacher, M., Schaer, S., Springer, T., Weber, R. in: Hugentobler, U., Schaer, S., Fridez, P. (Eds.), Bernese GPS Software Version 4.2. Astronomical Institute, University of Berne, 2001. Bevis, M., Businger, S., Herring, T., Rocken, C., Anthes, R., Ware, R. GPS meteorology: remote sensing of atmospheric water vapor using the global positioning system. J. Geophys. Res. 97 (14), 15787–15801, 1992.
M. Troller et al. / Advances in Space Research 37 (2006) 2211–2217 Doms, G., Scha¨ttler, U. (Eds.). Newsletter No. 3, February 2003. Technical Report, COSMO (Consortium for Small-scale Modelling), 2003. Ducrocq, V., Ricard, D., Lafore, J.-P., Orain, F. Storm-scale numerical rainfall prediction for five precipitating events over France: on the importance of the initial humidity field. Weather Forecast. 17, 1236–1256, 2002. Elgered, G. Atmospheric Remote Sensing by Microwave Radiometry, Chapter Tropospheric Radio-Path Delay from Ground-based Microwave Radiometry, Wiley, New York, pp. 215–258, 1993. Flores, A., Ruffini, G., Rius, A. 4d tropospheric tomography using GPS slant wet delays. Ann. Geophys. 18, 223–234, 2000. Hirahara, K. Local GPS tropospheric tomography. Earth Planets Space 52, 935–939, 2000. Hofmann-Wellenhof, B., Lichtenegger, H., Collins, J. GPS Theory and Practice. Springer, New York, 2001. Kruse, L. Spatial and temporal distribution of atmospheric water vapor using space geodetic techniques. Geoda¨tisch-geophysikalische Arbeiten in der Schweiz, vol, 61, Schweizerische Geoda¨tische Kommision, 2001. Mendes, V.B. Modeling the neutral-atmosphere propagation delay in radiometric space techniques. Ph.D. dissertation, University of New Brunswick, 1999.
2217
Mikhail, E.M. Observations and Least Squares. University Press of America, 1976. Park, S.K. Nonlinearity and predictability of convective rainfall associated with water vapour perturbations in a numerically simulated storm. J. Geophys. Res. (Atmospheres) 104, 31,575– 31,587, 1999. Saastamoinen, J. Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites. In: Henriksen, S.W., Mancini, A., Chovitz, B.H. (Eds.). The Use of Artificial Satellites for geodesy, Geophys. Monogr. Ser. AGU, 15, pp. 247–251, 1972. SMA. Charakteristiken der ANETZ-Daten. In: Beitra¨ge zum ANETZDaten Kolloquium vom 17. April 1985 in Zu¨rich, p. 56. Schweizerische Meteorologische Anstalt Zu¨rich, 1985. Troller, M. GPS based determination of the integrated and spatially distributed water vapor in the troposphere. Diss ETH No. 15513, 2004. Troller, M., Bu¨rki, B., Cocard, M., Geiger, A., Kahle, H.-G. 3d refractivity field from GPS double difference tomography. Geophys. Res. Lett. 29 (24), 2149–2152, 2002. Troller, M., Cocard, M., Geiger, A. Modellierung 4 dimensionaler Refraktionsfelder zur Berechnung von Wegla¨ngen-Korrekturen bei Satellitenmessungen. In: Simulation raumbezogener Prozesse: Methoden und Anwendungen, 9, IfGI prints, 2000.