Atmospheric anomalies over Mt . Etna using GPS signal delays and tomography of radio wave velocities

Due to the prominent topography of Mt. Etna, the use of satellite geodetic techniques may significantly suffer from atmospheric heterogeneities. This problem mainly affects the DInSAR technique. To overcome these drawbacks the present study attempts to make headway in measuring and interpreting atmospheric anomalies. We used the GAMIT software to obtain the ZTD (Zenith Total Delay) values for the GPS sessions performed in 1996-1997, during ERS-2 passes at Mt. Etna. GAMIT software also characterizes the statistical behaviour of the tropospheric effects, by using residuals for each station-satellite pair, and locates the atmospheric anomalies, present mostly at low altitudes. The attempt at using these results to produce a tomography of radio waves velocity of the troposphere suggests that the number of GPS stations used to investigate atmosphere is a critical point in such a study. The three stations are too few to invert any anomalies existing in the lower atmosphere. This result is a good starting point for better direct future study to verify the applicability of this tomographic technique to a geodetic network with a higher number of stations, with the aim of characterizing the lower atmosphere of Mt. Etna for a more accurate monitoring of ground deformations. Mailing address: Dr. Giuseppe Puglisi, Istituto Nazionale di Geofisica e Vulcanologia, Sezione Catania, P.zza Roma, 2, 95123 Catania, Italy; e-mail: puglisi-g@gt.ingv.it


Introduction
Monitoring of ground-deformations belongs to geophysical methodologies applied to the study of volcanic and seismogenetic areas for over a century.Satellite based techniques, such as GPS (Global Positioning System) or differential SAR interferometry (DInSAR), are used to meas-ure ground-deformations at Mt. Etna volcano.DInSAR measurements are based on the comparison between two SAR images, acquired at different times, during the ERS-1 or ERS-2 satellite passes over Mt.Etna.If on the surveyed areas a ground deformation phenomenon occurs between the two passes, differences in the path between the sensor and the ground in the two passes lead to a difference in the phase of the back-scattered radio wave.Based on this assumption, the interferogram resulting from the comparison between the two images gives a picture of the ground-deformation pattern, measured as displacement along the line of the sight of the sensor.
The main question on carrying out effective monitoring is due to the non-ideality of transmissive medium; electromagnetic waves, emitted from satellites, propagate in atmosphere at veloc-ity lower than the vacuum, depending on the local characteristics of the medium.The layers that most affect radio waves propagation are the ionosphere and neutral atmosphere (stratosphere and troposphere), which is the most important term in the DInSAR technique because of its large spatial and temporal variability.As the main effect is due to the troposphere (Hanssen, 1998), effects from both tropospheric and stratospheric regions have been indicated as «tropospheric» by geodesists (Gregorius and Blewitt, 1998), referring to the area of the atmosphere which is 0-50 km above the Earth's surface.
The differences in the path between the sensor and the ground may be, partially or totally, also due to differences in the tropospheric properties between the two satellite passes; furthermore, the prominent topography of the volcano, reaching an elevation of about 3000 m a.s.l., implies that the atmosphere condition might be dramatically different at the top of the volcano with respect to its base.Thus, the removal or the reduction of this unwelcome effect on the interferogram is a critical point in the DInSAR technique, because it could be confused with a deformation of volcanic edifice (Puglisi and Coltelli, 2001).To overcome this problem, it is possible to calibrate interferograms using measurements carried out with GPS (Williams et al., 1998), if on the surveyed area a GPS geodetic network exists for ground-deformation monitoring.
Such an approach has already been adopted in previous researches to analyse the same data set (Bonforte et al., 2001) or similar data set (Wadge et al., 2002;Webley et al., 2002).However, the present work represents progress for two aspects: the software used to analyze GPS data is different from that adopted in the previous studies and, moreover, the delay deriving from GPS data is here exploited in an attempt to obtain a tomography of the lower atmosphere.This analysis is not usual in atmospheric application of GPS and represents a new approach in geodetic applications at Mt. Etna.

