Data Systems

Water vapor plays an important role in meteorological applications; GeoForschungsZentrum (GFZ) therefore developed a tomographic system to derive 3-D distributions of the tropospheric water vapor above Germany using GPS data from about 300 ground stations. Input data for the tomographic reconstructions are generated by the Earth Parameter and Orbit determination System (EPOS) software of the GFZ, which provides zenith total delay (ZTD), integrated water vapor (IWV) and slant total delay (STD) data operationally with a temporal resolution of 2.5 min (STD) and 15 min (ZTD, IWV). The water vapor distribution in the atmosphere is derived by tomographic reconstruction techniques. The quality of the solution is dependent on many factors such as the spatial coverage of the atmosphere with slant paths, the spatial distribution of their intersections and the accuracy of the input observations. Independent observations are required to validate the tomographic reconstructions and to get precise information on the accuracy of the derived 3-D water vapor fields. To determine the quality of the GPS tomography, more than 8000 vertical water vapor profiles at 13 German radiosonde stations were used for the comparison. The radiosondes were launched twice a day (at 00:00 UTC and 12:00 UTC) in 2007. In this paper, parameters of the entire profiles such as the wet refractivity, and the zenith wet delay have been compared. Before th validation the temporal and spatial distribution of the slant paths, serving as a basis for tomographic reconstruction, as well as their angular distribution were studied. The mean wet refractivity differences between tomography and radiosonde data for all points vary from −1.3 to 0.3, and the root mean square is within the range of 6.5–9. About 32 % of 6803 profiles match well, 23 % match badly and 45 % are difficult to classify as they match only in parts.


Introduction
Water vapor plays a key role in various atmospheric processes; e.g., its long-term variations are related to climatological temperature trends.To precisely characterize the influence of water vapor on these processes, the knowledge of its temporal and spatial distribution is a key prerequisite.The poor resolution of the atmospheric water vapor observations in space and time limits the accuracy of short-term weather forecasts.In the recent years, ground-based GPS sounding systems were in operation and the ability to provide such spatiotemporally resolved information on the water vapor distribution was demonstrated.The provision of the integrated water vapor (IWV) and atmospheric zenith total delays (ZTDs) was established as a standard sounding technique.ZTD and IWV are currently operationally assimilated into regional weather forecast models.The derivation of 3-D water vapor distributions is currently under development.Since the first feasibility of water vapor tomography with GPS data in Bevis et al. (1992), many GPS-based water vapor retrieval systems have been set up (Flores et al., 2000;Troller et al., 2002;Champollion et al., 2009;Perler et al., 2011).
Radio signals emitted by the GPS satellites are bent and delayed on the way to the ground depending on atmospheric temperature, pressure and humidity.Based on the tropospheric delays from a network of GPS receivers, the Published by Copernicus Publications on behalf of the European Geosciences Union.
water vapor distribution within the troposphere can be determined by tomography techniques.The quality of the GPS tomography is affected by a number of factors such as the spatiotemporal distribution of the observations, the reconstruction method, the initial field, etc., and independent observations are required to validate its quality.Radiosonde data are often used to evaluate the capabilities of GPS tomography (Gradinarsky et al., 2004;Troller et al., 2006;Perler et al., 2011).However, only several radiosonde profiles from one or two stations were used for comparison in most of these studies.In this study, more than 8000 radiosonde profiles for the year 2007 provided by 13 German radiosonde stations at 00:00 UTC and 12:00 UTC are available.A set of criteria for selecting reliable profiles was defined before validation.All the chosen profiles are divided into well or poorly matched profiles with the help of defined parameters.This classification into three groups of profiles describes the quantitative and qualitative agreement between the tomographically reconstructed and radiosonde profiles.

Data sources
The current investigation involved the tropospheric products of the German ground GPS network and the radiosonde data from the German Weather Service (DWD) for the year 2007.In Fig. 1 the red dots show the 272 GPS stations used for tomography.The horizontal resolution of the generated tomography is about 50 km, which cannot be much better than the mean GPS interstation distance of 40 km.The white rhombs are the radiosonde stations in Germany, which provide two profiles per day at 00:00 UTC and 12:00 UTC for the year 2007.The horizontal distance between radiosonde stations is about 200 km on average, and the observations pressure, temperature and relative humidity are measured by radiosonde sensors.

