(electron density equal to zero at 80 km, but in the order of $10^{11} - 10^{12}e/m^3$ at 300 km), EM waves emitted by the radar are particularly sensitive to the zone where the rays are reflected, where the plasma frequency $f_p$ and the emission frequency $f_e$ are approximately the same. This defines the area of good coverage (independent of inversion method), which is the region between 100 and 300 km altitude and 500–2500 km in range (Figure 5.9).

<figure><img src="image_60.png" /><figcaption>Figure 5.9: Ratio of plasma $f_p$ to emission frequency $f_e$, calculated along the ray-path illustrating the sensitivity of the rays to the medium. The highest sensitivity is the zone where the ratio is 1.</figcaption></figure>

In summary, all grid tests show that the *v*- and the *v&r*-method generally enables the identification of perturbations above 100 km altitude. Furthermore, these tests indicate independence of the inversion results of the grid, or, in other words, the grid does not influence drastically the inversion results.

In general, for large scale perturbations, the *v&r*-method seems to perform better as compared to the *v*-method, which is illustrated in Figures 5.4 and 5.1 for the checkerboard and the localized perturbations. The checkerboard tests demonstrate that the *v&r*-method can reproduce more details of the target model, whereas the *v*-method is strongly damped and only reproduces the target up to 800 km horizontal distance.

As an even better result, both methods are able to recognize and separate small scale anomalies of 165 km × 20 km (Figure 5.6) and reproduce the correct sign of the perturbations. This works well for anomalies located at 200–300 km altitude, where
rays are reflected and therefore are most sensitive to the medium. The checkerboard perturbation tests (Figures 5.3 and 5.4) confirm this and in addition clearly indicate that anomalies at lower altitude are not well reconstructed, neither using the *v*-method nor with the *v&r*-method.

Figure 5.10 compares the model misfit for the *v* and the *v&r*-method for all performed tests, except for the second checkerboard test, as a function of the grid size for the *v* and the *v&r*-method.

<figure><img src="image_61.png" /><figcaption>Figure 5.10: Model misfit as a function of the grid size for the *v* and the *v&r*-method. The *x*-axis displays the grid size in distance × altitude in km, respectively.</figcaption></figure>

In four of the six performed tests, the model misfit is smaller for the *v&r*-method than for the *v*-method. Only for the radical checkerboard test and the perturba-
tion in one cell, the *v*-method performs better in terms of model misfit. This was
confirmed by the inversion results in Figure 5.5 and 5.7.
# Chapter 6

## Iterative approach

In order to improve the inversion results obtained with a single inversion, an iterative method was implemented. The solution found after each iteration can be used as a new *a priori* model for following inversions. Then, after a certain number of iterations, the synthetic target model should be obtained. In each iteration, rays are traced again in the new, updated model, for both the *v* and *v&r*-method described in Chapter 4, the tomography matrix is recomputed accordingly and the differences between observed (or, in the present case, synthetic) and computed propagation times are calculated and inverted.

Based on the previous synthetic tests, the grid was fixed to 25 × 20 km in distance and altitude, respectively, and the background ionosphere was generated for October at noon with a solar flux of 198 SFU. The following tests were performed for a localized electron density perturbation of 0.1 %. If not stated differently, data was simulated by tracing 1071 rays with elevation angles between 10° and 60° and in the frequency range 6 – 16 MHz. For evaluation of the difference between target and obtained solution, data and model misfits are calculated throughout iterations. The data misfit is the sum over the propagation time difference of all rays, given by

$$ \text{data misfit} = \sum_{j=1}^{\text{rays}} \left( \frac{\delta T_j}{T_j} \right)^2, \qquad (6.1) $$

where $\delta T$ is the difference between the propagation time in the target and the propagation time $T$ after each iteration in the updated model. The model misfit is given as the difference between target and solution summed over all cells that is

$$ \text{model misfit} = \sum_{i=1}^{\text{cells}} \left[ \left( \frac{\delta Ne}{Ne^0} \right)_i^{\text{target}} - \sum_{k=1}^{\text{iter}} \frac{\frac{\delta Ne}{Ne_{k,i}} \cdot Ne_{k,i}}{Ne_i^0} \right]^2. \qquad (6.2) $$

The model misfit can only be calculated for synthetic tests and not for real data. The
use of relative perturbations in this work means that in each iteration the electron
density perturbation $\delta Ne$ is linked to the previous, updated ionospheric model $Ne_k$.
For comparison with the target model that is with respect to the a priori ionospheric
model $Ne^0$, the obtained perturbation $\delta Ne$ has to be renormalized by $Ne^0$ after each
iteration.

6.1 Computational Costs

The iterative approach has a very large execution time (especially the v&r-method),
due to the calculation of the ray-path deflection that need a recalculation of the
ray-path after each perturbation of a cell. To reduce the execution time, it was
necessary to parallelize the code to take advantage of multi-core architectures of
modern computers. Therefore, the TDR FORTRAN code was modified using Open
Multi-Processing (OpenMP) to be able to run it in parallel on a multiprocessor
computer. In the following, the performances of the parallel code will be analyzed
quickly.

Figure 6.2 compares the time required to trace ~ $10^3$ rays on a grid with 1000 pixels for the v-method, using a Mac laptop with four cores and a Desktop computer with 32 cores. The Mac laptop, however, is able to execute eight threads in parallel, although it has only four cores. This feature is called hyper-threading, meaning that each core runs two threads at the same time. As expected, the OpenMP code is 32

<figure><img src="image_62.png" /><figcaption>Figure 6.1: Computational costs for numerical calculation of ray-path deflection in the v&r-method for ~ $10^3$ rays on a grid with 1000 pixels. The execution time in parallel (openMP) and on a single processor is compared using a Mac laptop with four cores and a 32 core desktop computer.</figcaption></figure>

times faster in parallel on the 32 cores desktop computer, but on the Mac laptop it
<figure><img src="image_63.png" /><figcaption>Figure 6.2: Computational costs for tracing ~$10^3$ rays on a grid with 1000 pixels for the *v*-method. The execution time in parallel (openMP) and on a single processor is compared using a Mac laptop with four and a 32 core desktop computer.</figcaption></figure>

is only around 3 × faster. This maybe associated with the hyper-threading, which
allows to run two threads on one processor and leads to slowing-down.

Figure 6.1 illustrates the time needed to calculate the ray-path deflection for ~$10^3$ rays on a grid with 1000 pixels. Certainly, the execution time is grid dependent, because the number of pixels perturbed varies with the pixel dimensions. Again, the expected gain of 32 is reached for the desktop computer, and for the Mac computer the gain is around 3. The above performed analysis showed that the *v*-method is much faster in a single iteration than the *v&r*-method in parallel on a grid with 1000 pixels and ~$10^3$ rays, with ~36s, compared to ~9 minutes on the 32 core desktop computer. That means, for one iteration with the *v&r*-method, 15 iterations with the *v*-method can be calculated in the same time.

6.2 Preliminary results

Figure 6.3 shows some preliminary results published by Roy et al. (2014) obtained with the iterative *v* and *v&r*-method applied to a localized perturbation of 0.1%. In the interest of computational speed, no separate *L*-curve analysis was conducted (which would involve several inversions) at each iteration, but the regularization parameter *λ* was fixed to the trace of the matrix **A**^T · **A** divided by 10 (*v*-method) or the trace of **M**^T · **M** divided by 100 (*v&r*-method). The effect of this regularization will be discussed more in detail in Sections 6.3.2 and 6.4.2. Additionally, after each iteration the solution was not re-injected entirely, but multiplied with a feedback value of 0.5. The interest of this value is explained in the following Section 6.2.1.

In both cases, data and model misfit decrease with iterations, indicating that the
<figure><img src="image_64.png" /><figcaption>Figure 6.3: Data misfit $\sum_{j=1}^{rays} (\frac{\delta T}{T})_j$ and model misfit as a function of iterations for the v-method (top) and the v&r-method (bottom). The inversion results for the minimum (red point) of the model misfit (iteration 5 for v-method and iteration 13 for v&r-method) are shown on the right panels.</figcaption></figure>

<table><caption>Table 6.1: Model misfit and corresponding data misfits for results in Figure 6.3.</caption><thead><tr><td>Method</td><td>iteration</td><td>optimum model misfit</td><td>corresponding data misfit</td></tr></thead><tbody><tr><td>v</td><td>5</td><td>~ 1.1 · 10<sup>-5</sup></td><td>~ 0.4 · 10<sup>-5</sup></td></tr><tr><td>v&r</td><td>13</td><td>~ 0.5 · 10<sup>-5</sup></td><td>~ 0.8 · 10<sup>-5</sup></td></tr></tbody></table>

model improves with iterations. The difference between target and solution is less than 30% for v&r at iteration 13 and less than 40% for the v-method at iteration 5 (Roy et al., 2014). The values for data and model misfit at the model misfit minima are summarized in table 6.4. The model misfit in the v&r at its minimum is the half of that for the v-method, on the other hand, the data misfit at iteration 13 for the v&r-method is the double compared to that of the v-method. However, the v&r-method reproduces better the shape of the target model and the smaller model misfit is a sign of the better performance of this method compared to v-method.

After a critical number of iterations, the discrepancy between solution and target model starts growing for both inversions, indicating that the system becomes very unstable with iterations. Nevertheless, the data misfit decreases further, and converges in the case of the v-method. The model misfit increase can be rationalized by two explanations.

First, the unknown endpoints of the ray allow the algorithm to fit the data by strongly deflecting the ray, or by strongly perturbing the electron density to fit the propagation time of the data. Any given observation can be explained by a combination of both negative and positive heterogeneity whose amplitude might
grow indefinitely as further iterations are performed. This problem is similar to that of selecting the regularization $\lambda$, and the following tests can serve to determine adequate values for such parameters, to be used in real inversions with the same data coverage. The second possible explanation is that after a critical number of iterations, the propagation time difference vector that is inverted is in the same order as the numerical noise introduced by discretization of the problem. Hence, the inversion just adds noise to the solution. To explore the reason for the instability of the system, the effect of regularization and feedback were checked with iterations for both methods in the following sections.

### 6.2.1 What is the feedback ?

Initial iteration tests published by Roy et al. (2014) showed that this inverse problem is very unstable and the solution diverges from the target model with iterations (Figure 6.3). One possible explanation is that the re-injected perturbation may strongly perturb rays when a large perturbation is added at the reflection points of the rays, where they are most sensitive to the medium (Chapter 4).

The present inverse problem can be considered as a closed loop control system (Åström and Murray, 2008) (Figure 6.4). The general objective in a control system is to induce a certain output behavior in the desired way by manipulating the input. This can be achieved by using a feedback. In addition, a feedback can stabilize an unstable system and also reduce the sensitivity with regard to external disturbances or to changing parameters in the system itself. For instance, in a car a controller sensing the current speed can maintain the current speed in the presence of disturbances such as hills and headwinds. The term feedback simply means that a certain part of the output is “fed back” as input into the process.

Here, the electron density perturbation $\delta Ne$, obtained from the inversion of the propagation time difference $\delta T$ between data and synthetics, is multiplied by a feedback value $\alpha$ before it is added to the electron density $Ne_{k-1}$. $Ne_{k-1}$ is the sum of previous electron density perturbation, namely

$$Ne_{k-1} = Ne^0 + \alpha \sum_{j=1}^{k-1} \delta Ne_j, \qquad (6.3)$$

with iteration $k$ and $Ne^0$ the electron density from the a priori ionospheric model.

The present problem, is, however, a nonlinear nested closed loop system, because the electron density perturbation $\delta Ne$ depends on the inversion of the propagation time difference $\delta T$, meaning $\delta Ne = \text{inv}(\delta T)$. $\delta T$ itself is a function of the raytracing in $\delta Ne$, meaning $\delta T = \text{TDR}(\delta Ne)$. Here, the feedback $\alpha$ applies to the electron density perturbation $\delta Ne$.
<figure><img src="image_65.png" /><figcaption>Figure 6.4: The inverse problem represented as closed loop control system. The propagation time difference $\delta T$ is inverted to obtain the electron density perturbation $\delta N_e$. $\delta N_e$ is multiplied by a feedback $\alpha$ and added to the electron density $N_{e-k-1}$. Then, rays are traced in the updated model $N_{e_k}$. The propagation time $T^{synth}$ is calculated, which is then compared with the observed propagation time $T^{obs}$. The parameter $\alpha$ is the feedback and $\Delta\tau$ means a delay of an iteration.</figcaption></figure>

## 6.3 Iterating the $v&r$-method

The previous preliminary results showed that the iterative $v&r$-method improves the inversion result with iterations. Data and model misfit decrease and are smaller than in the first iteration. However, the severe ray-path deflection destabilizes the system and after a critical number of iterations the solution diverges from the target model. Therefore, to reduce the effect of severe ray-path deflection in the $v&r$-method and to achieve convergence of the algorithm, the effect of feedback and regularization were analyzed. The results will be discussed more detailed in the following subsections.

### 6.3.1 Effect of the feedback on the iterative $v&r$-method

Different feedback values $\alpha$ (0.1, 0.3, 0.5, and 1) were tested for the inversion using the $v&r$-method. 30 iterations were performed each, except for the feedback of 1 that is already totally dominated by noise at the eleventh iterations and the solution does not agree with the target model in amplitude and shape at all (Figure 6.6). After each iteration, the obtained solution was multiplied by the feedback in order to stabilize the following inversion. In this way, the ionospheric model is slightly updated and not abruptly, so the rays are not as strongly perturbed as if they are traced in the model plus the entire perturbation. This approach should assure convergence.

Figure 6.5 shows the evolution of data and model misfits in the iterative $v&r$-method. Both misfits decrease and increase after a critical number of iterations. The increase depends on the feedback value used and starts later for smaller feedback values. The minima of data and model misfit do not occur after the same number of iterations. In the present case, the model misfit starts increasing after 15 iterations
<figure><img src="image_66.png" /><figcaption>Figure 6.5: Data misfit (left) and model misfit (right) as a function of iterations using different feedback values in the v&r-method. Diamonds mark minima.</figcaption></figure>

($\alpha = 0.1$), 5 iterations ($\alpha = 0.3$) and 3 iterations ($\alpha = 0.5$). The model misfit at its minima are for all four feedbacks nearly the same (Table 6.2 and Figure 6.5), as well as the data misfits at the corresponding iterations, expect the very large data misfit for a feedback of 1 at the first iteration.

Figure 6.5 illustrates that the increase of the model misfit occurs less rapid for the smallest feedback ($\alpha=0.1$) than for all others, indicating some stabilization effect of the feedback on the method, because the model is updated less radically. Nevertheless, the data misfit curve for the feedback of 0.1 suggests that the minimum is not yet reached and further iterations are necessary to reach the same level as for the other feedback values.

<table><caption>Table 6.2: Minimum model misfit and corresponding data misfit for results in Figure 6.5 using different feedback values $\alpha$ in the iterative v&r-method. The L-curve criterion was used to choose the best regularization parameter.</caption><thead><tr><td>Method</td><td>iteration</td><td>optimum model misfit</td><td>corresponding data misfit</td></tr></thead><tbody><tr><td>α = 0.1</td><td>15</td><td>0.63 · 10<sup>-5</sup></td><td>0.86 · 10<sup>-5</sup></td></tr><tr><td>α = 0.3</td><td>5</td><td>0.607 · 10<sup>-5</sup></td><td>0.96 · 10<sup>-5</sup></td></tr><tr><td>α = 0.5</td><td>3</td><td>0.605 · 10<sup>-5</sup></td><td>1 · 10<sup>-5</sup></td></tr><tr><td>α = 1</td><td>1</td><td>0.594 · 10<sup>-5</sup></td><td>2 · 10<sup>-5</sup></td></tr></tbody></table>
As mentioned above, the increasing model misfit can be explained by the trade-off between velocity heterogeneity and ray-path deflection (Section 4.6.2) or by the numerical noise introduced by the discretization of the medium. The data misfit decreases for a certain number of iterations, indicating that the data is fitted by the combination of both negative and positive heterogeneities.

<figure><img src="image_67.png" /><figcaption>Figure 6.6: Inversion results at the model misfit minima (top panel) and the data misfit minima (bottom panel) marked in Figure 6.5 by diamonds. Different feedback values ($\alpha$) were used and the L-curve as regularization was applied in the iterative v&r-method. The bottom line of each panel shows the target model.</figcaption></figure>

Subsequently, it increases again, followed by a large step that is ascribed to
the fact that rays with a larger $\delta T$ than in the first iteration were removed from
the data set for inversion. Since ray end points are not fixed and unknown in the
present problem, the ray-path can be greatly perturbed compared to Earth seismic
tomography.

Figure 6.6 illustrates the inversion results for the minima in the model (top panel) and the data misfit (bottom panel) marked by diamonds in Figure 6.5. All solutions at the model misfit minima display negative anomalies that are do not exist in the target model, but the position of maximum perturbation anomaly (at around 500 km horizontal distance and 200 km altitude) is reproduced. The shape of the perturbation does not change with iterations, only the amplitude of the solutions increases. These negative perturbation anomalies can be explained by the above mentioned nonlinear effects of ray-path deflections, which (endpoints not being fixed), can in principle, result in a faster propagation time even if the velocity perturbation is negative.

In summary, the inversion results are very similar, confirmed by the nearly same
values for the data and model misfit in Table 6.2. The inversion result after iteration
15 with a feedback $\alpha = 0.1$ resembles the same as after the first iteration using
$\alpha = 1$. Therefore, this parameter, did not lead to convergence in the model misfit,
but prevented the model misfit from increasing very fast. Furthermore, it slowed
down the model update and did not lead to significant improvement of the final
result.

The above tests on the feedback show that, for the *v&r*-method, a feedback value
of 0.1 is too small, because more iterations are needed to obtain the same inversion
result then with a larger value. On the other hand, re-injecting the entire solution
after each iteration seems to destabilize the system and leads to a rapid increase of
data and model misfit. In terms of iterations, the above tests showed that with the
*L*-curve, maximal 5 iterations are necessary to obtain a satisfying result.

6.3.2 Effect of the regularization on the iterative v&r-method

As discussed in Section 4.5.2, the present inverse problem is an under-determined
one and the matrix to be inverted ($\mathbf{A}^\mathrm{T}\mathbf{A}$) is sparse due to ray coverage. As a
result, the determinant of $\mathbf{A}^\mathrm{T}\mathbf{A}$ is close to zero and the matrix cannot be inverted.
However, the inverse problem can be solved by regularization.

In the iteration results presented before, the *L*-curve criterion was applied to
select the best regularization parameter. This criterion is rather arbitrary, because
there is no uniform definition of the *x*-axis of a *L*-curve plot, and it depends on the
user's choice. Additionally, the *x* and *y* axis have different dimensions. As a result,
the curvature of this curve is not uniquely defined, and neither is its corner. This
<figure><img src="image_68.png" /><figcaption>Figure 6.7: 50 largest eigenvalues $\lambda$ of matrix $M^T \cdot M$ in descending order and different values of regularization. $\lambda_{best}$ is the regularization parameter selected by the L-curve, trace means the trace of the matrix $M^T \cdot M$ and max Curvature the maximum curvature of the eigenvalue curve (dark blue).</figcaption></figure>

criterion is like a rule of thumb to obtain a regularization parameter.

Due to the limitations of the L-curve selection criterion described above, another tool to select a regularization parameter was tested. The ill-conditioning of the inverse problem manifests itself in a large dynamic of eigenvalues of the matrix. The trace of a matrix is equal to the sum of eigenvalues and allows for a quick estimation of those values that are well determined, which are the largest ones. Based on these results, the regularization can be adjusted (Press et al., 1992, Chapter 18.5).

Figure 6.7 shows the 50 largest eigenvalues of the matrix to invert for the v&r-method at the first iteration, as well as the levels of different regularization parameters compared to them. Clearly, the best regularization parameter selected by the L-curve criterion (red line) is too small compared to the eigenvalues. That means a lot of values in the matrix, corresponding to numerical noise with no useful information, are inverted and not suppressed.

On the other hand, the trace of the matrix $M^T \cdot M$ (green line) as regularization is too large, suppressing all useful information and strongly damping the solution. The value of the trace of the matrix divided by 100 (black line) lies exactly at the corner of the eigenvalue curve that separates the well estimated from the poorly estimated eigenvalues. The maximum curvature of the eigenvalue curve (purple
line) lies at a local maximum and has therefore no physical relevance.

Using a feedback of 0.5, different regularizations (trace divided by 10 and 100 and the maximum curvature of the eigenvalue curve) were tested in 30 iterations. The regularization parameter changes in iterations, because it is determined by the trace of the tomographic matrix that is recalculated in each iteration. The evolution of the data and model misfit is shown in Figure 6.8 and the inversion results for each method at the minimum of the model misfit and the data misfit minima are illustrated in Figure 6.9. The model misfit increases most quickly when the *L*-curve is used in iterations (after iteration 3), for the other methods it increases slowly after iterations 13 and 16 using the trace/100 and the maximum curvature, respectively. The larger the regularization, the later the increase of the model misfit.

The model misfits and corresponding data misfits are summarized in Table 6.3. As for the feedback test described above, the model misfit minima are very similar, as well as the data misfits at the corresponding iterations, except for the iterations with the trace of the matrix divided by 10. Consequently, the inversion results, shown in Figure 6.8 are very similar.

<figure><img src="image_69.png" /><figcaption>Figure 6.8: Data (left) and model misfit (right) evolution with iterations using the iterative *v&r*-method and different approaches to identify a regularization parameter. The feedback value was $\alpha = 0.5$.</figcaption></figure>