Atmospheric effects on GPS signals
GPS satellites transmit radio frequency signals, the L1 and the L2 carrier waves, that are modulated with pseudorandom noise codes.For geodetic purposes, carrier phase measurements are performed: phase difference between the signal received from every satellite and an equal signal produced from every receiver is observed.A common technique used to mitigate error sources is the double differencing of measurements between receivers and satellites.This effectively eliminates the receiver and satellite clock errors and reduces orbital errors and other spatially correlated errors (for details see e.g., Hofmann-Wellenhof et al., 1992).Denoting A and B receivers that simultaneously observe the same satellites k and j, the doubledifference equation is given as follows: (2.1) where λ is the wavelength of the signal, φ jk AB is a linear combination of the single difference carrier phase observations, ρ jk AB is an analogous linear combination of geometric ranges, N jk AB of integer ambiguity, ∆I jk AB of errors due to ionospheric refraction of radio signals, ∆T jk AB of errors due to tropospheric refraction, ∆m jk AB of error due to multipath and ε jk AB of error due to receiver.
Given that the ionosphere is a dispersive medium, it is possible to eliminate its effects on GPS radio signals through a dual frequency method, using an «ionospheric-free» Linear Combination (LC) of the L1 and L2 phase measurements (Bender and Larden, 1985;Bock et al., 1986;Dong and Bock, 1989).Instead, the neutral atmosphere is a nondispersive medium with respect to radio waves up to frequencies of 15 GHz, and thus the propagation is frequency independent.Consequently, tropospheric refraction cannot be eliminated through a dual frequency method.The refractive index of troposphere n is larger than that of free space (n=1) and varies with altitude, latitude and meteorological conditions; this causes a decrease of the light velocity in the medium below its vacuum value c (c≅3 ⋅10 8 m/s) and a bending of the signal path with respect to the geometric straight line path.Thus slowing and bending effects cause an excess path length or tropospher- where n(s) is the refractive index as a function of position s along the actual curved ray path, while ∆g is the difference S −Sg between the length of the actual curved path and the straight line path from satellite to receiver.This last difference is typically quite small S −S g < 0.1 m, except at low elevation angles θE (Saastamoinen, 1973).The error caused by neglecting of path curvature is less than 3 mm for θE ≥ 20°and 2 cm for θE=10° ( Spilker, 1996).As we fix cut-off angle to 15°, we can neglect this difference.Equation (2.2) is often formulated in terms of atmospheric refractivity N, defined by N= =10 6 (n−1), rather than index of refraction.The formula adopted by the International Association of Geodesy in 1963 for atmospheric refractivity, is that of Essen and Froome (Essen and Froome, 1951), and is valid supposing that air is, with good approximation, an ideal gas (2.3) where P is the total atmospheric pressure (in millibars), T is the atmospheric temperature (in degrees Kelvin) and e is the partial pressure of water vapour (in millibars), while the empirical constants k1, k2, k3, have the following numerical values: k1=77.624°K/mb;k2=−12.92°K/mb;k3= =371900°K 2 /mb.Equation (2.3) will be used in the tomographic analysis (see Section 6).Tropospheric delay can be divided into dry delay, which is due to dry air, and in particular to N2 and O2, and wet delay, which is a function of water vapour distribution.Dry delay represents about 90% of the total delay; it reaches values of about 2.3 m in the zenith direction and varies with temperature and pressure.Measuring pressure and supposing that atmosphere is in hydrostatic equilibrium, it is usually possible to model and remove the hydrostatic delay with an accuracy of a few millimetres or better (Bevis et al., 1992).The zenith wet delay can be as small as a few centimetres in arid regions and as large as 35 cm in humid regions.Although the wet delay is

#
much smaller than dry delay, it is usually more variable and more difficult to remove, and so it is not possible to accurately predict the wet delay from surface meteorological measurements (Bevis et al., 1992).Therefore, tropospheric delay ∆T can be expressed as where ∆zd and ∆zw are respectively the dry zenith delay and the wet zenith delay, while MAP d (θE) and MAPw(θE) are respectively the dry mapping function and the wet mapping function, which describe the «geometrical» dependence between zenith and slant delays (the delay along a path of arbitrary elevation angle θE).There are different models for atmospheric refraction; in this work we have used Saastamoinen model (Saastamoinen, 1972, 1973) that estimates the dry zenith delay and the wet zenith delay (in meters), through the following expressions: (2.5) (2.6) where P, e and T refer to the observation site.The numerical coefficients in eqs.((2.5), (2.6)) depend on local latitude and altitude of the station, while the factor f accounts for the variation in gravitational acceleration with latitude λ and the height h of the station above the ellipsoid (in kilometres) (Saastamoinen, 1973) .(2.7) Various forms have been proposed for the wet and the dry mapping functions; in this work we have used Niell mapping functions (Niell, 1996) that are expressed in terms of annual variation, with sinusoidal behaviour, of seasonal dependence on atmosphere, and in terms of latitude and altitude of receivers with respect to sea level.
Niell mapping functions have a less than 1 cm accuracy even with an elevation angle of 3° ( Niell, 1996).