GPS data processing
The IGS processing center located at GFZ operationally provides several geodetic products, e.g., orbits and clock corrections of all GPS satellites, which are estimated every three hours (Gendt et al., 2011).Based on these data a processing system for ground-based GPS atmosphere sounding is running operationally and estimates the ZTD, IWV and STDs on an hourly basis.STD data are available with a sampling rate of 2.5 min providing 6-12 GPS slant delays per station and epoch.Currently about 300 stations from Germany and neighboring countries are processed resulting in hourly data sets of about 50 000 STDs.
The GPS signal is bent and delayed by the atmosphere in comparison to undisturbed propagation in a vacuum.The GPS data processing needs to estimate this delay in order to obtain precise positions.The GPS atmosphere processing provides the most reliable atmospheric delays by utilizing the available geodetic information.Two mostly used processing strategies are network solutions providing double differences of carrier phase measurements and the precise point positioning (PPP) technique (Zumberge et al., 1997), which estimates the slant delays independently.At GFZ the PPP uses precise orbit and clock estimates, which are available from the IGS service center and can be further refined for the atmosphere processing, e.g., from 3 to 1 h.Each station is processed independently by PPP, and the processing of a large number of stations can easily be parallelized in order to obtain nearreal-time solutions.
STDs are provided by the GPS data processing software EPOS (Gendt et al., 2004) developed at GFZ using the PPP method.The STD can be expressed as a combination of different estimates (MacMillan, 1995;Bender et al., 2010): where ZHD and ZWD are the hydrostatic and wet zenith delay, m h and m w are the hydrostatic and wet mapping function (Niell, 1996), G N and G E are the delay gradient parameters in north and east directions, is the elevation, φ is the azimuth and δ is the postfit phase residual.
Real-time processing is important for several meteorological applications yet has the disadvantage that only observations arriving in time can be considered.Missing data lead to gaps in the time series, and stations without a fast data connection cannot be processed.For validation studies it would be more beneficial to use a reliable reprocessed data set with minimum gaps and a maximum number of observations.The currently largest and most carefully reprocessed data set available from the GFZ is the COPS reanalysis for 2007.During the Convective and Orographicallyinduced Precipitation Study (COPS, June to August 2007) and the Global Observation Period (GOP, January to December 2007), an extensive meteorological data set was taken and additional GPS stations were set up in eastern France and southern Germany (Wulfmeyer et al., 2008;Crewell et al., 2008).Not all stations had online data connection and it was necessary to collect all available GPS raw data and the necessary meteorological observations after the campaign.The reprocessing at GFZ was carried out in several steps.The quality of the available data and the remaining gaps were analyzed carefully after each step and attempts were made to obtain the missing data.

GPS tomography
The GPS tomography combines a large number of STD observations in order to reconstruct a spatially resolved refractivity field (Flores et al., 2000).In many cases it is more desirable to obtain the humidity distribution which is closely related to the wet refractivity and the slant wet delays (SWDs).The STD can be subdivided into a hydrostatic and a wet part, at which the latter is related to the water vapor.The SWD can be separated from the STD by estimating the slant hydrostatic delay (SHD) using the zenith hydrostatic delay (ZHD) (Saastamoinen, 1973;Bevis et al., 1992): where p 0 is the surface pressure in hPa, φ is the geographic latitude and h 0 is the station height above geoid in km.The ZHD can be mapped onto the slant path using the hydrostatic mapping function (Niell, 1996) and the SWD is given by the difference between the STD and the SHD: The SWD is related to the wet refractivity N w by the integral along the signal path S: This integral is discretized using a regular WGS84 grid and linearized by assuming that the signal path does not depend on N w within each cell: Here, s i is the slant subpath within the i-th grid cell.Equation (6) leads to a system of linear equations if all SWD observations available are combined to a vector m, the refractivity N w field is mapped on a vector n and the design matrix A is given by the slant subpaths in each voxel.
The wet refractivity can be obtained from Eq. ( 7) by means of an inverse reconstruction process.Equation ( 7) defines an illposed inverse problem which can be solved iteratively, e.g., by using the multiplicative algebraic reconstruction technique (MART) (Subbarao et al., 1997): where j denotes the grid cell, i the observation, k the iteration step and λ = 0.2 is the relaxation parameter.A i , n k is the dot product of the row vector A i with the the state vector n after the k-th iteration.The MART method starts with an initial field which is iteratively improved by applying small corrections to each grid cell.The initial field was obtained by using all the available synoptic data from the DWD and GPS meteo stations in the selected region.The stations observations were interpolated on the grid by means of a leastsquares collocation.It is a numerical approximation method with parameter estimation and statistical interpolation of observations in space and time (Hirter, 1996).Besides additional observtions, e.g.IWVs, and inter-voxel constraints were used to stabilize the tomographic reconstruction.The quality of the reconstructed N w field depends on the quality of the initial field, the number and quality of the input slant data, the resolution of the chosen grid and many parameters of the reconstruction algorithm.
The tomography system developed at GFZ (Bender et al., 2010) works with the iterative reconstruction technique MART with the GPS slant path delays from the German network as input data (Fig. 1).A good spatial coverage of the atmosphere by the slant paths is required for the tomography.Here, water vapor fields with a spatial resolution of about 50 km horizontally, 500 m vertically up to 10 km and a temporal resolution of 30 min were reconstructed at 00:00 UTC and 12:00 UTC from 1 January 2007 to 31 December 2007.The input slant data and synoptic data were in the period ±15 min of the radiosonde launch time.The reconstruction region is located from longitude 4.9 • to 15 • , and latitude 47.3 • to 54.8 • with the grid described in Fig. 1, and the height is from 0 to 10 000 m with 20 intervals.The interpolation of the synoptic data was used for the initial field with the method described above.With the above setting and applying the MART method, GPS water vapor tomography for 2007 has been reconstructed with 100 iterations.