The black and magenta curves in the model misfit are similar to that in Figure 6.5 using the *L*-curve and a feedback of 0.1 (green curve). Also data and model
<table><caption>Table 6.3: Minimum model misfit and corresponding data misfit for results in Figure 6.8 using different regularization methods and a feedback of $\alpha = 0.5$ in the iterative v&r-method.</caption><thead><tr><th>Method</th><th>iteration</th><th>optimum model misfit</th><th>corresponding data misfit</th></tr></thead><tbody><tr><td>L-curve</td><td>3</td><td>0.61 · 10<sup>-5</sup></td><td>1 · 10<sup>-5</sup></td></tr><tr><td>trace(M<sup>T</sup> · M)/100</td><td>13</td><td>0.6 · 10<sup>-5</sup></td><td>0.78 · 10<sup>-5</sup></td></tr><tr><td>max curv</td><td>16</td><td>0.61 · 10<sup>-5</sup></td><td>0.76 · 10<sup>-5</sup></td></tr><tr><td>trace(M<sup>T</sup> · M)/10</td><td>30</td><td>0.78 · 10<sup>-5</sup></td><td>1.2 · 10<sup>-5</sup></td></tr></tbody></table>

misfits are similar at there minima (Table 6.2 and 6.3). This indicates that feedback and regularization are linked together. The same inversion result as for the L-curve and a feedback of 0.1 can be obtained with a larger regularization and a larger feedback after nearly the same number of iterations.

Upon iterating with the trace of the matrix divided by 10, the regularization parameter is so large (Figure 6.7) that the solution is strongly damped and the model is not significantly updated after each iteration. That is the reason why the model misfit is still large and decreases at iteration 30. Probably more iterations are required to reach the same level of misfit as for the other regularization methods.

As a consequence, the inversion result in Figure 6.8 has a very small amplitude and does not match the target model. In terms of iterations, depending on the chosen regularization, maximal 16 iterations are necessary to reach the minimum of the model misfit, not considering the iterations of the trace of the matrix divided by 10.
<figure><img src="image_70.png" /><figcaption>Figure 6.9: Inversion results at the model misfit minima (top panel) and the data misfit minima (bottom panel) marked in Figure 6.8 by diamonds. The feedback was set to 0.5. The bottom line of each panel shows the target model.</figcaption></figure>

### 6.3.3 Iterative Damped least squares inversion with constraints

The noise in the damped least squares inversion, visible in particular at a low altitude of 100 km for all regularization methods (Figure 6.9), may strongly affect rays when passing from one cell into another. In order to avoid this and the resulting increasing noise in the solution, a smoothed solution re-injected after each iteration
tion could provide better results in the iterative approach. Thus, the damped least squares inversion with constraints (4.5.3) was applied to the iterative v&r-method. In this method, no L-curve criterion is used to choose a regularization parame- ter. Instead, the factors λ and λ₁ were initialized as described in Section 4.5.3 and are given by equation 4.41. The two other parameters were fixed to α₀ = 2 and α₁ = 14 based on the tests in Section 4.6.4. λ and λ₁ may change slightly with iterations, because the matrices (A and M) evolve with iterations. Nevertheless, the best regularization parameter chosen from the L-curve (λ = 0.63) is very close to α₀ · trace (M^T M) / trace(I) = 0.4. Consequently, the regularization is the same, but the solution is smoothed.

13 iterations were calculated with the constraint damped least squares inversion, where the feedback value was set to 0.5. The evolution of data and model misfits as well as the inversion result for the minima model misfit are compared to the thirteen iterations of the v&r-method using the L-curve and the same feedback value in Figure 6.10. In both cases the data misfit decreases, and the model misfit increases after a critical number of iterations, but much slower with the constraints. Additionally, the model and data misfit reach a smaller minimum with constraints than without the constraints matrices. The data misfit is, nevertheless, in both cases identically at iteration 3 and 4 (the model misfit minima).

Comparing the inversion results obtained for the iterations with and without constraints matrix (Figure 6.10), the effect of the smoothing is clearly visible. All the noise around the solution, in particular at low altitude (100 km) and for larger distances between 600 km disappeared, and the solution is smooth. As compared to the iterations with the L-curve only one additional iterations is necessary to obtain a smoothed solution with a smaller model misfit.

<table><caption>Table 6.4: Minimum model misfit and corresponding data misfit for results in Figure 6.10.</caption><thead><tr><td>Method</td><td>iteration</td><td>optimum model misfit</td><td>corresponding data misfit</td></tr></thead><tbody><tr><td>without constraints</td><td>3</td><td>0.61 · 10<sup>-5</sup></td><td>0.74 · 10<sup>-5</sup></td></tr><tr><td>with constraints</td><td>4</td><td>0.4 · 10<sup>-5</sup></td><td>0.74 · 10<sup>-5</sup></td></tr></tbody></table>
<figure><img src="image_71.png" /><figcaption>Figure 6.10: Iterative damped least squares inversion using the L-curve criterion as regularization vs. iterative damped least squares inversion with constraints. The top panel illustrates the evolution of data and model misfits as a function of iterations for both inversions and the bottom panel the inversion results at data misfit minima (left) and model misfit (left). The feedback was set to 0.5.</figcaption></figure>
## 6.4 Iterating the *v*-method

In the following, the results obtained with the iterative *v*-method will be discussed in detail. The method was applied to the same target model (localized perturbation), same grid and ionospheric conditions. Again, the effect of feedback and regularization were tested on the convergence of the solution. The smoothing was not tested, because the solutions after the first iteration in the *v*-methods are already smooth.

### 6.4.1 Effect of the feedback on iterative *v*-method

As for the *v&r*-method, 30 iterations where calculated for the *v*-method using feedback values of 1, 0.5, 0.3, and 0.1. In each iteration, the *L*-curve was computed again to choose the best regularization parameter. The evolution of data and model misfit is shown in Figure 6.11, as well as the inversion results at the data misfit minima.

For all feedback values, both misfits are decreasing and not increasing again like for the *v&r*-method. The decreasing rate depends, again, on the feedback $\alpha$: the smaller the feedback, the slower the misfits are decreasing. That means with a larger value, the solution is evolving more quickly. This is also visible in the inversion result that illustrate a horizontal movement of the maximum of electron density perturbation to the correct position of the target model with further development of the solution. While the maximum of electron density perturbation is located at around 350 km horizontal distance for $\alpha = 0.1$, it has already moved to approximately 400 km when re-injecting the entire solution after each iteration.

Nevertheless, with the large regularization imposed by the *L*-curve, it will take probably a lot of additional iterations to move the solution to the correct position. Interestingly, the minima of model and data misfit are for all feedback values at the 30. iteration, except for $\alpha = 0.5$, where the model misfit is minima at iteration 24. Maybe this indicates that the real misfit minima have not been reached yet and further iterations are necessary.

Compared to the iterative *v&r*-method described in Section 6.3.1, more iterations are necessary, regardless of which feedback value is used to reach a model misfit minimum. A test (not shown here) illustrates that a minimum is not even reached after 60 iterations for a feedback of 0.1. Nevertheless, 30 iterations for the *v*-method (18 minutes) are calculated faster than 3 iterations for the *v&r*-method (27 minutes).

The values for the model misfit minima and the corresponding data misfits are summarized in Table 6.5. Also here, the model misfit minima are similar for all feedback values, as well as the corresponding data misfits, except the result for the feedback of 0.1 with a larger data misfit. Figure 6.11 leads to the assumption that
further iterations are necessary in this case to reach the same level as for the other feedback values.

<figure><img src="image_72.png" /><figcaption>Figure 6.11: Effect of a feedback value $\alpha$ on the inversion results of the $v$-method with the $L$-curve as regularization. The top panel shows the evolution of data and model misfit, and the bottom panel the inversion results at the data misfits minima.</figcaption></figure>

<table><caption>Table 6.5: Model misfit minimum and corresponding data misfit for results in Figure 6.11 depending on different feedback values $\alpha$ in the iterative $v$-method.</caption><thead><tr><th>Method</th><th>iteration</th><th>optimum model misfit</th><th>corresponding data misfit</th></tr></thead><tbody><tr><td>$\alpha = 0.1$</td><td>30</td><td>$0.89 \cdot 10^{-5}$</td><td>$1.09 \cdot 10^{-5}$</td></tr><tr><td>$\alpha = 0.3$</td><td>30</td><td>$0.86 \cdot 10^{-5}$</td><td>$0.69 \cdot 10^{-5}$</td></tr><tr><td>$\alpha = 0.5$</td><td>24</td><td>$0.88 \cdot 10^{-5}$</td><td>$0.63 \cdot 10^{-5}$</td></tr><tr><td>$\alpha = 1$</td><td>30</td><td>$0.9 \cdot 10^{-5}$</td><td>$0.47 \cdot 10^{-5}$</td></tr></tbody></table>
<figure><img src="image_73.png" /><figcaption>Figure 6.12: Inversion results at the model misfit minima (top panel) and the data misfit minima (bottom panel) marked in Figure 6.5 by diamonds. Different feedback values ($\alpha$) were used and the L-curve as regularization was applied in the iterative v-method. The bottom line of each panel shows the target model.</figcaption></figure>
6.4.2 Effect of regularization on the iterative *v*-method

The above described tests showed that the *L*-curve is not a good tool to choose
a regularization parameter for the iterative *v*-method neither, because it strongly
damps the solution. Additonally, after 30 iterations the model misfit minimum
is still not reached, regardless which feedback value was used. Therefore, other
regularizations were tested as before for the iterative *v&r*-method in Section 6.3.1.

Figure 6.13 shows this regularization parameters in the first iteration compared
to the 50 largest eigenvalues of the matrix **A**^T . **A**. Clearly, the value selected by

<figure><img src="image_74.png" /><figcaption>Figure 6.13: 50 largest eigenvalues λ of matrix A^T . A in descending order and different values of regularization. λbest is the regularization parameter selected by the L-curve, trace means the trace of the matrix A^T . A and max Curvature the maximum curvature of the eigenvalue curve.</figcaption></figure>

the *L*-curve (red line) and the trace of the matrix (green line) are too large for
regularization. They are larger than the largest eigenvalue of the matrix, so they
strongly damped the solution and no useful information is added to the solution,
just white noise. This is the reason, why for the *L*-curve criterion, the feedback
is not destabilizing the system in iterations, because only very few information is
added.

On the other hand, the trace of the matrix divided by 100 (black line) is very
small and may introduce much noise into the inversion, but gives more freedom to
the solution to change. Here, a careful choice of the feedback value α is essential.
The smaller the regularization parameter, the smaller α needs to be chosen in order
to change slightly the solution and avoid the destabilization of the system with iterations. Test with iterations for different feedback values for the trace of the matrix divided by 10, 20, and 100 as regularization parameter proved that: The system becomes unstable for the smallest regularization parameter (trace divided by 100) already for $\alpha$ larger than 0.1, while for the trace of the matrix divided by 10 the iterations are stable for feedback of up to 0.5.

Figure 6.14 compares the evolution of data and model misfits for four regularization methods (L-curve, trace divided by 10, 20, and 100). In all methods, a feedback value of 0.1 was used. Since the L-curve criterion chooses the largest regularization parameter (Figure 6.13), the model misfit decreases very slowly with iterations and has not yet reached a minimum at iteration 30, while for the other regularization methods a minimum is reached because the solution evolves more quickly.

<figure><img src="image_75.png" /><figcaption>Figure 6.14: Data and model misfits evolution with iterations for four regularization pa- rameter applied to the *v*-method and a feedback value of 0.1.</figcaption></figure>

Model misfit minima and corresponding data misfits are summarized in Table 6.6. Model misfits for the regularizations without *L*-curve are similar, but the corre- sponding data misfits are quite different, in particular the small data misfit for the trace divided by 100 is remarkable.

Figure 6.15 illustrates the inversion results at data and model misfit minima. The results obtained after 30 iterations with the *L*-curve is less noisy as compared to the other inversion results after the same number of iterations. This is, as mentioned before, due to the large regularization, larger than the largest eigenvalue of the
<table><caption>Table 6.6: Data and model misfit for results in Figure 6.15 depending on different regularization methods and a feedback of $\alpha = 0.1$ in the iterative v-method.</caption><thead><tr><th>Method</th><th>iteration</th><th>optimum model misfit</th><th>corresponding data misfit</th></tr></thead><tbody><tr><td>L-curve</td><td>30</td><td>0.89 · 10<sup>-5</sup></td><td>1.09 · 10<sup>-5</sup></td></tr><tr><td>trace(M<sup>T</sup> · M)/100</td><td>13</td><td>1.1 · 10<sup>-5</sup></td><td>0.05 · 10<sup>-5</sup></td></tr><tr><td>trace(M<sup>T</sup> · M)/20</td><td>2</td><td>1 · 10<sup>-5</sup></td><td>1.3 · 10<sup>-5</sup></td></tr><tr><td>trace(M<sup>T</sup> · M)/10</td><td>5</td><td>0.92 · 10<sup>-5</sup></td><td>0.88 · 10<sup>-5</sup></td></tr></tbody></table>

matrix that suppresses all useful information in the matrix to invert and just adds white noise. Additionally, the solution is not on the right position of the target model. In contrast to that the solution for the three other methods is at the correct position, and, the lower the regularization, the higher the noise level in the solution. This is clearly confirmed by the model misfit curve in Figure 6.14.

Again, the solution moves slightly horizontally with iterations, which is clearly visible when comparing the inversion results at model misfit minima (Figure 6.15, bottom panel). While for the L-curve regularization the solution is still mislocated at iteration 30, it is at the right position with smaller regularization (trace divided by 100) already after iteration 13. After that, the algorithm only adds noise to the solution that is the reason why the model misfit starts increasing again (Figure 6.14). The trace of the matrix divided by 10 imposes a larger regularization to the solution, so the maximum of electron density perturbation has only slightly moved horizontally at iteration 5, but the difference to the 30. iteration with the L-curve is clearly visible.
<figure><img src="image_76.png" /><figcaption>Figure 6.15: Effect of regularization on the iterative *v*-method. Inversion results with the *v*-method at model misfit (bottom panel) and data misfit (top panel) minima in Figure 6.14. For all four methods, a feedback value of 0.1 was applied.</figcaption></figure>

## 6.5 Summary & Conclusions

In this chapter, both, the *v*-method and the *v&r*-method were iterated in order to improve the inversion results. The preliminarily results published by Roy *et al.* (2014) demonstrated that this is archived by the decreasing model and data misfit. At their minima, they are smaller compared to the first iteration. Additionally, they illustrated the better performance of the iterative *v&r*-method in terms of
model misfit (30% for *v&r* and 40% for the *v*-method). However, while the data misfit decreases for both methods, the model misfit increases after a critical number of iterations. This indicates that the system becomes unstable and the solution diverges from the target model. This instability can be explained by the numerical noise, introduced by the discretization of the model by blocs or the non-linearity of the ray-path deflection.

To deal with this problem and to reduce the instability in the iterative approach, the effect feedback and regularization were tested for both methods. The first parameter controls the amplitude, the second damps the re-injected solution. For both methods, feedback values of 0.1, 0.3, 0.5, and 1 were tested, as well as another regularization based on the eigenvalues of the matrices in addition to the *L*-curve criterion.

In the iterative *v&r*-method with the *L*-curve as regularization (Section 6.3.1), the solutions have similar model misfit minima and corresponding data misfit, but, depending on the used feedback value, more iterations are needed to obtain more or less the same inversion result. This means, the feedback slows down the convergence, but the test also showed a more slowly increase in the model misfit. Contrary, in the *v*-method, a minimum in the model misfit is not yet reached after 30 iterations and probably more iterations are necessary. In terms of execution time this means, regardless which feedback is used with the *L*-curve as regularization, the *v&r*-method will faster reach a model misfit minimum.

In a second test, the regularization method was tested in the iterative approach. The analysis of the eigenvalues of the matrices to invert compared to the regularization revealed that for the *v&r*-method, the *L*-curve returns a value that is too small and introduces noise in the inversion, but for the *v*-method the chosen value by the *L*-curve is to large, suppressing all useful information in the matrix. Regardless with regularization method is used in the *v* and the *v&r*-method, model and data misfit minima are similar and are reached later, the larger the regularization parameter. While the maximum of electron density perturbation is already at the right position in the first iteration of the *v&r*-method, it is mislocated in the *v*-method, but moves horizontally with further iterations, improving the result.

In conclusion, the *v&r*-method is computationally more expensive and numerically more complex, but, more rapid, because it reaches model misfit minima generally faster than the *v*-method. This is probably linked to the fact that the *v&r*-method takes into account the ray-path deflection starting from the first iteration, the *v*-method considers ray-path deflection from the second iteration on, because rays are traced in an updated, new ionospheric model. Contrary to the *v*-method, the *v&r*-method needs a smoothing in order to reduce the noise in the solution at low altitude.
# Chapter 7

## Real data inversion

In this chapter, the developed tomography method was applied to real data measured by the OTH radar Nostradamus. The synthetic data were replaced by real data. Section 7.1 clearly describe the preprocessing used to prepare the data to the inversion. Subsequently, the real data were inverted and vertical profiles of electron density calculated (Section 7.2). At the end of this chapter, some perspectives for further applications and improvements are discussed.

### 7.1 Data Preprocessing

The operation mode of the OTH radar Nostradamus is explained in Section 3.4. The radar transmits a chirp signal (linear frequency modulation) with a transmission time of 3 ms with a recurrence of 30 ms. After digitalization of the received signal, cross-correlation with the emitted signal, Doppler processing and digital beamforming for different elevation angles, a curve of elevation angle as a function of group path or group propagation time is obtained (e.g., Figure 7.1).

Figure 7.1 shows real data obtained on March 14, 2006 at 1855 UT measured by Nostradamus. Clearly, the data contain much scatter and reflections from different ionospheric layers are visible. It is clear that these data need some preprocessing before they can be inverted. The orange and red data points at elevation angles larger than 50° are probably multi-path data, meaning that the signal is reflected twice or more often at the ground. In the ray tracing methods, only rays with a single reflection are generally considered (Figure 7.6).

Some preprocessing is necessary to automatically remove multi-path data, sidelobes and noise peaks from the dataset. A good tool to accomplish this is the virtual height of reflection $h_v$ that can be easily calculated from the ground distance $d_g$ and
<figure><img src="image_77.png" /><figcaption>Figure 7.1: Real data (group propagation time as function of elevation angle) obtained from the OTH radar Nostradamus on March 14, 2006 at 1855 UT in azimuth 202°. The frequency is color coded.</figcaption></figure>

the elevation angle $\phi$, using the theorem of Martyn (Martyn, 1935)

$$h_v = \sqrt{\left(\frac{d_g}{2}\right)^2 + R_0^2 + 2R_0 \cdot \frac{d_g}{2} \cdot \sin\phi - R_0}, \quad (7.1)$$

with $R_0$ being the Earth radius. This theorem is a simple approximation for the propagation of radio waves in the ionosphere, allowing pre-filtering of points corresponding to the propagation in the same ionospheric layer (Figure 7.2). Multi-path rays normally have virtual heights of reflections that are twice the virtual height of single-path rays and that are often larger than 500 km. Consequently, the virtual height is a useful tool to exclude these data from the dataset. Even in the case of wrong virtual heights (because they are just an approximation), reflections corresponding to the same ionospheric layer must have approximately same virtual heights.

Figure 7.3 shows the relationship between virtual height of reflection, elevation angle, and group distance for rays with a frequency of 6 MHz and elevation angles between 10° and 20°, traced in the a priori ionospheric model NeQuick, generated for October at 12 UT. This epoch was chosen because it allows to see the reflection in different ionospheric layers, while for the data shown in Figure 7.1, only reflections in the F-layer are visible, and not in the E-layer. One distinguishes between high-
<figure><img src="image_78.png" /><figcaption>Figure 7.2: Explanation of Martyn’s Theorem. The path of the signal (red) is approximated by the green colored path. The virtual height is always larger than the true height of reflection.</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>Figure 7.3: Relationship between virtual height of reflection, group distance and elevation angle for rays with a frequency of 6 MHz and elevation angles between 10°-20°, traced in the *a priori* ionospheric model NeQuick, generated for October at 12 UT. yellow: high-angle ray. Red: low-angle ray.</figcaption></figure>
and low-angle rays. For rays contributing to the low-angle ray (red in Figure 7.3), the virtual height of reflection varies little, but in the high-angle ray (yellow), little variation in the elevation angle of the ray increases significantly the group distance, and consequently the virtual height of reflection (equation 7.1).

<figure><img src="image_80.png" /><figcaption>Figure 7.4: Virtual height of reflection as a function of signal frequency. The color denotes the number of ionospheric layers given in the preprocessing of the radar data to each data point (left). Data points corresponding to the same ionospheric layer are arranged in clusters (right).</figcaption></figure>

The virtual height of reflection is calculated already during the preprocessing of the radar data. It is included in the data set for each azimuth and one frequency, and the data points are assigned to ionospheric layers based on histogram analysis. Layers are numbered from 1 to 5. The number of layers does not correspond to the real number of ionospheric layers for the day of measurement (5 ionospheric layers do not exist), and is not comparable between frequencies (Figure 7.4, left) because it varies with frequency, and consequently the layer numbering too.