GPS measurements
The data set used in this work comes from GPS measurements that were carried out by Istituto Internazionale di Vulcanologia of CNR of Catania, during the ERS-2 passes over Mt.Etna volcano, from September 1996 to December 1997 (since 2001, IIV-CNR has been incorporated within Istituto Nazionale di Geofisica e Vulcanologia -INGV).The geodetic network used to carry out GPS measurements consists of three sites (fig.1a): the Istituto Internazionale di Vulcanologia roof (IIV) in Catania, about 50 m a.s.l., the Serra la Nave Astrophysics Observatory (SLN), at about 1700 m a.s.l. on the south-western side of Mt.Etna and the Etnean Volcanological Observatory (OSV), at about 2800 m a.s.l. on the north-eastern side (table I).The sites different heights allow any atmospheric anomalies present at different elevations to be investigated.
For stations in a local network, the elevation angles viewing a particular satellite will be nearly equal, producing high correlations among the estimated zenith delays (King and Bock, 2000).Therefore GPS data collected at Cagliari, Matera (belonging to the International GPS Service Network (IGS)) and Noto (belonging to European REFerence network (EUREF)) were taken into account during processing (fig.1b, table I), to enlarge the GPS local network used to study the troposphere.
Twenty-five 4 h long sessions were considered, starting 2 h before and ending 2 h after each ascending or descending ERS-2 pass over Mt.Etna (table II).Each station was equipped with double-frequency receivers (table II).

Estimation of zenith total delay
Tropospheric delays above zenith direction of each receiver, or Zenith Total Delays (ZTD), are estimated from the GPS measurements simultaneously with the station coordinates, using GAMIT software (King and Bock, 2000).Processing details are shown in table III.
Since meteorological parameters, such as temperature t, atmospheric pressure P and relative humidity H are not known at each station, they are calculated through the following expressions starting from default meteorological values at sea level (t 0= 20°C, P0=1013.25 mb, H0=50%): where T is the atmospheric temperature in degrees Kelvin and h is the station height in kilometres.Corrections to the zenith delay derived through atmospheric models, are estimated by considering it as a parameter to adjust in the least-square analysis.Although the contribution to the total delay coming from dry air should be very well modelled with respect to that deriving from the wet air, this is not assured in this case.The correction to ZTD calculated through atmospheric model indeed may be caused not only by the wet component but also by the dry one, because the default meteorological values at sea level we used in this work could be different from the «true values».
A variation in the zenith delay was allowed during the observation span of each GPS session by specifying a piecewise-linear model with stochastic constraints.The allowed variation between tabular points of the piecewise-linear function is defined by a Gauss-Markov process (Tralli and Lichten, 1990;King and Bock, 2000;Bock and Doerflinger, 2001).We fixed a «knot» spacing of 1 h, so we obtained a new ZTD value every hour for the 4 h long GPS sessions.Figure 2 shows the ZTD for each site and session at the time of ERS-2 passes over Mt.Etna (corresponding to the third «knot»).
As expected, tropospheric delay depends on station height: it diminishes with increasing station height, since the tropospheric layer traversed by electromagnetic waves is thinner.Furthermore, ZTD values corresponding to IIV, SLN and, when present, also to OSV site, have a similar trend during the investigated period.Finally, the results reported in fig. 2 are in agreement with those obtained previously (Bonforte et al., 2001), confirming the goodness of our processing.
Generally the ZTD estimations are larger during summer than winter time, although great variations are often present also during the same period.As regards SLN site, ZTD differences between different sessions are generally less than the same differences corresponding to IIV site; for example, the difference between ZTD corresponding to 17/09/1997 and that corresponding to 15/10/1997 is (6.2 ±1.5) cm (fig.3b), while for IIV site it is (17.3 ±1.8) cm (fig.3a).Although OSV site is present only in a few sessions, it has a similar behaviour to SLN site from 06/08/1997 session to 10/09/1997.This is in agreement with the greater stability of atmosphere at high altitudes than in proximity of sea level and this is more evident if we plot relative variation of ZTD with respect to a reference time (fig.4a,b).