Radiosonde profiles
A set of radiosonde (RS) profiles from 13 German RS stations (see Fig. 1) was available for the year 2007.The data set with the profiles of the 00:00 UTC and 12:00 UTC ascents contained information about the altitude, pressure, temperature, humidity, wind speed, etc. From these data the wet refractivity profiles can be estimated using the Smith and Weintraub (1953) or Thayer (1974) formula: Here, e is the partial pressure of the water vapor, T is the temperature and Z −1 w is the inverse compressibility factor of water vapor.The coefficients k 2 and k 3 can be found in Bevis et al. (1994).
However, the vertical distribution is rather variable and many profiles show large gaps.Most of RS profiles reach up to 20-35 km, but for the validation only RS data within the tomography grid up to 10 km are relevant.Most profiles consist of 30-50 observations by 10 km, which correspond to a mean vertical distance of the observations between 200 and 300 m.
For the sake of validating the quality of the tomographically reconstructed humidity fields the RS profiles were used as a reference.However, it should be considered that radiosondes also have their difficulties (Vömel et al., 2007;Miloshevich et al., 2009).The average accuracy of the widely used radiosonde Vaisala RS92 is about ±4 % of the measured relative humidity (RH) value for nighttime soundings and ±5 % for daytime sounding (Miloshevich et al., 2009).Compared to nighttime soundings there is a solar radiation error in daytime measurements caused by solar heating of the RH sensor.A time-lag error can be caused by slow sensor response to changing RH conditions at low temperatures (Miloshevich et al., 2004).Comparisons between RS humidity observations and GPS ZWDs can be found in Schneider (2010).

Validation of N w fields with radiosonde profiles
The GPS tomography provides N w fields on a rather coarse spatial grid and a reconstruction quality, which depends on the geometry of GPS satellite constellation and many other parameters.These fields have to be compared with a limited number of RS profiles with a high vertical resolution.The quality of the tomographically reconstructed profiles is limited firstly by the low spatial resolution of the tomography.The generated tomographic reconstruction has a horizontal resolution of about 50 km and vertical resolution of 500 m.Moreover the wet refractivity is assumed to be constant within a voxel.Therefore the interpolation is applied for the comparison with radiosonde data.Horizontally bilinear interpolation was used with four adjoining grid points whereas vertically cubic spline interpolation was performed (see Fig. 2) (de Boor, 2001).
Secondly, the quality of the tomographically reconstructed profiles is limited by the errors of observations.It is not only from estimated STDs in the EPOS but also from estimated SHD using Saastamoinen model.Thirdly, the solution is not unique and not stable due to the ill-posed inverse problem.GPS tomography usually generates a partly ill-posed problem with incomplete input data, in which a portion of the unknowns can be determined whereas the other part is underdetermined.
It is obvious that not all grid columns shown in Fig. 1 will have the same reconstruction quality.The spatial and temporal distribution of observations within the grid is highly variable, and it must be determined which parts might have been reconstructed with sufficient quality before the validation is applied.Otherwise, the validation will provide distorted results, e.g., by comparing regions which were not covered by slant observations.This is also an important question for further applications of the reconstructed humidity fields, where some sort of confidence level needs to be supplied together with the fields.Defining criteria for "reliable" columns which can later be applied without having reference data is therefore an important task.Another problem is the method to compare a gridded field with the profiles.There are several possibilities, e.g., comparing the data point by point or regarding the whole profiles.A statistically sound strategy has to be developed.