In a first step of real data analysis, a k-means cluster analysis (Forgy, 1965) was applied to the virtual height of reflection, in order to assign automatically points with the same virtual height of reflection to the same ionospheric layer. The aim of this method is to cluster the data in *k* portions, so that the sum of the squared difference between the cluster barycenters is minimal. For the present data example, Figure 7.4 shows the virtual height of reflection as function of frequency before and after clustering for the real data. Data points with the same virtual height have now
<figure><img src="image_81.png" /><figcaption>Figure 7.5: Raw data from March 14, 2006 in azimuth 202° (top) and assigned to different ionospheric layers using a clustering algorithm (bottom).</figcaption></figure>

been clustered together. In this plot, reflections from the F-layer are clearly visible (at 200 km), but not from the E-layer, because, this layer has already disappeared
at night (1855 UT). The cyan and purple colored points having twice the virtual
height of the black and red points may be multi-path reflections from the F-layer.

<figure><img src="image_82.png" /><figcaption>Figure 7.6: Group propagation time measured by the radar March 14, 2006 in azimuth 202° compared to the propagation time obtained by ray tracing (dots) in the *a priori* ionospheric model (top) and corresponding rays (bottom).</figcaption></figure>

Based on this analysis, the real data can be color coded as ionospheric layers, in
order to see which data should be excluded from the dataset (Figure 7.5). The data
has been assigned to 4 different layers, where the data with larger propagation time
that probably represents multi-path data, has been assigned to layers 3 and 4. In
this way, the multi-path data can be excluded from the data set.

With the signal frequencies and elevation angles Nostradamus takes measure-
ments, corresponding synthetic data can be computed and a propagation time differ-
ence between data and synthetics can be calculated. Figure 7.6 in the top compares
the group propagation time as a function of elevation for real data and synthetics.
The ray tracing was performed in the *a priori* ionospheric model NeQuick generated
for the epoch and solar flux of the data. The corresponding rays are shown in Figure
7.6 in the bottom plot.

All rays are reflected between 200 and 250 km altitude, so in the F region and not
in the E region that already has disappeared. The ray tracing does not reproduce
the multi-path data at elevation angles larger than 50°. Additionally, the radar is
emitting a beam in a certain azimuth with a specific frequency and elevation angle.
This signal is backscattered at several points, therefore for one frequency and one
elevation angle, multiple propagation time measurements are possible. Contrary,
the ray tracing only returns one propagation time. In general, Figure 7.6 illustrates
that the ray tracing is overestimating the data in terms of propagation time.

<figure><img src="image_83.png" /><figcaption>Figure 7.7: Explanation of path effects in the ionosphere. The signal propagates from the radar to the ground along a path in the E-layer (red), but along a path in the F-layer on its way back to the radar (black).</figcaption></figure>

This might been explained by the fact that the signal does not propagate along
the same path from the radar to the ground, and back to the radar (Figure 7.7).
While propagating in the E-layer in the forward run, the signal might propagate in the F-layer when returning to the radar. This effect is not taking into account in the ray tracing in the *a priori* model.

## 7.2 Inversion results

### 7.2.1 Spring data set

The tomography method was tested on real data using the data set of march 14, 2006 at 1855 UT. It contains measurements in eight azimuths of 22°, 67°, 112°, 157°, 202°, 247°, and 337°. All data were preprocessed in the described way before inversion. Then, for the remaining data ray tracing in the *a priori* ionospheric model NeQuick given for the appropriate day, time and solar flux was performed. The solar flux data is daily measured at local noon by the radio telescope in Penticon, Canada at 2800 MHz (10.7 cm wavelength) and the data is freely available on the NOAA FTP server. For March 14, 2006, the solar flux is 73.5 SFU.

Figure 7.8 shows the comparison between real data and synthetics used to compute the vector $\delta\mathbf{T} = T^{\text{real}} - T^{\text{synth}}$ of travel time perturbations for the azimuths 67°, 157°, 247° and 337°. The difference in the group propagation time is between +20 % and -40 %. After ray tracing for the calculation of the synthetic data, the

<figure><img src="image_84.png" /><figcaption>Figure 7.8: Difference in percent between the propagation time obtained by ray tracing in an *a priori* ionospheric model and real data of March 14, 2006 at 1855 UT for four different azimuths. The frequency is color coded.</figcaption></figure>

tomographic matrices for the *v* and the *v&r*-method were calculated, and the inver-
sions using the damped least squares inversion method were performed as described
in Section 4.5.

The inversion results of the real data after the first iteration in one selected azimuth (247°) for three different values of regularization for the *v*- and the *v&r*-method are shown in Figure 7.9. The data set consists of 428 rays. Both methods find electron density perturbations located between 200 and 300 km altitude and stretching up to 1500 km distance. It is not surprising that the inversion detects only a perturbation at that altitude, because at the time of measurement (1855 UT) the E-layer has nearly disappeared and the EM signals are reflected mainly in the F-layer. However, while the *v&r*-method reconstructs positive and negative

<figure><img src="image_85.png" /><figcaption>Figure 7.9: Real data inversion for data of March 14, 2006 at 1855 UT in azimuth 247°. In the *v&r*-method (left) and the *v*-method (right), different values of regularization λ were used. The grid size is 50 × 20 km in distance and altitude. The results for the best regularization parameter are shown in the bottom left and the middle right plot, respectively.</figcaption></figure>

perturbations, the *v*-method finds an entire positive and unrealistic perturbation of large amplitude. For the real data also applies that the regularization parameter λ does not change the inversion results for both methods, but larger regularization leads to a more damped solution, reducing the amplitude of the solution. For the largest regularization λ = 0.1, the perturbation is mainly concentrated above 200 km altitude, while for smaller regularization, also perturbations at 100 km are visible.

The best regularization parameters for both methods (λ = 0.1 and λ =0.02) were
again calculated from the *L*-curve, where the λ range obtained for the checkerboard
λ = [10⁻³, 10³] was used. Then, for any distance from the radar, a vertical profile
<figure><img src="image_86.png" /><figcaption>Figure 7.10: Vertical profile of electron density in 500 km distance from the radar calculated from real data inversion of OTH radar Nostradamus. The data was collected March 14, 2006 in azimuth 247°. The electron density perturbations $\delta Ne/Ne$ obtained by inversion are shown in the two subplots in the right. The electron density perturbation of one point in the left plot corresponds to the block starting from the altitude of the point in the right plot.</figcaption></figure>

of electron density can be calculated from an inversion result like in Figure 7.9 in order to correct a vertical profile given from an *a priori* model.

Figure 7.10 illustrates the vertical profiles of electron density in 500 km distance from the radar in the *a priori* ionospheric model NeQuick as well as the profiles corrected with the perturbation obtained from the *v&r* and the *v*-method. The vertical profile for NeQuick was calculated for the corresponding epoch and solar flux, in azimuthal direction of 247° in 500 km distance and was parametrized like the grid with a value every 20 km. It was then corrected with the electron density perturbations found by the inversion marked by vertical dashed lines in Figure 7.9. The resulting vertical profiles show a sensitivity to variations located between 180 km and 260 km altitude (Roy et al., 2014).

In both methods, there are no corrections of the NeQuick electron density profile for altitude lower than 160 km, because the inversion only finds perturbations in the F-layer and nothing for altitudes above 260 km. The electron density perturbations for the *v&r*-method are in average around 7%, for the *v*-method around 15% with a large perturbation between 220-240 km altitude.

Figure 7.11 illustrates the same three vertical profiles of electron density at 1000 km
<figure><img src="image_87.png" /><figcaption>Figure 7.11: Vertical profile of electron density in 1000 km distance from the radar calculated from real data inversion of OTH radar Nostradamus. The data was collected march 14, 2006 in azimuth 247°. The electron density perturbations $\delta Ne/Ne$ obtained by inversion are shown in the two subplots in the right. The electron density perturbation of one point in the left plot corresponds to the block starting from the altitude of the point in the right plot.</figcaption></figure>

distance from the radar (Roy et al., 2014). Here, both methods find entirely positive perturbations of electron density in the order of 40% in the altitude range 180-260 km. There is no much difference between the results of the two inversion methods ($v$ and $v&r$). Nevertheless, the vertical profiles show the strong variability of the ionosphere between two local positions (500 km and 1000 km).

In order to show the potential applications of the OTH radar Nostradamus in ionospheric tomography, the data for azimuths 67°, 157° and 337° were inverted, too for the same day and time. The inversion was performed with the damped least squares method using the $v&r$-method, and for each azimuth the best regularization parameter was chosen from the $L$-curve in the range $\lambda = [10^{-3}, 10^3]$ that was found for the checkerboard test as best range. The results are plotted in Figure 7.12 in the corresponding azimuthal direction over Europe. The positive and negative electron density perturbations are visible in all azimuthal directions in the altitude range 200 and 300 km and in the distance range up to 1500 km. As expected, the electron density perturbation is largest in the azimuth directions 247° and 337°, because of the position of the Sun at the time of measurement (1855 UT). This figure shows the potential of the developed ionospheric tomography method. With this method,
it is possible to calculate a vertical profile of electron density in any desired distance
in any azimuth direction where the radar took measurements. In the intermediate
azimuth, a vertical profile can be calculated by interpolation.

<figure><img src="image_88.png" /><figcaption>Figure 7.12: Electron density perturbation obtained by real data inversion. The data as collected at OTH radar Nostradamus March 14, 2006, at 1855 UT in four azimuth directions, nominally 67°, 157°, 247° and 337°. Source: Roy et al. (2014)</figcaption></figure>

Further perspectives, based on the potential of this method, for instance the
combination with GPS tomography are discussed in Section 7.3.2.

7.2.2 Winter data set

The tomography method was applied to another data set recorded in winter (De-
cember) 2008 after the above described preprocessing. In that year, the solar flux
is still small ~ 60-70 SFU, because in 2008 the Sun is undergoing a solar minimum.
Consequently, ionization in the ionosphere is smaller and less rays are reflected.
This time, a data set for noon was chosen, because ionization is highest at noon,
and it might be possible to see reflections from the E and the F-layer.

In Figure 7.13 the measured group propagation time in azimuth 67° and the
propagation time obtained by ray tracing in the *a priori* model are compared. Also
here, as for the spring data set, there a very few reflections from the E layer in
the synthetic data. This is probably due to the fact that the E layer is generally
not well estimated in ionospheric models. The propagation time measured by OTH
radar Nostradamus is larger than that from the ray tracing, probably due to the
path effects mentioned above.

Figure 7.14 shows the inversion results using the data measured December 13,
<figure><img src="image_89.png" /><figcaption>Figure 7.13: Group propagation time measured by OTH radar Nostradamus December 13, 2008 in azimuth 67° compared to the propagation time obtained by ray tracing (dots) in the *a priori* ionospheric model.</figcaption></figure>

2008 at noon with different regularization parameters $\lambda$. The inversion results with the best regularization parameter chosen from the *L*-curve are $\lambda = 0.1$ for *v&r* and $\lambda = 0.05$ for *v*. The perturbations are located between 100 and 300 km altitude and are in the order of 20%. They appear in distance range up to 1500 km for the *v&r*-method, while the *v*-method seems to be strongly damped. Additionally, such large scale entire negative perturbations are hardly realistic, indicating the better performance of the *v&r*-method.

In the *v*-method inversion results, perturbations are only visible up to 600 km. To go a step further, vertical profiles of electron density were calculated for 200 and 500 km distance from the radar in 67° azimuth. These profiles are shown in Figure 7.15 and 7.16. At 200 km distance from the radar, the perturbations found by the two inversion between 100 to 180 km altitude) are controversial. While the *v&r*-method finds positive perturbations, those for the *v*-method inversion are negative. Above 180 km altitude, both methods agree and reconstruct negative anomalies. In higher altitudes, no perturbations are found any longer because of the lack of reflected rays, and the vertical profiles coincide with the NeQuick model.

At 500 km distance from the radar (Figure 7.16), both inversion methods agree up to 160 km altitude in finding negative perturbations (Figure 7.16), where those reconstructed by the *v&r*-method are larger (20%) than for the *v*-method (<10%).
<figure><img src="image_90.png" /><figcaption>Figure 7.14: Real data inversion for data of December 13, 2008 at 1154 UT in azimuth 67°. In the v&r-method (left) and the v-method (right), different values of regularization $\lambda$ were used. The grid size is 25 × 50 km in distance and altitude, respectively.</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>Figure 7.15: Vertical profiles of electron density at 200 km distance from the radar using data from December 13, 2008 measured in 67° azimuth. The electron density perturbations are obtained from the inversion in Figure 7.14 for the best regularization parameter. The electron density perturbation of one point in the left plot corresponds to the block starting from the altitude of the point in the right plot.</figcaption></figure>
<figure><img src="image_92.png" /><figcaption>Figure 7.16: Vertical profiles of electron density at 500 km distance from the radar using data from December 13, 2008 measured in 67° azimuth. The electron density perturbation of one point in the left plot corresponds to the block starting from the altitude of the point in the right plot.</figcaption></figure>

Above 220 km altitude, no electron density perturbations are found, and the
vertical profiles for *v* and *v&r*-method fit the profile from the *a priori* model.

7.3 Perspectives

7.3.1 Comparison with ionosonde data

As showed in the sections above, this ionospheric tomography method allows to
calculate vertical profiles of electron density at any desired distance within reach of
the radar. Nevertheless, they need to be validated by other instruments that provide
vertical profiles, for instance ionosondes. A couple of them supply profiles of electron
density above the station after ionogram inversion and they can be downloaded
directly from the European DIgital upper Atmosphere Server (DIAS) homepage
page. Other stations provide only the 3 ionospheric parameters f_c, h_m, and y_m,
but based on them, a vertical profile can be calculated using for instance the Multi-
Quasi-Parabolic (MQP) modelisation (Croft and Hoogasian, 1968).

However, for the choice of an ionosonde for comparison with the inversion results
from Nostradamus, the distance of the ionosonde from the radar and the azimuth
need to be considered. Nostradamus data was collected in eight different azimuths,
starting with 22° and a step of 45°. This is the normal operation mode of this radar. The radar range is, depending on the frequency, in the order of 800 to 1600 km. Consequently, ionosondes that can be considered for comparison should be sited in one of the azimuths Nostradamus took measurements and between the half distance range that is 400 to 800 km. Data from Nostradamus is available between 2006 and 2012. Table 7.1 lists some ionosondes in Europe, as well as their distance to Nostradamus and their azimuth. The ionosondes that fulfill the condition in distance and azimuth mentioned above, are highlighted in Table 7.1. These are the stations in Chilton, Dourbes, Fairford, Pruhonice, and Ebre.

<table><caption>Table 7.1: Ionosonde stations in Europe, distance and azimuth to OTH radar Nostradamus. Blue highlighted stations can be considered for comparison with inversion results from Nostradamus. Source: http://car.uml.edu/common/DIDBFastStationList</caption><thead><tr><th>Ionosonde<br>Station</th><th>latitutde<br>[°]</th><th>longitude<br>[°]</th><th>Azimuth<br>[°]</th><th>distance<br>[km]</th></tr></thead><tbody><tr><td>Juliusruh, Germany</td><td>54.60</td><td>13.40</td><td>47.36</td><td>1075</td></tr><tr><td>Athens, Greece</td><td>38.00</td><td>23.50</td><td>114.97</td><td>2151</td></tr><tr><td>Chilton, U.K.</td><td>51.50</td><td>359.40</td><td>339.96</td><td>339.5</td></tr><tr><td>Dourbes, Belgium</td><td>50.10</td><td>4.6</td><td>56.20</td><td>301.82</td></tr><tr><td>El Arenosillo, Spain</td><td>37.1</td><td>353.30</td><td>209.03</td><td>1429</td></tr><tr><td>Fairford, U.K</td><td>51.70</td><td>358.50</td><td>332.59</td><td>386</td></tr><tr><td>Pruhonice, Czech Republic</td><td>50.00</td><td>14.60</td><td>76.14</td><td>989.47</td></tr><tr><td>Rome, Italy</td><td>41.8</td><td>12.50</td><td>126.21</td><td>1171</td></tr><tr><td>Ebre, Spain</td><td>40.8</td><td>0.50</td><td>183.2</td><td>872.65</td></tr><tr><td>San Vito, Italy</td><td>40.60</td><td>17.80</td><td>117.9</td><td>1591</td></tr></tbody></table>

The other ionosondes are too far away or lie in an azimuth that is too different to those of Nostradamus. In a preliminary test, the vertical profiles of electron density obtained by the ionospheric tomography with Nostradamus data were compared to the vertical profile of the ionosonde station Chilton. Chilton, sited at azimuth 340° and in 340 km distance (Table 7.1) can be compared with measurements of Nostradamus taken in 337° azimuth. The ionogram of Chilton at December 19, 2008 at noon is illustrated in Figure 7.17 and shows that there are only a E and a F2 layer visible. Based on the ionospheric parameters listed in the left column of the ionogram, the vertical profile (black in Figure 7.17) can be reconstructed using MQP modeling.

The inversion results of the Nostradamus data using the *v* and the *v&r*-method are shown in Figure 7.18. For both, the *L*-curve was used to choose the best regularization parameter $\lambda$. The electron density perturbation is located between 100 and 300 km altitude and are in the order of 50% for the *v&r* method and very small for the *v*-method (5%). This can be explained by the large regularization parameter
<figure><img src="image_93.png" /><figcaption>Figure 7.17: Ionogram measured at Chilton at December 19, 2008 at noon and inverted vertical profile of plasma frequency over the station (black line). Source: http://car.uml.edu/common/DIDBDayStationStatistic?ursiCode=RL052&year=2008&month=12&day=19</figcaption></figure>

$\lambda$ that strongly damps the solution. This is the reason why in the *v*-method, all electron density perturbations are located at 200 km altitude.

For comparison with the vertical profile of Chilton, first, a vertical profile for December 19, 2008 in the *a priori* ionospheric model NeQuick at the position of Chilton (Table 7.1) was computed. This was subsequently corrected by the perturbations found by the inversion with the *v*- and the *v&r*-method at the corresponding distance between 325 and 350 km marked by dotted lines in Figure 7.18.

The comparison of the vertical profiles of electron density is shown in Figure 7.19. The profile from the *a priori* model NeQuick is shifted upwards in altitude compared to the profile from Chilton. Additionally, the electron density from the model NeQuick is generally smaller than that from the vertical sounder. In the selected distance range (325 and 350 km), the *v&r*-method reconstructs positive and negative electron density perturbations that modify the *a priori* ionospheric model to even smaller perturbations at altitude 120 to 160 km (E-layer). Consequently, the corrected profile moves away from the Chilton profile. Contrary, above that altitude
<figure><img src="image_94.png" /><figcaption>Figure 7.18: Inversion results using real data measured by OTH radar Nostradamus December 19, 2008 at noon in 337° azimuth. *v&r*-method (top) and *v*-method (bottom). The perturbations found at distances marked by dashed lines were used to correct a vertical profile of electron density from the *a priori* ionospheric model and are shown in Figure 7.19.</figcaption></figure>

the profile corrected by the perturbation found by the *v&r*-method is closer to the profile at Chilton.

However, the profile of Chilton is not a direct measurement and obtained by inversion of the measured ionogram. The used inversion method is not known in details. Since the inversion result found by the *v*-method is strongly damped and the perturbations very small (5%), the profile corrected by the *v*-method coincides with the profile in the *a priori* model NeQuick. As the synthetic tests showed (Chapter ??), generally more iterations are necessary with the *v*-method to obtain a satisfying solution. This has to be investigated further and was not yet tested.
<figure><img src="image_95.png" /><figcaption>Figure 7.19: Comparison of the vertical profile of electron density obtained at the ionosonde station Chilton (blue) December 19, 2008 with the one obtained by ionospheric tomography using the *v*- and the *v&r*-method.</figcaption></figure>

### 7.3.2 Combination with GPS

As discussed in Chapter 4, ionospheric tomography based on GPS is limited to the region of maximum ionization in the ionosphere (F2-layer), because of the high signal frequency of the GPS. Contrary, ionospheric tomography with OTH radar allows a good reconstruction of the ionospheric plasma up to 300 km altitude within the radar range, as the tests with synthetics and real data showed. In order to have a complete reconstruction of ionospheric plasma over the entire altitude scale, GPS and OTH radar should be combined in one inversion.

A first synthetic test was conducted to test the feasibility of such an approach. Virtual GPS stations were put on the ground along a line in 89° azimuth and a satellite is assumed to move at 20 200 km altitude in that direction. Then, an electron density perturbation $\delta Ne$ in the form of a checkerboard test of 5% was added to the *a priori* model, and the vector with differential TEC for the GPS inversion and the vector of propagation time difference $\delta T/T$ for the OTH radar were calculated. The radar data was simulated for an ideal case by tracing 1071 rays with elevation angles between 10° and 60° and in a frequency range of 6–16 MHz. The differential TEC,
$dTEC$ is given by

$$dTEC = \int \delta Ne \cdot ds, \qquad (7.2)$$

with $ds$ being the ray-path between a satellite and a ground station. The problem that has then to be solved is given by equation 7.3. The two matrices for the radar and the GPS form one matrix that has to be inverted.