Analysis of residual observations
Atmospheric anomalies refer mostly to irregular distribution of atmospheric water vapour throughout the troposphere.They produce mismodelled effects which are contemplated by postfit residuals produced by GAMIT software, i.e. the difference between observed and com-puted values of observations, after parameter adjustments.The postfit residuals contain multipath effects and receiver error, besides atmospheric anomalies or fluctuations.The presence of these three errors bias the residuals (Parkinson, 1996), so that we performed a statistical analysis of postfit residuals to identify them in the experimental data.In order to do this we have analyzed the residuals corresponding to the 10/09/1997 session, in which data are available for the three Etnean stations (IIV, SLN, OSV) and ZTD values are higher than in the other sessions; consequently, the adjustments estimated by GAMIT software to ZTD computed through atmospheric models are the highest.Sky-Plot of available satellites for this session is reported in fig. 5.
Residuals show an increase in fluctuations when elevation angles decrease (fig.6a-f).This behaviour could be produced by either multipath effects or tropospheric anomalies because they become larger at lower elevation angles.ronment near the antenna, thus is site-depending and elevation-depending but not altitudedepending (Georgiadou and Kleusberg, 1988).Furthermore, antennas and receivers with the same characteristics were used at IIV, SLN and OSV stations for the analysed session (see table II), so differences on residuals cannot be due to different electronic noise levels.
Moreover, the amplitudes of fluctuations (i.e. the difference between maximum and minimum residual values) increase when altitude of site decreases (compare fig. 6a,d with fig.6b,c), independently from elevation of satellites.It is unlikely that multipath effects produce this peculiar behaviour, because the multipath depends primarily on the reflectivity of the envi-Fig.6a-f.Residuals for the satellites PRN5 and PRN7 and the stations IIV, SLN and OSV with respect to the epoch of measurement; in the secondary ordinate axes the elevation angles of satellites are reported.For brevity's sake, we have omitted plots corresponding to the other visible satellites (PRN9, PRN23, PRN26).Excluding multipath effects and instrumental differences, we may suppose that the increasing of the amplitude of fluctuations of residuals with decreasing altitude of sites is due to tropospheric anomalies.To analyse thoroughly the presence of this bias in the postfit residuals, we studied their distribution.The residuals have a normal distribution (with mean value r ⎯ G and standard deviation σR G ) in absence of mismodelled errors, so we applied a chisquare test to quantify the deviations of experimental residuals distributions from the normal one.The histograms in fig.7a-f report the number of residuals in each of the 10 intervals in which we divided the entire range of residual values.We compared these values with the the-oretical residual values in order to perform the chi-square analysis.Table IV shows the reduced χ 2 0 values and the probability of finding a χ 2 value greater or equal to χ 2 0 (P( χ 2 ≥ χ 2 0 )) for each visible satellite-receiver pair, assuming that measurements have a Gaussian distribution.By assuming a confidence level of 5%, we have a «highly significant agreement» with attended distribution when P( χ 2 ≥ χ 2 0 )≥ 5%.This agreement results for PRN5, PRN9 and PRN23 satellites at SLN and OSV stations and only for PRN23 at IIV (see table IV).This confirms that the atmosphere is more stable and homogeneous at higher altitude; in particular, atmospheric anomalies are present mostly beneath 1770 m a.s.l., during the analyzed session, although a ge-Fig.7a-f.Histograms representing the number of residuals contained in each of the 10 intervals: 1) rG< r odetic network with higher number of stations would provide more detailed information.