Estimating the quality of reconstructed profiles
About 50 000 slant delays per hour are used to reconstruct the humidity fields, but the distribution of slant paths between the receivers and satellites is always inhomogeneous in space and time, depending on the specific GPS satellite constellation.Consequently, the information available from the observations changes dynamically.Figure 3 shows an example of the spatial coverage of the atmosphere with slants paths.It shows that the slant paths are very dense in some areas while some other areas show almost no slant paths.During COPS a dense GPS receiver network was deployed in the Rhine Valley.Most stations were providing data from June 2007 until the end of the year.Therefore, the slants in the Black Forest region in southwestern Germany are very dense during this period.On the other hand there are regions in northern Germany and above the North Sea with very few observations.A reliable tomographic reconstruction requires a large number of intersecting slant paths from a wide angular range (Kak and Slaney, 1999;Champollion et al., 2005) within each grid cell.It can therefore be assumed that cells with no observations and cells that are pierced from the parallel paths will not be reconstructed very well.As more than 50 % of the cells remain empty and other cells with several GPS stations are pierced by a large number of parallel paths (Bender et al., 2009), these criteria can be used to identify regions with unreliable data within the grid.

Density of slant paths
The first step to estimate the quality of the tomographic reconstruction is to count the number of slants per grid cell and per column.In an ideal case only columns with a large number of evenly distributed slants would be selected, but due to the large number of "empty" cells, even the columns with one or more empty cells cannot be totally rejected.For comparison with the RS profiles, four neighboring columns are required to interpolate the gridded data on the profile (see Fig. 2).
Figure 4 shows a typical situation.The number of slants per voxel within four adjoining columns and at different altitudes can be significantly different (Fig. 4, right).Here, one column was not pierced by any slant path, and two columns show a rather small but evenly distributed number of slants which increases slightly with increasing height.This is typical for columns which do not contain any GPS stations but are surrounded by several stations.The last column shows a very high number of slants per cell and contains several stations.The bilinear interpolation of the N w data from these columns gives the highest weights to the least distant points without considering the number of observations.In Fig. 4 the column second closest to the RS profile gets a higher weight than the column with most observations; i.e., the interpolated value has a large contribution from the initial field.The reconstructed profile (Fig. 4, left) shows segments with a strong impact of the inversion, i.e., the STD observations, but also segments which are very close to the initial profile.A combination of such contradictory information often leads to strong artifacts in the reconstructed field.Such profiles were rejected for the validation.It is rather difficult to define general criteria for rejecting poor profiles as too-restrictive criteria remove too many profiles, with some of them being relatively good.By varying several parameters it was decided to reject all sets of four neighboring profiles where the maximum number of slants per cell in one of the four columns is smaller than two.In this case there are no intersecting slant paths and very often lots of cells without any data.This rather weak requirement already rejects ∼ 7 % of the profiles (see Table 1).

