# Three-dimensional ionospheric tomography of post-seismic perturbations produced by the Denali earthquake from GPS data

R. Garcia, F. Crespon, V. Ducic, P. Lognonné

► To cite this version:

R. Garcia, F. Crespon, V. Ducic, P. Lognonné. Three-dimensional ionospheric tomography of post-seismic perturbations produced by the Denali earthquake from GPS data. Geophysical Journal International, 2005, 163, pp.1049-1064. 10.1111/j.1365-246X.2005.02775.x. insu-03601099

HAL Id: insu-03601099

https://insu.hal.science/insu-03601099

Submitted on 8 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-dimensional ionospheric tomography of post-seismic perturbations produced by the Denali earthquake from GPS data

R. Garcia,<fnref n="1" /> F. Crespon,<fnref n="1" />,<fnref n="2" /> V. Ducic<fnref n="1" /> and P. Lognonné<fnref n="1" />

Accepted 2005 August 4. Received 2005 June 30; in original form 2004 July 29

## SUMMARY

The coupling between the solid earth and its atmosphere is responsible for vertically propagating infrasonic waves generated by seismic surface waves. These pressure waves are amplified as they propagate upward, and produce perturbations of ionospheric electron density. The electron density perturbations above California, due to the seismic surface waves generated by the Denali earthquake on 2002 November 3, have been imaged from GPS data by a tomographic method. The integrated electron content along GPS ray paths presents a noise level that is lower in the acoustic wave frequency band than in the gravity wave frequency band. Therefore, the filtered GPS data from Californian networks are inverted for a tomographic reconstruction of electron density perturbations in the acoustic frequency band. The inversion is properly resolved only in a small number of areas due to the geometry of GPS ray paths. In these areas, a wave propagating upward at 1.2 ± 0.3 km s⁻<fnref n="1" /> and horizontally at 4 ± 1 km s⁻<fnref n="1" /> is observed, with a timing consistent with an infrasonic wave generated by the path of seismic surface waves. The discrepancies between the observed electron density perturbation structure and the expected infrasonic wave can be explained by the poor resolution of the inverse problem or by a simple model of interactions between the neutral wave and the plasma. Future development of dense GPS networks and the advent of the Galileo system will overcome the resolution problems, and allow us to relate ionospheric perturbations to the seismic signal. Such a relation can be used to constrain the source and propagation of seismic waves as well as upper atmosphere characteristics.

**Key words:** acoustic waves, ionosphere, GPS, seismic waves, solid-atmosphere coupling, tomography.

## 1 INTRODUCTION

The Earth's ionosphere presents various kinds of electron density perturbations coming from above or from below. The most widely known long wavelength perturbations are the travelling ionospheric disturbances, which are associated to atmospheric internal gravity waves (Hooke 1968; Huang et al. 1998). However, infrasonic waves could also generate ionospheric perturbations. Such effects have been observed for 40 yr due to atmospheric explosions and to post-seismic atmospheric infrasounds (see Blanc 1985, for a review). These studies have demonstrated the mechanical coupling between the solid earth, the ocean and the atmosphere. A theory has been developed in order to incorporate the radiation in the atmosphere (Lognonné et al. 1998) and the atmospheric attenuation by viscosity effects (Artru et al. 2001). This coupling can also generate gravity waves by the interaction between the ocean and the atmosphere in the case of tsunami propagation (Peltier & Hines 1976) and was confirmed recently by observations (Artru et al. 2005). The phys-

ical mechanism producing the atmospheric infrasonic waves from seismic surface waves and their propagation are depicted in Fig. 1. The coupling at the Earth's surface generate infrasonic and gravity waves propagating upwards, amplified by the exponential density decrease of the atmosphere (Pokhotelov et al. 1995; Artru et al. 2001). Due to the viscosity of the atmosphere and non-adiabatic effects, these waves are attenuated (Pitteway & Hines 1963; Bass et al. 1984). Thus, the atmosphere acts as a low pass filter allowing only infrasonic waves with frequencies smaller than 0.05 Hz to reach the altitude of the *F* ionospheric layer (Farges et al. 2002). Below this frequency, the neutral waves interact with the ionosphere generating electron density perturbations.

The perturbations associated with these waves present a wide range of applications. First, these perturbations can be inverted to characterize the source of the waves. For seismic waves, the mapping of ionospheric perturbations associated to earthquakes can constrain the rupture process, or allow the reconstruction of the wavefield at the Earth's surface (Najita & Yuen 1979), including the

<footnote n="1">¹Département de Géophysique Spatiale et Planétaire, IPGP, CNRS UMR7096, 4 Ave de Neptune, 94107 St Maur des Fossés, France.</footnote>

<footnote>E-mail: garcia@ipgp.jussieu.fr</footnote>

<footnote n="2">²Noveltis, 2 Ave de l'Europe, 31520 Ramonville Saint Agne, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Diagram describing the coupling between the seismic waves, the infrasonic pressure waves and the ionosphere. The pressure waves are amplified through their vertical propagation in the atmosphere. Vs is the speed of sound in the atmosphere, and Vr is the velocity of seismic surface waves.</figcaption></figure>

measurement of the horizontal velocity of seismic surface waves with a precision compatible with earth tomographic models (Ducic et al. 2003). Moreover, by knowing the source mechanism, the atmospheric propagation of these waves can be studied to characterize the upper atmosphere and the ionosphere. For example, the study of the attenuation mechanism of these waves in the atmosphere can constrain the upper atmosphere viscosity and thermal conductivity (Artru et al. 2004).

Due to the development of new observation systems, the study of electron density perturbations of the ionosphere has recently gained some new interest. Two main observation systems have been used to map the ionospheric perturbations. First, Doppler ionospheric sounders have proven to be very efficient at observing the ionospheric perturbations associated with earthquakes at the base of the ionosphere (Blanc 1985; Chmyrev et al. 1997; Artru et al. 2004). Secondly, a new powerful ionospheric observation system is now given by dense networks of GPS receivers. This tool has been used to observe ionospheric perturbations associated with human explosions (Calais et al. 1998), atmospheric gravity waves (Calais & Minster 1995), seismic waves (Ducic et al. 2003) and tsunamis (Artru et al. 2005). These studies have detected the ionospheric perturbations and discussed the source mechanisms and the propagation characteristics. However, these analyses are subject to an important limitation related to their ability to reconstruct the 3-D structure of the electron density perturbation from the integrated total electron content (TEC) perturbation along satellite to receiver rays (Georges & Hooke 1970). In particular, the vertical resolution is important to separate infrasonic waves from gravity waves, the latter having a much smaller vertical velocity. The vertical resolution is also important for a better observation of seismic surface waves because the maximum of the energy spectrum will occur at decreasing frequencies with increasing altitudes.

First, the data acquisition system, the data analysis of the GPS observables to TEC and the data corrections for multipath effects are presented in the next section. Then, data filtering and the 3-D inversion method are presented. These algorithms are tested on a synthetic data set in order to investigate the intrinsic resolution of the inverse problem. The inversion process is then applied to GPS data from Californian networks in order to retrieve the structure of the ionospheric post-seismic perturbations associated with the Denali earthquake on 2002 November 3. Finally, our ability to characterize the source of the atmospheric infrasounds from these data, by taking into account the interaction between the neutral wave and the ionosphere, is discussed, as well as perspective with higher density networks at faster sampling rates, and limitations of the available 30 s sampling GPS data.

## 2 GPS DATA AND POST-SEISMIC IONOSPHERIC SIGNAL

### 2.1 GPS data and total electron content