Tomography of atmosphere
One of the aims of this study is to make headway with respect to previous research (Bonforte et al., 2001) in characterizing atmospheric anomalies.For this purpose, we attempted to carry out a three-dimensional tomography of electromagnetic wave propagation velocity fields relative to the investigated area, using atmospheric models and data from the processing of GPS measurements.The tomographic inversion aims at obtaining the spatial distribution of radio waves velocity and thus to identify anomalies with respect to values predicted by atmospheric models.To this end, we have used the Simulps12 software (Evans et al., 1994), devised for seismic tomography, which has been modified for atmospheric tomography in this work.
Since we know the satellite and receiver coordinates and the slant total delay, it is possible to determine the radio waves propagation time T R Sobs, which are input data for Simulps12 software (see eq. ( 6.7)).The propagation time of radio waves from a GPS satellite «S» to a receiver «R» can be expressed in analogy with the «ray theory» of Thurber (Thurber, 1993) (6.1)where v is the wave velocity and ds is an element of path length along the direction of wave propagation.Starting from an initial model of the radio waves velocity field (a priori informa- The residuals can be related to the model parameter changes through the following expression (Thurber, 1993): (6.3)where ∆v n are the perturbations to the velocity parameters and ∂T R S /∂vn are partial derivatives of propagation time with respect to velocity parameters.The complete system of simultaneous equations for the inversion can be written in the form (6.4) where r is the residual vector, M is the matrix of velocity parameter partial derivatives and ∆m is the vector of velocity parameter perturbations.The inversion technique used by Simulps12 is that of Thurber (Thurber, 1983), modified by the same Thurber (Thurber, 1993) and Eberhart-Phillips (Eberhart-Phillips, 1993); it leads to the following solution for the set of eq.(6.4): (6.5)where ε is the damping parameter; the matrix size ∆m is fixed by the number of velocity mod- el parameters.Thurber inversion is a «dampedleast-squares» inversion; this means that the norm of the model perturbations is weighted and combined with the squared data misfit and the combination is minimized at each iteration.In Simulps12, the distribution of radio waves velocity is represented through a threedimensional grid of nodes, in which velocity varies continuously in all directions, with linear interpolation among nodes (Evans et al., 1994).Nodes where velocity is computed are defined by the intersections of a set of planes in three orthogonal directions, one horizontal and two vertical.The a priori atmospheric model we used is characterized by a grid in which velocities vary only with respect to the altitude, while there are no variations along horizontal planes.Since radio waves velocity depends on the refraction index n (v=c/n, with c velocity of light in free space), we can express radio waves velocity in terms of refractivity N as (see Section 2) . (6.6) From the dependence of atmospheric parameters on height (eq.(4.1)) and by using the Essen and Froome eq. ( 2.3) for refractivity, we extracted the propagation velocity of radio waves for different heights from sea level (see table V).
The propagation times T R Sobs of radio waves, together with the above described a priori atmospheric model and co-ordinates of satellites and receivers, represent the input data to the Simulps12 software.They were obtained from the following expression: where ρ R S is the geometrical distance from satellite to receiver, SPD R S is the Slant Path Delay, calculated along the line of the sight of the satellite and finally (rG) R S is the post-fit residual, calculated from the GAMIT software, which allows unmodeled atmospheric effects to be considered.
We imagine ideally enclosing satellites, receivers and the portion of atmosphere that is to be analyzed inside a parallelepiped.As GPS satellites orbit at an altitude between 18000 and 22000 km and their elevation angles vary from 15°to 90°, the base of the parallelepiped is extremely large if compared to the Etnean area which is of the order of 10×30 km 2 .Since we are only interested in neutral atmosphere, we have decided to «reduce» the dimension of our system, shifting down the satellites ideally along the direction which connects them with GPS receivers.Thus we have «reduced» the height of satellites to 10 km a.s.l.obtaining a parallelepiped with a base of about 80×60 km 2 .Therefore the geometrical distance ρ R S between satellite S and receivers R is not the real one, but depends on new co-ordinates of satellite.The reduction of the height of GPS satellites led to suppose that tropospheric delay is due only to the first 10 km of atmosphere (troposphere) and not to all neutral atmosphere (troposphere and stratosphere) extending as far as about 50 km a.s.l.Although the troposphere contains about 90% of the air mass and 99% of water vapour of atmosphere, this represents a limit for this work.Atmospheric tomography was carried out for the 10/09/1997 session in correspondence with ERS-2 passes.As there are five satellites which are visible from stations IIV, SLN and OSV, we have only 15 observations for tomographic inversion.To increase the number of observations we considered also the positions of the same satellites 15 min before ERS-2 pass and 15 and 30 min after ERS-2 pass in the inversion.

Procedure for tomographic inversion
We chose a preliminary grid for inversion on the basis of the ray path distribution, so that the planes were closely spaced in the centre of the target volume, where there are the most rays, and more widely spaced in the periphery.Then, Thur- Height (km a.s.l.) 2 4 8 Velocity (km/s) 2.9972⋅10 5 2.9974⋅10 5 2.9976⋅10 5 ber's 3D ray tracing routine (Thurber, 1993) was applied using the preliminary grid structure in which velocities vary only with respect to altitude (as discussed in previous section).This allowed computing the Derivative Weight Sum (DWS) Fig. 8a-c.Derivative Weight Sum (DWS) at horizontal planes placed respectively at the height of 8, 4 and 2 km a.s.l.function which indicates the relative ray density in the volume of influence of the generic node, and the inversion Resolution Diagonal Element (RDE) matrix giving information on the influence by the ray travelling near the node on the velocity value estimated for the node itself.As these parameters indicate the level of resolution for the node itself, these were used for making changes to the grid.Finally, we chose for inversion a grid with a base of 80 ×60 km 2 , centred on the point of co-ordinates 14.92°E and 37.63°N, fixed as origin for the reference system of Simulps12.Vertical planes were placed at −12, −6, 0, 6, 12, 20 km with respect to origin, along x direction and to −16, −8, 0, 8 km, along y direction, while horizontal planes were placed at heights of 8, 4, 2 km with respect to sea level.DWS for each horizontal plane are shown in fig.8a-c.The highest values are obtained at 2 and 4 km suggesting that the tomographic results will be more reliable at these planes than at 8 km.
The outermost nodes were held fixed in the inversion.Moreover, as suggested by several authors (see for example Evans et al., 1994), we checked DWS during the inversion iterative process in order to identify the nodes for which sampling becomes poor at some step of the process: we set the condition that when DWS falls below 10 at a certain iteration in a given node, the velocity in that node at that iteration is fixed on the already existing value.We used a damping value of 50 in the least-squares inversion procedure.

Tomographic results
Figure 9a-c display the radio waves velocity distributions obtained through tomographic inversion.The obtained values for electromagnetic waves velocity seem to belong to a gas mixture with a higher density than air; in fact, regardless of atmospheric conditions, such as temperature, total pressure and mostly water vapour content, it is hardly probable that radio waves velocity (v i ± ± σi) resulting from the tomography had such small values as (2.9962±0.0001)⋅10 5 km/s, (2.9964±0.0002)⋅10 5 km/s and (2.9967±0.0002)⋅⋅10 5 km/s, respectively at the height of 2, 4 and 8 km a.s.l.If we consider the atmospheric parame-   ters at sea level (different from the standard one) which minimize the radio waves velocity values, consistently with the values recorded in the area of interest during 1996-1997, we find that the absolute value of difference between these values v mini and the smallest values of vi corresponding to each horizontal plane considered for tomographic analysis, is higher or equal to 4 times the sigma values σi.This is probably due to the reduction of the height of GPS satellites, with the implicit consequence that the greater part of air mass and water vapour of atmosphere is considered to be confined to the atmospheric layer beneath 10 km a.s.l.
Since the results of the tomography are greatly influenced by the number and distribution of the GPS stations, we tested the resolution of the images obtained by assuming a particular «feature» of the atmosphere: we considered an a priori model for radio waves velocity analogous to that used for tomographic inversion in the previous step, but supposed that there was an «atmospheric anomaly» where radio wave velocity was smaller than values calculated from the atmospheric model (Aloisi et al., 2002).The atmospheric anomaly was localized in the zone between 0 and 20 km with respect to the origin of Simulps12 reference system, along x direction, and for each plane perpendicular to y and z axes.We assumed that this value was 2.9965⋅10 5 km/s.The theoretical propagation times of radio waves were estimated in this velocity model for all the satellite-receiver pairs of the dataset used for 3D inversion.These propagation times, known as «synthetic times», were used later as input data for a tomographic inversion with an a priori model for radio waves velocity that is the same as the model used in the first inversion (the model without the «atmospheric anomaly»).A good reconstruction of the «atmospheric anomaly» should indicate a good resolution.
Figure 10a-c show that the «atmospheric anomaly» is not reproduced.However, the smallest value obtained at 2 km (2.9969 ±0.0001)⋅10 5 km/s) near IIV station and the values obtained at 4 km near the SLN and OSV stations (2.9970 ± ±0.0001) are relatively close to the expected one (2.9965±10 5 km/s), in agreement with the obtained DWS values that indicate ray densities at 2 and 4 km higher than at 8 km.These results confirm that the number of GPS stations is a critical point in the investigation of atmosphere by means of GPS signals delay.