Angular distribution of slant paths
The quality of a tomographic reconstruction depends not only on the number of slant paths but also on the way these paths are intersecting.Therefore, the crossing angles of any pair of slant paths within each voxel were investigated in the next step.In an ideal case an almost flat distribution of intersection angles ranging from 0 • to 180 • would be expected.To investigate a real situation the intersection angles and their distribution must be computed for each single grid cell.This was done for an example shown in Fig. 5.The two profiles shown were chosen from a region with a rather high density of slant paths, and all grid cells contain many observations.The RS profiles (dashed line) and the tomographically reconstructed profiles (solid line) are shown together with the N w data at the surrounding grid points.The humidity profile at Bergen, 31 January 2007, 12:00 UTC (Fig. 5, left), shows rather large deviations from the initial profile but could be reconstructed quite well.Contrarily, the humidity profile at Schleswig, 9 June 2007, 12:00 UTC, which was reconstructed from a comparably large number of STD observations, does not match the RS profile.Neither the fluctuations in the boundary layer nor the step at 3000 m can be seen in the tomography profile.The number of STDs around the RS station was sufficient and it could be assumed that the angular distribution of the slant paths is the reason for the poor quality of the reconstruction.This aspect was investigated in more detail.The N w data indicated by the black box, i.e., from the column closest to the interpolation value, are chosen as an example.The histograms of crossing angles of the two selected voxels at an altitude of 2000 m are shown in Fig. 6 and 2500 m in Fig. 7.
The crossing angle φ of every two slants r 1 and r 2 is as follows: In Fig. 6 (left) φ of Bergen are widely distributed and indicate intersecting slant paths from a wide angular range.Contrarily, the histogram for Schleswig (Fig. 7, left) shows a narrow distribution with slant paths pointing in two or three distinct directions.The number of slants in Fig. 5 (right) is 80, which is a little less than 87 in Fig. 5 (left) but should be sufficient for a reliable reconstruction.However, in the Schleswig case it is difficult to locate the information provided by the slant delays as all signals propagate almost in the same direction and there are no real intersections between the paths.This leads to the rather poor reconstruction as compared to Bergen.
Analyzing distributions as shown in Fig. 6 (left) and Fig. 7 (left) for each single grid cell would be rather difficult and cannot provide definite information on the reconstruction quality as the slant path distribution in the entire region has to be considered.For the sake of simplicity, a classification number is defined which is much easier to compute and to analyze: for each slant path the angles to the principal axes of the local Cartesian system are computed; i.e., φ z is the angle to the vertical axis (90 • < φ z < 180 • ), and φ x and φ y are the angles to the horizontal axes (0 • < φ x,y < 180 • ).The formula for the cutting angle between slant paths and the local axis vector is described by Eq. ( 11) (see Fig. 2): Here r is the slant vector and e i (i = x, y, z) are the Cartesian unit vectors.Figure 7 (right) shows the distribution of φ i for Bergen, Fig. 6 (right, for Schleswig).The number of entries in each histogram is much smaller than in Fig. 6 (left) as only one angle is computed for each slant path and each principal axis, while in Fig. 6 (left) any possible pair of slant paths is regarded.Again, a wide distribution of angles is found for the station Bergen (Fig. 6), while the angles for the station Schleswig are centered at some specific directions (Fig. 7).To quantify the width of the distribution, the differences between the maximum and minimum cutting angles are analyzed for each voxel: All the direction vectors of slant paths only were downward (90 • < φ z < 180 • , dφ z < 90 • ).It turned out that dφ z was not directly correlated with the reconstruction quality, while large numbers of dφ x and dφ y usually indicate a wide slant path distribution and a good reconstruction quality.To identify unreliable profiles the dφ x and dφ y values of all cells of the four grid columns surrounding the RS station were analyzed.The profile was rejected if the maximum of dφ x and dφ y were both smaller than 90 • in all of the four columns.Another 692 profiles were rejected by this criterion, regarding the profiles passing the criteria defined in the previous section (see Table 1).

Radiosonde data
A total of 8109 radiosonde profiles for the year 2007 were available for the study.However, not all of them could be used for the validation.Some of them show large vertical gaps or are truncated at a certain height.A total of 27 RS profiles with no data above 4000 m were rejected.In addition the wet refractivity at some start points of the radiosonde profiles is much larger than N w at the following points.According to Miloshevich et al. (2009), measurement errors tend to occur in the first 100 m or when the temperature gradient changes abruptly, due to the fact that the RH sensor possesses different thermal time constants at which the measured air temperature does not accurately represent the temperature from the RH sensor.Hence if the vertical height difference between the start point and the second point is smaller than 20 m and the wet refractivity difference is larger than 50, then the measured data of the start point is deleted.Altogether 22 start points were deleted for the whole year.Some radiosonde profiles have two values at the same height and sometimes the height even decreases at some points.For interpolation and integration the sequence of points with decreasing height was changed and 1cm was added at the points with the same height.
www  gaps.For 18 profiles the STD data were not sufficient to start the tomography.A total of 569 tomographically reconstructed profiles were excluded from the validation because the number of slant data in the vicinity of the RS station was insufficient, and 692 profiles were excluded because the angular coverage was poor.
Most of the profiles rejected from the tomography because of a insufficient number of slants belong to the station Schleswig.As can be seen in Fig. 1, there are only few GPS stations around Schleswig and they are located only on the western side.Certainly not all the profiles selected in this way provide reliable data, but due to the necessity to interpolate from the rather coarse tomography grid to the RS station four surrounding columns are required.It must be pointed out that the columns of the reconstructed field not matching these criteria do not necessarily contain poor data.Columns with few observations are dominated by the initial field, which was only modified due to the inter-voxel constraints.If the initial field was chosen well, then the "reconstruction" will match the real conditions, i.e., show a smooth exponential profile.However, this is not the result of the tomographic reconstruction and these data will not be considered further.