The GPS data used in this study come from 30 s sampling rate receivers of the South Californian Integrated GPS Networks (SCIGN), the Bay Area Regional Deformation (BARD) and the International GPS Service (IGS). Although the raw 1 s sampling GPS data is demonstrated below to be much superior in quality, our main data are the only data provided freely by the GPS dense networks at 30 s sampling rate.

These dense GPS networks, shown in Fig. 2, have been mainly developed for geodetic applications such as crustal deformation monitoring. Nevertheless, dense networks, in which the distance between GPS receivers is of the order of a few tens of kilometres, appear to be powerful observation systems for ionospheric remote sensing. By combining dual frequency GPS receiver data one can compute
<figure><img src="image_3.png" /><figcaption>Figure 2. location of the GPS receivers used in this study (networks BARD,SCIGN,IGS).</figcaption></figure>

the slant path delay due to the TEC of the ionosphere. The GPS data, for both frequencies $f_1$ (1575.42 MHz) and $f_2$ (1227.60 MHz), are respectively, the pseudo-ranges $P_1$ and $P_2$, and the phase measurements $L_1$ and $L_2$. Therefore, the slant TEC, expressed in TEC units, for a receiver-satellite couple at epoch $t$ along ray $i$ is given by the ionospheric combination (eq. 1).

$$d_i^t = K [L_{gf}^t - \langle L_{gf} + P_{gf} \rangle], \quad (1)$$

with

$$L_{gf} = \lambda_1 L_1 - \lambda_2 L_2, \quad (2)$$

$$P_{gf} = P_1 - P_2, \quad (3)$$

$$K = \frac{8\pi^2 m_e \epsilon_0}{e^2} \frac{f_1^2 f_2^2}{(f_1^2 - f_2^2)}, \quad (4)$$

where $\lambda_1, \lambda_2$ are the wavelength for $f_1$ and $f_2$, respectively, $e$ is the charge of one electron, $m_e$ is the mass of one electron and $\epsilon_0$ is the vacuum permittivity. The coefficient $K$ is derived from the second-order approximation of the refractive index of the ionosphere (Artru 2001). All non-dispersive effects on pseudo-ranges and phase measurements are avoided by the geometry-free linear combinations (eq. 2) and (eq. 3) (Mannucci et al. 1999). The averaged sum is subtracted to resolve the phase ambiguities. This data processing is biased by dispersive errors due to electronic systems. These biases are the interfrequency bias (IFB) and the transmitter group delay (TGD) for receivers and satellites, respectively (Sardon & Zarraoa 1997). Modelling of the slant TEC $d_i^t$ taking into account the electron density and the electronic systems biases is discussed in the next section.

The accuracy of TEC estimates is related to the $L_1$ and $L_2$ accuracy in RINEX format. Effectively, because our interest is focused on high-frequency TEC perturbations, the constant term $\langle L_{gf} + P_{gf} \rangle$ can be neglected, and high-frequency variations of $d_i^t$ have an accuracy proportional to the accuracy of $L_{gf}$. Despite a value written with three to five significant decimal digits in the RINEX file, $L_1$ and $L_2$ have an accuracy of only three significant decimal digits because the last digits are fixed for a receiver-satellite pair. Thus $L_1$ and $L_2$ have, respectively, an accuracy of $10^{-3}\lambda_1$ and $10^{-3}\lambda_2$, and the accuracy of $d_i^t$ is about $4.10^{-3}$ TEC units. So, the given $L_1$ and $L_2$ accuracy in RINEX format leads to a low effective sensitivity of GPS receivers to the slant TEC that is about 5 per cent of the peak-to-peak amplitude of the post-seismic ionospheric signal described below.

## 2.2 Post-seismic ionospheric signal

The Fig. 3 is presenting the slant TEC along the ray path from SBCC receiver to GPS satellite PRN08, high pass filtered above 3.3 mHz. The post-seismic ionospheric signal on this 1 s sampling GPS receiver is clearly visible around 22.8 hr GMT. Its peak-to-peak maximum amplitude is about 0.08 TEC units and its frequency content is in the range 4–15 mHz. The lowest frequency of this range is about the acoustic cut-off frequency of the atmosphere, below which acoustic waves could not propagate, and the higher frequency is explained by the atmospheric attenuation of high-frequency acoustic waves during their upward propagation (Artru et al. 2004). Panels (a), (b) and (c) of Fig. 4 give an example of 30 s sampling TEC data filtered in the seismic frequency band between 4.3 and 5.8 mHz, for all the receivers to the satellite PRN26 during three consecutive days. The filtered TEC amplitude data are represented by a colour code as a function of time and epicentral distance of the ionospheric piercing point, which is defined by the piercing point between the receiver-satellite ray and a thin ionospheric shell at 350 km altitude. On 2002 November 3 (panel b), a ionospheric perturbation is propagating horizontally at about 3.5 km s⁻¹ along great-circle paths passing through the epicentre, and it appears about 15 min after the arrival of seismic surface waves on the ground. This signal does not appear on the previous and following days. This perturbation has been described in detail by Ducic et al. (2003). The peak-to-peak amplitude of this signal is six times larger than the average noise level at $2\sigma$ in this frequency band. The structure and time evolution of this ionospheric perturbation is fully compatible with a production by infrasonic waves generated by Rayleigh seismic surface waves.

## 2.3 Noise estimate

In the framework of the Denali earthquake study we estimate the power spectral density of the GPS observations provided by eq. (1). We used more than 70 000 time-series sampled at 30 s, which are continuous $d_i^t$ of more than 100 time steps without cycle-slips, corresponding to the observations of the three following days: 2002 November 2, 3 and 4. Next, the logarithmic average of the power spectral density, and its standard deviation were computed and presented in Fig. 5. The power spectral density highlights the difference of power between the gravity and acoustic frequencies domains. The high standard deviation value reveals the permanent disturbed state of the ionosphere due to gravity and acoustic waves generated by various sources (Blanc 1985; Calais & Haase 2003). Finally, the noise level in the frequency band of the expected post-seismic acoustic waves is computed. The rms noise level for the frequency band
<figure><img src="image_4.png" /><figcaption>Figure 3. Slant TEC along the ray path from 1 s sampling GPS receiver SBCC to satellite PRN08 is presented after a high pass filtering above 3.3 mHz. Spectrogram and time-series are shown, respectively, on top and bottom of the figure.</figcaption></figure>

(4.3–5.8 mHz) used in the following section, is to $2 \times 10^{-3}$ TEC units and $5 \times 10^{-3}$ TEC units at one standard deviation, close to the GPS receivers accuracy.

## 2.4 Aliasing and multipath effects

In addition to electronic system biases, the slant TEC GPS data are subject to two clearly identified sources of noise: the aliasing of signals at periods lower than 60 s in the 30 s sampling data, and multipath effects that are predominant are low elevation angles.

The noise coming from the aliasing of signal at periods lower than 60 s has been analysed on a data set of five GPS receivers sampled at 1 s by comparing data resampled at 30 s with and without aliasing. Two examples of such comparisons are given in Fig. 6, where 1 s sampling slant TEC data of SBCC and DVSW GPS receivers are resampled at 30 s with and without applying an anti-aliasing filter. As seen on this figure, the peak-to-peak amplitude of the aliasing noise is highly variable from one GPS receiver to another and it is ranging from 0.01 TEC units to 0.08 TEC units. However, Fig. 7 is demonstrating that the amplitude of the post-seismic signal in the selected frequency range (4.3–5.8 mHz) is higher than the amplitude of the noise due to aliasing. Moreover, the aliasing noise is not coherent from one GPS receiver to another, and will be filtered out in the 3-D tomography inversion. Further studies will need 1 s sampling GPS data to correct for aliasing effects, or 30 s sampling GPS data already corrected for such effects.