Discussion and conclusive remarks
The aim of the present study was to make headway in measuring and interpreting atmospheric anomalies, following on from results of previous researches (e.g., Bonforte et al., 2001).ZTD values were obtained by means of the GAMIT software for the GPS sessions performed on 1996-1997 during ERS-2 passes at Mt. Etna volcano, by correcting the values estimated from atmospheric models.The results are in agreement with those obtained with Bernese software (Bonforte et al., 2001).GAMIT also provided residuals for each station-satellite pair allowing to characterize the statistical behaviour of the tropospheric effects and to locate the atmospheric anomalies, mostly present at low altitudes.Furthermore the analysis of residuals has shown that the analysed signals did not suffer high multipath effects and that the noise level was low.The attempt at using these results to produce a tomography of the troposphere suggests that the number of GPS stations used to investigate atmosphere is a critical point in such a study.Indeed, the three stations used in this work seem too few to invert any existing anomalies in the lower atmosphere.Furthermore they are insufficient to fix a priori grid fine enough to obtain the spatial resolution needed to correct DInSAR images.This result however represents a useful starting point to better direct future studies to verify the applicability of this tomographic technique by exploiting data provided by a geodetic network with higher number of stations.This will characterize the lower atmosphere of Mt.Etna with higher spatial resolution, providing suitable information to correct DInSAR images.
technical support and to Dr. R.W. King (Massachusetts Institute of Technology, Cambridge) for stimulating discussions and suggestions in this study.
stations used in this work.Columns from left to right respectively display station name, latitude, longitude, ellipsoid height and the network to which the station belongs (INGV = Istituto Nazionale di Geofisica e Vulcanologia; IGS = International GPS Service Network; EUREF = European REFerence network).