$$\begin{pmatrix} \frac{\delta T}{T} & 1 \\ \frac{\delta T}{T} & 2 \\ \vdots \\ \frac{\delta T}{T} & i \\ \hline dTEC_1 \\ dTEC_2 \\ \vdots \\ dTEC_k \end{pmatrix} = \begin{pmatrix} A_{11}^{OTH} & A_{12}^{OTH} & \dots & A_{1j}^{OTH} \\ A_{21}^{OTH} & A_{22}^{OTH} & \dots & A_{2j}^{OTH} \\ \vdots & \vdots & \vdots & \vdots \\ A_{i1}^{OTH} & A_{i2}^{OTH} & \dots & A_{ij}^{OTH} \\ \hline A_{11}^{GPS} & A_{12}^{GPS} & \dots & A_{1i}^{GPS} \\ A_{21}^{GPS} & A_{22}^{GPS} & \dots & A_{2i}^{GPS} \\ \vdots & \vdots & \vdots & \vdots \\ A_{k1}^{GPS} & A_{k2}^{GPS} & \dots & A_{kj}^{GPS} \end{pmatrix} \cdot \begin{pmatrix} \frac{\delta Ne}{Ne} & 1 \\ \frac{\delta Ne}{Ne} & 2 \\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \frac{\delta Ne}{Ne} & j \end{pmatrix} \quad (7.3)$$

where $A^{OTH}$ is the matrix $A + A'$ for the v&r-method (Section 4.5.4). The matrix for the GPS is given by

$$A_{ij}^{GPS} = Ne_i^0 \cdot ds_{ij}, \qquad (7.4)$$

with $Ne^0$ being the electron density from the a priori ionospheric model.

The two data sets have very different values, the values in the GPS $dTEC$ vector are larger (~$10^{15}$) than in the OTH propagation time difference vector (~0.25). Additionally, due to the chosen GPS station configuration, the amount of $dTEC$ measurements is double compared to the radar data. In order to do a joint inversion that considers both data sets and is not masking one of them, the two datasets need to be weighted, so that

$$|d_{OTH}| \sim \text{weight}^2 \cdot |d_{GPS}|. \qquad (7.5)$$

$|d_{OTH}|$ is the $L_2$-norm of the vector with travel time differences $\frac{\delta T}{T}$ and $|d_{GPS}|$ is the norm of the differential TEC values. The weight is applied to the latter and the matrix for the GPS.

Figure 7.20 illustrates the joint inversion result for the simulated checkerboard perturbation of 5% on a grid of 50 × 20 km in distance and altitude, respectively. The best regularization parameter was chosen in the way described in Section 4.6.3.
The reconstruction of the target model is quite well between 300 and 400 km altitude that is due to the GPS and the negative and positive anomalies can be clearly separated. Between 100 and 300 km altitude, the OTH radar contributes to reconstruction of the target model up to 800 km distance and the negative as well as the positive electron density perturbations are indicated. Probably, this inversion result can be improved further by applying a damped least squares inversion with constraints (Section 4.5.3). For larger distances, also the GPS finds no perturbation at low altitude, because of the high frequency of the GPS signal. However, the high altitude electron density irregularities are a little bit enlarged to the top boundary of the grid and do not end clearly before as in the target model.

<figure><img src="image_96.png" /><figcaption>Figure 7.20: Joint inversion of OTH radar and GPS (top) using simulated data for a checkerboard perturbation of 5%. Perturbation gaps in the target model are due to missing ray-paths. The grid size is 50 × 20 km in distance and altitude, respectively.</figcaption></figure>

The above described test showed that a joint inversion of GPS and OTH radar can improve the inversion result. However, this test is based on an unrealistic GPS station distribution along a line in the azimuth direction of the radar.

## 7.4 Conclusions

In this chapter, the ionospheric tomography method, previously tested for synthetic data, was applied to real data measured by OTH radar Nostradamus. The performance was analyzed on two data sets, for early evening in spring 2006 and at noon in winter 2008. Contrary to the synthetic data, inversion of real data needs some preprocessing of the data, in order to consider multi-path and path effects in the
real ionosphere compared to the *a priori* ionospheric model. A good tool to accomplish that was found in the virtual height of reflection, easily calculated by Martyns theorem.

For both data sets, the *v* and the *v&r*-method were used for inversion, and their performances were compared after one iteration. Generally, for the best regularization parameter chosen from the *L*-curve, the *v* method inversion results are strongly damped, their amplitude reduced and the electron density perturbations very localized. This can be explained, as discussed in Chapter 6 with the best value of regularization being larger than the eigenvalues in the matrix, suppressing all useful information. With the *v&r*-method and the same technique to choose a regularization parameter, anomalies are found within the range of the radar up to 1500 km distance. To go a step further, the electron density profiles from the *a priori* ionospheric model NeQuick, generated for the corresponding epoch and solar flux, were corrected with the electron density perturbations found by inversion.

However, the real data suffers from the reduced number of rays, compared to the previous synthetic tests, making the problem more non unique and ill-posed. Additionally, the Sun underwent a deep solar minimum in 2006 and still in 2008, consequently the solar flux and ionization in the ionosphere are small. This reduces the number of reflected rays and the ray coverage. This becomes evident in the inversions, where electron density anomalies are mainly reconstructed in the F-layer, even at noon where also an E-layer is present. Therefore, *a priori* vertical profiles of electron density are mainly corrected in the altitude range between 180 and 240 km and apart from that coincide with the vertical profile from the *a priori* ionospheric model.

As discussed in the perspectives, the vertical profiles of electron density can be compared with those from ionosonde stations in Europe. In a preliminary test this was demonstrated to be feasible with the station in Chilton and the *v&r*-method showed a better performance. Nevertheless, this needs to be developed further, with more iterations, since the *v*-method generally needs more iterations than the *v&r*-method. Additionally, the methodology used to invert the Chilton data needs to be understood deeply for a comparison to minimize the effect of the methodical differences on the results. The raw data from the Chilton ionosonde could be included into our general join inversion using several instruments (i.e., GPS, OTH radar, Ionosondes).

With the inversion of real data in 4 different azimuths, the impressive potential of the tomography method was illustrated, showing that a 3D tomography over Europe would be generally possible. As discussed and demonstrated in a synthetic test in the last section of this chapter, this ionospheric tomography method could be combined with the tomography done by GPS. This would improve the missing sensitivity of
the GPS to the lower ionosphere at 100 km altitude and the two methods would complement one another. Since the OTH radar Nostradamus is not available at the moment, only data from the previous years can be used, measured in the 8 azimuths in the normal operation mode. For this data, corresponding GPS orbits need to be found. Another option is to use radar data from one of the other OTH radar in the U.S, Australia or the SuperDARN. In order to apply the tomography method to these other radars, they need to provide measurements of the elevation angle and the frequency, or the tomography method needs to be adapted at their measurements.
# Chapter 8

# Conclusions and Perspectives

In the first part of this chapter the conclusions of this research will be summarized. There are several areas of this study that can be further explored or in which further developments can be made based on these results. These perspectives will be summarized in the second part of the chapter.

## 8.1 Conclusion

In this thesis, a method for ionospheric tomography for monostatic OTH radar (receiver and transmitter are at the same place) was developed for the first time. All previous established ionospheric tomography methods are mainly based on GPS. During the last decade, many efforts were made to improve these methods and to overcome limitations due to ray-coverage and the non-uniqueness of the inverse problem (Section 4.2). However, due to the high frequency of the GPS signal, the attempts at an ionospheric tomography were only sensitive to the region of maximum ionization (the F2-layer).

### Validation of the method

The developed tomography method for OTH radar was validated on a number of different synthetic benchmark tests (Chapter 4) that consists of checkerboard and localized perturbations. The *v*-method, considering only the effect of electron density perturbation on the propagation time of the EM waves, was found to be insufficient in a first iteration, because it does not reconstruct well the target models. Therefore, the *v&r*-method was established, taking into account the ray-path deflection introduced by the unknown endpoints of the rays.

The concept of ray-path deflection was validated by the idealized case of *frozen rays*, where the ray-path is not perturbed, as if the ray end-points were known. In this case, synthetic target models are well reconstructed in the zone of high ray
coverage within the range of the radar.

The *v&r*-method can be improved further, reducing the noise by smoothing with a damped-least squares inversion with constraints (Section 4.6.4). The *v*-method does not need a smoothing in the first iteration.

**Resolution analysis**

The *v* and the *v&r*-method were applied to a bunch of different synthetic benchmark tests of different size and amplitude. This detailed sensitivity and accuracy study allowed to investigate the resolution of the tomography method (Chapter 5).

In these tests, the influence of the grid, cell dimensions, regularization parameter, as well as smallest perturbation amplitude where analyzed. The main conclusions of these resolution tests are that both inversion results are independent of the grid size when using checkerboard or localized perturbations. Even better, two small localized perturbations of 165 × 20 km can be identified and separated by both methods.

However, in the first iteration, the *v&r*-method seems to perform better for checkerboard perturbations and localized perturbation because it can identify the correct position of the maximum anomaly and is able to reconstruct more details of the checkerboard perturbation in larger distance. The *v*-method mislocates the perturbation, is not able to reproduce its shape and in the case of the checkerboard it is strongly damped.

In both inversion results (*v* and *v&r*), large-scale negative anomalies are visible for the localized perturbation that are not present in the target model and do not appear in the *frozen ray* inversion. These artifacts can be explained with nonlinear effects of ray-path deflections, which (endpoints not being fixed and unknown) can in principle result in a faster propagation time even if the velocity perturbation is negative.

All synthetic tests illustrated the dependence of the solution on the ray-coverage and to the zone of sensitivity of rays to the medium, where the plasma frequency is close to the signal frequency. The checkerboard perturbation demonstrated clearly these dependences, because perturbations are not well resolved neither at large distances and zone with few rays, nor at low altitude.

**Development of an iterative methodology**

The inversion results after the first iteration for a localized perturbation showed that for the *v&r*-method the amplitude of the solution is too small, and the *v*-method mislocates the perturbation. In order to improve these inversion results, an iterative methodology was developed for the *v* and the *v&r*-method. The preliminary results for an iterative approach published by Roy *et al.* (2014) showed a decrease in data
and model misfit with iterations and a subsequent increase in the model misfit after
a critical number of iterations. At their model misfit minima the difference between
target and solution is less than 30% for *v&r* at iteration 13 and less than 40% for
the *v*-method at iteration 5 (Roy *et al.*, 2014). This lead to the conclusion that
the iteration approach generally improves the inversion results, because data and
model misfit are reduced compared to the first iteration, but the system becomes
unstable after a critical number of iterations. To deal with this problem and to
reduce the instability, the effect of feedback and regularization were tested for both
methods. The first parameter controls the amplitude, the second damps the re-
injected solution.

Both parameters were found to slow down convergence in both inversion meth-
ods, but they do not reduce the instability. However, the *v&r*-method reaches model
misfit minima generally faster than the *v*-method. This may be explained by the
fact that the *v&r*-method takes into account the ray-path deflection starting from
the first iteration, while the *v*-method seems to consider it indirectly with the ray
tracing in an updated model after the first iteration. This became obvious with
the horizontal movement of the electron density perturbation to the correct position
with iterations in the *v*-method. The model misfit for the *v* and the *v&r*-method
is in the same order at their minima. The instability in the inversion, more severe
in the *v&r* method than in the *v*-method, may be explained by the numerical noise
induced by the cells discretization, or by the non-linearity of the ray-path deflection.

Since the *v&r*-method is computationally very expensive, the code was par-
allelized to take advantage of multi-core architectures of modern computers. In
general, for one iteration with the *v&r*-method, 15 iterations with the *v*-method can
be calculated in the same time. However, often the *v*-method did not reach a model
misfit minimum after 30 iterations yet.

In conclusion, the *v&r*-method is computationally more expensive, numerically
more complex and needs a smoothing in order to reduce the noise in the solution at
low altitude, but, on the other hand, is faster compared to the *v*-method.

Real data inversion

The *v* and the *v&r*-method were applied successfully to real data collected with the OTH radar Nostradamus. Two data sets were tested, one of spring 2006 in the early evening, and another of winter 2008 at noon. The inversion provide a map of electron density perturbation in the selected azimuth for the given day and time. In the first data set, only electron density perturbations in the F-layer are visible, in the second data set an E-layer is also present.

Based on such maps, electron density profiles from an *a priori* ionospheric model
can be corrected in any desired distance within the range of the radar. This is a huge
improvement compared to, for instance, vertical profiles calculated on the base of
backscattered ionogram inversions, because this method allows only the calculation
of the electron density at the middle of the ray-path and it is assumed to be valid
for all other distances or it is interpolated.

For the first time a reconstruction of the real electron density distribution in
the ionosphere was obtained from OTH radar tomography, revealing vertical irreg-
ularities in the region between 200 km and 300 km height. Performing the inversion
in different azimuths showed that these irregularities are visible in other directions,
too at the same day. Additionally, this first test illustrated the impressive potential
of the ionospheric tomography method to reconstruct the electron density over the
entity of Europe with Nostradamus.

In a preliminary synthetic test, vertical profiles of electron density obtained after
the first iteration were compared with those from the ionosonde in Chilton, U.K in
order to validate the inversion results. The *v&r*-method corrects the profile from the
*a priori* ionospheric model towards the profile at the ionosonde, but more iterations
are necessary to reduce the difference between them.

However, the Sun underwent a solar minimum in 2006 and 2008 and as a con-
sequence, solar flux and ionization in the ionosphere are weak. Therefore, less rays
are reflected in the *a priori* ionospheric model in the ray tracing, and the num-
ber of rays is reduced additionally by using the frequency and elevation angles of
Nostradamus measurements. That makes the inverse problem more ill-posed and
non-unique compared to the inversion of synthetic data.

A second preliminary synthetic test showed that the tomography by OTH radar
can be combined with the tomography by GPS and that these two methods can
complement each other. The tomography with GPS is mainly sensitive to the region
of maximum ionization (F2 layer) because of the high frequency signal frequency.
On the other hand, the OTH radar can reconstruct quite well the electron density
in the lower ionosphere up to 300 km altitude.

8.2 Perspectives