Multipath effects are seen in Fig. 4 in the form of ‘ghosts’ of ionospheric perturbations, which are characterized by daily repetition. These effects have strong amplitudes at low elevation angles, because of reflections of GPS signals on structures surrounding

GPS receivers. The post-seismic ionospheric perturbation is better retrieved by TEC data at low elevation angles because of the structure of this perturbation. So, multipath effects must be corrected on filtered TEC data before inversion. In order to correct for the multipath effects, which may reach the same amplitude as expected signals, we developed a method derived from those proposed by Bock et al. (2000) and Nikolaidis et al. (2001). Thanks to the repeatability of GPS satellite positions from one orbit to the other, the noise resulting from multipath effects is a signal at one orbit period. The orbits of GPS satellites are close to sidereal day orbits, but their periods are changing from one satellite to the other with a difference to the sidereal day period that can be as large as 80 s (Choi et al. 2004). The orbit period repeatability allows us to estimate the multipath signal for each filtered time-series as the correlated information over a 3-day period. In fact, for the day of interest j, the corrected time-series $s_j^c$ are defined by eq. (5).

$$s_j^c = s_j - (\alpha s_{j-1} + \beta s_{j+1}). \qquad (5)$$

The coefficients $\alpha$ and $\beta$ are optimally computed by a least-square method, that is, minimizing $\|s_j - (\alpha s_{j-1} + \beta s_{j+1})\|^2$ for each time-series. Finally, in order to evaluate the multipath correction, we compute the variance reduction $VR$ and select the corrected time-series if variance reduction exceeds 70 per cent. This method is summarized by the following steps done for each filtered time-series:

(1) Readjusting the time-series to the orbit period by adding a time lag to $s_{j-1}$ (respectively, $s_{j+1}$) which is obtained from the maximum of cross-correlation between the time-series $s_j$ and $s_{j-1}$ (respectively, $s_{j+1}$) within a ±4 min time window around the sidereal day period,
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Filtered TEC data (colour bar) for all the receiver time-series to satellite PRN26, as function of time (in hours GMT) and epicentral distance range (in km) between the ionospheric piercing point at 350 km altitude along the receiver-satellite ray and the earthquake location. (a), (b) and (c) filtered TEC data for, respectively, 2002 November 2, 3 and 4 and (d) filtered TEC data on November 3 after multipath correction. On panels (b) and (d) the dashed line indicates the arrival time on the ground of the seismic surface wave. The average elevation angle of this satellite is about 36° at the time of the post-seismic signal.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Logarithmic average of the TEC power spectral density (PSD in $UTEC^2.s$) and its standard deviation (STD) as function of frequency (in Hz) for the GPS observations of 2002 November 2, 3 and 4. Dashed lines indicate from left to right: the Brunt-Väisälä frequency and the acoustic cut-off frequency of the atmosphere, and the corner frequencies of the selected frequency band.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** High pass filtered (above 3.3 mHz) slant TEC GPS data of SBCC (a) and DVSW (b) GPS receivers looking at satellite PRN08. The 1 s sampling GPS data are resampled at 30 s without any aliasing correction (top) and with an anti-aliasing filter correction (middle). The aliasing noise is estimated by a difference between the two traces (bottom).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Power spectral density of post-seismic ionospheric signal in 30 s slant TEC GPS data corrected for aliasing effects (plain lines), and corresponding aliasing noise on the same GPS receivers (dashed lines). Vertical bars indicate the frequency range (4.3–5.8 mHz) used for filtering of 30 s GPS data.</figcaption></figure>

(2) Optimally estimating the coefficients $\alpha$ and $\beta$,

(3) Computing the eq. (5) to get the corrected filtered time-series,

(4) Assessing the correction of multipath effects by computing the variance reduction between $s_j$ and $s_j^c$,

(5) Replacing $s_j$ by $s_j^c$ if the reduction variance is greater than 70 per cent.

This method has been applied to the observations of satellite PRN26 and the results are presented in Fig. 4(d). Comparing panels (b) and (d) of Fig. 4 shows that few time-series have been corrected. Fig. 8 gives an example of successful multipath correction on a single time-series. Although this simple method corrected 12 per cent of the time-series for satellite PRN26, only 6.5 per cent of all the time-series have been corrected. The low number of corrected time-series can be explained by the strong constraint imposed by the variance reduction criterion, and by a majority of satellites with high elevation

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Sample of multipath correction: filtered TEC data for the ray between the receiver AZRY and satellite PRN31 (dashed line, $s_j$), and same data after multipath correction (plain line, $s_j^c$). VR indicates the variance reduction.</figcaption></figure>

angles that are less affected by multipath effects. Moreover, the efficiency of the method is greatly improved by using 1 s sampling GPS data with precise orbit periods (Choi *et al.* 2004). The data processing and noise level estimates, presented above, provide the input and a set of parameters for the inversion described below.

# 3 INVERSION

This section presents the data filtering, the model parameters, the inversion problem and its solution and finally some inversion tests on synthetic data.

## 3.1 GPS data and filtering

The slant TEC data extracted from the processing of GPS observables can be described by the following equation at each time
step t:

$$d_i^t = \int_{\text{ray}(i)} \rho(t, \mathbf{r}) \, \mathrm{d}\mathbf{r} + TGD_k^t + IFB_i^t, \qquad (6)$$

where $d_i^t$ is the slant TEC along the GPS ray from satellite k to receiver l expressed in TEC units, $\rho(t, \mathbf{r})$ is the ionospheric electron density at position r and time step t, and $TGD_k^t$ and $IFB_i^t$ are the biases of satellite k and of receiver l at time step t, respectively. All the model parameters ($\rho(t, \mathbf{r})$, $TGD_k^t$ and $IFB_i^t$) are normalized in order to respect the units defined by $d_i^t$.

The slant TEC $d_i^t$ along the GPS ray from satellite k to receiver l is bandpass filtered between corner frequencies $F_1$ and $F_2$ ($F_1 < F_2$). The filtered signal is noted $\tilde{d}_i^t$. For acoustic post-seismic waves, $F_2$ is taken below the upper frequency of infrasonic waves imposed by the atmospheric attenuation, and $F_1$ is chosen close to the acoustic cut-off frequency of the atmosphere ($\approx 3.8$ mHz, Lognonné et al. 1998).

This procedure is implemented in order to enhance the ionospheric perturbations with frequencies between $F_1$ and $F_2$. However, because the GPS satellites are moving, the filtering effect will depend on the receiver-satellite geometry. Therefore, the apparent wave frequency of the wave will be different from one ray to another. In the worst case, a ionospheric piercing point along the receiversatellite ray at 350 km of altitude, corresponding to the average maximum of electron density, will move at about 0.5 km s⁻¹ relative to the wave front. This information is directly extracted from the slope of GPS time-series presented in Fig. 4. If we are dealing with ionospheric perturbations generated by Rayleigh waves, their horizontal speed is about 3.5 km s⁻¹. Thus, a satellite following the wave will sample an apparent frequency about 86 per cent of the wave frequency, and a satellite crossing the wave will sample a frequency about 14 per cent higher than the wave frequency. The frequency band can be enlarged in order to take into account this effect. However, this effect is responsible for an overlapping of the apparent frequency bands of acoustic and gravity waves. The data are filtered in a very narrow frequency band around 200 s period ($F_1 = 4.3$ mHz and $F_2 = 5.8$ mHz) in order to limit the contamination by gravity wave signals.

Because the receiver and satellite biases variations are slow enough, their contribution in this frequency band is neglected (Sardon & Zarraoa 1997; Mannucci et al. 1999). So, in this frequency band, eq. (6) reduces to:

$$\hat{d}_i^t = \int_{\text{ray}(i)} \Delta\rho(t, \mathbf{r}) \, \mathrm{d}\mathbf{r}, \qquad (7)$$

where $\Delta\rho(t, \mathbf{r})$ is the high-frequency perturbation of the electron density. Because this variable is a perturbation, it could be either positive or negative, and it's probability distribution is very likely Gaussian.

## 3.2 Model parameters and inverse problem

The eq. (7) defines the forward problem, and has been mainly inverted for 2-D ionospheric models, see Lognonné et al. (2005) for such an application. This equation can, however, be used for inverting a 3-D structure, sampled with a 3-D grid of nodes covering the ionospheric volume studied. The electron density perturbation $\Delta\rho(\mathbf{r})$ is discretized in a model parametrized by the electron density perturbation at each node, and interpolated between the nodes by a function ensuring the continuity of the model at the cell boundaries. The nodes are situated on a 3-D regular grid along the latitude, longitude and radius coordinates. The grid is shown in Fig. 9. The latitudinal, longitudinal and radial sizes of the grid cells are $\Delta\theta = 0.7^\circ$, $\Delta\phi = 1.2^\circ$ and $\Delta r = 20$ km. A total number of 9702 nodes constitute the grid. Between the nodes, the electron density perturbation is interpolated in spherical coordinates by the following function (Thurber 1983):

$$\Delta\rho(\mathbf{r}) = \Delta\rho(\theta, \phi, r) \qquad (8)$$

$$= \sum_{\alpha\beta\gamma} \Delta\rho_{\alpha\beta\gamma} \left(1 - \frac{|\theta - \theta_{\alpha}|}{\Delta\theta}\right) \left(1 - \frac{|\phi - \phi_{\beta}|}{\Delta\phi}\right) \left(1 - \frac{|r - r_{\gamma}|}{\Delta r}\right) \qquad (9)$$

$$= \sum_{\alpha\beta\gamma} \Delta\rho_{\alpha\beta\gamma} f_{\alpha\beta\gamma}(\mathbf{r}), \qquad (10)$$

<figure><img src="image_10.png" /><figcaption>Figure 9. Grid on which model parameters are inverted in this study.</figcaption></figure>
where $\Delta\rho_{\alpha\beta\gamma} = \Delta\rho(\theta_{\alpha}, \phi_{\beta}, r_{\gamma})$ is the model parameter giving the electron density perturbation at the node, $\theta$ is the latitude, $\phi$ is the longitude, and $r$ is the radius in the WGS84 reference frame. The function $f_{\alpha\beta\gamma}(\mathbf{r})$ expressing the interpolation between the nodes is non-zero only for $\theta_{\alpha-1} < \theta < \theta_{\alpha+1}$, and $\phi_{\beta-1} < \phi < \phi_{\beta+1}$, and $r_{\gamma-1} < r_{\gamma} < r_{\gamma+1}$. For clarity the three indices $\alpha\beta\gamma$ will be replaced by the index $j$ in the rest of the text. So, with this parametrization, eq. (7) is reformulated as:

$$ \hat{d}_i = \int_{\text{ray}(i)} \sum_j \Delta\rho_j f_j(\mathbf{r}) \, d\mathbf{r} = \sum_j \left( \int_{\text{ray}(i)} f_j(\mathbf{r}) \, d\mathbf{r} \right) \Delta\rho_j. \quad (11) $$

The inverse problem now appears as a linear problem, which can be summarized in matrix notation:

$$ \hat{\mathbf{d}} = \mathbf{G}\Delta\rho \quad (12) $$

where $\mathbf{G}$ is a matrix of the form $G_{ij} = \int_{\text{ray}(i)} f_j(\mathbf{r}) d\mathbf{r}$, in which the integral is evaluated numerically along a straight line ray connecting the receiver to the satellite. Because not all the parameters of the grid are sampled by the data, the matrix $\mathbf{G}$ includes a lot of zero values. In order to reduce the computation time, the matrix $\mathbf{G}$ is rebuilt into a new matrix including only the parameters sampled by the rays.

### 3.3 A priori information

The inverse problem is solved at each time step by a least-squares estimate of model parameters. The a priori model parameters are set to zero. The a priori covariance matrix of the model parameters at time step $t$ ($C_p^t$) is computed assuming that the electron density perturbations are proportional to the IRI2001 model at this time step (Bilitza 2001). This strong constraint is a crude approximation of the interaction between the neutral acoustic wave and the plasma (Davies & Archambeau 1998). But, in the absence of a proper modelling of such interactions, this hypothesis appears as the most reasonable. Therefore, each model parameter will have an a priori standard deviation of 0.1 per cent relative to the IRI2001 model. Moreover, a horizontal correlation length of 100 km is introduced between the parameters $\Delta\rho_j$. These two a priori constraints perform a strong damping of the underdetermined inverse problem.

The a priori covariance matrix $C_d$ of the data vector $\hat{\mathbf{d}}$ is a diagonal matrix, and for each slant TEC data, the standard deviation is set to $\sigma_i = 0.01$ TEC units, corresponding to two times the noise level

obtained in Section 2. This a priori error does not take into account possible contamination by multipath effects at low elevation angles, because these effects are supposed to be corrected by the processing in Section 2.

### 3.4 Least-squares inversion

The least-squares inversion at each time step is performed following the Newton method (Tarantola 1987). The filtered data $\hat{\mathbf{d}}$ are inverted following one of the two procedures described below. These two procedures are mathematically identical, but their computation time depends on the number of data ($n$) and the number of inverted parameters ($m$) at each time step.

The first procedure is summarized by the following formulas:

$$ \Delta\rho = (\mathbf{G}^T \mathbf{C}_d^{-1} \mathbf{G} + \mathbf{C}_p^{-1})^{-1} \mathbf{G}^T \mathbf{C}_d^{-1} \hat{\mathbf{d}}, \quad (13) $$

$$ \mathbf{C}_m = (\mathbf{G}^T \mathbf{C}_d^{-1} \mathbf{G} + \mathbf{C}_p^{-1})^{-1}, \quad (14) $$

$$ \mathbf{R} = (\mathbf{G}^T \mathbf{C}_d^{-1} \mathbf{G} + \mathbf{C}_p^{-1})^{-1} \mathbf{G}^T \mathbf{C}_d^{-1} \mathbf{G}, \quad (15) $$

where $\mathbf{C}_m$ is the a posteriori covariance matrix of the model parameters and $\mathbf{R}$ is the resolution matrix. This procedure involves the computation of the inverse of two squared matrices of size $m$ by $m$. Because the a priori covariance matrix of the data vector $\mathbf{C}_d$ is diagonal, the computation of its inverse is trivial.

The second procedure is summarized by the following formulae:

$$ \Delta\rho = \mathbf{C}_p \mathbf{G}^T (\mathbf{C}_d + \mathbf{G} \mathbf{C}_p \mathbf{G}^T)^{-1} \hat{\mathbf{d}}, \quad (16) $$

$$ \mathbf{C}_m = \mathbf{C}_p - \mathbf{C}_p \mathbf{G}^T (\mathbf{C}_d + \mathbf{G} \mathbf{C}_p \mathbf{G}^T)^{-1} \mathbf{G} \mathbf{C}_p, \quad (17) $$

$$ \mathbf{R} = \mathbf{C}_p \mathbf{G}^T (\mathbf{C}_d + \mathbf{G} \mathbf{C}_p \mathbf{G}^T)^{-1} \mathbf{G}. \quad (18) $$