Comparison of GPS tomography with radiosonde data
After cleaning the available RS profiles and interpolating the reliable parts of the tomographically reconstructed N w field on the individual points of the RS profiles, two sets of N w profiles defined on the same points in space are available and can be compared.Many methods can be used for the comparison such as a point-by-point validation or the validation of entire profiles.

Point-by-point validation
In a first step the data are compared point by point; i.e., the mean differences and their rms are compared for each station without regarding the height or observation time of the data with the number of data a.The absolute differences might be misleading for the almost-exponential vertical N w profiles, and the corresponding relative differences were also taken into account: Table 2. Statistics of the differences between radiosonde and GPS tomography; a is the number of observations.N w and rms are the average difference of N w between radiosonde and tomography and its root mean square.N r w and rms r are the average and root mean square of the relative differences.

RS
The results for one year of data are shown in Table 2.For each station the number of observations a is given.There were many points at different refractivities between 0 and 100 and different value scale from 0 to 100, and as a result the corresponding root mean squares were quite variable.
Table 2 shows the overview of the result.The N w vary from −1.3 to 0.3 and an overall rms of about 6.5-9 has been reached, depending on different radiosonde stations (see Table 2.) Most of the N r w are positive because of many small radiosonde observations at high altitude (N RS w,i ) as denominator.At high altitude the measured values of RS were always very small, and therefore the corresponding N r w in this region have always very large values in spite of small N w .These results are rather unspecific and provide little information about the quality of the tomographic reconstruction.It cannot be seen how many profiles show a good representation of the real vertical atmospheric structure.
More meaningful than the mean of all observations is the mean deviation at a certain altitude.For this purpose the differences of RS and tomography data were classified in 20 vertical layers defined for the tomographic reconstruction, and the means as well as their root mean square were computed using Eqs.( 13) and ( 14) for each individual layer.N w is dominated by deviations in the boundary layer, where N w is large and leads to small relative deviations.N r w is dominated by deviations near the tropopause, where N w is small and small deviations can lead to very large relative errors up to several 100 %.The N w observations between 0 and 5 and the large relative uncertainty of humidity data at higher altitudes do not necessarily correspond to poorly matching profiles.Data from three RS stations (Lindenberg, Essen and Meiningen) were chosen as examples.The difference of the wet refractivity distributed respectively in large, middle and small rms (Table 2).The result for 20 equidistant layers are shown in Fig. 8.It can be seen that plots for the different stations are following the same pattern.The bias of different RS stations shows rather large fluctuation in the lower layers.The rms is rather similar for all stations and decreases with increasing altitude.The relative bias of different RS stations shows large fluctuations in the upper layers.rms r increases with increasing altitude for all stations.Below 2000 m the mean difference is very variable because of an insufficient number of GPS observations at low altitudes.Above 8000 m the wet refractivity of radiosonde data is rather small and has a large effect on the relative value.The presented discrepancies can be affected by the quality of GPS tomography and radiosonde data.