Some remaining issues retain open and are interesting subjects of forthcoming work.
The presented ionospheric tomography method was developed for monostatic OTH
radar. These kind of radars can be found also in the USA (ROTHR, Australia
(JORN and the North-Pole (SuperDARN) and the method can be applied to the
previously collected data of the french radar Nostradamus. If the resolution of the
other OTH radars is sufficient, and they can measure the polarization of a wave,
the effect of the magnetic field need to be included in the ray tracing to distinguish
between ordinary and extraordinary wave.

At the same time, this tomography method is not limited to monostatic OTH
radars, and can be easily adapted to bistatic OTH radars. Although here the ray-
path deflection needs to be considered, because the signal is still backscattered at
the ground, at an unknown point, and returns to the receiver along another way
through the ionosphere.

This tomography approach can be improved further by including other iono-
spheric sounding techniques, in particular TEC measured by GPS at the ground
with dense arrays or by occultation with on-boarded GPS receivers. That latter
idea is interesting, because the two methods complement each other: While OTH
ionospheric tomography is sensitive to the lower ionosphere (up to 300 km altitude)
during daytime, GPS ionospheric tomography is only sensitive to the region of maxi-
mum ionization (~300 km), because of the high frequency of the signal. In addition,
TEC GPS data can compensate for the lack of reflected HF rays during night. The
feasibility of such a joint inversion was demonstrated in a preliminary synthetic test
(Section 7.3.2).

GPS ionospheric tomography was already combined with other sounding tech-
niques in previous work, for instance with ionosonde data and backscattered iono-
gram inversions (see Chapter 4), but never with an OTH radar tomography method.
This could open exciting new perspectives, improving the vertical resolution and
supplying more information on electron density at the lower ionosphere.

In this thesis the developed ionospheric tomography method was applied success-
fully to real data from OTH radar Nostradamus (Chapter 7), showing the potential
to obtain 3D maps of electron density over entire Europe. This tomographic ap-
proach can be developed further to four dimensions, including the time and providing
the evolution of electron density during day and night, or over an entire year.

The obtained vertical profiles of electron density need to be validated by ionosonde
data in Europe or other inversion methods. The validation with the data of the ver-
tical sounder in Chilton has been tested, but needs further development with more
iterations.

The developed iterative approach showed the general improvement of the solu-
tion, but also that the system becomes unstable after a critical number of iterations.
This could be explained by the numerical noise introduced by the discretization of
the model and can be addressed by a more non-linear inversion method, for instance
the CGM (Section 4.6.5).
Acknowledgments

In the following, I want to thank some people for supervision and support, both technically or morally, through the many ups and downs but nevertheless exciting and informative last three years. First of all, I thank the **ONERA** for giving me the possibility and means to realize this thesis as well as the **Planetary and Space Sciences group** at the IPGP.

I express my sincere gratitude to Giovanni Occhipinti for giving me the chance to continue his work such an exciting field. Our lively discussions have always led to new insights and ideas and I acknowledge his always positive and encouraging attitude. Additionally, I thank him for the possibility to participate in numerous conference, as well as his support regarding invited seminars and their preparation in form of trial presentations.

I thank Gilbert Auffray for his warm welcome in his research team at the ONERA, his care and help to bring my thesis to a good direction. Many thanks to all my ONERA colleges of the RBF team. Thanks to them, I enjoyed working at the ONERA. Without their friendship, help with my thesis, in administrative or technical problems, and their humor in everyday life, this thesis would have been much more difficult to finalize. In this context, it was a pleasure for me sharing an office with Nicolas, Jean-François and later with Marie-José, giving the distraction I needed every now and then.

I deeply thank Jean-Philippe Molinié for supervising me continuously during the three years. His ideas and suggested solutions have significantly contributed to the success of this work and I appreciate the informative and, due to his immense knowledge, very helpful discussions. I also appreciate his availability to listen to any of my personal problems, his patience especially with my explanations in french, and advice.

Moreover, I thank Mark Wieczorek for having been my official thesis director until the end, as well as for his support, advice and helpful suggestions during my thesis.

A special thank goes to Lapo Boschi for helping me in the development of the conjugate gradient method and his endless and patient explanations. He always took the time to answer my questions, even outside office-hours, and his e-mails and
photo explanations helped me a lot.

Further thanks go to **Véronique Rannou** for fruitful discussions on real data inversions and all issues concerning the ionosphere. I thank **Pierdavide Coïsson** for his explanations of the ionosphere, the discussions he had, and that he provided me a (genius !!) LaTeX thesis template.

Last but not least, I am grateful to everybody who generously offered me help during my quest of Ph.D degree and to all people who spent time in reading my manuscript for improvement and error search, in particular my **mother, Veronique, Pierdavide, and Valérie**.

Finally, I thank my **parents and sisters** for their encouragement during bad times, and their support in any situation I needed them. Without them, I would never have finished this thesis. They gave me the strength to go further when I doubted in my possibilities.
List of Figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Densité des ions dans atmosphère</td>
      <td>3</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Carte global du TEC (Contenu Electronique Total)</td>
      <td>7</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>La tomographie de l'ionosphère par GPS</td>
      <td>8</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Tracé de rayon</td>
      <td>9</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Résultats d'inversion après la première itération pour un test de damier (bas) et une perturbation localisée de 0.1 % (haut) pour différents paramètres de régularisation λ en utilisant l'approche des rayons figés, la méthode v et la méthode v&r.</td>
      <td>12</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>Résultats d'inversion pour trois modèles cibles et différentes grilles en utilisant la méthode v et la méthode v&r</td>
      <td>13</td>
    </tr>
    <tr>
      <td>1.7</td>
      <td>Ecarts en temps de propagation et des écarts par rapport au modèle cible avec les itérations pour la méthode v et v&r</td>
      <td>15</td>
    </tr>
    <tr>
      <td>1.8</td>
      <td>L'inversion des données réelles mesuré par Nostradamus</td>
      <td>16</td>
    </tr>
    <tr>
      <td>1.9</td>
      <td>L'inversion conjoint des données simulées du radar transhorizon et du GPS</td>
      <td>18</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Daytime ion densities in the atmosphere</td>
      <td>23</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Ionospheric layers</td>
      <td>24</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Global TEC maps</td>
      <td>25</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Examples for vertical electrons density profiles</td>
      <td>27</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>NeQuick2 vs. IRI</td>
      <td>31</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Example of an ionogram measured at the ionosonde Station Juliusruh, Germany</td>
      <td>37</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Worldwide distribution of ionosondes created on the basis of the station map from NOAA as at June 2014.</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>An elevation scan backscatter ionogram from OTH radar Nostradamus.</td>
      <td>39</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>Inversion result after Fridman (1998) over Scandinavia for days in summer and winter during night using backscatter ionogram inversion</td>
      <td>40</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Principle of GPS occultation measurement</td>
      <td>42</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Schematic representation of OTH and a classical radar</td>
      <td>44</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Fields of view of currently operational SuperDARN radars in northern (left) and southern hemisphere (right)</td>
      <td>47</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Example of an antenna of Nostradamus.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Simulations for the capacities of OTH radar Nostradamus</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Panoramic sounding of the ionosphere with Nostradamus</td>
      <td>49</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>4.1</td>
      <td>Principle of GPS ionospheric tomography</td>
      <td>53</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Example of traced rays in the ionospheric model NeQuick</td>
      <td>61</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Validation of the linearization of the refractive index</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Schematic representation of the parametrization of the ionospheric model and the effect of ray-path deflection</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Sketch for creating eight constraint matrices</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Proof of the linearity of the inverse problem</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Inversion results for a checkerboard benchmark test using analytically and numerically calculated frozen rays</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.8</td>
      <td>Inversion results for a localized perturbation using analytically and numerically calculated frozen rays</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.9</td>
      <td>Ratio of plasma to signal frequency along the rays</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.10</td>
      <td>v-method vs v&r-method for a localized perturbation of 0.1 %</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.11</td>
      <td>v-method vs v&r-method for a checkerboard perturbation of 0.1 %</td>
      <td>77</td>
    </tr>
    <tr>
      <td>4.12</td>
      <td>L-curves for the localized and the checkerboard benchmark test</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.13</td>
      <td>Dependence of the inversion results on the range of regularization parameters used for the calculation of the L-curve</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.14</td>
      <td>Error curves for the determination of the best combination of parameters α₀ and α₁</td>
      <td>80</td>
    </tr>
    <tr>
      <td>4.15</td>
      <td>Damped least squares inversion vs. damped least squares with constraints for a localized perturbation</td>
      <td>82</td>
    </tr>
    <tr>
      <td>4.16</td>
      <td>Damped least squares inversion vs. damped least squares with constraints for a checkerboard perturbation</td>
      <td>82</td>
    </tr>
    <tr>
      <td>4.17</td>
      <td>Conjugate Gradient vs. steepest descent Method</td>
      <td>84</td>
    </tr>
    <tr>
      <td>4.18</td>
      <td>Model misfit as a function of iterations for different model covariances and a fixed data covariance of 8.3 · 10⁻⁶s</td>
      <td>86</td>
    </tr>
    <tr>
      <td>4.19</td>
      <td>Data misfit as a function of iterations for different model covariances and a fixed data covariance of 8.3 · 10⁻⁶s</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.20</td>
      <td>Solution of the Conjugate Gradient Method applied to a localized perturbation</td>
      <td>88</td>
    </tr>
    <tr>
      <td>4.21</td>
      <td>Data and model misfit as a function of iterations for the Conjugate Gradient Method</td>
      <td>88</td>
    </tr>
    <tr>
      <td>4.22</td>
      <td>Solution of the Conjugate Gradient Method applied to a checkerboard perturbation</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4.23</td>
      <td>Data and model misfit as a function of iterations for the Conjugate Gradient Method applied to a checkerboard perturbation</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Effects of increasing cell dimensions on the inversion results for a localized perturbation of 0.1%</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Effects of increasing cell dimensions in altitude on the inversion results for a localized perturbation in one cell of 0.1%</td>
      <td>93</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Effects of increasing cell dimensions on the inversion results for a checkerboard perturbation of 0.1%</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Effects of larger cell dimensions on the inversion results for another checkerboard perturbation of 0.1%</td>
      <td>96</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>5.5</td>
      <td>Effects of increasing cell dimensions on the inversion results for a radical checkerboard test of 0.1%</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Effects of larger cell dimensions in distance on the inversion results for two localized perturbations of 0.1%</td>
      <td>99</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Effects of increasing cell dimensions on the inversion results for a localized perturbation in one cell of 0.1%</td>
      <td>101</td>
    </tr>
    <tr>
      <td>5.8</td>
      <td>Effect of the amplitude of perturbation on the inversion result</td>
      <td>103</td>
    </tr>
    <tr>
      <td>5.9</td>
      <td>Ratio of plasma to emission frequency</td>
      <td>104</td>
    </tr>
    <tr>
      <td>5.10</td>
      <td>Model misfit as a function of the grid size</td>
      <td>105</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Computational costs for numerical calculation of ray-path deflection in the v&r-method</td>
      <td>108</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Computational costs for ray tracing in the v-method</td>
      <td>109</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Data and model misfit as a function of iterations for the v and v&r-method</td>
      <td>110</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>The inverse problem represented as closed loop control system</td>
      <td>112</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Data misfit and model misfit as function of iterations using different feedback values in the v&r-method</td>
      <td>113</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Inversion results after iterations for different feedback values using the v&r-method and the L-curve as regularization criterion.</td>
      <td>114</td>
    </tr>
    <tr>
      <td>6.7</td>
      <td>50 largest eigenvalues λ of matrix M<sup>T</sup>·M in descending order</td>
      <td>116</td>
    </tr>
    <tr>
      <td>6.8</td>
      <td>Data and model misfit evolution in the iterative v&r-method</td>
      <td>117</td>
    </tr>
    <tr>
      <td>6.9</td>
      <td>Inversion results for different regularization methods in the v&r-method</td>
      <td>119</td>
    </tr>
    <tr>
      <td>6.10</td>
      <td>Iterative damped least squares inversion vs. iterative damped least squares inversion with constraints</td>
      <td>121</td>
    </tr>
    <tr>
      <td>6.11</td>
      <td>Effect of a feedback value on the inversion results of the v-method with the L-curve as regularization</td>
      <td>123</td>
    </tr>
    <tr>
      <td>6.12</td>
      <td>Inversion results after iterations for different feedback values using the v-method and the L-curve as regularization criterion</td>
      <td>124</td>
    </tr>
    <tr>
      <td>6.13</td>
      <td>50 largest eigenvalues λ of the matrix A<sup>T</sup>·A in descending order</td>
      <td>125</td>
    </tr>
    <tr>
      <td>6.14</td>
      <td>Effect of regularization on the iterative v-method - Data and model misfits</td>
      <td>126</td>
    </tr>
    <tr>
      <td>6.15</td>
      <td>Effect of regularization on the iterative v-method - Inversion results</td>
      <td>128</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Real data obtained from the OTH radar Nostradamus on March 14, 2006 at 1855 UT</td>
      <td>132</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Explanation of Martyn's Theorem</td>
      <td>133</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Relationship between virtual height of reflection, group distance and elevation angle</td>
      <td>133</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Clustering</td>
      <td>134</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Raw data from March 14, 2006 in azimuth 202° and assigned to different ionospheric layers using a clustering algorithm (right)</td>
      <td>135</td>
    </tr>
    <tr>
      <td>7.6</td>
      <td>Raw data from March 14, 2006 in azimuth direction 202° compared to raytracing</td>
      <td>136</td>
    </tr>
    <tr>
      <td>7.7</td>
      <td>Explanation of path effects in the ionosphere</td>
      <td>137</td>
    </tr>
    <tr>
      <td>7.8</td>
      <td>Difference in the propagation time between ray tracing in an a priori ionospheric model and real data in percent</td>
      <td>138</td>
    </tr>
  </tbody>
</table>
List of Figures

<table>
  <tbody>
    <tr>
      <td>7.9</td>
      <td>Real data inversion with the v&r and the v-method</td>
      <td>139</td>
    </tr>
    <tr>
      <td>7.10</td>
      <td>Vertical profile of electron density in 500 km distance from the radar calculated from real data inversion</td>
      <td>140</td>
    </tr>
    <tr>
      <td>7.11</td>
      <td>Vertical profile of electron density in 1000 km distance from the radar calculated from real data inversion</td>
      <td>141</td>
    </tr>
    <tr>
      <td>7.12</td>
      <td>Electron density perturbation over Europe for March 14, 2006 at 1855 UT</td>
      <td>142</td>
    </tr>
    <tr>
      <td>7.13</td>
      <td>Raw data from December 13, 2008 in azimuth direction 67° compared to ray tracing in the a priori ionospheric model</td>
      <td>143</td>
    </tr>
    <tr>
      <td>7.14</td>
      <td>Inversion results using the data measured December 13, 2008 at 1154 UT</td>
      <td>144</td>
    </tr>
    <tr>
      <td>7.15</td>
      <td>Vertical profiles of electron density at 200 km distance from the radar Nostradamus</td>
      <td>144</td>
    </tr>
    <tr>
      <td>7.16</td>
      <td>Vertical profiles of electron density at 500 km distance from the radar Nostradamus</td>
      <td>145</td>
    </tr>
    <tr>
      <td>7.17</td>
      <td>Ionogram measured at Chilton at December 19, 2008 at noon</td>
      <td>147</td>
    </tr>
    <tr>
      <td>7.18</td>
      <td>Inversion results using real data measured by OTH radar Nostradamus December 19, 2008 at noon in 337° azimuth</td>
      <td>148</td>
    </tr>
    <tr>
      <td>7.19</td>
      <td>Comparison of the vertical profile of electron density obtained at the ionosonde station Chilton with the one obtained by ionospheric tomography</td>
      <td>149</td>
    </tr>
    <tr>
      <td>7.20</td>
      <td>Joint inversion of OTH radar and GPS</td>
      <td>151</td>
    </tr>
  </tbody>
</table>
List of Tables

<table>
  <tbody>
    <tr>
      <td>6.1</td>
      <td>Model misfit and corresponding data misfits for results in Figure 6.3.</td>
      <td>110</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Minimum model misfit and corresponding data misfit for results in Figure 6.5 using different feedback values α in the iterative v&r-method. The L-curve criterion was used to choose the best regularization parameter.</td>
      <td>113</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Minimum model misfit and corresponding data misfit for results in Figure 6.8 using different regularization methods and a feedback of α = 0.5 in the iterative v&r-method.</td>
      <td>118</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Minimum model misfit and corresponding data misfit for results in Figure 6.10.</td>
      <td>120</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Model misfit minimum and corresponding data misfit for results in Figure 6.11 depending on different feedback values α in the iterative v-method.</td>
      <td>123</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Data and model misfit for results in Figure 6.15 depending on different regularization methods and a feedback of α = 0.1 in the iterative v-method.</td>
      <td>127</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Ionosonde stations in Europe, distance and azimuth to OTH radar Nostradamus</td>
      <td>146</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_97.png" /></figure>
# Appendices
<figure><img src="image_98.png" /></figure>
# Appendix A

## Article: New constraints on the 3D shear wave velocity structure of the upper mantle underneath Southern Scandinavia revealed from non-linear tomography

I strongly contributed to this work as a student assistant at the Geophysical Institute at the Karlsruhe Institute of Technology (KIT) in Germany. The collaboration continued during my Ph.D thesis at the IPGP, and the article was finished and published in January 2013.

### Abstract

In this work, the 3D shear-wave velocity underneath Southern Scandinavia was revealed from seismic tomography, using the MAntle investigations of Norwegian Uplift Structures (MAGNUS) seismic station network in Norway. The travel time residuals are corrected for the known crustal structure of Southern Norway and weighted to account for data quality and pick uncertainties.

The resulting residual pattern of subvertically incident waves is very uniform and simple. It shows delayed arrivals underneath Southern Norway compared to fast arrivals underneath the Oslo Graben and the Baltic Shield. The 3D upper mantle $v_s$ structure underneath the station network is determined by performing non-linear travel time tomography. As expected from the residual pattern the resulting tomographic model shows a simple and continuous $v_s$ perturbation pattern: a negative $v_s$ anomaly is visible underneath Southern Norway relative to the Baltic Shield in the east with a contrast of up to 4 % $v_s$ and a sharp W-E dipping transition zone. Reconstruction tests reveal besides vertical smearing a good lateral reconstruction of the dipping $v_s$ transition zone and suggest that a deep-seated anomaly at depth is real and not an inversion artifact.

The upper part of the reduced $v_s$ anomaly underneath Southern Norway (down to 250 km depth) might be due to an increase in lithospheric thickness from the Caledonian Southern Scandins in the west towards the Proterozoic Baltic Shield in Sweden in the east. The deeper-seated negative $v_s$ anomaly (300 - 410 km depth) could be caused by a temperature anomaly possibly combined with effects due to fluids or hydrous minerals. The determined simple 3D $v_s$ structure underneath Southern Scandinavia indicates that mantle processes might influence and contribute to a Neogene uplift of Southern Norway
<figure><img src="image_99.png" /></figure>
# New constraints on the 3D shear wave velocity structure of the upper mantle underneath Southern Scandinavia revealed from non-linear tomography

B. Wawerzinek *, J.R.R. Ritter, C. Roy <fnref n="1" />

Karlsruhe Institute of Technology, Geophysical Institute, Hertzstr. 16, 76187 Karlsruhe, Germany

## ARTICLE INFO

**Article history:**

Received 25 May 2012

Received in revised form 2 December 2012

Accepted 23 December 2012

Available online 12 January 2013

**Keywords:**

Mantle processes

Seismic tomography

Shear waves

Relative travel time residuals

Southern Scandinavia

## ABSTRACT

We analyse travel times of shear waves, which were recorded at the MAGNUS network, to determine the 3D shear wave velocity ($v_S$) structure underneath Southern Scandinavia. The travel time residuals are corrected for the known crustal structure of Southern Norway and weighted to account for data quality and pick uncertainties. The resulting residual pattern of subvertically incident waves is very uniform and simple. It shows delayed arrivals underneath Southern Norway compared to fast arrivals underneath the Oslo Graben and the Baltic Shield. The 3D upper mantle $v_S$ structure underneath the station network is determined by performing non-linear travel time tomography. As expected from the residual pattern the resulting tomographic model shows a simple and continuous $v_S$ perturbation pattern: a negative $v_S$ anomaly is visible underneath Southern Norway relative to the Baltic Shield in the east with a contrast of up to 4% $v_S$ and a sharp W-E dipping transition zone. Reconstruction tests reveal besides vertical smearing a good lateral reconstruction of the dipping $v_S$ transition zone and suggest that a deep-seated anomaly at 330–410 km depth is real and not an inversion artefact. The upper part of the reduced $v_S$ anomaly underneath Southern Norway (down to 250 km depth) might be due to an increase in lithospheric thickness from the Caledonian Southern Scandes in the west towards the Proterozoic Baltic Shield in Sweden in the east. The deeper-seated negative $v_S$ anomaly (330–410 km depth) could be caused by a temperature anomaly possibly combined with effects due to fluids or hydrous minerals. The determined simple 3D $v_S$ structure underneath Southern Scandinavia indicates that mantle processes might influence and contribute to a Neogene uplift of Southern Norway.

© 2013 Elsevier B.V. All rights reserved.

## 1. Introduction

At the western edge of the Baltic Shield, Northern Europe, the Scandinavian Mountains (Scandes) form an about 1400 km long mountain range (Fig. 1). The Scandes are the second largest mountain range in Europe and extend from 59°N to 63°N. The present high topography of the Scandes is still under debate (Lidmar-Bergström and Bonow, 2009; Nielsen et al., 2009). The main mountain building phase occurred during the Scandinavian Caledonian orogeny 440–420 Ma ago when Baltica and Laurentia collided (Torsvik and Cocks, 2005). Today the mountains still reach maximum heights of 2.5 km and the Southern Scandes in Norway have roughly a dome-like shape (Fig. 1). This dome has a high-level low-relief landscape with a base at about 1000 m altitude on which the mountain massifs are sited (Lidmar-Bergström and Bonow, 2009; Lidmar-Bergström et al., 2000). There are different hypotheses which try to explain the present high topography of the Scandes far

away from current plate boundaries and in the absence of active compressional tectonics (Cloetingh et al., 2007). Some studies come to the conclusion that there were later uplift phases after the Caledonian orogeny, for example during the Neogene (e.g. Lidmar-Bergström et al., 2000; Rohrman and van der Beek, 1996; Smelror et al., 2007). Japsen and Chalmers (2000) and Anell et al. (2009) summarise widespread high topography areas, Cenozoic uplift and their geodynamic driving forces around the North Atlantic. Another point of view is that the Caledonides were formed during one main orogenic phase and that the present topography is a result of the interaction between erosion, climate and ice (Nielsen et al., 2009). However, this hypothesis is in dispute (Gabrielsen et al., 2010; Lidmar-Bergström and Bonow, 2009 and references therein).

The debate on this controversial issue is mainly due to sparse and insufficient data which limit our knowledge on the deep structure of the Scandes and thus prevent the determination of a consistent geodynamic model. However, recent seismological experiments improved this situation. Using receiver functions to map the Moho, Svenningsen et al. (2007) found a 10–12 km thick crustal root underneath the Southern Scandes and claimed this root as Airy-type compensation for the mountains. Stratford and Thybo (2011a) reviewed existing crustal models, presented the results from recently measured seismic refraction lines (MAGNUS-REX project) and compiled a Moho

<footnote>* Corresponding author at: Leibniz Institute for Applied Geophysics, Stilleweg 2, 30655 Hannover, Germany. Tel.: +49 511 643 3389.</footnote>

<footnote>E-mail address: britta.wawerzinek@liag-hannover.de (B. Wawerzinek).</footnote>

<footnote n="1">¹ Now at: Institut de Physique du Globe de Paris, Université Paris Diderot, 5 rue Thomas Mann, 75205 Paris Cedex 13, France.</footnote>
<figure><img src="image_100.png" /><figcaption>**Fig. 1.** Map of the study region and the MAGNUS station network. Numbers correspond to mobile stations (circles, the leading NWG station code is omitted for clarity), while permanent stations (triangles) are labelled with complete station code. The elevation is colour-coded and the 800 m contour line is plotted to visualise the dome-like shape of the Southern Scandes. For global context of the study area see Fig. 2.</figcaption></figure>

map. Their model has a thinner crustal root and its 40 km thick crust is regarded as typical shield-type crust. Stratford and Thybo (2011a) interpreted the observed 60 km lateral offset between the deepest Moho topography and the highest surface topography by a possible lithospheric flexure. The increased seismic velocities (up to 5%) around the Palaeozoic Oslo Graben, which were determined from the MAGNUS-REX data, are explained by mafic intrusions (Stratford and Thybo, 2011b).

Below the 22 NORSAR (Fig. 1) subarrays, Aki et al. (1977) found a relatively low seismic velocity region to the west compared to a higher seismic velocity region to the east. This difference in seismic P-wave velocity ($v_p$) is interpreted as Baltic Shield material in the east, which is undisturbed by the Caledonian orogeny, and the rifting of the Oslo Graben or its related volcanism. Medhus et al. (2009) determined travel time residuals of teleseismic P-waves which were recorded during field experiments and at permanent stations in western Scandinavia. Their residual pattern of late arrivals indicates reduced seismic velocity in the upper mantle below the Southern Scandes and the Norwegian-Danish Basin, whereas fast arrivals are observed towards east at the Oslo Graben and into Sweden (Medhus et al., 2009). Two teleseismic tomography studies by Medhus et al. (2012) present relative and absolute P-wave velocity perturbations in the upper mantle below the Southern Scandes and their surroundings: low $v_p$ is recovered underneath the Southern Scandes and the Danish and German basins. The contrast to the faster upper mantle $v_p$ under the Baltic Shield is up to 5% and its western boundary is clearly mapped (Medhus et al., 2012). Reduced S-wave velocity ($v_s$) is also known below the Scandes from surface wave analyses: Weidle and Maupin (2008) determined a low $v_s$ anomaly reaching from Iceland across the

eastern North Atlantic below southern Norway at about 70–150 km depth. Maupin (2011) confirmed the low $v_s$ anomaly in the mantle underneath Southern Norway. The influence of these low seismic velocity anomalies on the topography on the Scandes is still unclear. Although dynamic topography may play a role for generation of high elevation (for a review see Braun, 2010), its role in Scandinavia is disputed. Pascal and Olesen (2009) demonstrated with integrated gravity and thermal modelling that Cenozoic uplift of the Scandes cannot be explained only by an asthenospheric diapir.

In the above context we want to address the following questions: How does the mantle structure look like at the transition between the offshore and onshore Norwegian lithosphere including the Southern Scandes? How does the mantle structure change further east towards the Baltic Shield? Are there structural elements at depth which may help to describe the mountain building processes in the region and what is the lateral and depth extension of possible anomalies? The presence of seismic anomalies and interpretations concerning their origin help to explain the building process at depth. For example, a low velocity anomaly may be caused by increased temperature and relatively shallow asthenospheric material. Such knowledge can indicate that a fraction of the present high topography is due to dynamic topography (Braun, 2010). If no velocity anomalies are present underneath the Scandes a completely different mountain building process without active mantle processes would be the cause for the present topography. The idea of this work is to measure travel times of teleseismic shear waves recorded at the stations of the MAGNUS network and to determine travel time residuals relative to a laterally homogeneous standard Earth model, e.g. iasp91 (Kennett and Engdahl, 1991). These residuals are inverted into 3D shear
wave velocity perturbations underneath the station network. Our inversion technique is very similar to the famous first seismic tomography by Aki et al. (1977) which was also conducted in Southern Norway.

## 2. The MAGNUS experiment

The MAGNUS experiment (MAntle investiGations of Norwegian Uplift Structures) was a broadband seismological field experiment to record teleseismic and local earthquakes across Southern Norway (Fig. 1). Thirty-one mobile broadband recording stations of the KArlsruhe BroadBand Array (KABBA) were deployed in and around the Southern Scandes between September 2006 and June 2008. The KABBA stations together with the permanent broadband stations KONO (Kongsberg), BER (Bergen), HFC2 (Hagfors) and 7 stations of the Norwegian Seismic Array (NORSAR) covered a nearly equally spaced network across the Southern Scandes. The MAGNUS dataset was integrated into the EUROCORE programme TOPO-EUROPE (Cloetingh et al., 2007) of the European Science Foundation and therein specifically into the Collaborative Research Project *TopoScandiaDeep* (http://www.mn.uio.no/geo/english/research/projects/toposcandiadeep/, state 29th November 2012) which concentrates on the deep structure of the Scandes. For details on MAGNUS see Weidle et al. (2010).

A relatively good azimuthal coverage of teleseismic earthquakes was achieved due to continuous recording with very few station failures. Seismic phases from 179 earthquakes were studied in order to determine travel time residuals across the MAGNUS network. 128 of these events (Fig. 2) with signal-to-noise ratios larger than 2 were suitable for the shear wave tomography, because seismic phases of these events could be clearly seen at most stations. The moment

magnitude $M_w$ of the earthquakes was mostly above 6. However, weaker events underneath Africa were specifically examined to fill observational gaps in the south.

## 3. Analysis of travel time residuals

### 3.1. Determination of relative, weighted, crust corrected travel time residuals

At first the raw waveforms are pre-processed including a removal of possible offsets and trends in the recordings. The individual instrument response is removed by deconvolution to obtain true amplitude velocity-proportional seismograms. Finally all seismograms are bandpass filtered from 8 s to 20 s (0.05–0.125 Hz). This frequency band is identified as most suitable, because the seismic phases can be seen best at most stations in this band. Generally, the signal-to-noise ratio (SNR) is better at inland stations compared to stations close to the Atlantic and North Sea coasts. In order to enhance the polarised shear waves a coordinate rotation is done which transforms the horizontal N–S and E–W recording components into transverse (T) and radial (R) components. The theoretical backazimuth to the epicentres is used for the rotation.

We determine travel time residuals for 128 teleseismic events which we use as input data for the tomography. Following Evans and Achauer (1993), we pick the first clearly visible phase (up or down in a seismogram) of a shear wave front across the whole network. The arrival time picking was done by hand at magnified waveforms (Douglas et al., 1997) of the direct S or SKS wave fronts. Both R and T component waveforms could be used from 79 events (direct S-waves), only T component waveforms from 28 events (direct S-waves) and only R component waveforms from 21 events (mainly SKS-waves). Altogether 8348 pick times are determined

<figure><img src="image_101.png" /><figcaption>**Fig. 2.** Distribution of the analysed 128 earthquakes (diamond symbols) in dependence of their epicentral distance (in degree; radial axis) and backazimuth. The location of the MAGNUS network is indicated by a triangle in the centre. The size of the symbols indicates the moment magnitude $M_w$ of the events. Hypocentral parameters are taken from the International Seismological Centre and the U.S. Geological Survey.</figcaption></figure>
(Table 1). Depending on the shape of the S or SKS wavelet and its SNR, a quality factor is attributed to each pick. High-quality picks with a high SNR and symmetric waveform are “A”-class, less clear picks are “B”-class, acceptable picks are “C”-class and doubtful picks are “X”-class (see supplement by Kirschner et al., 2011 for data picking and quality issues). “X”-class picks are removed from the further analysis. In this way 6196 travel time data remain for the residual analysis, including 2859 T component picks and 3337 R component picks. The quality distribution is 942 “A”, 3109 “B” and 2145 “C” class picks. If both component recordings, R and T, can be used for picking, only the picks of the component recording with the better data quality are kept and the picks of the other component recording are removed. The component selection (R or T) is done event wise by determining the better pick quality, e.g. the larger number of “A” picks, of the component recordings. This sorting reduces our dataset to 3999 high-quality travel time residuals including 757 “A”, 2041 “B” and 1201 “C” class picks.

To account for station elevation and known variations of the gross crustal structure, we apply a 1D crustal travel time correction following Martin et al. (2005). This correction takes into account the travel time between sea level and true station height, sediment effects (here not necessary due to missing sediments below the MAGNUS stations) and Moho topography relative to the *iasp91* reference Earth model (Moho at 35 km depth; Kennett and Engdahl, 1991). Ray tracing is done with the theoretical slowness of the wave. The 1D crustal travel time correction is isotropic, so that travel times on R and T components are affected in the same way. For Norway the Moho topography was taken from Stratford et al. (2009) and it varies between 28 km depth (station NWG22 at the coast) and 40 km depth (stations NWG05 and NWG08 in the NE of the MAGNUS network). In Fig. 3 the crustal travel time corrections are shown for each station site based on an average slowness of 7 s/° for teleseismic shear waves. Along the coast teleseismic shear waves arrive about 0–0.5 s earlier relative to *iasp91* (35 km Moho depth; Kennett and Engdahl, 1991) due to the thinner crust (about 28–30 km). Inside the Southern Scandes the arrival time is delayed by up 0.4 s due to a thin crustal root and the high station elevation. Compared to the residuals described below, the crustal correction is about 20%–30% and is needed to avoid smearing of crustal travel time anomalies into the underlying mantle during the inversion.

The calculation of the residuals is done in the following standardised way (see also Evans and Achauer, 1993 or Kirschner et al., 2011): Hypocentral parameters including origin time were mostly taken from the International Seismological Centre, preferably updated Engdahl et al. (1998) solutions. The theoretical travel time is determined based on the *iasp91* reference Earth model (Kennett and Engdahl, 1991). The measured travel time is corrected for first order crustal travel time effects due to surface and Moho topography (Fig. 3). The difference between the measured and the theoretical travel times is weighted according to the quality classes of the picks (weights: “A” 1.0, “B” 0.5 and “C” 0.25), see Evans and Achauer (1993) or Kirschner et al. (2011). Finally, the average residual value of each event phase is subtracted to eliminate far-field effects such as hypocentre uncertainties and source side structural anomalies.

<caption>**Table 1**</caption>

Overview on number of determined S- and SKS-wave residuals; the overall analysed residuals are given on the left side; the residuals used for the tomographic inversion are given on the right side.

<table><thead><tr><th>Phase</th><th>Comp.</th><th>No. overall<br>A / B / C / X quality</th><th>No. for inversion<br>A / B / C quality</th></tr></thead><tbody><tr><td>S</td><td>R</td><td>242 / 945 / 614 / 685</td><td>126 / 362 / 198</td></tr><tr><td>S</td><td>T</td><td>594 / 1698 / 1045 / 1002</td><td>556 / 1365 / 712</td></tr><tr><td>SKS, SKKS</td><td>R</td><td>106 / 466 / 486 / 465<br>Σ 942 / 3109 / 2145 / 2152</td><td>75 / 314 / 291<br>Σ 757 / 2041 / 1201</td></tr><tr><td>Sum</td><td></td><td>8348</td><td>3999</td></tr></tbody></table>

<figure><img src="image_102.png" /><figcaption>**Fig. 3.** Crustal travel time correction terms for an average teleseismic shear wave slowness of 7 s/°. The correction time is relative to the *iasp91* reference Earth model (Kennett and Engdahl, 1991) and includes station elevation (Weidle et al., 2010) and Moho topography (Stratford et al., 2009) effects on the travel time.</figcaption></figure>

In doing so, the knowledge of the absolute velocities at depth is removed whereby only relative velocity changes will be relevant in the following. This procedure results in first-order crust corrected, weighted, relative travel time residuals, just called residuals in the following. These residuals are caused by seismic velocity anomalies underneath the station network and we use them as input data to invert for a 3D image of the shear wave velocity contrasts.

A comparison of the R and T component residual pattern at the stations shows that there is hardly any difference between both. Both components display a quite similar distribution of early and late arrival times (Fig. 4). E.g. at station NWG28, which is located in the southern part of the MAGNUS network (Fig. 1), S-waves from NW, N, NE, E and SE and SKS-waves are delayed (reddish diamonds) relative to S-waves from SW and W (bluish diamonds) both on R and T component recordings (Fig. 4). The only exception is found at stations NWG18, NWG32, NC602, NAO01 and NB201 for events in a backazimuthal range of 200° to 270° and with about 11–13 s/° slowness (epicentres at the South Atlantic ridge, see Fig. 2). These stations are all placed just north of the Oslo Graben (Fig. 1). There the 36 R component arrivals are faster compared to the 45 T component arrivals. This is the only observation of a dependence on polarisation on the residuals and hence a very weak indication of anisotropy in our entire residual dataset. However, all other measurements (approx. 98% of the dataset) do not show any difference between R and T component residuals. This overwhelming coincidence is astonishing, because a shear wave splitting analysis of SKS phases revealed clear shear wave splitting below the Southern Scandes (Roy and Ritter, in press). The SKS splitting, which is depth- and backazimuth dependent, must be related to a complex anisotropic structure at depth and seems to average out in our residual measurements. The reason for this difference between residuals measurements and SKS splitting is not understood yet, but also observed in the Eifel (Keyser et al., 2002). However, we infer that an isotropic travel time
<figure><img src="image_103.png" /><figcaption>**Fig. 4.** Radial (left column) and transverse (right column) component residuals at selected stations of the MAGNUS network (Fig. 1). The azimuthal direction depends on the backazimuth of the events (Fig. 2). The radial axis represents slowness in 2 s/° intervals. The travel time residuals are relative to their average across the network, weighted by quality and crust corrected (see Section 3).</figcaption></figure>

inversion for a 3D tomography can be justified in our case due to the prevailing isotropic residuals.

## 3.2. Residual pattern beneath the MAGNUS network

For the tomography we choose for each event phase the residuals of either the R or T component recordings depending on the better data quality which is estimated from the weighting procedure (see above). In this way only the most reliable residual measurements are used for the tomography. In Fig. 5 we summarise all these 3999 residuals (Table 1) in a spatial view. These are averaged per station and colour-coded per station (circles) as well as smoothly interpolated between the station sites. This residual plot displays very obvious large-scale trends in the residual pattern and indicates their

related shear wave velocity anomalies ($\Delta v_s$) at depth: waves arriving from the NE or SE (epicentres in the NW Pacific or Asia, Fig. 2) arrive earlier at stations east of ~10°E on the Baltic Shield (Fig. 5a and b) whereas waves from SW (South American epicentres) arrive mostly delayed at the western Baltic Shield (Fig. 5c). The four stations NWG19, NWG26, NWG34 and HFC2 in the SE of the MAGNUS network show only fast arrivals, with the exception of NE incident waves that propagate underneath the Southern Scandes.

In the centre of the MAGNUS network, inside the Southern Scandes, mainly delayed arrivals are measured (Fig. 5). Since first order travel time effects due to surface and Moho topography are already excluded, mantle anomalies are the main reason for the displayed residuals assuming that errors in the crustal model as well as lateral heterogeneities in the crust cause only second order
<figure><img src="image_104.png" /><figcaption>**Fig. 5.** Spatial distribution of direction-dependent residuals. The travel time residuals are relative to their average across the network, weighted by quality and crust corrected (see Section 3). a) Average residuals for events from NE backazimuth (see arrow in lower left corner), b) average residuals for events from SE backazimuth, c) average residuals for events from SW backazimuth, d) average residuals for events from NW backazimuth and e) average residuals for events with steeply incident core phases.</figcaption></figure>