This procedure implies the computation of the inverse of one squared matrix of size $n$ by $n$. The computations of the inverse squared matrix and the matrix multiplications are performed by LAPACK and BLAS mathematical libraries. The computation cost of the inverse of an $n$ by $n$ squared matrix is proportional to $n^{3/8}$. Therefore, if $n^3 > 2m^3$ the first procedure will be chosen. If not, the second procedure is applied.

<figure><img src="image_11.png" /><figcaption>Figure 10. Isosurfaces of resolution for the synthetic inversion test seen from top (on the left) and from south west (on the right). The planes delimited by yellow and red lines are the vertical cut planes used in the presentation of the results of synthetic and real data inversion, respectively.</figcaption></figure>
## 3.5 Synthetic tests and tomographic resolution

The inversion method gives an estimate of the *a posteriori* error on the electron density model and the resolution of the inverted parameters. The inversion method has been tested on synthetic TEC data computed in a model presenting electron density perturbations similar to what is expected for a post-seismic ionospheric perturbation associated with Rayleigh waves. In order to represent this electron density perturbation, a Gaussian wavelet propagating horizontally at 3.5 km s⁻¹ and vertically at 0.6 km s⁻¹ with a dominant period of 200 s has been computed by the following formula:

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between the input and output models of the synthetic test. The colour scales for the input and output electron density perturbation models are on the right and on the left of the first picture, respectively. The time evolution of the ionospheric perturbation is followed from left to right and top to bottom. t0 is the entry time of the perturbation in the inverted grid, corresponding to the GPS time 22:37:00 on 2002 November 3.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Variance reduction (in per cent) of the inverse problem as a function of time (in hours).</figcaption></figure>

$$ \Delta\rho(\mathbf{r}) = 0.01\rho(\mathbf{r}) \exp\left[-\left(\frac{r-r_{\max}}{2h_{\max}}\right)^2\right] \times \exp\left[-\left(\frac{\Phi(t, \theta, r)}{8}\right)^2\right] \cos[\Phi(t, \theta, r)], \quad (19) $$

where $\Phi(t, \theta, r) = \omega t - k_{\theta} \theta - k_{r} r$ is the phase of the wave, and $\rho(\mathbf{r})$ is the electron density predicted by the IRI2001 model. The first exponential term reduces the vertical extent of the wave in order to simulate the amplification factor of the infrasonic waves, which has a maximum at some altitude ($r_{\max} = 300$ km) and decreases on both sides of this maximum ($h_{\max} = 30$ km) due to amplification and attenuation effects (Pitteway & Hines 1963; Bass et al. 1984; Farges et al. 2002). The second exponential term reduces the horizontal extent of the wave to about two wavelengths. The cosine term is the phase variation of the wave. The angular frequency is $\omega = \frac{2\pi}{T}$, corresponding to a period $T$ of 200 s. The wave vector parameters are $k_{\theta} = -\frac{2\pi}{V_{\theta}T}$ and $k_{r} = \frac{2\pi}{V_{r}T}$ with $V_{\theta} = 0.0312^{\circ}/s$ and $V_{r} = 0.6$ km s$^{-1}$. Thus, the synthetic ionospheric perturbation is a 1 per cent perturbation propagating from north to south at the Rayleigh wave speed and from bottom to top at infrasonic wave speed. This perturbation mimics the expected post-seismic ionospheric perturbation by respecting symmetries and amplitudes, but it is not a modelling of the interactions between the neutral infrasonic wave and the ionosphere.

Synthetic TEC data have been computed by integration of the electron density perturbation along satellite to receiver rays with the actual data geometry of the Californian Continuous GPS Network on 2002 November 3. Next, these data were inverted with the inversion method described in the previous section. Fig. 10 presents isosurfaces of the resolution of the model parameters. The maximum resolution is low ($\approx 0.15$) because of the underdetermination of the inverse problem (2400 data and 7500 parameters at each time step), and the strong damping constraint imposed on the model parameters through the small apriori standard deviation (0.1 per cent of background electron density). As can be seen in the figure, the resolution is not uniform in the model space. Regions of crossing GPS rays can be resolved, whereas others not. The best-resolved regions are located between 200 and 400 km altitude, and laterally in one main region centred around 36° latitude north and 239° longitude east. We have chosen to present the results by vertical cuts perpendicular to the wave fronts and centred on this region. The output ionospheric perturbations presented in Fig. 11 are one order

of magnitude lower than the input ones. As clearly demonstrated by Fig. 10, this is due to the low a posteriori resolution of the parameters. In Fig. 11, the perturbation is clearly detected, but the output model is smoothed along the vertical and horizontal directions due to the geometry of GPS rays and to the horizontal correlation length. The poor vertical resolution does not allow the recovery of the vertical wavelength of the perturbation and the altitude of the maximum of the perturbation, but the horizontal wavelength is properly recovered. The amplitude of the signal is enhanced in the regions with the best resolution; and reduced at the southern boundary due to a lack of resolution. The time evolution of the output model is in phase with the input, allowing the recovery of the propagation speed in the horizontal direction. Thus, this synthetic example demonstrates the ability of our inversion to detect small ionospheric perturbations with the correct temporal evolution and with a crude estimate of the horizontal wavelength despite an overall smoothing, particularly in the vertical direction.

# 4 POST-SEISMIC IONOSPHERIC PERTURBATION

## 4.1 Ionospheric perturbation associated to the Denali earthquake