Validation of entire profiles
In a second step attempts were made to compare entire profiles and to quantify the degree of consistency between RS and tomography profiles.The visual validation of 6803 in-dividual profiles one by one requires a lot of work.Therefore, defining general criteria, which can be checked automatically, is necessary to determine the agreement between radiosonde and tomography profiles in a consistent way.To check if two profiles are almost similar it is beneficial to compare some integrated quantities, e.g., the Zenith Wet Delay (ZWD) of the profile up to the top grid layer.
The ZWD should be almost the same for both profiles; i.e., the difference D should be small.the area K between the two N w (h) functions should also be considered: The refractivity changes considerably with the weather situation, and rather small differences on a warm humid summer day can lead to K data much larger than it was found for poorly matching profiles at a dry and cold winter night.Therefore, the normalized area k is used for comparing the quality in a more general way.
Most profiles look very similar if d and k are sufficiently small; however, there are some profiles with outliers which have a small impact on the integrated quantities d and k but still need to be rejected.Such outliers can be due to artifacts in the tomographic reconstruction but are also present in the RS data.The maximum differences m is used to identify such outliers: By visual inspection certain quantities for m, d and k were set.They characterize well-matching profiles and poorly matching profiles (Table 3).Profiles belonging to neither group were regarded as indifferent.It turned out that even the relative quantities m, d and k depend to some degree on the atmospheric humidity, and different thresholds were defined for different weather situations.More relaxed criteria can be used in case of very humid situations, while rather rigid settings are required in dry cases.The ZWD is used to quantify the total amount of atmospheric humidity, and four classes  with different thresholds for d, k, and m as given in Table 3 are defined to identify the quality of the tomographically reconstructed profiles.For very dry situations (ZWD < 60 mm) it is important to reject outliers (small m), but rather large values of k can be tolerated.With an increasing amount of water vapor the outliers become less important, but the normalized area k between the N w (h) curves must be rather small.A small value of k implies that d is also small and no specific threshold for d needs to be given for the well-matching profiles.However for rejecting poorly matching profiles d should be specified.
To determine the thresholds of different parameters we use a coarse value above or below mean value for the poor and good criteria at first.Then a sensitive adjustment with smaller and larger value was made to test the results.Consequently the thresholds for the criteria were empirically determined and could be changed for different situations.
In Fig. 10 four characteristic profiles are given together with the corresponding quantities of the ZWD, m, k and d.Example (a) shows two profiles with almost the same ZWD (small d) and outliers below the threshold for poor profiles but with a large area between the N w (h) functions (k = 59.4 % > 40 %).The situation in Fig. 10, (b) shows one fundamental problem of the GPS tomography: the reconstructed profile is almost identical to the RS profile for altitudes above 1500 m but absolutely wrong at lower altitudes because of the sparse GPS observations in the boundary layer.The profile is regarded as poorly matching because of the large m = 25.93.In case (c) the total amount of humidity is wrong and the profile was rejected because d = 78.14% is too large.Figure 10, (d) shows a well-reconstructed profile.

Results of the validation
The results of the validation study are summarized in Table 4. Regarding only the tomographically reconstructed profiles with sufficient slant data (6803), about 32 % of the profiles match very well and about 23 % do not match at all.For almost half of the profiles (45 %) the classification was rather difficult as they show well-matching parts together with some discrepancies.However, it must be pointed out that there are very few observing systems which provide spatially resolved humidity fields and that most humidity observations come with a rather large error.A fraction of 77 % profiles with no major discrepancies to the RS profiles and 32 % of wellmatching profiles is therefore a good result for the tomography.
There are several reasons for the large number of indifferent profiles which do not match very well with the RS profiles but also show no major discrepancies.One fundamental problem is the different nature of the data: RS take point measurements which represent the atmospheric state only in close vicinity to the RS sensor, while the GPS tomography provides voxel means which cover 50 km × 50 km horizontally and several hundred meters vertically.It cannot be expected that both observations are almost identical.Another reason are the errors of the observation systems.RS profiles can sometimes be rather wrong and it is difficult to find a statistical representation for these errors.The GPS tomography has to deal with several difficulties.The STD data provided by the GPS processing systems are the result of a rather sophisticated analysis process and the error is difficult to estimate (Deng et al., 2011).These data enter an ill-posed inverse reconstruction process which potentially amplifies the errors of the input data.Furthermore, the STD data are incomplete and do not cover the whole region, which leads to a rather low resolution of the tomographic reconstruction (Bender et al., 2011).Due to the low vertical resolution and the necessity to interpolate from the spatial grid on the RS profile, it is not possible to resolve details of the RS profile.Especially thin layers of increased or reduced humidity cannot be resolved, as long as their vertical dimensions are below the vertical grid spacing.