travel time effects. In the south of the network, residuals are delayed if the waves propagated from N and NE and they arrive earlier for waves approaching from SW. The stations in the NW of the network (NWG03, NWG04, NWG07 and NWG09) have mostly delayed arrivals with the biggest delays for waves propagating from SE (Fig. 5b).

The overall residual pattern in Fig. 5 has one major common feature: waves propagating underneath the Southern Scandes are always delayed whereas waves travelling underneath the Baltic Shield arrive early. Along the west coast of Southern Norway waves arriving from the Atlantic side are faster at the southern stations compared to the northern stations. The backazimuthal direction has a greater influence on the residuals compared to the slowness (or incidence angle) as can be seen for example in Fig. 4.

The transition between slow seismic velocity at depth in the west and fast seismic velocity at depth in the east can be visualised best by plotting only residuals of steeply arriving SKS phases (Fig. 5e). The maximum residual contrast reaches nearly 4 s. The change between fast and delayed arrivals occurs at 10°E longitude in Southern Norway. This change is not very sharp but appears as transition between about 9°E and 11°E in a 100 km wide zone. The biggest SKS delays do not coincide with the highest surface topography of the Southern Scandes. These are found towards SW Norway at the North Sea coast. The area around the Oslo Graben is characterised by fast arrivals relative to the Southern Scandes. On the Baltic Shield the SKS phases arrive increasingly earlier towards east (Fig. 5e).

The scatter of the picked arrival times can be assessed in Fig. 4. The scatter is mainly due to picking uncertainties and timing problems at the station. Recordings with known timing problems due to missing GPS time synchronisation had been excluded from this analysis. The remaining data scatter can be assessed by comparing residuals from

epicentre clusters, e.g. Sumatra, Japan or South America. The residuals from such clusters show similar values in Fig. 4 and their scatter is mostly less than 0.1 s. Larger residuals variations may be related to identified picking uncertainties which are treated by weighting the data quality (see Sections 3.1 and 4.1). Thus picks with possible larger uncertainties are down-weighted during the inversion.

Based on the residual pattern we expect the following shear wave velocity anomalies $\Delta v_s$ in the upper mantle: positive $\Delta v_s$ in the east of the MAGNUS network and hence faster $v_s$ underneath the Baltic Shield in contrast to negative $\Delta v_s$ in the centre of the network and hence a pronounced low $v_s$ anomaly underneath the Southern Scandes. Towards south and southwest of the network negative $\Delta v_s$ are expected. Compared to previous residual datasets that we analysed in Europe and Africa our residuals at the MAGNUS network in Norway display a very simple low-velocity anomaly underneath its centre.

### 3.3. Residual pattern of former studies in Norway and adjacent areas

Our S-wave residuals can be compared to published P-wave residuals. Using teleseismic P-waves, Berteussen (1975) found relative residuals of $\pm 0.8$ s across the NORSAR stations and a general trend for late arrivals in the west and early arrivals in the east. He also demonstrates that a variation in Moho depth alone can only explain about 20% of these observations and would require an increase of more than 10 km in Moho depth towards west (that is thicker crust below the south-eastern part of the Southern Scandes relative to the Baltic Shield in the border region between Norway and Sweden). Thus Berteussen (1975) realised that more complicated seismic models were necessary to explain the residuals and Aki et al. (1977)
presented such a 3D model in their pioneering work on seismic tomography.

The ±0.8 s P-wave residual by Berteussen (1975) corresponds to roughly ±1.4 s for S-wave residuals based on a $v_p/v_s$ of about 1.73. We find weaker residual changes (<1 s) across the NORSAR stations (Fig. 5) which might be caused by our computation of weighted residuals or due to a reduced $v_p/v_s$ in Southern Norway. Waves propagating from SE to NW and nearly vertically incident SKS-waves show residual changes of only 0.6–0.75 s across the NORSAR stations. We observe later arrivals at the western NORSAR stations compared to the eastern NORSAR stations for vertically incident waves (Fig. 5e) as well as delayed arrivals at the north-western NORSAR stations relative to the south-eastern NORSAR stations for waves propagating from SE to NW (Fig. 5b). In contrast, waves with incidence from NW to SE show faster arrivals (approx. 1 s) at northern and western NORSAR stations compared to the southern and eastern NORSAR stations (Fig. 5d). However, no significant residual variations (<0.3 s) between the NORSAR stations are found from waves propagating from NE to SW and vice versa (Fig. 5a, c). The residual variations of nearly vertically incident waves, which are caused by upper mantle velocity perturbations directly underneath the station network, show the same residual distribution pattern as Berteussen (1975): late arrivals in the west and early arrivals in the east of the NORSAR network. However the residual variation across the NORSAR stations is smaller (0.6 s) compared to the observation of Berteussen (1975).

Analysing teleseismic P-waves, Medhus et al. (2009) found delayed arrivals underneath Southern Norway compared to the Baltic Shield in the east. The sharp transition between delayed and faster arrivals appears in the area of the Oslo Graben and the residual variations are smaller than 1 s. Medhus et al. (2009) assumed that a change in lithospheric thickness might cause the observed residual pattern. According to this idea, an increase in lithospheric thickness towards the Baltic Shield and therefore faster velocities underneath the Baltic Shield compared Southern Norway might lead to faster arrivals in the area of the Baltic Shield and delayed arrivals in Southern Norway. The residual pattern of Medhus et al. (2009) is in good agreement with our S-wave residual pattern (Fig. 5), although the transition zone between faster and delayed arrivals is located further west in our S-wave residual pattern compared to the P-wave residual pattern of Medhus et al. (2009).

A comparison to other travel time analyses in the region shows related relative residual patterns: Shomali et al. (2006) analysed P- and S-wave travel time data of the TOR experiment and found a residual pattern linked to the Trans European Suture Zone (TESZ): South of TESZ the residuals are positive (late arrivals) whereas negative residuals (early arrivals) are observed north of the TESZ. Eken et al. (2007) analysed P-wave travel times which were recorded at the Swedish National Seismological Network (SNSN). Their pattern of averaged relative residuals illustrates early wave arrivals in the southern part of the network compared to later wave arrivals in the northern part of the network (Eken et al., 2007) and their resulting mean P-wave velocity pattern illustrates higher seismic velocity in the southern part of Sweden compared to northern Sweden. Eken et al. (2007) presume that the change of the sign of the travel time residuals coincides with the Archean–Proterozoic boundary. Although these travel time residuals are not directly comparable to our study (different networks and mean values), they show that early wave arrivals are observed on the Baltic Shield in southern Sweden.

# 4. Travel time tomography

## 4.1. Inversion method

To compute the teleseismic travel time tomography model and hence determine the 3D upper mantle $v_s$ structure, we use the non-linear JI-3D inversion method (Jordan, 2003). JI-3D is a further

development of the ACH-type inversion scheme (Aki et al., 1977). For a detailed description of the JI-3D inversion method see Jordan (2003) or Wawerzinek et al. (2008).

The model space is composed of layers which are subdivided into variably wide blocks to take into account the station distribution, the ray density and criss-crossing rays (Jordan, 2003). In the centre of the model the block size corresponds roughly to the station distance whereas the block size is enlarged at the model edges where only few criss-crossing rays are observed (Fig. 6a). Thereby, the values of the diagonal elements of the resolution matrix are almost the same which leads to a stable inversion result (Jordan, 2003). The ray paths are computed iteratively by using a 3D ray tracer (Steck and Prothero, 1991). To include a priori information, e.g. uncertainties of the residuals or the allowed variation of the individual model parameters, a Bayes algorithm is applied (Zeyen and Achauer, 1997). To account for uncertainties of onset picking and therefore residuals picks of quality class "C" are weighted as 0.25, "B"-picks as 0.5 and "A"-picks as 1. Comparable to model damping the admitted model parameter variation is set variably for each block (Jordan, 2003). For the starting model the admitted parameter variation is determined by performing a trade-off analysis, during iterative inversion the admitted model parameter variation is computed by the inversion code (Jordan, 2003). To account for crustal 3D travel time effects, which are not eliminated by applying the 1D crustal correction (see Section 3), the model parameters of the crustal blocks are not totally fixed but only small variations are allowed (Martin et al., 2005).

The crust and upper mantle parameterisation underneath our study area (Table 2) is subdivided into 9 layers down to 500 km depth (Fig. 6b) and 272 blocks with a maximum horizontal model space extension of 1200 km N–S and 1200 km E–W. The vertical model extension is chosen in dependence of the network extension and the observed criss-crossing rays. In addition, the relative P-wave tomography study of Medhus et al. (2012) image no significant seismic velocity contrast underneath 400 km depth. Therefore, the vertical model extension was set to 410 km depth to image the complete upper mantle velocity structure. The upper first and second layers contain the crust and account for station site effects (first layer) and not yet corrected 3D crustal effects (second layer). The deepest layer (410–500 km) is required to absorb effects of underlying shear wave velocity anomalies or heterogeneities (e.g. Evans and Achauer, 1993). Thereby the absorbing effect is dependent on the spatial resolution in the overlying layers. However, these three layers (the two crustal layers and the deepest layer) will not be interpreted later on, because they are mainly used to absorb inversion artefacts. The interpretation will only include the upper mantle structure from 35 km to 410 km depth (layers 3–8, Table 2).

The starting model comprises the average $v_s$ of the iasp91 Earth reference model (Kennett and Engdahl, 1991) as background velocities. To account for non-linear effects the inversion is performed iteratively and we achieve a stable result by computing two iterations which leads to a residual variance reduction of almost 63% which means that the data are underfitted. The offset and averaging technique (Evans and Zucca, 1988) is performed to reduce the influence of the block parameterisation. Thereby the parameter grid is shifted along the block edges with an offset of 1/3 of the smallest block size and for each of these nine parameter grids an individual inversion is performed. To obtain the final inversion result the individual nine inversion results are averaged afterwards (Evans and Achauer, 1993).

## 4.2. Inversion results

Our inversion results are presented as horizontal (Fig. 7) and vertical (Figs. 8 and 9) cross sections through the model with $\Delta v_s$ anomalies. Following the description of the model we discuss its significance based on reconstruction tests. As we use relative residuals across the station network with a mean value of zero, no knowledge
<figure><img src="image_105.png" /><figcaption>**Fig. 6.** Ray distribution and parameterisation of the tomography model. a) Horizontal cross section through the model layer at 180–250 km depth with ray path segments and the block parameterisation (Wawerzinek, 2012). b) Vertical W–E cross section between 60°N and 61°N with ray geometry to the seismic stations and the vertical layer parameterisation.</figcaption></figure>

of absolute velocities can be given, therefore only velocity contrasts between the anomalies are meaningful. Hence, absolute velocity perturbations cannot be determined but these are mostly in the same percentage range when compared with absolute velocity models. Five major $\Delta v_s$ anomalies can be identified in the tomography model:

I. The shallow mantle (<180 km depth) along the Norwegian coast has negative $\Delta v_s$ of about –2% compared to the continental region (35–180 km depth slices in Fig. 7 and W–E cross sections in Fig. 9).

II. At shallow depth (35–75 km depth slice in Fig. 7) the mantle underneath the Oslo Graben area has positive $\Delta v_s$ of about 1.5% relative to the surrounding lithosphere.

III. The mantle below the Baltic Shield has positive $\Delta v_s$ compared to the mantle underneath the Southern Scandes; this boundary is not vertical but dips underneath the Baltic Shield between 35 km and 180 km depth (W–E cross sections in Fig. 9).