A seen in Fig. 4, the TEC data filtered between 4.3 and 5.8 mHz presents perturbations of peak-to-peak amplitudes between 0.03 and 0.06 TEC units, after the Denali earthquake. This signal is inverted following the procedure described above. The a priori constraints imposed by the covariance matrix allow us to reduce the amplitude of the output signal in the poorly determined regions and to enhance the output signal in the properly resolved areas. As shown in Fig. 12, the variance reduction at each time step varies from 10 to 20 per cent before and after the ionospheric perturbation crossing to 60–70 per cent when the ionospheric perturbation is passing above the GPS network. The inversion results are presented in Fig. 13 on a vertical cut perpendicular to the seismic wave front. The post-seismic ionospheric wave is clearly visible in the best-resolved region. The wave propagates from north to south and from bottom to top (movie available on http://ganymede.ipgp.jussieu.fr/~garcia). By following the maxima and minima of the wave we estimate the horizontal velocity of the wave to $4 \pm 1$ km s$^{-1}$, and the vertical velocity to $1.2 \pm 0.3$ km s$^{-1}$. These estimates are limited by the horizontal and vertical resolutions of the model grid and by the time sampling of GPS data (30 s), but they are in good agreement with a perturbation propagating horizontally at the speed of seismic surface waves and vertically at sound speed (Fig 14c from Hedin 1991). The vertical wavelength of the wave is smoothed out by the inversion process. However, the horizontal wavelength is about $750 \pm 100$ km, in agreement with surface Rayleigh waves of 200 s period propagating at $4$ km s$^{-1}$. The peak-to-peak amplitude of electron density perturbation is about $10^9$ e/m$^3$ in the best-resolved regions with about 0.1 resolution. The synthetic test and the low output resolution indicate that the true electron density perturbation could be one order of magnitude higher than the output value, reaching $10^{10}$ e/m$^3$. When compared to the ionospheric model predicted by the IRI2001 model, plotted in Fig. 14, the maximum electron density perturbation is of the order of 0.5 per cent of the background electron density.

Fig. 15 presents the Rayleigh wave vertical displacement, decimated at 30 s sampling and filtered between 4.3 and 5.8 mHz, at the broad-band seismic station CMB located at 239.6°E and 38°N. Two wavelengths are clearly seen on the ionospheric perturbation and
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Vertical cut of the output electron density perturbation model on a plane perpendicular to the seismic surface wave front (red plane on Fig. 7) as a function of time. This cut is seen from west-southwest. The electron density perturbation associated with infrasonic post-seismic waves is moving horizontally at ≈4 km s⁻¹ and vertically at ≈1.2 km s⁻¹.</figcaption></figure>

could be explained by the two wavelengths seen on the ground displacement. The peak-to-peak amplitude of the displacement source signal is ≈1 mm. Therefore, we can conclude that a Rayleigh wave seismic signal at a 200 s period with ≈1 mm peak-to-peak amplitude is able to produce a ≈0.5 per cent electron density perturbation in the ionosphere.

## 4.2 Discussion

The electron density perturbation associated with a Rayleigh wave can be imaged by TEC data from a dense network of continuous GPS receivers. Even if the picture is smoothed by the inversion process, the main features of the infrasound wave are retrieved in
<figure><img src="image_15.png" /><figcaption>**Figure 14.** From left to right: (a) electron density (e m⁻³), (b) total mass density (kg m⁻³) and (c) sound speed (m s⁻¹) in logarithmic scale as a function of altitude (in km) predicted by the IRI2001 model at 38°N and 238°E the 2002 November 3 at 22h44 GPS time, and by the MSISE-90 atmospheric model (Hedin 1991) for (b) and (c), respectively.</figcaption></figure>

the electron density perturbation structure and evolution. The most important discrepancy between the expected infrasound wave structure and the electron density perturbation structure is the inclination of the wave front. Effectively, for an infrasound shock wave propagating at 4 km s⁻¹ horizontally and 700 m s⁻¹ vertically, the inclination of the wave front relative to the ground is predicted to be around 10°. But, the inclination of the electron density perturbation wave front in Fig. 13 is about 18° ± 3°. This discrepancy is mainly due to the lack of vertical resolution caused by the geometry of the GPS ray sampling. However, a second explanation can be formulated in terms of plasma dynamics. As described by Hooke (1968) for gravity waves, the displacement of ions is constrained to be along magnetic field lines with a velocity proportional to the neutral wave velocity. Hooke's demonstration applies here for low-frequency acoustic waves because their spatial wavelength and frequency are close to those of gravity waves at the Brunt-Väisälä frequency, even if their horizontal speed is much higher. Above the GPS network, the horizontal projection of the dipolar magnetic field is approximately perpendicular to the Rayleigh wave front (within 18°), and the magnetic field inclination is about 60°. So the combined effects of the magnetic field orientation and the increase of

neutral wave velocity with altitude tend to align the wave front perpendicular to the magnetic field direction, consequently increasing its inclination.

This simple interpretation of magnetic field effects demonstrates that a full modelling of the interaction between neutrals and ions must be performed before interpreting the ionospheric electron density perturbation at these altitudes in terms of acoustic wave structure. Therefore, more studies must be performed before retrieving Rayleigh wave characteristics from GPS TEC data measurements.

## 5 CONCLUSION

A method of 3-D ionospheric tomography from TEC GPS measurements has been developed and applied to the detection of an ionospheric perturbation generated by seismic Rayleigh waves. The inversion presents only a small number of areas with a reasonable resolution due to crossing ray paths from different satellites to the receiver network. The inversion results in these regions show an electron density perturbation similar to what is expected for an infrasound wave generated by Rayleigh surface waves: in terms of timing, structure and temporal evolution. From this observation, we
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Vertical displacement on the ground at the broad-band seismic station CMB and vertical profile of electron density perturbation at the closest grid node as a function of time. CMB station is located at 38°N and 239.6°E. The seismic trace has been decimated to 30 s time sampling and filtered between 4.3 and 5.8 mHz. Scaling factors are 0.03 mm km⁻¹ and 2 × 10⁷ (e/m³)/km for the ground displacement and for the electron density perturbation, respectively.</figcaption></figure>

have deduced an amplitude scaling relating a 1 mm Rayleigh wave displacement at 200 s period to a 0.5 per cent electron density perturbation in the ionosphere. The discrepancies between the retrieved electron density perturbation structure and the infrasound wave are mainly related to the lack of vertical resolution of the inverse problem, but they can also be explained by the effect of the magnetic field on the movement of ions.

The study of post-seismic ionospheric perturbations is justified by different potential applications to the study of the structure of the solid earth and the atmosphere. For signals retrieved above lands, the ionospheric signal above source regions can be used to infer the long-period strong motion signal that cannot be obtained from seismometers, due to their saturation, nor with accelerometers, due their short-period instrument response. For ionospheric signals above the oceans, these could be used to retrieve the long-period vibrations of the Earth where the seismometers are lacking. For this application, a specific space-borne imaging system must be developed. However, by using dense GPS networks in coastal areas, we can already image the lithospheric structure at the margins. In either case, these ionospheric perturbations are able to characterize the atmospheric properties (viscosity, winds) from the propagation of atmospheric infrasounds.

Despite their great potential, these studies must face some important technical problems. First, the low sampling rate (30 s) of GPS receivers is hiding the high-frequency atmospheric dynamics, and it is creating significant noise due to aliasing of high-frequency signals. These limitations will be overcome with the availability of high sampling rate GPS data or with 30 s data sampled after applying an anti-aliasing filter. Then, GPS TEC data inversions are facing a lack of resolution due to the low number of GPS satellites visible at the same time. This second problem will be reduced by using receivers of Global Navigation Satellite Systems combining NAVS-TAR and GLONASS systems, and will certainly be solved with the advent of the Galileo system. Finally, the modelling of neutrals-ions interactions at ionospheric heights, which has been thoroughly investigated in the 1960's (Hooke 1968; Nelson 1968; Georges & Hooke 1970), would need a reinvestigation with modern numerical modelling methods.

## ACKNOWLEDGMENTS

We thank Raul Madariaga, Kristine Larson and an anonymous reviewer for improving the manuscript by their reviews. Many thanks also to Juliette Artru and Jeanine Gagnepain Beyneix for helpful discussions, and to the ionospheric team at Noveltis company. Scripps Orbit and Permanent Array Center (SOPAC) has provided GPS data, International GPS Service (IGS) has provided satellite orbits and seismic data of CMB station come from Berkeley Digital Seismograph Network. This study has been partly supported by the 'Space Weather Network' from ESA and the 'Réseau Terre et Espace' from CNES (French space agency) and French Ministry of Research. IPGP contribution number 2085.

## REFERENCES

* Artru, J., 2001. Observations au sol ou par satellite et modelisation des signaux ionosphériques post-sismiques, *PhD thesis*, Institut de Physique du Globe de Paris, Paris.

* Artru, J., Lognonné, P. & Blanc, E., 2001. Normal modes modelling of post-seismic ionospheric oscillations, *Geophys. Res. Lett.*, **28**, 697–700.

* Artru, J., Farges, T. & Lognonné, P., 2004. Acoustic waves generated from seismic surface waves: propagation properties determined from doppler sounding observation and normal-modes modeling, *Geophys. J. Int.*, **158**, 1067–1077.

* Artru, J., Ducic, V., Kanamori, H., Lognonné, P. & Murakami, M., 2005. Ionospheric detection of gravity waves induced by tsunamis, *Geophys. J. Int.*, **160**, 840–848.

* Bass, H., Evans, L., Piercy, J. & Sutherland, L., 1984. Absorption of sound by the atmosphere, in *Physical acoustic*, Vol.XVII, pp. 145–232, Orlando, FL, USA, Academic Press, Inc.

* Bilitza, D., 2001. International reference ionosphere 2000, *Radio Science*, **36**, 261–275.

* Blanc, E., 1985. Observations in the upper atmosphere of infrasonic waves from natural or artificial sources: a summary, *Ann. Geophysicae*, **3**, 673–688.

* Bock, Y., Nikolaidis, R. & de Jonge, P., 2000. Instantaneous geodetic positioning at medium distances with global positioning system, *J. geophys. Res.*, **105**, 28 223–28 253.

* Calais, E. & Haase, J., 2003. Detection of ionospheric perturbations using a dense GPS array in southern California, *Geophys. Res. Lett.*, **30**, 1628.

* Calais, E. & Minster, B., 1995. GPS detection of ionospheric perturbations following the January 17, 1994, Northridge earthquake, *Geophys. Res. Lett.*, **22**, 1045–1048.
Calais, E., Minster, B., Hofton, M. & Hedlin, M., 1998. Ionospheric signature of surface mine blasts from global positioning systems measurements, *Geophys. J. Int.*, **132**, 191–202.

Chmyrev, V., Isaev, N., Serebryakova, O., Sorokin, V. & Sobolev, Y., 1997. Small-scale plasma inhomogeneities and correlated elf emissions in the ionosphere over an earthquake region, *J. Atm. Sol. Terr. Phys.*, **59**, 967–974.

Choi, K., Bilich, A., Larson, K.M. & Axelrad, P., 2004. Modified sidereal filtering: Implication for high-rate GPS positioning, *Geophys. Res. Lett.*, **31**, L22608, doi:10.1029/2004GL021621

Davies, J. & Archambeau, C., 1998. Modelling of atmospheric and ionospheric disturbances from shallow seismic sources, *Phys. Earth planet. Inter.*, **105**, 183–199.

Ducic, V., Artru, J. & Lognonné, P., 2003. Ionospheric sensing of the Denali earthquake Rayleigh surface waves, *Geophys. Res. Lett.*, **30**, doi: 10.1029/2003GL017812.

Farges, T., Artru, J., Lognonné, P. & Le Pichon, A., 2002. Effet des séismes sur l'ionosphère, *Chocs*, **26**, 7–18.

Georges, T. & Hooke, W., 1970. Wave-induced fluctuations in ionospheric electron content: a model indicating some observational biases, *J. geophys. Res.*, **75**, 6295–6308.

Hedin, A., 1991. Extension of the MSIS thermosphere model into the middle and lower atmosphere, *J. geophys. Res.*, **96**, 1159–1171.

Hooke, W., 1968. Ionospheric irregularities produced by internal atmospheric gravity waves, *J. Atmos. Terr. Phys.*, **30**, 795–823.

Huang, C.-S., Sofko, G. & Kelley, M., 1998. Numerical simulations of mid-latitude ionospheric perturbations produced by gravity waves, *J. geophys. Res.*, **103**, 6977–6989.

Lognonné, P., Clévédé, E. & Kanamori, H., 1998. Normal mode summation of seismograms and barograms in a spherical earth with realistic atmosphere, *Geophys. J. Int.*, **135**, 388–406.

Lognonné, P. et al., 2005. Ground based GPS tomography of ionospheric post-seismic signal during demeter: the SPECTRE project, *Planet. Space Science*, in press.

Mannucci, A., Iijima, B., Lindqwister, U., Pi, X., Sparks, L. & Wilson, B., 1999. 25. GPS and ionosphere, in *Review of Radio Science 1996–1999*, pp. 625–665, eds Stone, W., Wiley-IEEE Press.

Najita, K. & Yuen, P., 1979. Long-period Rayleigh wave group velocity dispersion curve from HF doppler sounding of the ionosphere, *J. geophys. Res.*, **84**, 1253–1260.

Nelson, R.A., 1968. Response of the ionosphere to the passage of neutral atmospheric waves, *J. Atm. Terr. Phys.*, **30**, 825–835.

Nikolaidis, R.M., Bock, Y., de Jonge, P. & Shearer, P., 2001. Seismic wave observations with global positioning system, *J. geophys. Res.*, **106**, 21897–21916.

Peltier, W. & Hines, C., 1976. On the possible detection of tsunamis by a monitoring of the ionosphere, *J. geophys. Res.*, **81**, 1995–2000.

Pitteway, L. & Hines, C., 1963. The viscous damping of atmospheric gravity waves, *Can. J. Phys.*, **41**, 1935–1948.

Pokhotelov, O., Parrot, M., Fedorov, E., Pilipenko, V., Surkov, V. & Gladychev, V., 1995. Response of the ionosphere to natural and man-made acoustic sources, *Ann. Geophysicae*, **13**, 1197–1210.

Ronchi, C., Iacono, R. & Paolucci, P., 1996. The cubed sphere: a new method for the solution of partial differential equations in spherical geometry, *J. Comp. Phys.*, **124**, 93–114.

Sardon, E. & Zarraoa, N., 1997. Estimation of total electron content using GPS data: how stable are the differential satellite and receiver instrumental biases?, *Radio Science*, **32**, 1899–1910.

Tarantola, A., 1987. *Inverse Problem Theory*, Elsevier Science Publishers B. V., Amsterdam.

Thurber, C., 1983. Earthquake locations and three-dimensional crustal structure in the Coyote lake area, central California, *J. geophys. Res.*, **88**, 8226–8236.

## APPENDIX A: RAY INTEGRATION IN CUBED SPHERE COORDINATES

We present here an analytic method to integrate a field variable defined on a 3-D grid along a straight line, and subsequently construct the matrix G relating slant TEC data to the electron density in the ionosphere. The 3-D variable is interpolated between the grid nodes by two different methods: a constant value in each cell, and an interpolation between the grid nodes ensuring first order continuity at the grid cell boundaries. The grid is defined using the cubed sphere coordinates (Ronchi et al. 1996). In this coordinate system a cube is mapped on a sphere. A coordinate system is defined on each face of the 'cubed sphere'. Because we are interested in this paper in regional networks of dense GPS data, we will restrict our attention to only one face of the cubed sphere. Extension of our results to the whole cubed sphere by using the appropriate coordinate systems on each face is straightforward as long as only Cartesian and cubed sphere coordinate systems are used. Moreover, the GPS Network is supposed to be centred on the surface point of Cartesian coordinates ($R_E$, 0, 0) with $R_E$ the Earth's radius at this point. Such a situation is obtained by rotating the WGS84 Cartesian reference frame.

The cubed sphere coordinate system is defined on the face centred on ($R_E$, 0, 0) by two angles $\xi$ and $\eta$, and the radius $r$. The points of constant $\xi$ lie on a plane defined by the poles and a point of Cartesian coordinates (cos ($\xi$), sin ($\xi$), 0). The points of constant $\eta$ lie on a plane defined by the three points with the following Cartesian coordinates: (0, -1, 0), (0, 1, 0), and (cos ($\eta$), 0, sin ($\eta$)). $\xi$ and $\eta$ are spanning the range [-45°, 45°]. The relations between the cubed sphere and the Cartesian coordinate systems are defined by:

$$ \frac{y}{x} = \tan(\xi), \qquad (A1) $$

$$ \frac{z}{x} = \tan(\eta), \qquad (A2) $$

$$ r = \sqrt{x^2 + y^2 + z^2}. \qquad (A3) $$

A grid is defined by regular sampling of the variables $\xi$, $\eta$ and $r$. Grid nodes are defined by the coordinates ($\xi_\alpha$, $\eta_\beta$, $r_\gamma$), with

$$ \xi_{\alpha} = \alpha * \Delta\xi, \qquad (A4) $$

$$ \eta_{\beta} = \beta * \Delta\eta, \qquad (A5) $$

$$ r_{\gamma} = r_{\alpha} + \gamma * \Delta r, \qquad (A6) $$
where $\alpha, \beta$ and $\gamma$ are integers. The first advantage of the cubed sphere coordinate system is that a regular sampling of the coordinates produces grid cells of approximately the same size. The straight line (i.e. ray), along which the variable (i.e. electron density) is integrated, is parametrized by:

$$M(l) = (x, y, z) = M_0 + l * e = (x_o, y_o, z_o) + l * (e_x, e_y, e_z). \quad (A7)$$

With $M_0$ the receiver position, $e$ a unit vector in the receiver to satellite direction, and $l$ the distance between the receiver and a point along the line. Another important advantage of the cubed sphere coordinate system is that the borders of the cells are defined by planes. This property facilitates the determination of the intersection point between the ray line and the border of the cell. So, the ray line will intersect the border of the cell:

(i) defined by $\xi = \xi_\alpha$ at $l = \frac{x_o \tan(\xi_\alpha) - y_o}{e_y - e_x \tan(\xi_\alpha)}$,

(ii) defined by $\eta = \eta_\beta$ at $l = \frac{x_o \tan(\eta_\beta) - z_o}{e_z - e_x \tan(\eta_\beta)}$,

(iii) defined by $r = r_\gamma$ at $l = -(x_o e_x + y_o e_y + z_o e_z) \pm \sqrt{(x_o e_x + y_o e_y + z_o e_z)^2 - (x_o^2 + y_o^2 + z_o^2 - r_\gamma^2)}$.

So, the entry and exit points of the ray line in each cell could be easily computed, allowing integration of the field variable for a particular interpolation within the cell. For example, the values of the parameter $l$ corresponding to the intersection of the ray line with all the cell's borders could be tabulated and ordered in a list of increasing values {$l_c$}.

First, for a parametrization with a constant value of the variable in each cell ($\rho_j$ in the cell number $j \equiv c$), the matrix formulation of the integral defined by $d_i = \int_{\text{ray}(i)} \rho(l) dl$ is easily retrieved by setting:

$$d_i = \sum_j G_{ij} \rho_j, \quad (A8)$$

with $G_{ij} = l_{c+1} - l_c = l_{j+1} - l_j$ the length of the ray number $i$ in the cell number $j \equiv c$.

A second parametrization is defined by using a variation of the formula (10):

$$\rho(\mathbf{r}) = \rho(\xi, \eta, r) \quad (A9)$$

$$= \sum_{\alpha\beta\gamma} \rho_{\alpha\beta\gamma} \left( 1 - \frac{|\tan(\xi - \xi_\alpha)|}{\tan(\Delta\xi)} \right) \left( 1 - \frac{|\tan(\eta - \eta_\beta)|}{\tan(\Delta\eta)} \right) \left( 1 - \frac{|r - r_\gamma|}{\Delta r} \right) \quad (A10)$$

$$= \sum_{\alpha\beta\gamma} \rho_{\alpha\beta\gamma} \left( 1 - \left| \frac{(y - \tan(\xi_\alpha)x)}{\tan(\Delta\xi)(x + \tan(\xi_\alpha)y)} \right| \right) \left( 1 - \left| \frac{(z - \tan(\eta_\beta)x)}{\tan(\Delta\eta)(x + \tan(\eta_\beta)z)} \right| \right) * \quad (A11)$$

$$\left(1 - \frac{|\sqrt{x^2 + y^2 + z^2} - r_\gamma|}{\Delta r}\right) \quad (A12)$$

$$= \sum_{\alpha\beta\gamma} \rho_{\alpha\beta\gamma} f_{\alpha\beta\gamma}(\mathbf{r}), \quad (A13)$$

where $\rho_{\alpha\beta\gamma}$ is the inverted parameter at node $\alpha\beta\gamma \equiv j$, and $f_{\alpha\beta\gamma}(\mathbf{r})$ is a positive function of the Cartesian coordinates of the point, defined for $\xi_{\alpha-1} < \xi < \xi_{\alpha+1}$, $\eta_{\beta-1} < \eta < \eta_{\beta+1}$, and $r_{\gamma-1} < r_{\gamma} < r_{\gamma+1}$. Along the ray line, the Cartesian coordinates of the integration point $M(l)$ are defined by $x = x_o + l * e_x$, $y = y_o + l * e_y$, and $z = z_o + l * e_z$. If these values are replaced in the function $f_{\alpha\beta\gamma}(\mathbf{r})$, we obtain:

$$f_{\alpha\beta\gamma}(\mathbf{r}) = f_{\alpha\beta\gamma}(l) \quad (A14)$$

$$= \left( 1 - \frac{1}{\tan(\Delta\xi)} \left| \frac{(y_o - \tan(\xi_\alpha)x_o) + (e_y - \tan(\xi_\alpha)e_x) * l}{(x_o + \tan(\xi_\alpha)y_o) + (e_x + \tan(\xi_\alpha)e_y) * l} \right| \right) * \quad (A15)$$

$$\left( 1 - \frac{1}{\tan(\Delta\eta)} \left| \frac{(z_o - \tan(\eta_\beta)x_o) + (e_z - \tan(\eta_\beta)e_x) * l}{(x_o + \tan(\eta_\beta)z_o) + (e_x + \tan(\eta_\beta)e_z) * l} \right| \right) * \quad (A16)$$

$$\left( 1 - \frac{|\sqrt{(x_o^2 + y_o^2 + z_o^2) + 2(x_o e_x + y_o e_y + z_o e_z) * l + l^2} - r_\gamma|}{\Delta r} \right). \quad (A17)$$

Because the expressions in the absolute values keep the same sign in one cell, the absolute values can be removed if the expression is computed in one particular cell. The integral formulation $d_i = \int_{\text{ray}(i)} \rho(l) dl$ gives a matrix term of the form:

$$G_{ij} = \sum_c \int_{\text{ray}(i) \text{ in cell } c} f_{\alpha\beta\gamma}(l) dl, \quad (A18)$$
where the sum is over the cells adjacent to the node $\alpha\beta\gamma \equiv j$. So, for each cell $c$ involving the node $\rho_{\alpha\beta\gamma}$, along the ray $i$, we compute an integral of the form:

$$ \int_{l_c}^{l_{c+1}} \left( \frac{a_1 + b_1 * l}{c_1 + d_1 * l} \right) \left( \frac{a_2 + b_2 * l}{c_2 + d_2 * l} \right) \left( \frac{\sqrt{a_3 + b_3 * l + c_3 * l^2 + d_3}}{\Delta r} \right) dl, \qquad (A19) $$

with $(a_k, b_k, c_k, d_k, k = 1, 2, 3)$ constant terms deduced from eq. (A17) for each node $\rho_{\alpha\beta\gamma}$ and each adjacent cell $c$. This integral presents an analytical solution and is evaluated between $l_c$ and $l_{c+1}$ for each cell adjacent to the grid node $\alpha\beta\gamma \equiv j$. The sum of the integrals involving the node $\rho_{\alpha\beta\gamma}$ along the ray $i$ gives the matrix term $G_{ij}$.

This parametrization is much more complex than the first one, even if an analytical solution is obtained for the computation of the matrix term $G_{ij}$. Its use is justified by two advantages. First, its intrinsic continuity ensures the complete absence of discontinuities in the model, which is useful if ray tracing or amplitude computations must be performed by using this 3-D model. Second, this parametrization produces a kernel matrix $G$ with much more non-zero elements than the first one. By the way performing a natural correlation between the model parameters over one cell size.

This Appendix demonstrates that the cubed sphere coordinate system is able to give analytical estimates of the forward problem of integration of electron density models along straight line rays.