Fig.
Fig. 1a,b.a) The Mt. Etna volcano and the Etnean GPS measurement sites are shown; b) the geographic location of the investigated area, in Southern Italy and the three ITRF stations taken into account during processing.

Fig. 2 .
Fig. 2. ZTD values corresponding to each Etnean site and session, at the time of ERS-2 passes over Mt.Etna (21:16 and 09:41 respectively for ascending and descending passes).

Fig
Fig. 4a,b.Variations of ZTD for the ascending (at 21:16) (a) and descending (at 09:41) (b) ERS-2 passes over Mt.Etna with respect to a reference time, that is conventionally assumed as the session carried out in a) 25/09/1996 and in b) 02/10/1996.The error bars are the errors on the variations of ZTD.
software calculates propagation times T R Scal from eq. (6.1).The misfit between observed propagation times T R Sobs and propagation times calculated by Simulps12, gives the residuals r R S (they are different from residuals (rG) R S obtained by GAMIT software) (6.2) Fig.9a-c.Radio waves velocity model obtained through tomographic inversion at horizontal planes placed respectively at the height of 8, 4 and 2 km a.s.l.

Fig
Fig.10a-c.Results of the resolution test for radio wave velocity model obtained through tomographic inversion at horizontal planes placed respectively at the height of 8, 4 and 2 km a.s.l.Vertical black lines indicate the zone where we have supposed there is the «atmospheric anomaly».

Table II .
List of the GPS sessions and the SAR passes over Mt.Etna considered in this work.For ascending passes, SAR satellite passes over Mt.

Table III .
Processing parameters for the GAMIT software.The parameter «Station error» indicate the variance of GPS measurements for each receiversatellite pair.

Table IV .
The obtained reduced χ 2 0 values and the probability of finding a χ 2 value greater or equal to χ 2 0 (P(χ 2 ≥ ≥ χ 2 0 )) for each satellite-receiver pair, if we assume that measurements have a Gaussian distribution.

Table V .
The