IV. Below the Southern Scandes there is no specific $v_s$ anomaly in the upper 180 km depth (35–180 km in Fig. 7 and cross sections in Fig. 9) besides the general contrast to the Baltic Shield lithosphere which is characterised by higher $v_s$.

V. In the lowermost upper mantle, at about 250–410 km depth, there is a clear low $v_s$ anomaly with a high contrast in $v_s$ of maximum 3% relative to the mantle underneath the Atlantic towards SW and the Baltic Shield in the east.

<table><caption>**Table 2**

Parameterisation of the inversion model; the used background shear wave velocities $v_s$ are taken from iasp91 (Kennett and Engdahl, 1991).</caption><thead><tr><th>Layer</th><th>Depth range in km</th><th>Node depth in km</th><th>v<sub>s</sub> background in km/s</th></tr></thead><tbody><tr><td>1</td><td>0–2</td><td>0</td><td>3.36</td></tr><tr><td>2</td><td>2–35</td><td>34</td><td>3.75</td></tr><tr><td>3</td><td>35–75</td><td>36</td><td>4.47</td></tr><tr><td>4</td><td>75–120</td><td>95</td><td>4.49</td></tr><tr><td>5</td><td>120–180</td><td>150</td><td>4.51</td></tr><tr><td>6</td><td>180–250</td><td>215</td><td>4.53</td></tr><tr><td>7</td><td>250–330</td><td>290</td><td>4.66</td></tr><tr><td>8</td><td>330–410</td><td>409</td><td>4.87</td></tr><tr><td>9</td><td>410–500</td><td>411</td><td>5.07</td></tr><tr><td>Half space</td><td></td><td></td><td>5.26</td></tr></tbody></table>

### 4.3. Reconstruction of the inversion

In order to access the reliability of the inversion results we first compare the spatial distribution of the measured residuals with the anomalies in the inversion model. Then we evaluate numerical reconstruction tests in order to test the confidence and limitations of our 3D model as well as the smearing of anomalies along the ray paths.

A comparison of the inversion result (Fig. 7) with the distribution of residuals of nearly vertically incident waves (Fig. 5e) reveals a very good correlation at 120–410 km depth: The inversion result shows negative $\Delta v_s$ (reduced $v_s$) in the area of delayed arrivals (Southern Norway) relative to positive $\Delta v_s$ (increased $v_s$) in the area of fast arrivals (Baltic Shield). The transition between delayed and fast arrivals coincides well with the transition zone between slower and faster velocity at 180–250 km depth. Thus, the inversion result visualises a simple $v_s$ distribution as expected by analysing the distribution of the travel time residuals.

Four numerical reconstruction tests are performed in the following way: first a 3D model with synthetic $\Delta v_s$ anomalies is constructed. 3D ray tracing is used to calculate travel times in these models with the same ray distribution as for the real MAGNUS tomography. Afterwards synthetic residuals are determined. Then the synthetic residuals are inverted with the same parameterisation as the real model in order to estimate the resolution properties of the input anomalies. We display the results of a checkerboard test (CBT in Fig. 10) and three anomaly tests (AT1, AT2 and AT3 in Figs. 11, 12 and Supp. Fig. S4, respectively). Each reconstruction test displays the inversion result after applying the offset and averaging technique (Evans and Zucca, 1988) as interpolated velocity pattern. Additionally, the inversion results of the nine individual inversion models are displayed as circles in the depth slices in Figs. 10–12.

The input velocity model of the checkerboard test (CBT, Fig. 10) comprises alternating anomalies with $v_s$ amplitudes of +1.5% (blue outlines) and –1.5% (red outlines) at 75–125 km and 180–250 km depth. The synthetic input anomalies have extensions of 100 km by 100 km at 75–120 km depth and 200 km by 200 km at 180–250 km depth. The inversion result of the CBT (Fig. 10) recovers the input velocity anomalies but also illustrates the vertical smearing both upwards and downwards. At first, we note that in some model areas the single inversion results (circles) reproduce the synthetic anomalies
<figure><img src="image_106.png" /><figcaption>**Fig. 7.** Horizontal slices through the inversion result after applying the offset and averaging technique. The S-wave velocity contrast is colour-coded and the model parameterisation is indicated by block boundaries and node positions (circles). The S-wave velocity is reduced underneath Southern Norway compared to its surroundings especially to the Baltic Shield in the east. The sharp transition zone between slower and faster $v_S$ is displaced towards east with increasing depth.</figcaption></figure>

better than the final offset and averaging inversion result. However, in other areas the image of the finally averaged inversion model better represents the synthetic anomalies. This outcome demonstrates the influence of parameterisation on the inversion result. Since sharp vertical boundaries are not expected in the upper mantle and the influence of the parameterisation on the inversion result has to be mitigated, we apply offset and averaging technique (Evans and Achauer, 1993; Evans and Zucca, 1988). Second, the mapping of the synthetic anomalies is good in the centre of the model space, but there is reduced resolution towards the model edges. This variation

in resolution is caused by the ray distribution (see Fig. 6) which has a decreasing ray density and few criss-crossing rays towards the model edges. Third, the synthetic anomalies at 180–250 km depth are better resolved than those at 75–120 km depth which may be caused by the larger extension of the input anomalies at 180–250 km depth and therefore more traversed criss-crossing rays inside the parameterisation blocks. Fourth, the amplitudes of the synthetic anomalies are underestimated. The CBT indicates that the vertical resolution is reduced due to strong vertical smearing (Fig. 10) whereas the lateral resolution is reduced by applying the
<figure><img src="image_107.png" /><figcaption>**Fig. 8.** Location map of the vertical cross sections (A–H) through the inversion result which are shown in Fig. 9. The intersection points of the profiles are labelled X1–X9.</figcaption></figure>

offset and averaging technique. However, large anomalies (approx. 200 km) are well reproduced in the centre of the model, but there is a decrease in resolution towards the model edges.

In reconstruction test AT1 (Fig. 11) synthetic anomalies are pre-set which are similar to those of the actual inversion result (Fig. 7) with amplitudes of +1.5% $v_S$ (blue outlines) and -1.5% $v_S$ (red outlines) at 75–250 km depth and a W–E dipping velocity transition. Reconstruction test AT1 (Fig. 11 and Supp. Fig. S2) shows moderate vertical smearing (upwards and downwards) across two layers (e.g. no smearing into the uppermost two layers). The synthetic anomalies (75–180 km depth) are recovered very well just as the W–E transition zone between slower and faster velocity anomalies (see vertical cross sections in Supp. Fig. S2), taking into account the lateral smearing due to the application of the offset and averaging technique. The amplitudes of the synthetic anomalies are underestimated. Based on test AT1 we infer that the deep negative $v_S$ anomaly in Fig. 7 is not caused by vertical smearing, because the amplitude contrast in AT1 at 330–410 km depth (Fig. 11) is much smaller than the $v_S$ contrast in the actual model (Fig. 7).

In order to learn more about the resolution of the deep anomaly, we perform an additional reconstruction test AT2 (Fig. 12 and Supp. Fig. S3) which includes a velocity anomaly with an amplitude of -2.5% $v_S$ (red outline) at 330–410 km depth. This reconstruction test AT2 recovers the negative $v_S$ anomaly at 330–410 km depth, but its perturbation amplitude (about -1.5% $v_S$) is underestimated. Instead, an increase in amplitude (approx. ±1.5% $v_S$) is observed at 75–250 km depth due to an increase in upward vertical smearing compared to test AT1 (Fig. 11). The reconstruction tests AT1 (Fig. 11 and Supp. Fig. S2) and AT2 (Fig. 12 and Supp. Fig. S3) show that besides the W–E $v_S$ contrast a deep additional anomaly is required to explain the negative $v_S$ anomaly at 330–410 km depth and that the transition zone between faster velocity in the east and slower velocity in the west can be resolved well with our dataset.

To show that the upper $v_S$ anomaly is not caused by vertical smearing of a strong deep-seated low velocity anomaly another reconstruction test AT3 (Supp. Fig. S4) is performed. A low velocity anomaly is pre-set at 330–410 km depth with an amplitude of -2.5% $v_S$. The inversion result of this reconstruction test illustrates a well reproduced low velocity anomaly at 330–410 km depth whose amplitude is underestimated and vertical upward smearing across the model is found. However, no similarities between the

inversion result of AT3 and the actual model in Fig. 7 are observed at 35–330 km depth. This confirms our assumption that the low velocity anomaly at 75–250 km depth may not be explained as inversion artefact due to a strong deep-seated anomaly.

From the four reconstruction tests (CBT, AT1, AT2 and AT3) we conclude: The ray distribution (ray density and criss-crossing rays) causes a decrease in resolution from the centre towards the edges of the model which leads to a reduced lateral resolution close to the model edges. By applying the offset and averaging technique weak lateral smearing occurs but the influence of the parameterisation is reduced. However, large anomalies (~200 km) are reproduced well. Due to strong vertical smearing the vertical resolution is low. Nevertheless, the reconstruction test AT1 (Fig. 11) and AT2 (Fig. 12) show that larger anomalies and their transitions (Supp. Figs. S2 and S3) can be mapped well and that an additional deep-seated negative $v_S$ anomaly can be recovered. Such a negative $v_S$ anomaly at 330–410 km depth explains the observed residuals at the MAGNUS network. Both, the upper and deeper-seated anomaly can be resolved independently (AT1 and AT3). The deeper-seated anomaly might continue below 410 km depth in the mantle but this is beyond our resolution. The positive anomaly at 250–410 km depth in the southwestern part of the inversion model may partly be due to the condition for a mean value of 0% per layer. This assumption may be supported by a comparison of the reconstruction tests AT1 and AT2 (Figs. 11 and 12). A positive anomaly appears in the south-western part of the model by including a strong negative anomaly in the centre of the model at 330–410 km depth. This effect demonstrates that only velocity contrasts are meaningful in regional tomography models based on normalised residuals. In addition, the anomaly is located at the edge of our inversion volume where the resolution is low caused due to only few criss-crossing rays.

## 5. Discussion and interpretation

A teleseismic travel time tomography is performed to determine the 3D $v_S$ structure beneath Southern Scandinavia. As relative travel time residuals are used as input data, $v_S$ contrasts can be interpreted in our tomographic images of the upper mantle (Figs. 7 and 9). The observed $v_S$ structure contains a stable $v_S$ pattern at 35–410 km depth: $v_S$ is reduced underneath Southern Norway compared to its surroundings, especially the Baltic Shield in the east. The transition between slower and faster $v_S$ is sharp and it dips towards east with increasing depth. The velocity contrast is maximum 4% $v_S$ and on average up to 3% $v_S$. In Fig. 13 a schematic illustration summarises the main features which are recovered by our teleseismic shear wave analysis. The reconstruction tests (AT1 and AT2) imply that the W–E dipping transition zone as well as the deep anomaly down to 410 km depth are well resolved.