Conclusions and outlook
As outlined in the introduction, our aim is to evaluate the degree of agreement between radiosonde and GPS tomography based on one whole year of RS data.A GPS tomography system at developed at GFZ uses a large data set of slant delay, which is a quantity integrated along the signal path through the atmosphere.The humidity field is reconstructed with tomographic techniques.The GPS tomography uses a grid model with a vertical resolution of 500 m and horizontal resolution of about 50 km, which cannot be much better than the mean GPS interstation distance 40 km.The grid was initialized by interpolating the synoptic data onto the grid nodes.For comparison with the RS profiles the gridded data were interpolated on the profile.
The limitations of GPS tomographic technique are due to errors in the slant wet delays and their poor geometric distribution.The GPS contribution to the tomography is highly variable in space and time due to the variation of GPS satellite constellation, and a uniform quality of the reconstructed fields can therefore not be expected.Therefore, these studies have set a minimum number of slant paths and minimum difference cutting angles dφ y and dφ x for the used tomographic profiles.In this study a total of 6803 profiles, which have been considered with sufficient slant data using the described criteria, were validated with different methods.The results of point-by-point validation show a mean difference of wet refractivity within −1.3-0.3 and an overall root mean square of about 6.5-9 compared to radiosonde data, depending on different radiosonde stations.Thereof the comparison of tomographically reconstructed humidity profiles with radiosonde profiles shows about 32 % of profiles being good, but there is a non-negligible fraction of profiles with artifacts, especially in the lower atmosphere, and there are regions with insufficient GPS observations.Profiles with large discrepancies are considered poor-agreement profiles (23 %).This result provides a qualitative but aslo quantitative agreement in percent or in profile numbers.The change of different input parameter for tomography changes the result at the same time.This classification could be helpful to test the suitable parameters for reconstruction algorithms.
Apparently, the validation results depend on both the radiosonde and tomographic data.The quality of GPS SWDs is impaired not only by the experimental error but also by several assumptions and approximations made during the GPS data analysis.The quality of radiosondes is mainly affected by sensor error and uncertainty in time.In addition differences in atmospheric conditions sampled at different locations and times should be noted too.That means the error caused by the time and location interpolation.Especially when atmospheric situations change rapidly from heavy rain to dry periods and contrariwise, the differences between GPS tomography and radiosonde are larger.Furthermore, radiosondes provide in situ observations, while the GPS tomography estimates voxel means for discrete periods of time.
The number of GNSS slant observations will increase in the near future.GLONASS is already fully operational and slant data are basically available.Galileo and COMPASS observations will become available within the next years.Together with the continuously increasing number of GNSS stations this will lead to a considerable improvement of the www.ann-geophys.net/31/1491/2013/Ann.Geophys., 31, 1491-1505, 2013 reconstruction quality.In parallel the reconstruction algorithms will be improved and better error estimates for the tomographically reconstructed fields will be provided.

Fig. 1 .
Fig. 1.Location of the 272 GPS receiver sites (red dots) inside the tomography grid (black) and the German radiosonde stations (white rhombs).There were no data for station Altenstadt in the year 2007.A 14 × 18 × 20 cell grid was chosen for the reconstruction.

Fig. 2 .
Fig. 2. Bilinear/spline interpolation: the rhomb points are the center of the voxels; the round points are the horizontal locations of radiosonde stations in different layers.In the horizontal the bilinear form with four adjoining points and in the vertical cubic splines are used for the interpolation.

Fig. 3 .
Fig. 3. Spatial coverage of the atmosphere by GPS slant paths in Germany.The flag shows the location of COPS investigation area.The sampling rate of STD is 2.5 min.

Fig. 4 .
Fig. 4. Vertical profiles at station Kümmersbruck on 28 March 2007, 12:00 UTC (left panel) and the number of slant paths per voxel in four adjoining columns for interpolation (right panel).

Fig. 5 .
Fig. 5. Profiles on 31 January 2007, 12:00 UTC, at station Bergen (left) and on 9 June 2007, 12:00 UTC, at Schleswig (right).The grid cells indicated by the black boxes were chosen as examples as these cells are closest to the RS stations.

Fig. 6 .Fig. 7 .
Fig. 6.The histograms of the crossing angles of every two slants in the chosen voxels in Fig. 5 (left) at station Bergen.The corresponding cutting angles with axis vectors e x , e y , e z and a spread of dφ x = 142.75 • , dφ y = 129.37 • , dφ z = 53.81• are shown on the right side.

Fig. 8 .
Fig. 8. Mean difference of N w and their rms with all (6803) tomography and radiosonde profiles in 20 layers (left panel); relative mean differences of N r w and their rms r at the right side.

)Fig. 9 .
Fig. 9. Example profiles from the station Kümmersbruck on 3 January 2007, 12:00 UTC.Used parameter for the validation.m: maximal absolute difference of N w ; ZWD: zenith wet delay (left panel); D: difference ZWD; K: integral of absolute difference wet refractivity (right panel).

Fig. 10 .
Fig. 10.Example profiles for the different quality classes, poor (a, b, c) and good (d).

Table 1 .
Number of unused profiles according to different criteria.

Table 3 .
Criteria for unmatched profiles (left) and well-matched profiles (right) between radiosonde and GPS tomography.

Table 4 .
Profiles were divided into 3 classes.