There are also global and regional tomography studies which are based on body waves (P) and surface waves and cover the area of Southern Scandinavia. Weidle and Maupin (2008) analysed group velocities of surface waves (Rayleigh and Love waves) to derive the 3D $v_S$ distribution underneath Northern Europe. They revealed a negative velocity anomaly beneath Iceland that stretches out towards Southern Norway. In their model $v_S$ is reduced up to -10% of $v_{S\_ak135}$ at 70–150 km depth (Weidle and Maupin, 2008). Across Southern Norway the velocity contrast is about 5%–7% of $v_{S\_ak135}$ (Weidle and Maupin, 2008) what is significantly larger than our modelled $v_S$ reduction of -3% to -4% of $v_{S\_iasp91}$. This difference may be caused by performing different tomography methods with diverse sensitivities and resolutions. As we determine relative travel time residuals, a regionally broad $v_S$ anomaly would be removed and only variations within the anomaly would be left over. Thereby our model can only resolve these more local $v_S$ contrasts underneath the MAGNUS network and may underestimate a broader anomaly.
<figure><img src="image_108.png" /><figcaption>**Fig. 9.** Vertical slices through the inversion result after applying the offset and averaging technique. The S-wave velocity contrast is colour-coded with the same colour scale as in Fig. 7. Underneath Southern Norway the S-wave velocity is reduced compared to the Baltic Shield in the east. In west–east direction (profiles DD', EE', FF', GG' and HH') a sharp transition zone between slower and faster $v_S$ is visible which is displaced towards east with increasing depth.</figcaption></figure>

By analysing the phase velocity of Rayleigh waves at the MAGNUS stations, Maupin (2011) determined an average 1D $v_S(z)$ profile underneath Southern Norway. In this model the upper mantle $v_S$ is reduced compared to the reference Earth model ak135 (Kennett et al., 1995). At 150 km depth $v_S$ is 2% slower than $v_{S\_ak135}$ (Maupin, 2011) and approximately 0.3 km/s slower than $v_S$ expected in cratonic areas (Pedersen et al., 2009) as found underneath Sweden (Maupin, 2011). Maupin's (2011) Rayleigh wave model is in good agreement with our tomographic body wave model which shows a velocity change of –2% $v_{S\_iasp91}$ below Southern Norway and a velocity contrast of up to 4% $v_{S\_iasp91}$ between the Southern Norway in the west and the Baltic Shield in the east.

P-wave tomography models (Bijwaard and Spakman, 1999; Spakman published in Jones et al., 2010; Medhus et al., 2012) also show a reduction of mantle $v_P$ underneath Southern Norway compared to the Baltic Shield in the east. The global model of Bijwaard

and Spakman (1999) contains a $v_P$ reduction of –2% $v_{P\_ak135}$ at 300 km depth, whereas a more recent study of Spakman (published in Jones et al., 2010) image a $v_P$ reduction of 0.2–0.3 km/s relative to the Baltic Shield. In both studies $v_P$ increases towards the Baltic Shield and the transition between faster and slower velocity is located in the area of the Oslo Graben (Bijwaard and Spakman, 1999; Spakman published in Jones et al., 2010). Medhus et al. (2012) performed an absolute and a relative $v_P$ travel time tomography whose $v_P$ distributions show a strong reduction of $v_P$ underneath Southern Norway compared to the Baltic Shield. At a depth range of 100–300 km the $v_P$ contrast is up to 3%–4% of $v_{P\_ak135}$ in the absolute velocity model and also 3%–4% of $v_{P\_iasp91}$ in the relative velocity model (Medhus et al., 2012).

In contrast to the above mentioned studies, Ritsema et al. (2004) and Megnin and Romanowicz (2000) found no indication of a velocity reduction beneath Southern Norway compared to its surroundings.
<figure><img src="image_109.png" /><figcaption>**Fig. 10.** Reconstruction test CBT (checkerboard pattern test). The pre-set synthetic anomalies (blue outline: +1.5% $v_s$, red outline: -1.5% $v_s$) are placed in the velocity model at 75–120 km and 180–250 km depth. The interpolated S-wave velocity perturbations visualise the corresponding inversion result after applying the offset and averaging technique. Additionally, the circles show the block-wise inversion results of the individual nine offset inversion models before averaging. There is strong vertical smearing which leads to reduced vertical resolution within the model space.</figcaption></figure>

This might be due to a lower resolution in these global tomography studies caused by the sparse covering with permanent stations in Scandinavia.

The conclude, the seismic velocity ($v_p$ and $v_s$) studies, which are based on dense regional station networks, show a major common

feature: the reduced seismic velocity underneath the Scandinavian Mountains and relatively increased seismic velocity underneath the Baltic Shield (this study; Bijwaard and Spakman, 1999; Maupin, 2011; Weidle and Maupin, 2008; Spakman published in Jones et al., 2010; Medhus et al., 2012). However, the seismic velocity contrast
<figure><img src="image_110.png" /><figcaption>**Fig. 11.** Reconstruction test AT1 (anomaly pattern test 1). The pre-set synthetic anomalies are placed in the velocity model at 75–250 km depth (blue outline: +1.5% vS, red outline: -1.5% vS) and these represent a dipping W-E contrast as found in the real model (Fig. 7). The interpolated S-wave velocity perturbations visualise the corresponding inversion result after applying the offset and averaging technique (vertical cross sections see Supp. Fig. S2). Additionally, the circles show the inversion results of the individual nine offset inversion models before averaging. This reconstruction test illustrates a good reproduction (position and amplitude) of the synthetic anomalies, although weak vertical smearing is visible across two layers.</figcaption></figure>

and especially the amplitude of the velocity reduction vary as well as the exact location of the anomaly. Differences between these $v_P$ and $v_S$ structures are related mainly to the perturbation amplitude which might be caused by analysing P- and S-waves in different period

ranges, by different ray paths or propagation directions, different resolution properties and different model space parameterisations. Additionally, a variation in the $v_P/v_S$ ratio might affect differences between $v_P$ and $v_S$ tomographic structures.
<figure><img src="image_111.png" /><figcaption>**Fig. 12.** Reconstruction test AT2 (anomaly pattern test 2). The pre-set synthetic anomalies are placed in the velocity model at 75–250 km depth (blue outline: + 1.5% vS, red outline: – 1.5% vS) in addition to a stronger synthetic anomaly (red outline: – 2.5% vS) at 330–410 km depth. These anomalies represent vS contrasts as found in the real model (Fig. 7). The interpolated S-wave velocity distribution visualises the corresponding inversion result after applying the offset and averaging technique (vertical cross sections see Supp. Fig. S3). Additionally, the circles show the inversion results of the individual nine offset inversion models before averaging. This reconstruction test illustrates a good reproduction (position and amplitude) of the synthetic anomalies, although some vertical smearing is visible.</figcaption></figure>

In Scandinavia other regional tomography studies also imaged deep reaching velocity contrasts in the upper mantle across major tectonic boundary zones. E.g. the Tornquist Zone, the suture between the Proterozoic NE Europe and the Phanerozoic Central Europe, is

characterised by up to 5% P-wave velocity ($v_p$) contrast by non-linear tomography (Shomali et al., 2002). Sandoval et al. (2004) recovered an inclined velocity contrast at the suture of the Archean Karelian terrane and the Proterozoic Svecofennian terrane in Finland
<figure><img src="image_112.png" /><figcaption>**Fig. 13.** Vertical cross section through a schematic model of the shear wave structure underneath Southern Scandinavia at 60.5°N (Wawerzinek, 2012). At 35–410 km depth, $v_S$ is continuously reduced underneath Southern Norway (west of 11°E, $-Δv_S$, light grey) compared to the Baltic Shield ($+Δv_S$, dark grey). In addition, the schematic model contains discontinuities which were recovered by common conversion point stacks of S-receiver functions (Wawerzinek, 2012). Several transitions to low velocity zones are visible in a depth range of 50–250 km (dashed zone), but a unique identification of the lithosphere–asthenosphere-boundary is not possible. Observed depth variations of the 660 km discontinuity (lower grey dashed line) are marked by arrows which are not scaled correctly. The 660 km discontinuity is deeper underneath Southern Norway compared to the Baltic Shield. In the area of Southern Norway, the 410 km discontinuity is not detected in our S-RF (upper grey dashed line) which might be caused by the observed $v_S$ anomaly at 330–410 km depth.</figcaption></figure>

and Russia. Further in the west, in Sweden, Eken et al. (2007) determined a $v_P$ contrast of 2%–3% across the Proterozoic–Archean lithospheric domains. Their tomographic image shows a slab-like shape. In Eken et al. (2008) P- and S-wave velocity contrasts are displayed for the Proterozoic–Archean suture as well as for the transition between the Transcandinavian Igneouos Belt and the Svecofennian domain. These examples and our new results indicate that very old collisional processes from Proterozoic and even Archean times are still preserved in the upper mantle.

In the following, we discuss several reasons that may cause the velocity reduction beneath Southern Norway relative to the Baltic Shield: a variation in lithospheric thickness, a variation in mantle composition, a temperature variation, the presence of hydrated rocks or partial melts or a combination of some or all these possible causes. For petrophysical impacts on seismic velocity we refer to e.g. Cammarano et al. (2003), Goes et al. (2004) or Ritter (2007). The following interpretations may be valid for our $v_S$ mantle model:

I. We infer an increase in lithospheric thickness from the Scandes towards the Baltic Shield as presented e.g. by Plomerová et al. (2008). They compiled regional body-wave data to create a model of the lithospheric thickness of Fennoscandia. In their model the lithospheric thickness increases from the western coast of Southern Norway (<90 km) towards the Baltic Shield in the east (>130 km; Plomerová et al., 2008). Shear wave velocity is slightly faster in the lower lithosphere compared to the asthenosphere (Eaton et al., 2009) and slightly faster inside older shield lithosphere compared to younger Phanerozoic orogens (Jones et al., 2010). Thus a lithospheric thickness variation might explain the lateral velocity variations at a depth range of 75–250 km (Fig. 9). Thereby, the observed slower $v_S$ underneath Southern Norway could be interpreted as $v_S$ in a shallower asthenosphere relative to the faster $v_S$ in a deeper lithospheric root underneath the Baltic Shield. This lithosphere–asthenosphere transition seems not to be a steep boundary, but it is inclined with a dip towards east (Figs. 9 and 13).

II. The composition of the upper mantle certainly differs between the Phanerozoic Southern Norwegian region and the Baltic Shield in Southern Sweden, because the past geodynamic evolution of Southern Norway was strongly influenced by the Caledonian orogeny, the Cenozoic opening of the Northern Atlantic and the Neogene basin building in the North Sea (Ramberg et al., 2008; Torsvik and Cocks, 2005). In the mantle a related variation in mineral composition, fluid and melt content might affect the seismic velocity structure beneath Southern Scandinavia. However, we expect only small seismic

velocity changes due to variations in the mantle composition (Cammarano et al., 2003; Sobolev et al., 1997).

III. Since temperature variations cause larger seismic velocity perturbations than compositional variations in an anhydrous mantle (Goes et al., 2004; Sobolev et al., 1997), temperature differences could perhaps better explain our velocity perturbations in the upper mantle. A positive temperature anomaly, which may involve partial melts, leads to a reduction of the elastic moduli and therefore to a reduction of the seismic velocity (Berckhemer et al., 1982). If the $v_S$ reduction is caused by a pure temperature anomaly underneath Southern Norway, a temperature increase of up to 200–330 K (Cammarano et al., 2003) is required to explain a velocity contrast of up to 4% $v_S$ at 75–250 km depth based on an average mineral composition of the mantle. A deeper anomaly at 330–410 km depth, which has a velocity contrast of about 4% $v_S$, may be explained by a pure temperature anomaly of about 400 K (Cammarano et al., 2003). This relatively high temperature anomaly might be reduced by taking into account hydrous minerals at 330–410 km depth. Thus a smaller temperature anomaly of less than 300 K is enough to explain the observed contrast of about 4% $v_S$ (Ritter, 2007). However, the presence of fluids, partial melts or hydrous minerals inside the deep anomaly below Southern Norway cannot be confirmed at the moment.

In summary, our S-wave tomography study contributes the 3D shear wave velocity structure in the upper mantle to the on-going discussion on the origin of the Scandinavian Mountains. Especially our results should help to define the possible dynamic contribution of the mantle structure. The improved knowledge of the $v_S$ structure of the upper mantle is an important proxy for the determination of the driving forces for the uplift and high topography. If one accepts increased temperature and the presence of shallow asthenospheric material underneath the mountains as reason for the seismic velocity reduction, then recently induced dynamic topography (Braun, 2010) may be an option to explain at least a fraction of the present topography in Southern Norway. However, the mantle domain with reduced seismic velocities (S-wave model in Fig. 7 and P-wave model in Medhus et al., 2012) is wider than the area of high topography, thus this low velocity mantle domain may affect the lithosphere in Southern Scandinavia across a wider region (Medhus et al., 2012). The uplift mechanism of Southern Norway can be better tested now by performing geodynamic modelling to explore different hypothesis, mechanisms and processes based on proxies from seismology. To verify or exclude different processes the observed shear wave velocity structure has to be taken into account as boundary condition. In the second part of *TopoScandiaDeep*, such a geodynamic modelling will be
performed to determine the cause of uplift in Southern Norway (Weidle
et al., 2010).

**6. Conclusions**

To determine the 3D shear wave velocity structure underneath Southern Scandinavia, we analyse S-waves of 128 teleseismic events which were recorded at the MAGNUS network in Southern Norway. In the first step, we compute relative travel time residuals (Section 3) by taking into account 1D travel time effects due to the known crustal structure. These residuals are caused by velocity perturbations underneath the station network within the upper mantle. The resulting residual pattern of nearly vertically incident waves shows a simple residual distribution: late arrivals are visible in Southern Norway and early arrivals are observed in the area of the Oslo Graben and the Baltic Shield. The travel time difference is maximum 4 s across the MAGNUS network. S-wave residuals with flat incident ray paths depend significantly on backzimuth: the area of late arrivals is shifted systematically in the propagation direction of the teleseismic waves. This leads to the prediction that a low shear velocity zone underneath Southern Norway is the cause of the observed residual pattern.

In the second step, we compute a non-linear travel time tomography (Section 4) by using the JI-3D inversion code (Jordan, 2003). Taking into account pick uncertainties and data quality, weighted residuals are used as input data. The resulting *v*<sub>S</sub> structure shows a continuous velocity pattern: slower *v*<sub>S</sub> in the west compared to the east. The transition between fast and slow velocity dips towards east with increasing depth (Fig. 13). Below 180 km depth the area of slower velocity comprises Southern Norway compared to the area of faster velocity underneath the Baltic Shield. The low *v*<sub>S</sub> anomaly reaches as deep as the mantle transition zone at 410 km depth. Because relative travel time residuals are used as input data for the tomography, only the velocity contrast can be interpreted which is up to 4% of *v*<sub>S</sub> and on average 3% of *v*<sub>S</sub> across Southern Scandinavia. Four reconstruction tests (Section 4.3) indicate a good resolution in the centre of the model space, e.g. the area of the lateral velocity transition. Furthermore, the reconstruction tests suggest that the observed low velocity anomaly at 330–410 km depth can be recovered and that it is not an inversion artefact.

Finally, the observed 3D v_S structure can be explained by laterally
varying mantle properties. An increase in lithospheric thickness
towards east can be related with different lithologies and/or lateral
ambient temperature conditions. A shallower asthenosphere beneath
Southern Norway would cause a positive temperature anomaly rela-
tive to colder lithosphere of the Baltic Shield in the east. In the lower-
most upper mantle the deep-seated low v_S anomaly might be due to a
combination of a temperature anomaly with fluids, partial melts
and/or hydrous minerals at 330–410 km depth and might be due
to deeper reaching geodynamic processes in the mantle transition
zone or even the lower mantle.

The variable upper mantle structure indicates that mantle pro-
cesses might influence or even support a Neogene uplift of Southern
Norway. A more detailed geodynamic modelling analysis will be
conducted including our upper mantle v_S structure. This analysis of
possible uplift mechanisms will provide more insights on the regional
uplift processes such as dynamic topography in Southern Norway.

**Acknowledgements**

MAGNUS waveforms were recorded with the mobile KArlsruhe
BroadBand Array of the Karlsruhe Institute of Technology, Germany
as well as with permanent stations of the NORSAR array and the
Norwegian National Seismological Network. Financial support for
the MAGNUS experiment was provided by the Universities of Aarhus,
Copenhagen, Karlsruhe and Oslo as well as NORSAR. This work has
been done in the framework of the ESF EUROCORES TOPO-EUROPE

Programme 07-TOPO-EUROPE-FP-014: The Scandinavian mountain chain: deep processes (TopoScandiaDeep). We acknowledge financial support from the Deutsche Forschungsgemeinschaft (RI1133/8). We thank the many helping hands during the experiment (station deployment and station service) as well as data preparation. We thank the TopoScandiaDeep working group for helpful discussions. Figures were created with Generic Mapping Tools (Wessel and Smith, 1998) and in MATLAB®. Data pre-processing was done in MATLAB® using KaSP-Toolbox (Groos, 2010), S onset picking was done in SeismicHandler (Stammler, 1993) and the inversion was done by using JI-3D inversion code (Jordan, 2003). We thank the reviewers W. Spakman and T. Eken for helpful comments to improve the manuscript.

**Appendix A. Supplementary data**

Supplementary data to this article can be found online at http://dx.doi.org/10.1016/j.tecto.2012.12.033.

**References**

Aki, K., Christoffersson, A., Husebye, E.S., 1977. Determination of the three-dimensional seismic structure of the lithosphere. Journal of Geophysical Research 82, 277–296.

Anell, I., Thybo, H., Artemieva, I.M., 2009. Cenozoic uplift and subsidence in the North Atlantic region: geological evidence revisited. Tectonophysics 474, 78–105.

Berckhemer, H., Kampfmann, W., Aulbach, E., Schmeling, H., 1982. Shear modulus and Q of forsterite and dunite near partial melting from forced oscillation experiments. Physics of the Earth and Planetary Interiors 29, 30–41.

Berteussen, K.A., 1975. Crustal structure and P-wave travel time anomalies at NORSAR. Journal of Geophysics 41, 71–84.

Bijwaard, H., Spakman, W., 1999. Tomographic evidence for a narrow whole mantle plume below Iceland. Earth and Planetary Science Letters 166, 121–126.

Braun, J., 2010. The many surface expressions of mantle dynamics. Nature Geoscience 3, 825–833.

Cammarano, F., Goes, S., Vacher, P., Giardini, D., 2003. Inferring upper-mantle temperatures from seismic velocities. Physics of the Earth and Planetary Interiors 138, 197–222.

Cloetingh, S.A.P.L., Ziegler, P.A., Boogaard, P.J.F., Andriessen, P.A.M., Artemieva, I.M., Bada, G., van Baalen, R.T., Breekman, F., Ben-Avraham, Z., Brun, J.-P., Bunge, H.P., Burov, E.B., Carbonell, R., Facenna, C., Friedrich, A., Gallert, J., Green, A.G., Heidbach, O., Jones, A.G., Matenco, L., Mosar, J., Oncken, O., Pascal, C., Peters, G., Sliaupa, S., Soesoo, A., Spakman, W., Stephenson, R.A., Thybo, H., Torsvik, T., de Vicente, G., Wenzel, F., Wortel, M.J.R., TOPO-EUROPE Working Group, 2007. TOPO-EUROPE: the geoscience of coupled deep Earth-surface processes. Global and Planetary Change 58, 1–118.

Douglas, A., Bowers, D., Young, J.B., 1997. On the onset of P seismograms. Geophysical Journal International 129, 681–690.

Eaton, D.W., Darbyshire, F., Evans, R.L., Grütter, H., Jones, A.G., Yuan, X., 2009. The elusive lithosphere–asthenosphere boundary (LAB) beneath cratons. Lithos 109, 1–22.

Eken, T., Shomali, Z.H., Roberts, R., Bödvarsson, R., 2007. Upper-mantle structure of the Baltic Shield below the Swedish National Seismological Network (SNSN) resolved by teleseismic tomography. Geophysical Journal International 169, 617–630.

Eken, T., Shomali, Z.H., Roberts, R., Hieronymus, C.F., Bödvarsson, R., 2008. S and P velocity heterogeneities within the upper-mantle below the Baltic Shield. Tectonophysics 462, 109–124.

Engdahl, E.R., van der Hilst, R., Buland, R., 1998. Global teleseismic earthquake relocation with improved traveltimes and procedures for depth determination. Bulletin of the Seismological Society of America 88, 722–743.

Evans, J.R., Achauer, U., 1993. Teleseismic velocity tomography using the ACH method: theory and application to continental-scale studies. In: Iyer, H.M., Hirahara, K. (Eds.), Seismic Tomography: Theory and Practice. Chapman and Hall, London, pp. 319–360.

Evans, J.R., Zucca, J., 1988. Active high-resolution seismic tomography of compressional wave velocity and attenuation structure at Medicine Lake volcano, northern California Cascade range. Journal of Geophysical Research 93, 15015–15036.

Gabrielsen, R.H., Faleide, J.I., Pascal, C., Braathen, A., Nystuen, J.P., Etzelmuller, B., O'Donnell, S., 2010. Reply to discussion of Gabrielsen et al. (2010) by Nielsen et al. (this volume): latest Caledonian to present tectonomorphological development of southern Norway. Marine and Petroleum Geology 27, 709–723.

Goes, S., Cammarano, F., Hansen, U., 2004. Synthetic seismic signature of thermal mantle plumes. Earth and Planetary Science Letters 218, 403–419.

Groos, J., 2010. Broadband seismic noise: Classification and Green's function estimation. Dissertation, Geophysical Institute, Karlsruhe Institute of Technology (KIT).

Japsen, P., Chalmers, J.A., 2000. Neogene uplift and tectonics around the north Atlantic: overview. Global and Planetary Change 24, 165–173.

Jones, A., Plomerova, J., Korja, T., Sodoudi, F., Spakman, W., 2010. Europe from bottom up: a statistical examination of the central and northern European lithosphere–asthenosphere boundary from comparing seismological and elec- tromagnetic observations. Lithos 120, 14–29.

Jordan, M., 2003. JI-3D. A new approach to high resolution regional seismic tomography: Theory and applications, PhD thesis. Universität Göttingen, Göttingen, Germany.

Kennett, B.L.N., Engdahl, E.R., 1991. Traveltimes for global earthquake location and phase identification. Geophysical Journal International 105, 429–465.
Kennett, B.L.N., England, E., Buland, R., 1995. Constraints on seismic velocities in the Earth from travel times. Geophysical Journal International 122, 108–124.

Keyser, M., Ritter, J.R.R., Jordan, M., 2002. 3D shear wave velocity structure of the Eifel plume, Germany. Earth and Planetary Science Letters 203, 59–82.

Kirschner, S., Ritter, J., Wawerzinek, B., 2011. Teleseismic wave front anomalies at a Continental Rift: no mantle anomaly below the central Upper Rhine Graben. Geophysical Journal International 186, 447–462.

Lidmar-Bergström, K., Bonow, J.M., 2009. Hypotheses and observations on the origin of the landscape of southern Norway — a comment regarding the isostasy-climate-erosion hypothesis by Nielsen et al. 2008. Journal of Geodynamics 48, 95–100.

Lidmar-Bergström, K., Ollier, C.D., Sulebak, J.R., 2000. Landforms and uplift history of southern Norway. Global and Planetary Change 24, 211–231.

Martin, M., Ritter, J.R.R., the CALIXTO working group, 2005. High-resolution teleseismic body-wave tomography beneath SE Romania. — I. Implications for three-dimensional versus one-dimensional crustal correction strategies with a new crustal velocity model. Geophysical Journal International 162, 448–460.

Maupin, V., 2011. Upper-mantle structure in southern Norway from beamforming of Rayleigh wave data presenting multipathing. Geophysical Journal International 185, 985–1002.

Medhus, A.B., Balling, N., Holm Jacobsen, B., Kind, R., England, R.W., 2009. Deep-structural differences in southwestern Scandinavia revealed by P-wave travel time residuals. Norwegian Journal of Geology 89, 203–214.

Medhus, A.B., Balling, N., Jacobsen, B.H., Weidle, C., England, R.W., Kind, R., Thybo, H., Voss, P., 2012. Upper-mantle structure beneath the Southern Scandes Mountains and the Northern Tornquist Zone revealed by P-wave traveltime tomography. Geophysical Journal International 189, 1315–1334.

Megnin, C., Romanowicz, B., 2000. The three-dimensional velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. Geophysical Journal International 143, 709–728.

Nielsen, S.B., Gallagher, K., Leighton, C., Balling, N., Svenningsen, L., Jacobsen, B.H., Thomsen, E., Nielsen, O.B., Heilmann-Clausen, C.H., Egholm, D.L., Summerfield, M., Clausen, O.R., Piotrowski, J.A., Thorsen, M.R., Huuse, M., Abrahamsen, N., King, C., Lykke-Andersen, H., 2009. The evolution of western Scandinavian topography: a review of Neogene uplift versus the ICE (isostasy-climate-erosion) hypothesis. Journal of Geodynamics 47, 72–79.

Pascal, C., Olesen, O., 2009. Are the Norwegian mountains compensated by a mantle thermal anomaly at depth? Tectonophysics 475, 160–168.

Pedersen, H., Fishwick, S., Snyder, D., 2009. A comparison of cratonic roots through consistent analysis of seismic surface waves. Lithos 109, 81–95.

Plomerová, J., Babuška, V., Kozlovskaya, E., Hyvönen, L.V.L., 2008. Seismic anisotropy — A key to resolve fabrics of mantle lithosphere of Fennoscandia. Tectonophysics 462, 125–136.

Ramberg, I.B., Bryhni, I., Nottvedt, A., Rangnes, K., 2008. The Making of a Land — Geology of Norway. Norsk Geologisk Forening, Trondheim.

Ritsema, J., van Heijst, H.J., Woodhouse, J.H., 2004. Global transition zone tomography. Journal of Geophysical Research 109, B02302.

Ritter, J.R.R., 2007. The seismic signature of the Eifel Plume. In: Ritter, J.R.R., Achauer, U. (Eds.), Mantle Plumes: A Multidisciplinary Approach. Springer, Heidelberg, pp. 379–404.

Rohrman, M., van der Beek, P., 1996. Cenozoic postrift domal uplift of North Atlantic margins: an asthenospheric diapirism model. Geology 24, 901–904.

Roy, C., Ritter, J.R.R., in press. Complex deep seismic anisotropy below the Scandinavian Mountains. Journal of Seismology. http://dx.doi.org/10.1007/s10950-012-9325-4.

Sandoval, S., Kissling, E., Ansorge, J., SVEKALAPKO Seismic Tomography Working Group, 2004. High-resolution body wave tomography beneath the SVEKALAPKO array – II. Anomalous upper mantle structure beneath the central Baltic Shield. Geophysical Journal International 157, 200–214.

Shomali, Z.H., Roberts, R.G., TOR Working Group, 2002. Non-linear body wave teleseismic tomography along the TOR array. Geophysical Journal International 148, 562–574.

Shomali, Z.H., Roberts, R.G., Pedersen, L.B., the TOR Working Group, 2006. Lithospheric structure of the Tornquist Zone resolved by nonlinear P and S teleseismic tomography along the TOR array. Tectonophysics 416, 133–149.

Smelror, M., Dehls, J., Ebbing, J., Larsen, E., Lundin, E.R., Nordgulen, O., Osmundsen, P.T., Olesen, O., Ottesen, D., Pascal, C., Redfield, T.F., Rise, L., 2007. Towards a 4D topographic view of the Norwegian sea margin. Global and Planetary Change 58, 382–410.

Sobolev, S.V., Zeyen, H., Granet, M., Achauer, U., Bauer, C., Werling, F., Altherr, R., Fuchs, K., 1997. Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations. Tectonophysics 275, 143–164.

Stammler, K., 1993. SeismicHandler — Programmable multichannel data handler for interactive and automatic processing of seismological analyses. Computers and Geosciences 19, 135–140.

Steck, L.K., Prothero, W.A., 1991. A 3D raytracer for teleseismic bodywave arrival times. Bulletin of the Seismological Society of America 81, 1332–1339.

Stratford, W., Thybo, H., 2011a. Seismic structure and composition of the crust beneath the southern Scandes, Norway. Tectonophysics 502, 364–382.

Stratford, W., Thybo, H., 2011b. Crustal structure and composition of the Oslo Graben, Norway. Earth and Planetary Science Letters 304, 441–442.

Stratford, W., Thybo, H., Faleide, J.L., Olesen, O., Trygvasson, A., 2009. New Moho map of onshore southern Norway. Geophysical Journal International 178, 1755–1765.

Svenningsen, L., Balling, N., Jacobsen, B., Kind, R., Wylegalla, K., Schweitzer, J., 2007. Crustal root beneath the highlands of southern Norway resolved by teleseismic receiver functions. Geophysical Journal International 170, 1129–1138.

Torsvik, T.H., Cocks, L.R.M., 2005. Norway in space and time: a centennial cavalcade. Norwegian Journal of Geology 85, 73–86.

Wawerzinek, B., 2012. Untersuchung der elastischen Scherwellenstruktur unter dem Südkandinavischen Gebirge. Dissertation, Karlsruhe Institute of Technology, (in German).

Wawerzinek, B., Ritter, J.R.R., Jordan, M., Landes, M., 2008. An upper-mantle upwelling underneath Ireland revealed from non-linear tomography. Geophysical Journal International 175, 253–268.

Weidle, C., Maupin, V., 2008. An upper-mantle S-wave velocity model for Northern Europe from Love and Rayleigh group velocities. Geophysical Journal International 175, 1154–1168.

Weidle, C., Maupin, V., Ritter, J., Kværna, T., Schweitzer, J., Balling, N., Thybo, H., Faleide, J.I., Wenzel, F., 2010. MAGNUS — a seismological broadband experiment to resolve crustal and upper mantle structure beneath the Southern Scandes Mountains in Norway. Seismological Research Letters 81 (1), 76–84.

Wessel, P., Smith, W.H.F., 1998. New, improved version of the Generic Mapping Tools released. EOS, Transactions, American Geophysical Union 79, 579.

Zeyen, H., Achauer, U., 1997. Joint inversion of teleseismic delay times and gravity anomaly for regional structures. In: Fuchs, K. (Ed.), Upper Mantle Heterogeneities from Active and Passive Seismology. Kluwer Academic Publishers, Dordrecht, pp. 155–169.
# Appendix B

## Article: Complex deep seismic anisotropy below the Scandinavian Mountains

This article contains mainly my work of my diploma thesis at the Geophysical Institute (GPI) at the Karlsruhe Institute of Technology (KIT) in Germany. Additionally, results of my research on "Application of array processing to measure teleseismic shear wave splitting" during a one month research visit in 2011 at the Norwegian Seismic Array (NORSAR), supported by a NERA grant, are included. The collaboration with the GPI continued during my Ph.D thesis at the IPGP and the article was finished and published in August 2012.

### Abstract

To study seismic anisotropy the birefringence of SKS and SKKS phases were analyzed. These phases, which should be polarized radially, are split into an additional transverse component if they propagate through an anisotropic medium. The results are directions $\phi$ of the apparent fast shear wave polarization and delay times $\delta t$ between the split phases.

For station KONO in Southern Norway, $\phi$ and $\delta t$ are frequency-dependent $\phi$ and $\delta t$ values, indicating a depth-dependent anisotropy. Additionally, $\phi$ and $\delta t$ values vary with epicentre backazimuths in Norway, indicating a complex anisotropic structure in the crust and upper mantle.

Stacking of the SKS/SKKS waveforms improves the signal-to-noise ratio along one station line and allows us to better determine the splitting parameters. A unique and complete model of the complex anisotropy cannot be obtained due to the limited observed backazimuth range. Near-surface tectonic structures correlate with the splitting pattern and thus the crust is one anisotropic layer in the region. Partly preferred orientations in the rock fabric at the surface can be correlated with $\phi$.
<figure><img src="image_113.png" /></figure>