# Absolute earthquake locations using 3-D versus 1-D velocity models below a local seismic network: example from the Pyrenees

T. Theunissen, S. Chevrot, M. Sylvander, V. Monteiller, M. Calvet, A.

Villaseñor, S. Benahmed, H. Pauchet, F. Grimaud

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

## To cite this version:

T. Theunissen, S. Chevrot, M. Sylvander, V. Monteiller, M. Calvet, et al., Absolute earthquake locations using 3-D versus 1-D velocity models below a local seismic network: example from the Pyrenees. Geophysical Journal International, 2018, 212, pp.1806-1828. 10.1093/gji/ggx472. insu-04831808

HAL Id: insu-04831808

https://insu.hal.science/insu-04831808v1

Submitted on 12 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Absolute earthquake locations using 3-D versus 1-D velocity models below a local seismic network: example from the Pyrenees

T. Theunissen,<fnref n="1" /> S. Chevrot,<fnref n="2" />,<fnref n="3" /> M. Sylvander,<fnref n="4" />,<fnref n="5" /> V. Monteiller,<fnref n="6" /> M. Calvet,<fnref n="4" />,<fnref n="5" />
A. Villaseñor,<fnref n="7" /> S. Benahmed,<fnref n="4" />,<fnref n="5" /> H. Pauchet<fnref n="4" />,<fnref n="5" /> and F. Grimaud<fnref n="4" />,<fnref n="5" />

Accepted 2017 October 30. Received 2017 October 27; in original form 2016 September 5

## SUMMARY

Local seismic networks are usually designed so that earthquakes are located inside them (primary azimuthal gap <<180°) and close to the seismic stations (0–100 km). With these local or near-regional networks (0°–5°), many seismological observatories still routinely locate earthquakes using 1-D velocity models. Moving towards 3-D location algorithms requires robust 3-D velocity models. This work takes advantage of seismic monitoring spanning more than 30 yr in the Pyrenean region. We investigate the influence of a well-designed 3-D model with station corrections including basins structure and the geometry of the Mohorovicic discontinuity on earthquake locations. In the most favourable cases (GAP < 180° and distance to the first station lower than 15 km), results using 1-D velocity models are very similar to 3-D results. The horizontal accuracy in the 1-D case can be higher than in the 3-D case if lateral variations in the structure are not properly resolved. Depth is systematically better resolved in the 3-D model even on the boundaries of the seismic network (GAP > 180° and distance to the first station higher than 15 km). Errors on velocity models and accuracy of absolute earthquake locations are assessed based on a reference data set made of active seismic, quarry blasts and passive temporary experiments. Solutions and uncertainties are estimated using the probabilistic approach of the NonLinLoc (NLLoc) software based on Equal Differential Time. Some updates have been added to NLLoc to better focus on the final solution (outlier exclusion, multiscale grid search, S-phases weighting). Errors in the probabilistic approach are defined to take into account errors on velocity models and on arrival times. The seismicity in the final 3-D catalogue is located with a horizontal uncertainty of about 2.0 ± 1.9 km and a vertical uncertainty of about 3.0 ± 2.0 km.

**Key words:** Earthquake hazards; Seismic tomography; Crustal structure.

## 1 INTRODUCTION

Absolute earthquake location is probably the oldest inverse problem in seismology (Aki & Richards 1980). It is performed routinely in regional seismological observatories worldwide in order to build seismicity catalogues and/or to perform real-time earthquake location. The HYPO-71 or related software (Lee & Lahr 1972; Lahr 1989), based on the Geiger method (Geiger 1910), are still used in many seismic observatories. In this approach, traveltimes are usually computed in 1-D velocity models, sometimes considering static station delays for the contribution of topography and

shallow subsurface. Several studies have suggested that the use of a 3-D velocity model can significantly improve absolute earthquake locations (Font et al. 2004; Flanagan et al. 2007; Simmons et al. 2012; Theunissen et al. 2012b). However, these studies mostly focused on cases where seismic events are outside regional or global networks, that is, located far from the stations and with a significant azimuthal gap. Few studies have tried to quantify the benefit of using a 3-D rather than a 1-D velocity model (Husen et al. 2003). Indeed, it is widely believed that a large number of seismic stations around the event guarantees a well constrained hypocentre location (Bondar et al. 2004). However, in complex tectonic environments

<footnote n="1">¹Department of Earth Science, University of Bergen, Postboks 7803, N-5020, Bergen, Norway. E-mail: theunissen.thomas@gmail.com</footnote>

<footnote n="2">²Université de Toulouse; UPS-OMP; GET; F-31062 Toulouse, France</footnote>

<footnote n="3">³CNRS; GET; 14, avenue Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote n="4">⁴Université de Toulouse; UPS-OMP; IRAP; F-31062 Toulouse, France</footnote>

<footnote n="5">⁵CNRS; IRAP; 14, avenue Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote n="6">⁶LMA, CNRS, UPR7051, Aix-Marseille University, Centrale Marseille, F-13402 Marseille Cedex 20, France</footnote>

<footnote n="7">⁷Instituto de Ciencias de la Tierra Jaume Almera (ICTJA) – CSIC, E-08028 Barcelona, Spain</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Seismic network and structural map of the Pyrenees: (a) Seismic stations network used in this study. (b) Main geological and structural units of the study region modified from Chevrot *et al.* (2014). NPF: North Pyrenean Fault, NPFT: North Pyrenean Front Thrust, SPFT: South Pyrenean Front Thrust, SH: Sillon Houiller (sometimes also referred to as the ‘Toulouse Fault’), NPZ: North Pyrenean Zone, and SPZ: South Pyrenean Zone.</figcaption></figure>

such as continental orogens, characterized by complex 3-D struc-
tures, it has been shown that the use of several subregional 1-D mod-
els can improve absolute earthquake locations (Husen et al. 2011).
Several studies of 3-D earthquake location at a regional scale (e.g.
Satriano et al. 2006; Huang & Zhao 2012; Mostaccio et al. 2013)
have shown that traveltime residuals and clustering are clearly im-
proved but that the improvement in terms of accuracy has not been
clearly established. Active sources (quarry blast, active seismic ex-
periments) are sometimes used as Ground Truth (GT) reference
events to assess the accuracy of absolute earthquake locations at
local scale (Husen et al. 1999, 2003; Lin et al. 2006) as well as the
quality of 3-D velocity models (for example Satriano et al. 2006;
Lin et al. 2011). An important limitation in hypocentre determina-
tion of shallow events comes from the poorly constrained shallow
subsurface. Because of the non-linearity of the earthquake loca-
tion process, a well-constrained 1-D velocity model (large data set
or/and based on controlled seismics) may produce more accurate
earthquake locations than a poorly constrained 3-D model.

We quantify several aspects of the 3-D absolute earthquake lo-
cation process below a local or regional (0°-5°) seismic network,
compared to 1-D location procedures. In particular, we will focus
on traveltime estimation differences and station corrections, from
one velocity model to another. We use reference sources such as
shots, quarry blasts and reference earthquakes to estimate location
accuracy. The aim of this work is (1) to locate seismicity using a

simple and stable 3-D routine, (2) quantify errors of earthquake
location parameters and (3) estimate the accuracy of hypocentres'
locations.

Comparisons are made between 1-D and 3-D results and we
demonstrate the importance of using average layered 3-D models
including sedimentary basin geometry and Moho topography with
station corrections to obtain the best hypocentre localization.

2 DATA

## 2.1 The RSSP network

We will exploit the data from the Pyrenean seismic moni-
toring network (RSSP, Réseau de Surveillance Sismique des
Pyrénées http://rssp.irap.omp.eu/ in French) operated by the Ob-
servatoire Midi-Pyrenees (OMP). This seismic network has been
operating for more than 30 yr (e.g. Souriau & Granet 1995;
Souriau & Pauchet 1998; Souriau et al. 2001) (Details in sup-
plementary material). This study combines all available picks
(P- and S-wave arrival times) from earthquakes recorded in
the period 1978–2013 by the RSSP network, complemented
by other regional seismic networks in Spain (IGN, Instituto
Geografico Nacional http://www.ign.es/ign/layout/sismo.do) and in
Catalunya (ICGC, Institut Cartographic i Geologic de Catalunya
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Epicentral distribution and overview some quality criteria of the whole data set (Section 2.3). All calculations are made using the procedure and the 1-D velocity model from OMP (Gallart *et al.* 1981; Daignières *et al.* 1982; Pauchet *et al.* 1999) with associated station corrections (correction from topography) and based on the flat earth approximation. (a) Epicentral distribution of the compiled seismicity catalogue. Colour according to the depth. (b) Distribution of primary azimuthal gap (GAP I). (c) Distribution of the distance to the first station. (d) Root Mean Square (RMS) versus GAP I plot. (e) horizontal and vertical uncertainties, respectively ERH and ERZ, are estimated in that process by the average of differences between each result from HYPO-71 with different initial depth (commonly used in the observatory).</figcaption></figure>

http://www.igc.cat/web/en/terratremols.php) and several temporary experiments (Fig. 1 and Supporting Information Table S1).

## 2.2 Seismicity in the Pyrenees

The Pyrenees are an E-W striking orogen flanked by two foreland basins. The active deformation in the Pyrenees is low and mainly characterized by extension (Chevrot *et al.* 2011; Asensio *et al.* 2012; Rigo *et al.* 2015). Pyrenean seismicity is low to moderate with occasional earthquakes of magnitude up to 5.5 (e.g. Souriau & Granet 1995; Souriau & Pauchet 1998; Souriau *et al.* 2001). The largest historical earthquake is the 1660 Bagneres-de-Bigorre earthquake (21/06/1660) with a moment magnitude of $M_W$ 6.1 ± 0.1 (Cara *et al.* 2008) or $M_W$ 6.4 ± 0.3 (Stucchi *et al.* 2013; Locati *et al.* 2014). The seismicity is located inside the upper crust at depths shallower than 15 km (Fig. 2). On the west side, the seismicity is concentrated along the North Pyrenean Fault while its distribution is more diffuse on the east side.

Several events can be considered as Ground Truth (GT) Reference Events, since these events have very precise hypocentre locations (Supporting Information Fig. S1). The RSSP has identified and indexed about 2000 quarry blasts recorded between 1997 and 2011.

These blasts are usually recorded by less than 6 seismic stations, resulting in inaccurate locations. Upon request, we were provided information about the quantity of TNT and the exact location (but no information about origin time of blasts) at the Luzenac talc quarry located in the eastern part of the Pyrenean belt. Two quarry blasts from Luzenac have been selected in addition to a quarry blast from the quarry of Bustince located in the western part of the Pyrenean belt.

To define GT Reference earthquakes, we selected earthquakes located by dense temporary networks and therefore recorded at short distance and with small azimuthal gap. Four sequences of earthquakes have been well studied in the past using temporary seismic stations. These seismic sequences are used to test the accuracy of absolute location:

(i) the Argelès aftershocks sequence following the Argelès-Gazost earthquake ($M_L$ = 5.0) from 2006/11/17 18h19 to the end of December 2006. Aftershocks, relocated in a 3-D velocity model, are well distributed along an E-W trending normal fault (Sylvander *et al.* 2008).

(ii) the Cauterets sequence following two moderate earthquakes ($M_L$ = 4.6 and 4.3) between 2002 May 16 and 2002
May 20 is distributed in two small very clustered swarms that were relocated by a double-difference relative location process (Dubos et al. 2004).

(iii) The Agly sequence from 1996 February 18 to February 24 following the Saint-Paul-de-Fenouillet earthquake ($M_L = 5.2$) (Rigo et al. 1997; Pauchet et al. 1999; Sylvander et al. 2007) on the east side of the Pyrenees

(iv) The Lacq induced seismicity, related to hydrocarbon extraction. This seismicity for the period 1974–1997 has been relocated with a local 3-D velocity model (Guyoton et al. 1992; Bardainne et al. 2008).

## 2.3 Overview of *P*- and *S*-wave arrival times data set

### 2.3.1 OMP Bulletin

The current procedure to locate earthquakes and build the seismicity catalogue at OMP is based on HYPO-71 (Lee & Lahr 1972; Lahr 1989). Traveltime picks are weighted according to the epicentral distance from 1 for stations closer than 30 km from the epicentre to 0 for stations further than 120 km. A $V_P/V_S$ ratio of 1.75 is applied on the 1-D *P*-wave velocity model deduced from seismic refraction profiles (Gallart et al. 1981; Daignières et al. 1982; Pauchet et al. 1999). Station delays correspond to topography corrections. Several inversions are run with different initial depth and the final solution is determined by a statistical analysis of the results obtained in different inversion. The OMP catalogue contains about 20 500 events recorded by a total of 47 seismic stations (Supporting Information Table S1).

### 2.3.2 The final data set

All available bulletins summarized in Supporting Information Table S1 have been combined with the OMP bulletin. For this study, we also added manual picks from 26 accelerometers over the period 2001–2010 for 151 earthquakes and from 9 broad-band seismic stations as well as on the 133 temporary seismic stations from the PYROPE (Chevrot & Sylvander 2017)/IBERARRAY (Díaz et al. 2009) experiments for 313 earthquakes over the period 2010/09-2013/02. The distribution of epicentres is shown in Fig. 2(a). The analysis of the picks in the entire data set ($P_g$, $P_n$) provides an estimate pick uncertainty of 0.24 s for *P* waves and 0.41 s for *S* waves. The scatter of traveltime picks strongly depends on quality. It is respectively of about 0.13, 0.44, 0.37 and 0.59 s for qualities from 0 (the best) to 3 (the worst) for *P*-phases and about 0.18, 0.20, 0.33 and 0.56 s for qualities from 0 to 3 for *S*-phases.

The final data set contains 24 500 earthquakes from October 1975 to February 2013 (Fig. 2) for a total of 235 261 *P* phases and 235 172 *S* phases recorded by 377 different seismic stations. 80 per cent of this data set has been recorded after 1996 (Fig. 2).

The average primary azimuthal gap (GAP I) is 168 ± 64°. The OMP procedure gives a depth uncertainty of 2.7 km while the horizontal uncertainty is on average 1.3 km.

## 3 METHODOLOGY

To deal with the trade-off between earthquake location and *P*- and *S*-waves velocity structures (Crosson 1976) and to limit the impact of our imperfect knowledge of the 3-D structure (Chang et al. 1983; Myers & Schultz 2000; Bondar & McLaughlin 2009; Thurber 1992;

Eberhart-Phillips & Michael 1993; Kissling et al. 1995), we use the following procedure:

(i) Computation of an initial 1-D velocity model and associated stations delays and construction of a 3-D *a-priori* velocity model based on available geophysical and geological data

(ii) Absolute earthquake location with these initial velocity models and re-computations of station delays

(iii) 3-D local seismic tomography using initial earthquake location and the 3-D *a-priori* velocity model as initial state

(iv) Results of the inversion are used to compute a 3-D average velocity model which is used for the final 3-D absolute earthquake location.

### 3.1 Traveltime computation

Traveltimes are computed with a 3-D finite difference algorithm (Podvin & Lecomte 1991), within a Cartesian grid of homogeneous cubic cells. Velocity models are described by cells of $4 \times 4 \times 2 \text{ km}^3$. These velocities are interpolated on a grid composed by cells of $1 \times 1 \times 1 \text{ km}^3$ for traveltime computation.

A posteriori ray tracing is then performed to determine the rays geometry and to reduce the error on the traveltime estimation by a factor 10 (Monteiller et al. 2005). Errors on traveltimes are about $10^{-3}$ s for a cell size of $1 \text{ km}^3$ using such a procedure (Monteiller 2005). Since this approach is computationally demanding it is only applied in the local earthquake tomography (LET). Indeed, in 3-D absolute earthquake location, we need to compute traveltimes from all cells to all stations and for both *P* and *S* phases. The error is about $10^{-2}$ s for a cell size of 20 m. We use this resolution to calculate the reference traveltimes. We then use these reference traveltimes to calculate the errors made using coarser grid sizes. Supporting Information Fig. S3 shows that the error with a cell size of $1 \text{ km}^3$ is lower than 0.08 s for a layered 1-D velocity model and lower than 0.2 s in a complex 3-D model with velocity gradients. Since our final 3-D model is described by $4 \times 4 \times 2 \text{ km}^3$ cells, we estimate that errors on traveltimes are between 0.1 s and 0.2 s. The impact of these errors on earthquake location is discussed in Section 6.1.

### 3.2 Local tomography

We use a traveltime tomography method to simultaneously invert for the velocity model and the hypocentre parameters (Aki & Lee 1976; Spencer & Gubbins 1980; Spakman & Nolet 1988; Thurber 1992; Benz et al. 1996). The code developed by Monteiller (2005) has been already applied in several studies (Latorre et al. 2004; Vanorio et al. 2005; Gautier et al. 2006; Theunissen et al. 2012b). It is a linearized, iterative and damped least-square approach. The inverse problem is solved with the LSQR algorithm (Paige & Saunders 1982).

The inversion grid is composed of $153 \times 130 \times 40$ (respectively along the *x*, *y* and *z* dimensions), that is, 795 600, regular cells of $4 \times 4 \times 2 \text{ km}^3$. The computation of traveltime residuals is performed in a grid with cubic cells of $1 \text{ km}^3$. We invert for $V_P$ and $V_P/V_S$. The total number of unknowns in the inversion is about 1 600 000, given by the number of cells $\times 2$ ($V_P$ and $V_S$) + number of earthquakes $\times 4$ (hypocentre location and origin time $t_0$). The regularization of
the problem involves a 3-D Laplacian smoothing matrix (L) and a diagonal damping matrix (DP):

$$ \begin{pmatrix} D_1 \cdot A \cdot D_2 \\ L \\ DP \end{pmatrix} \cdot m = \begin{pmatrix} D_1 \cdot b \\ 0 \\ 0 \end{pmatrix} \iff A_R \cdot m = b_R \quad (1) $$

where A is the slowness partial derivatives matrix, m is the model vector and b the traveltime residual vector. As proposed by some authors (Spakman & Nolet 1988; Le Meur et al. 1997), normalization and scaling of the derivative matrix through diagonal matrices $D_1$ (weighting matrix) and $D_2$ (normalization matrix), is performed for better reconstruction of the different parameters. This operation will remove the influence of parameter units and will also take into account the data sensitivity to each class of parameters. The parameters used for the regularization (weighting, smoothing and damping) are dimensionless. The weighting is defined by three parameters: pick quality, amplitude of traveltime residual and length of rays. The smoothing is defined along the three spacial directions and a different coefficient can be applied to $1/V_p$ and $V_p/V_s$. All these parameters are fixed during the inversion and are defined by trial-and-error in order to obtain the best compromise between variance reduction, smoothing and amplitude anomalies.

The data set is sorted in order to remove the outliers identified from the hodochrones, the epicentral distribution of traveltime residuals and the Wadati diagram (Supporting Information Fig. S2). The initial velocity model is the 3-D a-priori model (see Section 4) and the best available earthquake locations are used as input for hypocentres' location. The cell size ($4 \times 4 \times 2 \text{ km}^3$) is a compromise between forward calculation efficiency and maximum expected resolution. The inversions have been performed in two steps: first by inverting $P_g$ and $S_g$ phases, and second by adding $P_N$ and $S_N$ phases, keeping fixed hypocentres location and origin times.

## 3.3 Absolute earthquake location computation

Our choice of the NLLoc software (NonLinLoc from Lomax et al. 2000) is motivated by three reasons: (1) the solution can be described by a probability density function (PDF; e.g. Tarantola & Valette 1982; Moser et al. 1992; Wittlinger et al. 1993) based on Equal Differential Time that is more robust in presence of outliers and which is independent of the origin time (Zhou 1994; Font et al. 2004; Lomax 2005; Theunissen et al. 2012b), (2) the weighting scheme is well designed and (3) the search algorithm, Oct-Tree, is efficient and permits an efficient browsing of the entire structure leading also to a good estimation of the location PDF.

We consider a traveltime-dependent Gaussian model error of 3 per cent for both P- and S-waves. This error is between 0.05 s and 0.5 for P-waves and between 0.125 and 1.0 s for S-waves. An initial hypocentre solution is used to compute a distance weighting that puts a weight of 1 to stations closer than 90 km to the source and a weight that linearly decreases to 0 for stations further than 120 km. Phase pick qualities 0,1,2,3,4 are used to define time uncertainties 0.05, 0.2, 0.3, 0.5 and 4.0 s respectively based on picks analysis. The quality given to a pick is higher for S-phase than P-phase for a given quality. A factor of 2.5 is thus applied to S-phase time uncertainties.

The Oct-tree algorithm is efficient for defining the PDF when grid search parameters are well defined. However, it is quite difficult to correctly define a grid search for the whole seismic network in order to keep the same degree of resolution for all final solutions. To avoid this, we have defined an optimal grid search 60 km × 60 km × 59 km horizontally centred on the solution. To focus on

this optimal grid search, an iterative multiscale approach inspired from Theunissen et al. (2012a) is used (Supporting Information Fig. S3). The iterative approach consists of first browsing the parameter space (x, y, z) using the first level of Oct-Tree search, that is, only the initial grid is browsed, for an increasingly small volume and then a refined grid is used. The final complete Oct-Tree search is done with an optimal grid search using an initial sampling of 20 × 20 × 24 cells leading to a starting resolution of 3 km × 3 km × 2.5 km respectively along x, y and z-directions (Supporting Information Fig. S3). The efficiency of the Oct-tree algorithm is reduced in the presence of outliers by lowering spatial resolution. In order to take this into account, we consider that if an absolute traveltime residual is higher than a threshold, the phase/station couple is temporarily removed from the PDF calculation of the centred solution of the Oct-Tree subdivision. This threshold is defined as the median of all residuals used in the process, in order to consider worse solutions mainly when subdivisions are far away from the final solution, plus a maximum time error. The maximum time error corresponds to the traveltime-dependent Gaussian model error of the considered phase/station couple.

## 4 CONSTRUCTION OF VELOCITY MODELS

Four velocity models are constructed for this study. The first model is a 1-D velocity model. This model is compared to the 1-D model from OMP (called 1-D flat OMP). It is derived for two different geographic references, flat earth (called 1-D flat) and geocentric (called 1-D geoc). The reader can find details about geographic references in the supplementary material. The second velocity model is a 3-D a-priori velocity model that includes all geophysical and geological constraints. This model is only used as an initial velocity model for the inversion procedure. The third velocity model is from an LET. This model is computed using a geocentric reference frame (called 3-D geoc tomographic). The fourth velocity model is directly derived by averaging the previous LET model. It is dedicated to absolute earthquake location (called 3-D geoc average). Table 1 summarizes all earthquake location procedures and their velocity models.

### 4.1 1-D models (called 1D flat and 1D geoc)

We derived a 1-D velocity model with its station corrections using the VELEST 3.1 program (Kissling 1988; Kissling et al. 1994, 1995). We performed a combined inversion of the velocity model, hypocentre parameters (location and origin time) and station delays using P and S arrival-times provided by a selected data set (2051 earthquakes with at least 12 P and 6 S phase observations with azimuthal gap <180°). This inversion is repeated as a trial-and-error process with different initial velocity models (derived from the OMP observatory routine 1-D model and the 1-D model from Souriau & Granet (1995)) and hypocentral parameters (from OMP bulletins), and for different damping values in order to get the so-called best Minimum 1-D velocity model (Fig. 3).

$P_g$ and $S_g$ are first used to jointly determine the velocities in the upper-crustal layers with hypocentre location and station corrections. Then refracted phases are added in order to constrain lower-crustal layers and determine station corrections.

Station corrections capture the contribution of shallow geological layers (basement and meso-cenozoic sedimentary basins) while those of remote seismic stations, north of the Aquitania basin and southwest of the Ebro basin, are also affected by the Moho
<table><caption>**Table 1.** List of absolute earthquake location processes used in this study.</caption><thead><tr><th>Alias</th><th>Process num.</th><th>Algorithm</th><th>P- and S-wave models</th><th>Geographic reference</th><th>Station delays</th></tr></thead><tbody><tr><td colspan="6"><strong>1-D structure</strong></td></tr><tr><td>1-D flat OMP routine</td><td>1</td><td>OMP procedure: multi-runs of HYPO71 with different initial depth</td><td>Approximative 1-D P-model based on refraction (Pauchet et al. 1999; Gallart et al. 1981; Daignières et al. 1982) and constant V<sub>P</sub>/V<sub>S</sub> ratio of 1.75 V<sub>P</sub> and V<sub>S</sub> 1-D VELEST models</td><td>Flat earth (simple conversion)</td><td>Yes, based on station elevation</td></tr><tr><td rowspan="2">1-D flat NLLoc</td><td>2</td><td rowspan="2">NLLoc</td><td rowspan="2">V<sub>P</sub> and V<sub>S</sub> 1-D VELEST models</td><td rowspan="2">Flat earth approximation (Lambert III projection)</td><td>No</td></tr><tr><td>3</td><td>Yes, jointly inverted into VELEST and corrected from mean residuals*</td></tr><tr><td rowspan="2">1-D geoc NLLoc</td><td>4</td><td rowspan="2">NLLoc</td><td rowspan="2">V<sub>P</sub> and V<sub>S</sub> 1-D VELEST models</td><td rowspan="2">Ellipsoidal earth approximation (geocentric conversion)</td><td>No</td></tr><tr><td>5</td><td>Yes, jointly inverted into VELEST and corrected from mean residuals*</td></tr><tr><td colspan="6"><strong>3-D structure</strong></td></tr><tr><td rowspan="2">3-D geoc average NLLoc</td><td>6</td><td rowspan="2">NLLoc</td><td rowspan="2">V<sub>P</sub> and V<sub>S</sub> 3-D a-priori models based on averaging 3-D tomographic models</td><td rowspan="2">Ellipsoidal earth approximation (geocentric conversion)</td><td>No</td></tr><tr><td>7</td><td>Yes, mean residuals</td></tr><tr><td>3-D geoc tomo NLLoc</td><td>8</td><td>NLLoc</td><td>V<sub>P</sub> and V<sub>S</sub> 3-D tomographic models</td><td>Ellipsoidal earth approximation (geocentric conversion)</td><td>No</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** 1-D velocity model (blue and red lines) and associated station delays calculated from multiple runs (grey lines) of VELEST (Kissling 1988; Kissling et al. 1994, 1995) (see the text, Section 4). 2051 earthquakes in the 1-D joint inversion with at least 12 *P* and 6 *S* phase observations with azimuthal gap < 180° have been used. (a) 1-D *P*- and *S*-wave velocities. (b) *P* phases station delays. (c) *S* phases station delays. Values for velocity model and station delays can be found in Supporting Information Tables S2 and S3.</figcaption></figure>

depth variations and the earth's sphericity that are not considered (Fig. 3).

The 1-D model from VELEST is interpolated on a Cartesian grid composed of 4 × 4 × 2 km cells similarly to other models. The aforementioned geographic conversion from flat-earth to a geocentric system is applied at this stage (see the Supporting Information for details). Station corrections are computed by running NLloc on the entire catalogue (Frohlich 1979; Pujol 1988), taking the average residual when the number of phases is larger than 400, or the median when the number of phases is smaller, in order to limit the effect of outliers. In the flat earth approximation, VELEST station delays differ by −0.02 ± 0.09 s on average compared to these new station delays. In the geocentric reference, VELEST station delays differ by 0.01 ± 0.32 s.

Figs 4(a) and (b) shows that with a 1-D velocity model, the flat earth approximation leads to an over-estimate of *P*-traveltimes of about 0.24 s at a distance of 250 km for an hypocentre located between 0 and 30 km depth. For epicentral distances of less than 120 km, traveltimes are underestimated by up to about 0.1 s. We may note that near field differences and the apparent asymmetry are enhanced (to greatly varying degree) by averaging 1-D layers into a Cartesian grid with an ellipsoidal shape. These observations explain the increase scatter of station delays differences from VELEST compared to those calculated by averaging residuals in the geocentric reference. Readers may also refer to the paper of Snoke & Lahr (2001) for a discussion of using spherical approximation versus flat approximation in earthquake location at regional distance.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** *P* traveltimes differences between velocity models. (1) GEOC1D (*1D geoc NLloc*) and FLAT (*1D flat NLloc*) to evaluate the effect of earth sphericity (a,b). (2) GEOC3D (*3D geoc average NLloc*) and GEOC1D effect of 3-D structure (c,d) with a seismic station located on the basement (e,f) with a seismic station located on a sedimentary basin. (3) 2-D forward velocity model (Line 8 from Pedreira *et al.* (2003)) and GEOC3D. PYLU seismic station is located at the centre of the seismic network.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Geological and geophysical inputs for the construction of the 3-D *a-priori* velocity model. (a) Thickness of Ebro and Aquitania sedimentary basins. Data come from a compilation of isobaths maps (Chantraine et al. 1996; Serrano et al. 2006; Filleaudeau 2011; Lozano 2012), geological cross-sections (Biteau et al. 2006; Serrano et al. 2006; Filleaudeau 2011) and well data (Filleaudeau 2011). (b) Smooth Mohorovicic discontinuity built from compilation of Pedreira et al. (2007); Diaz & Gallart (2009); Diaz et al. (2012) and receiver function calculation at PYROPE, TOPO-IBERIA and RLPB (Réseau Large Bande Permanent) (Chevrot et al. 2014, 2015) (see also Supporting Information Fig. S4). (c) smooth gradient used for $V_P$ within the basement (i.e. the crust minus sedimentary basins). It is built from a selection of reference velocity models (lines 1, 6, 7 and 8 of Pedreira et al. (2003)). (d) Wadati diagram for a selection of about 11 000 well located earthquakes with the 1-D velocity model with station delays computed with a flat earth approximation.</figcaption></figure>

## 4.2 3-D velocity models

The 3-D *P*- and *S*-velocity models are built in two steps. In the first step, we exploit all the available geophysical and geological data. The second step consists of constraining average velocity fields within the crust using the result of a 3-D LET.

### 4.2.1 Construction of a 3-D *a-priori* model

We build a 3-D model (Fig. 5) that includes:

(i) **Relief** topography and bathymetry, taken from ETOPO1 (Amante & Eakins 2009).

(ii) **Basement top** that corresponds to the bottom of the Meso-Cenozoic Ebro and Aquitania basins, taken from a compilation of isobaths maps (Chantraine et al. 1996; Serrano et al. 2006; Filleaudeau 2011; Lozano 2012), geological cross-sections (Biteau

et al. 2006; Serrano et al. 2006; Filleaudeau 2011) and well data (Filleaudeau 2011).

(iii) **Crust** We use a smooth velocity field throughout the crust. To constrain the crustal *P*-wave velocity field and to quantify the effect on traveltime estimation, we use the 3-D *P* velocity model of Guyoton et al. (1992) as a reference for sedimentary basins and the forward 2-D seismic refraction lines from Pedreira et al. (2003) for the crust. We compile data from selected reference velocity models (lines 1, 6, 7 and 8 from Pedreira et al. (2003)) to apply a smooth gradient fitting average features of the crustal velocity field within the basement (Fig. 5) in contrast to the sedimentary basin for which we applied a constant average velocity of 5.0 km s⁻¹ for *P*-waves velocity based on a 3-D velocity model of (Supporting Information Fig. S4). For *S*-waves, we apply a constant $V_P/V_S$ ratio. The Wadati diagram, based on a data set of well-located earthquakes determined within the 1-D minimum velocity model with station corrections, gives an average $V_P/V_S$ ratio of 1.732 ± 0.032.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The 3-D geocentric average velocity model constructed from the 3-D local earthquake tomography based on a 3-D *a-priori* initial velocity model. (a–c) Selection of trusted parts of the inversion results used for the averaging. (a) maps view of selected cells. (b,c) SSW–NNE cross-section across the central Pyrenean belt showing the location of selected cells on $V_P$ and $V_P/V_S$ inverted parameters, respectively. (d,e) $V_P$ and $V_P/V_S$ fields respectively for the resulting 3-D average velocity model dedicated to absolute earthquake location purposes. The inset shows velocity profiles for sedimentary basins (green), crust (red) and mantle (purple) that result from the horizontal average of selected nodes. Continue and dashed lines represents $V_P$ and $V_P/V_S$ respectively. (f,g) station delays associated with the 3-D average velocity model.</figcaption></figure>

(iv) **Crust/Mantle boundary** We use data from deep seismic sounding (Pedreira *et al.* 2007; Diaz & Gallart 2009; Diaz *et al.* 2012) and receiver function studies compiled in Chevrot *et al.* (2014, 2015) (Fig. 5 and Supporting Information Fig. S5)

### 4.2.2 Construction of a 3-D model dedicated to absolute earthquake location (called 3D geoc average)

We use the result of the 3-D local tomography inversion in order to average a well constrained layered velocity structure for the crust. The 3-D *a-priori* model is used as an initial velocity during the inversion. Best resolved cells are used to compute the average (Fig. 6). Velocities (*P* or *S*) from a cell of the inversion grid (4 × 4 × 2 km³) are only used if the minimum cumulative length of rays reaches 5 km. This selection is limited to a region around the Pyrenean belt where the maximum resolution is achieved. In the end, we use

horizontal averages for the crust and mantle resulting in a smooth vertical gradient. For sedimentary basins, lateral variations are so strong that the average exhibits an artificial discontinuity (inset of Figs 6d and e). We therefore keep an average velocity of 5 km s⁻¹ for basins. We use 1.780, 1.708 and 1.752 respectively for $V_P/V_S$ ratio of basins, crust and mantle. Station delays are then computed after applying NLloc with the entire catalogue following the same procedure used with the 1-D model.

## 5 RESULTS

Table 1 gives the different sets of absolute earthquake location used in this study. 22 784 earthquakes are relocated. In the following, we distinguish shallow (0–6 km) from deep upper-crustal events (6–15 km) among GT reference seismic events.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Probability Density Function sampling and uncertainties used in this study, example from the $M_L$5.0 2006/11/17 18h19 Argelès earthquake.</figcaption></figure>

## 5.1 Earthquake location comparison and uncertainties

Supporting Information Figs S6 and S7 present an overview of the quality of the results. The quality is described by the RMS (root mean square) of phase residuals used during the process and by vertical and horizontal uncertainties. OMP locations obtained with HYPO-71 are shown in Fig. 2. Uncertainties ERH and ERZ (horizontal and vertical uncertainties) given by the OMP bulletin are quantified from the average of the different results obtained using different initial depths. Regarding results from the NLloc process, uncertainties are defined from the PDF (Fig. 7). NLloc provides Gaussian error estimates such as the expectation hypocentre location which represents the isobarycentre of sampling dots defining the PDF, and the 68 per cent confidence ellipsoid (Lomax et al. 2000). These estimates reflect the errors obtained by lin-

earized location algorithms (Lomax et al. 2000) such as HYPO-71 (Lee & Lahr 1972) or HYPOELLIPSE (Lahr 1989). Large differences between the expectation hypocentre location and the maximum likelihood hypocentre location as estimated from the PDF can be assumed to result from an ill-conditioned location problem since the Gaussian estimates are only adequate uncertainty estimates if the PDF shows a clear, single, global minimum (Lomax et al. 2000). To be able to statistically describe uncertainties, we thus defined SEH and SEZ, respectively for horizontal and vertical uncertainty, that corresponds to the maximum horizontal uncertainty deduced from location of semi-major principal axis of the 68 per cent joint-confidence ellipsoid divided by 1.87 times. This definition is the same as the one in the HYPOELLIPSE software (Lahr 1999).
In the end, results are then statistically compared based on seven quality parameters: RMS, SEH, SEZ, edd (Euclidean distance between the solution and the expectation hypocentre), edh (same but horizontal distance on sphere), edz (same but vertical difference $Z_{\text{expectation}} - Z_{\text{solution}}$), updf (average distance between the solution and all nodes describing the PDF).

All results show the increase of these seven quality parameters with the azimuthal gap. The median SEH linearly increases close to zero when azimuthal gap is lower than 90° to about 5 km for a gap of 325°. Statistical analysis of quality parameters over the entire catalogue for each earthquake location process, are not really significant since differences are lower than 100 m (Supporting Information Fig. S7). However, note that locations without station delays systematically exhibit high location uncertainties. Note also that the locations of the expectation solution (based on edz) are systematically deeper than the final solution. This pattern can be explained by the trade-off (1) between the depth and the origin time as well as between (2) the depth, the seismic network configuration and the velocity structure resulting in a vertically distributed PDF. On average the RMS is 0.12 ± 0.06, the SEH 2.05 ± 1.92 km, the SEZ 3.00 ± 2.00 km, the edd 2.1 ± 2.6 km, the edh 0.86 ± 1.38 km, the edz 1.24 ± 2.60 km and at last the updf 1.88 ± 1.82 km.

Supporting Information Fig. S8 shows the earthquake location differences using station delays. Most of the earthquakes are located within the range of uncertainties. On average the horizontal misfit is close to 1.0 ± 1.2 km. The horizontal misfit is largest with results from OMP bulletin compared to those using the 3-D average model (2.0 ± 1.7 km). On average the vertical misfit is about 0.6 ± 3.0 km. Finally, comparing all results together, 95 per cent (~2σ) of earthquakes seem to be located within a volume extending to a maximum of 5 km horizontally and 10 km vertically.

## 5.2 Accuracy for shallow events (0–6 km)

### 5.2.1 Quarry blasts

Fig. 8 shows the absolute locations obtained for the three quarry blasts and Supporting Information Fig. S9 the PDF of the first Luzenac quarry blast. Supporting Information Fig. S10 shows the quality parameters for the three quarry blasts.

We first describe results obtained with station delays and the 3-D tomographic model without station delays. In most cases, the horizontal accuracy is less than 1 km (mainly between 100 m and 800 m). In the case of the Luzenac quarry blasts, the accuracy using the 3-D velocity models (3-D average and 3-D tomographic models) is 1.5 km for blast *n* 1 and 1.3 km for the *n* 2. The vertical accuracy is poorer for the HYPO-71 OMP routine (<7 km), quite good for the two NLLoc 1-D processes (<2 km), perfect for the 3-D average model (<100 m) and unstable using the 3-D tomographic model (<100 m for two blasts and 14 km for the third one). The best misfits and the more stable solutions are obtained with the 3-D average model. The 1-D locations show very good horizontal accuracy. The RMS using 3-D models is higher than 1-D models for the Luzenac quarry and the PDF for the Luzenac quarry blast *n* 1 (Supporting Information Fig. S9) shows that solution is less well constrained. The edd, edh, edz and updf (Supporting Information Fig. S10) show that this is the case for the three quarry blasts.

Results without station delays have poorer accuracy except for the case of the Bustince quarry. In the case of the Luzenac quarry

blast *n* 2, using the 3-D average model without station delays leads to a poor depth estimation. Using the 3-D tomographic model does not improve this depth determination, which suggests that the inversion from arrival time tomography does not lead to a better estimation of traveltimes in the shallow part of this region.

### 5.2.2 Lacq gas field seismicity

Figs 9 and 10 show the absolute earthquake locations of the induced Lacq gas field seismicity. The reference is based on a NLLoc 3-D absolute earthquake location process (Bardainne et al. 2008) using a flat earth approximation with simple geographic conversion to kilometres. A high resolution 3-D P velocity model (Guyoton et al. 1992) and a constant $V_p/V_s$ ratio of 1.73 has been used. Over the period 1974–1997, seismic events are recorded on average by about 9–10 stations. Only few S-waves have been picked because most seismometers were single-component. As a result, 50 per cent of events have no S-phases and others have less than 8 S-phases. The azimuthal gap is often much lower than 180°. All the closest stations were at distance shorter than 8 km.

To first order, clustering based on our NLLoc procedure is better than the HYPO-71 OMP routine whatever the velocity model and with or without applying station delays. Compared to the reference, the seismicity is more clustered. Moreover, several earthquakes are located below 10 km depth while this is not the case in our results using the 3-D velocity models except for two earthquakes. Some earthquakes are located close to the free surface in our processes. These features are amplified using a 1-D velocity model and limited using 3-D models.

Regarding the seismicity recorded during 1998–2013 (Fig. 10), since the Lacq gas field network has been removed, results are more contrasted. Azimuthal gaps are larger than 200° and the distance to the first recording seismic station is systematically higher than 28 km. While all 1-D processes give worst results, using the 3-D average velocity model (with or without station delays) gives a blurred image of the Lacq gas field seismicity meaning that the accuracy is quite good both horizontally (<2 km) and vertically (<3 km). Only 16 earthquakes out of 282 are located too deep (>10 km). The HYPO-71 OMP routine gives only deep earthquakes between 8 and 16 km depth. In contrast, our process using a 1-D velocity model leads to a very shallow seismicity close to the surface. Using the 3-D tomographic model does not improve the hypocentre determination.

## 5.3 Accuracy for deep upper-crustal events

### 5.3.1 Argelès-Gazost sequence 2006

Fig. 11 shows the results of absolute earthquake location for the Argelès-Gazost sequence. 263 earthquakes are compared to the 3-D GT reference locations of Sylvander et al. (2008). In this case, azimuthal gaps are systematically lower than 180°, and the distance to the first station smaller than 10 km. On average, there are about 10 stations recording each earthquake for which *P* and *S*-waves picks are available. In such an optimal configuration, the pattern and clustering along the normal fault is well retrieved with all location processes with or without applying station delays. However, using station delays increases the accuracy. Compared to the reference, all NLloc processes, whatever the velocity model or the geographic reference, show a systematic shift of about 2 km northwards of the seismic cluster. The worst clustering is obtained
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Quarry blast absolute locations comparison between the eight processes analysed in this study (Table 1).</figcaption></figure>

with the HYPO-71 OMP routine. Moreover, the accuracy of the 3-D average model seems lower than using 1-D velocity models. Similar to the localization of the Luzenac quarry blasts, this sequence of seismicity is located at the boundary between the basement and the sedimentary basins (Fig. 11 and Supporting Information Fig. S1). This probably explains the slight horizontal shift of some earthquakes to the north (towards the basin) resulting in a more diffuse clustering.

### 5.3.2 Saint-Paul-de-Fenouillet sequence 1996

Fig. 12 shows the results for the Saint-Paul-de-Fenouillet sequence that occurred between 1996/02/18 and 1996/02/23 (5 d). 334 GT

reference earthquakes have been located by Pauchet *et al.* (1999) and Rigo *et al.* (1997). Most of the events have an azimuthal gap lower than 180°. The distance to the first recording seismic station is always less than 10 km. The HYPO-71 OMP routine gives similar results to the reference. The cluster is vertically extended between 9 km and 4 km depth and between 11 km and 6.5 along the N-S direction in the reference case. Moreover, two branches seem to emerge at the top of the cluster. Some of the earthquakes are located very close to the surface when using the 3-D average model. This feature is more limited with other processes. Only few earthquakes are located below 10 km depth but we may note that using 1-D velocity models the main shock (gap = 134°, distance to the first station = 16 km, not recorded by temporary network) seems
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Clustering analysis of the Lacq gas field seismicity over the period 1974–1997. Two local permanent local seismic networks were recording ground motions during this period. The 3-D absolute earthquake location using NLLoc in a high resolution 3-D P velocity model from Bardainne *et al.* (2008) has been used as reference for comparison. 1030 earthquakes are compared.</figcaption></figure>

badly located between 1 and 2 km below the seismic cluster and between 3.5 and 5 km below with other estimations.

### 5.3.3 Cauterets and Pic-du-Midi sequences 2002

Figs 13 and 14 show respectively the first sequence of 12 earthquakes located close to Cauterets (11 earthquakes are compared to the reference) and the second sequence of 3 earthquakes located close to the Pic-du-Midi. The relative earthquake location based on HYPO-DD using cross-correlation waveforms from Dubos *et al.* (2004) are used as reference for the first sequence and a HYPO71 process for the second reference. In both case, the azimuthal gap ranges between 95° and 210° and the distance to the first recording seismic station is about 8 km in all cases. The first cluster of seismicity is located between 7.5 and 10 km depth while the second cluster is located around 12 km depth. The first sequence is located

in the centre of the Palaeozoic area while the second sequence is nearest to the sedimentary basin.

In the first sequence, the NLLoc process using the 3-D tomographic model gives the best accuracy horizontally and vertically (accuracy ≪1 km for both). The worst cases are both obtained when using 1-D velocity models. The HYPO-71 OMP routine gives worst horizontal location (accuracy of about 2–3 km). For other 1-D processes, accuracy is about 1 km horizontally and up to 5 km in depth. Regarding the 3-D average model, horizontal accuracy is lower than 1 km while the vertical accuracy is about 1–2 km. Shallowest events are those with the highest azimuthal gap around 200° and have a vertical uncertainty of about 2.3 km, in agreement with accuracy.

In the second sequence, the NLLoc process using the 3-D average model gives the best accuracy horizontally and vertically (accuracy ≪1 km for both). In that case, the three events are also well aligned
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Clustering analysis of the Lacq gas field seismicity over the period 1998–2013. The previous IPGS (Institut de Physique du Globe de Strasbourg) local seismic network has been removed in 1998. 282 earthquakes have been recorded during this period and compared here with the reference described in Fig. 9 (bold grey line).</figcaption></figure>

along the preferred nodal plane. Using the 3-D tomographic model leads to the worst depth accuracy (between 1 km and 4 km) and about 1 km horizontally. The 1-D HYPO-71 OMP routine gives an horizontal accuracy of about 1 km and 2 km in depth. Others 1-D processes give a good horizontal accuracy (≪1 km) and a depth accuracy of about 1 km.

# 6 DISCUSSION

## 6.1 Modelling errors

Even if a probabilistic earthquake location is more reliable, the quality of the estimated absolute hypocentre locations depends essentially on the quality of the velocity model (e.g. Husen &

Hardebeck 2010, and references therein). It is difficult to quantify errors on the imperfect knowledge of the velocity structure. Currently, the only way to do this is to use active seismic data for which location and origin time of sources are known (Satriano et al. 2006; Lin et al. 2011; Gesret et al. 2015). In most cases the shallow velocities are not well constrained. As shown in this study, considering the geometry of sedimentary basins in addition to station delays greatly improves the accuracy of the final computed traveltime. Station delays are close to simulated arrival-time residuals computed in reference velocity models (Fig. 15). The new station delays associated with the 3-D velocity model dedicated to absolute earthquake location (Fig. 6) are centred on zero and well explain the simulated residuals observed along the reference seismic refraction profile (Figs 4 and 15). Their absolute values are
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Clustering analysis of the Argelès-Gazost sequence following the *M*<sub>L</sub> 5.0 2006/11/17 18h19 main shock. The 3-D earthquake location using SIMULPS and a 3-D velocity model (Dubos *et al.* 2004) from Sylvander *et al.* (2008) has been used as reference. 263 earthquakes have been compared here based on Sylvander *et al.* (2008).</figcaption></figure>

lower than 0.8 s. These delays are necessary to improve the absolute earthquake location and they mainly correct for sedimentary structure, that is, the location of the basement and lateral variation of the velocity structure related to the sediment deposition history.

Using only station delays without considering sedimentary basin geometry (1-D case) leads to badly computed arrival-times for refracted waves below sedimentary basins. This is well illustrated by the Lacq seismicity relocation (Fig. 10). However, this work would greatly benefit from accurate mapping of basement geometry and stratigraphy since they can introduce additional errors which have a great influence in the horizontal accuracy (example from

Luzenac quarry blasts and Argelès-Gazost sequence). Deep structures at the scale of the full network can be constrained by passive seismic tomography but generally the resolved domains are sparse. In some regions where the seismicity is dense and homogeneously distributed below a well-designed permanent seismic network, it is possible to use the result of 3-D seismic tomography to compute the absolute earthquake location (Lin *et al.* 2007). But in most cases, amplitude of velocity anomalies is not well resolved and the resolution of the inversion is not homogeneous in space. The resolved volumes of the crust depend very much on the ray-distribution (*P*<sub>g</sub> and *P*<sub>n</sub>), the aperture of the network and the quality of the
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Clustering analysis of the Saint-Paul-du-Fenouillet sequence following the *M*<sub>L</sub> 5.2 1996/02/18 01h45 main shock. The 1-D HYPO71 earthquake location from Pauchet *et al.* (1999) and Rigo *et al.* (1997) has been used as reference. 334 earthquakes have been compared.</figcaption></figure>

data. Globally in most areas the 3-D data coverage of a LET can be much better than the density of 2-D seismic lines. The computation of 3-D models from 2-D lines therefore always includes interpolation and/or extrapolation, which explains why the use of an averaging approach based on a 3-D LET is useful to improve the velocity structure. The layered crustal velocity field could therefore be used to avoid additional random errors from the velocity structure.

In addition to the velocity model accuracy, the precision of the traveltime computation introduces an additional modelling error. Precision of the finite difference algorithm of Podvin & Lecomte (1991) to calculate traveltimes depends on grid-spacing. In our case (cell size of 1 km³), we have shown that errors are mainly between -0.1 s and 0 s and it may increase up to -0.2 s meaning that the error on calculated traveltimes can be significant (slowing down the velocity structure). These errors are not homogeneous in space (Supporting Information Fig. S11) and are 2 to 5 times, or more, lower (Fig. 15) than errors from 3-D model inaccuracy. Indeed, misestimation of the traveltimes due to the 3-D model inaccuracy is only partially corrected by including station delays. As demonstrated by the estimated accuracy obtained with the Pyrenean seismic network, we may assume that the error in traveltime computation, in addition to measurement errors, do not significantly affect the results of the absolute earthquake location. Further tests should be performed with better algorithms to fully consider this technical issue which is a direct consequence of the finite difference approach (see Noble *et al.* 2014).

## 6.2 Accuracy and uncertainties

Recurring questions from structural geologists, geodynamicists or people from industry when using earthquake catalogues are related to quality parameters and uncertainties. They want to know how far solutions are from the true location, in other words, how far uncertainties are from accuracy. The largest contribution to location and uncertainty errors is due to the fact that the velocity model errors are usually not correctly taken into account. At the reservoir scale and based on active seismic sources, Gesret *et al.* (2015) use a Bayesian approach to estimate the uncertainties on velocity structures and then propagate these uncertainties in the earthquake location process. Except in Taiwan where large active experiments have been done (Lin *et al.* 2011), such an approach is difficult to apply in 3-D at the scale of a local or regional seismic network where only sparse independent active data are usually available. Limiting modelling errors related to relief, spherical earth shape, sedimentary basin, basement geometry and Moho geometry is needed to better approximate ray geometry within the earth. Assuming that the 3-D model and ray geometry are close enough to reality and considering Gaussian model error of 3 per cent allow us to propagate velocity model errors into the earthquake location process. Full PDF should provide fair estimates of the accuracy following such a deterministic approach. The analysis of some GT reference earthquakes shows that the PDF can vary a lot depending on the velocity model, the data set and the depth of the considered earthquake. Two cases must be considered:
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Clustering analysis of the 2002 May earthquake sequence in the central Pyrenees close to Cauterets.</figcaption></figure>

(i) *good network configuration:* when the primary azimuthal gap is lower than 180° and the distance to the first station is lower than 15 km (~distance to earthquake depth), a 1-D solution can be better than 3-D for shallow events (0–6 km). In this case, horizontal location of the hypocentre is slightly less well determined in 3-D than in 1-D even if the sedimentary basins are taken into account because of the imperfect knowledge of velocity variations within them. However, in all cases 3-D approaches provide better estimates of the hypocentre depth.

(ii) *poor network configuration:* when the distance to the first seismic station is larger than the hypocentre depth and/or the primary azimuthal gap is higher than 180°, 1-D solution should be regarded with suspicion, since focal depths can be systematically wrong. The 3-D model still gives sufficiently good results in this case. The best example is given by the Lacq seismicity for earthquakes recorded after 1997 (Fig. 10). In that case, considering the sedimentary basin geometry results in more realistic ray paths within the crust (i.e. refracted waves below the sedimentary basins for more distant seismic stations) and improves the estimations of hypocentre depth. A similar observation has previously been made when locating earthquakes offshore along the east coast of Taiwan where Moho geometry and velocity structure of the upper-mantle are necessary to estimate the hypocentre depth of distant earthquakes using refracted first arrivals (Theunissen et al. 2012b; Lallemand et al. 2013).

If we consider that modelling errors are small in the 3-D model, we can use uncertainties rather than network criteria to estimate the quality of the solution. Some GT reference earthquakes are

either located too close to the surface or too deep (see Lacq gas field seismicity or Saint-Paul-de-Fenouillet sequence). Analysis of uncertainties (SEH, SEZ, edd, edz, edh or updf) shows that most of these earthquakes exhibit large EDD or/and SEZ, about 13.4 ± 11.1 km and 7.7 ± 4.8 km respectively. In other words, the problem is ill-conditioned and the PDF does not show a clear, single, global minimum. Despite of this perfect correlation between inaccurate location and uncertainties, some earthquakes seem correctly located but are characterized by large edd or/and SEZ.

To summarize:

(i) solution close to the GT reference with small uncertainties: mostly the case within the 3-D model

(ii) solution close to the GT reference and large uncertainties: these solutions are not well enough constrained. In a quality selection process these earthquakes would be excluded as being potentially far from the true solution.

(iii) solution far from the GT reference and small uncertainties: only a few of these are observed within the 3-D model. Possibly, these earthquakes do not belong to the seismic cluster used as reference.

(iv) solution far from the GT reference and large uncertainties: badly located events mostly belong to this category as already discussed.

In cases where independent geological and geophysical data are not enough to constrain the velocity structure and station delays in a given region, a fully Bayesian approach considering a range of possible velocity structure for each seismic station and 1-D velocity models should be considered.
<figure><img src="image_17.png" /><figcaption>Figure 14. Clustering analysis of the 2002 May earthquake sequence in the central Pyrenees close to Pic du Midi.</figcaption></figure>

## 6.3 The final catalogue

We use results obtained with the dedicated 3-D *a-priori* velocity model and associated station delays as the final catalogue. The low resolution of the 3-D velocity model in particular in the shallowest crust and within the sedimentary basins leads to a slight deterioration in horizontal location (up to ~1 km) of the solutions for shallow events. Despite this observation, hypocentre depth, one of the most important parameters, is well determined in most cases even when the distance to the first seismic stations is larger than the hypocentre depth or when the azimuthal gap is important (>180°). Moreover, the 3-D velocity model takes advantage of previous geophysical and geological studies as well as the seismic tomography. Consequently modelling errors are small enough and we expect the uncertainties to be close to the optimal estimation of the accuracy. In addition to removing extreme outliers, modifications made to the NLLoc Octree algorithm enable homogeneous sampling around the solution. When doing this, quality criteria (SEH, SEZ, edh, edz, edd) are all comparable from one event to another. We use these criteria to select best located earthquakes.

A qualitative comparison of absolute earthquake locations between 1-D flat NLloc and 3-D geoc average NLloc processes (Fig. 16) for the entire catalogue reveals that they are similar as expected. In details, the seismicity distribution is different on boundaries of the seismic network. To the east, the seismicity is on average between 5 and 10 km depth (one standard deviation) while it is between the surface and 15 km depth based on the 1-D location. A similar observation can be made to the west of the Pyrenees where average seismicity is between 15 and 35 km depth while it is between the surface and 25 km depth based on the 1-D location.

In order to evaluate the value of the new catalogue, a comparison with the 1-D solutions is made (Fig. 16 and Supporting Information Fig. S12). It reveals some new important features. Seismic structures are better resolved in particular close to boundaries of the seismic network. The best example concerns a north dipping structure in the western part of the Pyrenees that extends further west which is well imaged using the 3-D approach. A large number of earthquakes concentrate on a plane from 5–10 km depth to 30 km depth with a width of at least 40 km in cross-sections 36 to 26 (Fig. 16). This plane is only visible on one cross-section on the OMP catalogue (Supporting Information Fig. S12). Some seismicity clusters are identified better with 3-D results, for example the deep seismicity cluster at 30 km depth on cross-section n26. Some other clusters localized in the central part of the seismic network are also identified better or also sometimes with a different geometry and location from the 1-D solution (for example the vertical seismic cluster on cross-section 67 (Fig. 16)). Finally, the statistical distribution of the seismicity from east to west (Fig. 16b) shows that the seismicity is more dispersed vertically on average using the 3-D approach. The maximum depth reached by the seismicity is relatively constant, about 15–20 km depth, in the central and eastern part of the Pyrenees. While the 1-D OMP results show a peak of seismicity at 30 km depth at 200 km (following the horizontal axis in Supporting Information Fig. S12), a totally different feature is seen in the 3-D results where a transition to deeper depth occurs between 200 and 150 km (east of Pampeluna) to reach a deep seismicity domain up to 30–35 km depth where a north dipping structure can be identified (Fig. 16b).
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Errors estimation based on simulated observed residuals: (observed traveltime) minus (calculated traveltime) and compared to calculated station delays and RMS (based on all runs realized in this study). We simulate observed residuals in perfect cases where the origin of error is known. Two types of usual errors are compared: (1) measurement errors of seismic arrival times due to picking errors, Finite Impulse Response (FIR) effect (Scherbaum & Bouin 1997), phase misidentification and clock error at a seismic station, and (2) modelling errors due to sphericity of the earth (see the Supporting Information and Fig. 4), the use of a 1-D model (Fig. 4), the use of a 3-D *a-priori* model (Fig. 4), and due to the traveltimes calculation precision (see Section 3.1 and Supporting Information Fig. S11). Ranges of simulated observed residuals are based on few examples shown in Fig. 4 and Supporting Information Fig. S11 and considered as representative. Abbreviations: 1Dflat, 1-D velocity model in flat earth reference; 1Dgeoc, in geocentric reference; 3Dgeoc, 3-D average velocity model dedicated to absolute earthquake location; 2D line 8, reference 2-D seismic refraction profile from Pedreira *et al.* (2003).</figcaption></figure>

# 7 CONCLUSIONS

Local and regional absolute earthquake location (0°–5°) has been performed on the data of the Pyrenean seismic network. The velocity structure has been described by a grid composed by 4 × 4 × 2 km³ cells. Results of 1-D versus 3-D approaches have been compared in order to estimate the contribution of using a 3-D model in source location.

Successful application of a 3-D approach (rather than 1-D) is possible if modelling errors (ray geometry and travel time estimation) are small enough. Using the results of a local tomography leads to inhomogeneous resolution in space and deteriorates absolute earthquake location in some cases. Consequently, a dedicated 3-D *a-priori* model with well estimated station delays constitutes a

better approach. In the case of the local Pyrenean seismic network, the 3-D velocity model takes advantage of a long list of previous geophysical and geological studies and several years of seismic recording.

In most favourable cases (gap<180° and distance to the first station <15 km ~ hypocentre depth), the 1-D solutions can be better than 3-D for shallow events (0–6 km). In this case, horizontal location of the hypocentre is slightly less well determined in 3-D than in 1-D, owing to lateral variations in sedimentary basins which are not captured in the 3-D model. However, in all cases the 3-D approach leads to better estimates of hypocentre depth. The accuracy is lower than 2 km vertically and horizontally.

In a poor network configuration, when the distance to the first seismic station is larger than the hypocentre depth and/or the
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Comparison of seismicity distribution between the final catalogue using the 3-D dedicated velocity model (3Dgeoc average NLloc process) and the 1-D flat NLloc process. Error bars correspond to ERH and ERZ. For the 3-D process, sorting is for SEH < 5 km, SEZ < 5 km and edd < 5 km (respectively 16724/24427 and 15126/24400 earthquakes for 3-D and 1-D selection). Error bars correspond to SEH and SEZ. (a) Examples of differences (yellow triangle) along some SWS-NEN cross-sections. (b) statistical distribution of the seismicity along the E-W direction based on sorted catalogues. The average is framed by one standard deviation (vertical bar) and by the shallowest and deepest earthquakes respectively min and max on the figure. Cross-sections numbering is given every 4 km.</figcaption></figure>

primary azimuthal gap is larger than 180°, 1-D solution should be regarded with suspicion. A 3-D model can provide more reliable locations in such case, provided that the data set of phase arrivals is sufficient and the velocity model is well constrained.

The new catalogue based on the 3-D approach reveals new features mainly because of improved hypocentre depth. The final seismicity is located with a horizontal uncertainty of about 2 ± 2 km and 3 ± 2 km vertically on average. In the western part of the seismic network we find a clear and continuous north-dipping plane between 10 and 30 km depth and 40 km width. The general seismicity distribution along the belt changes from a continuous crustal seismicity between 0 and 15 km (with some events between 15 and 20 km depths) in the eastern part of the Pyrenees to a domain of transition (~50 km width east of Pampeluna) to reach a deep seismicity region (30–35 km) to the west.

In order to further reduce modelling errors, we need a better velocity model of the shallowest part of the crust, that is, sedimentary basins, as well as an improved 3-D finite difference traveltime

computation code. Finally, our new location algorithms could be applied routinely by the Réseau de Surveillance Sismique des Pyrénées (RSSP) managed by the Observatoire Midi-Pyrénées (OMP).

## ACKNOWLEDGEMENTS

This study was initiated during an ATER (attaché temporaire d'enseignement et de recherche) position funded by the University Paul Sabatier III (Toulouse, France) at the Institut de Recherche en Astrophysique et Planétologie (IRAP, UMR5277) part of the Observatoire Midi-Pyrénées (OMP). This work was funded by EDF (Electricité de France) under the SIGMA (Seismic Ground Motion Assessment) project. This work was included in the work package 1 (WP1), *The seismic source in hazard models*. We acknowledge Kevin Manchuel, leader of WP1. We also acknowledge people involved in the TOPO-IBERIA and PYROPE projects for providing us with the high quality data from these temporary
seismic deployments. We want also to acknowledge the institutions responsible for the maintenance of the different permanent seismic networks contributing data to this study: Instituto Geográfico Nacional of Spain (IGN), Catalanian Geological Institute (IGC), Commissariat à l'Energie Atomique (CEA), French Broad-band network (RLBP) and RAP French strong motion network provided by RESIF - Réseau Sismologique et géodésique Français. Seismic Network doi:10.15778/RESIF.FR. RESIF is a national Research Infrastructure, recognised as such by the French Ministry of Higher Education and Research. RESIF is managed by the RESIF Consortium, composed of 18 Research Institutions and Universities in France. RESIF additionally supported by a public grant overseen by the French National Research Agency (ANR) as part of the 'Investissements d'Avenir' program (reference: ANR-11-EQPX-0040) and the French Ministry of Ecology, Sustainable Development and Energy. We acknowledge Thomas Bardainne (CGG - Microseismic) and Mario Ruiz Fernandez (Institute of Earth Sciences Jaume Almera of the Spanish Scientific Research Council (ICTJA-CSIC)) for providing us access to traveltime picks and velocity models from temporary seismic experiments. We acknowledge Robert Sarda (IMERYS) for giving information about blasts at the Luzenac quarry. Thanks to Leo Zijerveld for his suggestions and help with the English. Authors thank Prof Egill Hauksson, editor and two anonymous reviewers for their constructive reviews and careful handling of the manuscript.

## REFERENCES

Aki, K. & Lee, W., 1976. Determination of 3-Dimensional velocity anomalies under a seismic array using 1ST-P arrival times from local earthquakes: 1. A homogeneous initial model, *J. geophys. Res.*, **81**(23), 4381-4399.

Aki, K. & Richards, P.G., 1980. *Quantitative Seismology: Theory and Methods*, vols I and II, University Science Books, pp. 331-349.

Amante, C. & Eakins, B.W., 2009. Etopol1 arc-minute global relief model: Procedures, data sources and analysis, Tech. Rep. 24, NOAA.

Asensio, E., Khazaradze, G., Echeverria, A., King, R.W. & Vilajosana, I., 2012. GPS studies of active deformation in the Pyrenees, *Geophys. J. Int.*, **190**(2), 913-921.

Bardainne, T., Dubos-Sallée, N., Sénéchal, G., Gaillot, P. & Perroud, H., 2008. Analysis of the induced seismicity of the lacq gas field (southeastern France) and model of deformation, *Geophys. J. Int.*, **172**(3), 1151-1162.

Benz, H., Chouet, B., Dawson, P., Lahr, J., Page, R. & Hole, J., 1996. Three-dimensional P and S wave velocity structure of Redoubt Volcano, Alaska, *J. geophys. Res.*, **101**(B4), 8111–8128.

Biteau, J.-J., Le Marrec, A., Le Vot, M. & Masset, J.-M., 2006. The Aquitaine Basin, *Pet. Geosci.*, **12**(3), 247-273.

Bondar, I. & McLaughlin, K., 2009. Seismic location bias and uncertainty in the presence of correlated and non-Gaussian travel-time errors, *Bull. seism. Soc. Am.*, **99**(1), 172-193.

Bondar, I., Myers, S., Engdahl, E. & Bergman, E., 2004. Epicentre accuracy based on seismic network criteria, *Geophys. J. Int.*, **156**(3), 483-496.

Cara, M., Alasset, P.J. & Sira, C., 2008. Magnitude of historical earthquakes, from macroseismic data to seismic waveform modelling: application to the Pyrenees and a 1905 earthquake in the Alps, in *Historical Seismology: Interdisciplinary Studies of Past and Recent Earthquakes*, vol. 2 of *Modern Approaches in Solid Earth Sciences*, pp. 369-384, Springer.

Chang, A.C., Shumway, R.H., Blandford, R.R. & Barker, B.W., 1983. Two methods to improve location estimates-preliminary results, *Bull. seism. Soc. Am.*, **73**, 281-295.

Chantraine, J. *et al.*, 1996. *Carte géologique de la France 1:1000000*, 6e edn.

Chevrot, S., Sylvander, M. & DelLouis, B., 2011. A preliminary catalog of moment tensors for the Pyrenees, *Tectonophysics*, **510**(1-2), 239-251.

Chevrot, S. et al., 2014. High-resolution imaging of the Pyrenees and massif central from the data of the pyrope and iberarray portable array deployments, *J. geophys. Res.*, **119**(8), 6399-6420.

Chevrot, S., Sylvander, M., Diaz, J., Ruiz, M., Paul, A. & the PYROPE Working Group, 2015. The Pyrenean architecture as revealed by teleseismic P-to-S converted waves recorded along two dense transects, *Geophys. J. Int.*, **200**(2), 1096-1107.

Chevrot, S., Sylvander, M. & RESIF, 2017. Seismic network X7:PYROPE PYRenean Observational Portable Experiment (RESIFSISMOB); RESIF - Réseau Sismologique et géodésique Français. doi:10.15778/RESIF.X72010.

Crosson, R., 1976. Crustal structure modeling of earthquake data: I. Simultaneous least-squares estimation of hypocenter and velocity parameters, *J. geophys. Res.*, **81**(17), 3036-3046.

Daignières, M., Gallart, J., Banda, E. & Hirn, A., 1982. Implications of the seismic structure for the orogenic evolution of the Pyrenean range, *Earth planet. Sci. Lett.*, **57**(1), 88-100.

Diaz, J. & Gallart, J., 2009. Crustal structure beneath the Iberian Peninsula and surrounding waters: a new compilation of deep seismic sounding results, *Phys. Earth planet. Inter.*, **173**(1-2), 181-190.

Díaz, J. *et al.*, 2009. The IBERARRAY broadband seismic network: a new tool to investigate the deep structure beneath Iberia, *ORFEUS Newsletter*, **8**(2), 1-6.

Díaz, J., Pedreira, D., Ruiz, M., Pulgar, J.A. & Gallart, J., 2012. Mapping the indentation between the Iberian and Eurasian plates beneath the Western Pyrenees/Eastern Cantabrian Mountains from receiver function analysis, *Tectonophysics*, **570**, 114-122.

Dubos, N., Sylvander, M., Souriau, A., Ponsolles, C., Chevrot, S., Fels, J. & Benahmed, E., 2004. Analysis of the 2002 May earthquake sequence in the central Pyrenees, consequences for the evaluation of the seismic risk at Lourdes, France, *Geophys. J. Int.*, **156**(3), 527-540.

Eberhart-Phillips, D. & Michael, A.J., 1993. Three-dimensional velocity structure, seismicity, and fault structure in the Parkfield region, central California, *J. geophys. Res.*, **98**(15), 737-758.

Filleaudeau, P.-Y., 2011. Croissance et dénudation des Pyrénées du Crétacé Supérieur au Paléogène: Apports de l'analyse de bassin et thermochronométrie détritique, *PhD thesis*, Université Pierre et Marie Curie.

Flanagan, M.P., Myers, S.C. & Koper, K.D., 2007. Regional travel-time uncertainty and seismic location improvement using a three-dimensional a priori velocity model, *Bull. seism. Soc. Am.*, **97**(3), 804-825.

Font, Y., Kao, H., Lallemand, S., Liu, C. & Chiao, L., 2004. Hypocentre determination offshore of eastern Taiwan using the Maximum Intersection method, *Geophys. J. Int.*, **158**(2), 655-675.

Frohlich, C., 1979. Efficient method for joint hypocenter Determination for large groups of earthquakes, *Comput. Geosci.*, **5**(3-4), 387-389.

Gallart, J., Banda, E. & Daignieres, M., 1981. Crustal structure of the Paleozoic axial zone of the Pyrenees and transition to the North-Pyrenean zone, *Ann. Geophys.*, **37**(3), 457-480.

Gautier, S., Latorre, D., Virieux, J., Deschamps, A., Skarpelos, C., Sotiriou, A., Serpetsidaki, A. & Tselentis, A., 2006. A new passive tomography of the aigion area (Gulf of Corinth, Greece) from the 2002 data set, *Pure appl. Geophys.*, **163**(2-3), 431-453.

Geiger, L., 1910. Herdbestimmung bei Erdbeben aus den Ankunftszeiten, *Nachrichten von der Koniglichen Gesellschaft der Wissenschaften zu Göttingen*, **4**, 557pp and 373pp.

Gesret, A., Desassis, N., Noble, M., Romary, T. & Maisons, C., 2015. Propagation of the velocity model uncertainties to the seismic event location, *Geophys. J. Int.*, **200**(1), 52-66.

Guyoton, F., Grasso, J.-R. & Volant, P., 1992. Interrelation between induced seismic instabilities and complex geological structure, *Geophys. Res. Lett.*, **19**(7), 705-708.

Huang, Z. & Zhao, D., 2012. Relocating the 2011 Tohoku-oki earthquakes (m 6.0-9.0), *Tectonophysics*, **586**, 35-45.

Husen, S. & Hardebeck, J., 2010. Earthquake location accuracy, Community Online Resource for Statistical Seismicity Analysis, p. 35pp. doi:10.5078/corssa-55815573. Available at: http://www.corssa.org.
Husen, S., Kissling, E., Flueh, E. & Asch, G., 1999. Accurate hypocentre determination in the seismogenic zone of the subducting Nazca Plate in northern Chile using a combined on-/offshore network, *Geophys. J. Int.*, **138**(3), 687–701.

Husen, S., Kissling, E., Deichmann, N., Wiemer, S., Giardini, D. & Baer, M., 2003. Probabilistic earthquake location in complex three-dimensional velocity models: application to Switzerland, *J. geophys. Res.*, **108**(B2), 2077, doi:10.1029/2002JB001778.

Husen, S., Kissling, E. & Clinton, J.F., 2011. Local and regional minimum 1D models for earthquake location and data quality assessment in complex tectonic regions: application to Switzerland, *Swiss J. Geosci.*, **104**(3), 455–469.

Kissling, E., 1988. Geotomography with local earthquake data, *Rev. Geophys.*, **26**(4), 659–698.

Kissling, E., Ellsworth, W., Eberhart-Phillips, D. & Kradolfer, U., 1994. Initial reference models in local earthquake tomography, *J. geophys. Res.*, **99**(B10), 19 635–19 646.

Kissling, E., Solarino, S. & Cattaneo, M., 1995. Improved seismic velocity reference model from local earthquake data in northwestern Italy, *Terra Nova*, 7(5), 528–534.

Lahr, J., 1999. revised 2012, HYPOELLIPSE: A computer program for determining local earthquake hypocentral parameters, magnitude, and first motion pattern, no. 99-23 in *Reports-Open file series*, U.S. Geological Survey.

Lahr, J.C., 1989. Hypoellipse/version 2.0: A computer program for determining local earthquake hypocentral parameters, magnitude and first motion pattern, U.S. Geol. Surv., 92 pp.

Lallemand, S., Theunissen, T., Schnürle, P., Lee, C.-S., Liu, C.-S. & Font, Y., 2013. Indentation of the Philippine sea plate by the Eurasia plate in Taiwan: details from recent marine seismological experiments, *Tectonophysics*, **594**, 60–79.

Latorre, D., Virieux, J., Monfret, T., Monteiller, V., Vanorio, T., Got, J. & Lyon-Caen, H., 2004. A new seismic tomography of Aigion area (Gulf of Corinth, Greece) from the 1991 data set, *Geophys. J. Int.*, **159**(3), 1013–1031.

Lee, W.H.K. & Lahr, J.C., 1972. Hypo-71: A computer program for determining hypocenter, magnitude and first motion pattern of local earthquakes, U.S. Geol. Surv.

Le Meur, E., Virieux, J. & Podvin, P., 1997. Seismic tomography of the Gulf of Corinth: a comparison of methods, *Ann. Geophys.*, **40**(1), doi:10.4401/ag-3931.

Lin, G., Shearer, P. & Fialko, Y., 2006. Obtaining absolute locations for quarry seismicity using remote sensing data, *Bull. seism. Soc. Am.*, **96**(2), 722–728.

Lin, G., Shearer, P.M., Hauksson, E. & Thurber, C.H., 2007. A three-dimensional crustal seismic velocity model for southern California from a composite event method, *J. geophys. Res.*, **112**, B11306, doi:10.1029/2007JB004977.

Lin, Y.-P., Zhao, L. & Hung, S.-H., 2011. Assessment of tomography models of Taiwan using first-arrival times from the TAIGER active-source experiment, *Bull. seism. Soc. Am.*, **101**(2), 866–880.

Locati, M., Rovida, A., Albini, P. & Stucchi, M., 2014. The ahead portal: a gateway to European historical earthquake data, *Seism. Res. Lett.*, **85**(3), 727–734.

Lomax, A., 2005. A reanalysis of the hypocentral location and related observations for the great 1906 California earthquake, *Bull. seism. Soc. Am.*, **95**(3), 861–877.

Lomax, A., Virieux, J., Volant, P. & Berge-Thierry, C., 2000. Probabilistic earthquake location in 3D and layered models - Introduction of a Metropolis-Gibbs method and comparison with linear locations, in *Advances in Seismic Event Location*, pp. 101–134, *Shalheveth Freier Ctr Peace Sci & Technol.*, Springer.

Lozano, J.F., 2012. Cenozoic deformation of Iberia, A model for intraplate mountain building and basin development based on analogue modelling, *PhD thesis*, Universidad Complutense de Madrid.

Monteiller, V., 2005. Tomographie à l'aide de décalages temporels d'ondes sismiques P : développements méthodologiques et ap-

plications, *PhD thesis*, Savoie University, 154 pp., https://tel.archives-ouvertes.fr/tel-00653731.

Monteiller, V., Got, J., Virieux, J. & Okubo, P., 2005. An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the Kilauea volcano, *J. geophys. Res.*, **110**, B12306, doi:10.1029/2004JB003466.

Moser, T., Vaneck, T. & Nolet, G., 1992. Hypocenter determination in strongly heterogeneous earth models using the shortest-path method, *J. geophys. Res.*, **97**(B5), 6563–6572.

Mostaccio, A., Tuve, T., Patane, D., Barberi, G. & Zuccarello, L., 2013. Improving Seismic Surveillance at Mt. Etna Volcano by Probabilistic Earthquake Location in a 3D Model, *Bull. seism. Soc. Am.*, **103**(4), 2447–2459.

Myers, S. & Schultz, C., 2000. Improving sparse network seismic location with Bayesian kriging and teleseismically constrained calibration events, *Bull. seism. Soc. Am.*, **90**(1), 199–211.

Noble, M., Gesret, A. & Belayouni, N., 2014. Accurate 3-D finite difference computation of traveltimes in strongly heterogeneous media, *Geophys. J. Int.*, **199**(3), 1572–1585.

Paige, C. & Saunders, M., 1982. Lsqr - an algorithm for sparse linear-equations and sparse least-squares, *ACM Trans. Math. Softw.*, **8**(1), 43–71.

Pauchet, H., Rigo, A., Rivera, L. & Souriau, A., 1999. A detailed analysis of the February 1996 aftershock sequence in the eastern Pyrenees, France, *Geophys. J. Int.*, **137**(1), 107–127.

Pedreira, D., Pulgar, J., Gallart, J. & Diaz, J., 2003. Seismic evidence of Alpine crustal thickening and wedging from the western Pyrenees to the Cantabrian Mountains (north Iberia), *J. geophys. Res.*, **108**(B4), 2204, doi:10.1029/2001JB001667.

Pedreira, D., Pulgar, J.A., Gallart, J. & Torne, M., 2007. Three-dimensional gravity and magnetic modeling of crustal indentation and wedging in the western Pyrenees-Cantabrian Mountains, *J. geophys. Res.*, **112**, B12405, doi:10.1029/2007JB005021.

Podvin, P. & Lecomte, I., 1991. Finite-difference computation of traveltimes in very contrasted velocity models - a massively parallel approach and it associated tools, *Geophys. J. Int.*, **105**(1), 271–284.

Pujol, J., 1988. Comments on the joint determination of hypocenters and station corrections, *Bull. seism. Soc. Am.*, **78**(3), 1179–1189.

Rigo, A., Pauchet, H., Souriau, A., Gresillaud, A., Nicolas, M., Olivera, C. & Figueras, S., 1997. The February 1996 earthquake sequence in the eastern Pyrenees: first results, *J. Seismol.*, **1**(1), 3–14.

Rigo, A. *et al.*, 2015. Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011, *Geophys. J. Int.*, **201**(2), 947–964.

Satriano, C., Zollo, A., Capuano, P., Russo, G., Vanorio, T., Caielli, G., Lovisa, L. & Moretti, M., 2006. A 3D velocity model for earthquake location in Campi Flegrei area: application to the 1982–84 uplift event, in *Geophysical Exploration of the Campi Flegrei (Southern Italy) Caldera' Interiors: Data, Methods and Results*, pp. 38–49, eds Zollo, A., Capuano, P. & Corciulo, M., Doppiavoce.

Scherbaum, F. & Bouin, M., 1997. FIR filter effects and nucleation phases, *Geophys. J. Int.*, **130**(3), 661–668.

Serrano, O., Delmas, J., Hanot, F., Vially, R., Herbin, J., Huel, P. & Tourlière, B., 2006. *Le Bassin d'Aquitaine: valorisation des données sismiques, cartographie structurale et potentiel pétrolier*, Rapports régionaux d'évaluation pétrolière de l'IFP, BRGM.

Simmons, N.A., Myers, S.C., Johannesson, G. & Matzel, E., 2012. LLNL-G3Dv3: Global P wave tomography model for improved regional and teleseismic travel time prediction, *J. geophys. Res.*, **117**, B10302, doi:10.1029/2012JB009525.

Snoke, J.A. & Lahr, J.C., 2001. Locating earthquakes: At what distance can the earth no longer be treated as flat?, *Seismol. Res. Lett.*, **72**(5), 538–541.

Souriau, A. & Granet, M., 1995. A tomographic study of the lithosphere beneath the PYRENEES from local and teleseismic data, *J. geophys. Res.*, **100**(B9), 18 117–18 134.

Souriau, A. & Pauchet, H., 1998. A new synthesis of Pyrenean seismicity and its tectonic implications, *Tectonophysics*, **290**(3–4), 221–244.
Souriau, A., Sylvander, M., Rigo, A., Fels, J., Douchain, J. & Ponsolles, C., 2001. Pyrenean Tectonics: main seismological constraints, *Bull. soc. Géol. France*, **172**(1), 25–39.

Spakman, W. & Nolet, G., 1988. Imaging algorithms, accuracy and resolution in delay time tomography, in *Mathematical Geophysics*, pp. 155–188, eds Vlaar, N.J., Nolet, G., Wortel, M.J.R. & Cloetingh, S.A.P.L., Reidel.

Spencer, C. & Gubbins, D., 1980. Travel-time inversion for simultaneous earthquake location and velocity structure determination in laterally varying media, *Geophys. J. R. astr. Soc.*, **63**(1), 95–116.

Stucchi, M. et al., 2013. The SHARE European Earthquake Catalogue (SHEEC) 1000–1899, *J. Seismol.*, **17**(2), 523–544.

Sylvander, M., Monod, B., Souriau, A. & Rigo, A., 2007. Analysis of an earthquake swarm (May 2004) in the French eastern Pyrenees: towards a new tectonic interpretation of the Saint-Paul-de-Fenouillet earthquake (1996), *C. R. Geosci.*, **339**(1), 75–84.

Sylvander, M., Souriau, A., Rigo, A., Tocheport, A., Toutain, J.-P., Ponsolles, C. & Benahmed, S., 2008. The 2006 November, M(L) = 5.0 earthquake near Lourdes (France): new evidence for NS extension across the Pyrenees, *Geophys. J. Int.*, **175**(2), 649–664.

Tarantola, A. & Valette, B., 1982. Generalized non-linear inverse problems solved using the least-squares criterion, *Rev. Geophys.*, **20**(2), 219–232.

Theunissen, T., Font, Y., Lallemand, S. & Gautier, S., 2012a. Improvements of the maximum intersection method for 3D absolute earthquake locations, *Bull. seism. Soc. Am.*, **102**(4), 1764–1785.

Theunissen, T., Lallemand, S., Font, Y., Gautier, S., Lee, C.-S., Liang, W.-T., Wu, F. & Berthet, T., 2012b. Crustal deformation at the southernmost part of the Ryukyu subduction (East Taiwan) as revealed by new marine seismic experiments, *Tectonophysics*, **578**(SI), 10–30.

Thurber, C., 1992. Hypocenter velocity structure coupling in local earthquake tomography, *Phys. Earth planet. Inter.*, **75**(1–3), 55–62.

Vanorio, T., Virieux, J., Capuano, P. & Russo, G., 2005. Three-dimensional seismic tomography from P wave and S wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera, *J. geophys. Res.*, **110**, B03201, doi:10.1029/2004JB003102.

Wittlinger, G., Herquel, G. & Nakache, T., 1993. Earthquake location in strongly heterogeneous media, *Geophys. J. Int.*, **115**(3), 759–777.

Zhou, H., 1994. Rapid 3-dimensional hypocentral determination using a master station method, *J. geophys. Res.*, **99**(B8), 15 439–15 455.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Table S1.** Data used in this study. OMP: Observatoire Midi-Pyrenees; ReNass: Réseau National de Surveillance Sismique; CEA: Commissariat à l'Energie Atomique; RAP: Réseau Accélérométrique Permanent; BRGM: Bureau de Recherche Géologique et Minière; IGN: Instituto Geographico Nacional; IGC: Institut Geologic de Catalunya; RLBP: Réseau Large Bande Permanent; PYROPE: http://w3.dtp.obsmip.fr/RSSP/PYROPE/; IBERARRAY: http://www.ictja.csic.es/gt/rc/LSD/PRJ/index TOPOIBERIA.html.

**Table S2.** 1-D seismic velocity model obtained following a usual approach based on VELEST 3.1 program [Kissling 1988; Kissling et al. 1994, 1995].

**Table S3.** Seismic stations compiled and used in this study.

**Figure S1.** Seismic events used as references. Up: Three quarry blasts. Down: Three seismic sequences for which temporary seismic stations have been deployed.

**Figure S2.** Raw data set overview and main selection based on hodochrones (lines 1 and 2), residuals distribution (third line) and wadati diagram (on the right). Calculations are made based on a 1D NonLinLOc process with flat earth approximation and station corrections.

**Figure S3.** Iterative multiscale approach for reaching the optimal grid search on which the complete Oct-Tree search must start. numxinit, numyinit and numzinit: initial number of cells in each direction given by the user (10 × 10 × 6 here). *dx*, *dy* and *dz*: cell size. The final Oct-Tree search is done in a grid of 60 × 60 × 59 km³ initially subdivided by 20 × 20 × 24 = 9600 cells, respectively in *x*,*y* and *z*. The Oct-Tree search is running until reach a maximum number of 200000 nodes browsed or until reach a cell smaller than 30 m.

**Figure S4.** $V_p$ median profile below Lacq (Adour-Arzacq Basin) based on the 3-D velocity model of Guyoton et al. (1992).

**Figure S5.** Compiled data set used for the Mohorovicic discontinuity construction.

**Figure S6.** RMS and uncertainties from NLLoc earthquake location processes using station corrections.

**Figure S7.** RMS and uncertainties from NLLoc earthquake location results.

**Figure S8.** Earthquake location differences between the 4 results M1 (flat earth, 1-D velocity models, station delays), M2 (ellipsoidal earth, 1-D velocity models, station delays), M3 (ellipsoidal earth, 3-D average models, stations delays) and M4 (ellipsoidal earth, 3-D tomographic models, NO station delays).

**Figure S9.** Representation of the Probability Density Function (PDF) for the Luzenac quarry blast n°1. Comparison is made between processes that use station delays.

**Figure S10.** Uncertainties analysis and comparison between the seven processes analysed in this study for the three quarry blasts.

**Figure S11.** Error on travel times estimation for two 2-D cases. A: left: 1D velocity model with constant velocity layers from Velest [Kissling 1988; Kissling et al. 1994, 1995] (see Section 4 in the paper). Right: difference in travel times estimation between a grid composed of large cubic cells of 1 km³ and a grid composed of small cubic cells of 20 m³. The initial maximum error of such small size is lower than $10^{-2}$ s, this last is thus used as reference to describe the error on travel times calculation. B: left: 2-D velocity model from seismic refraction forward modelling (line 8) [Pedreira et al. 2003]. This model is initially defined in a grid composed with a cell size of 500 m × 50 m. Right: difference in travel times estimation between a grid composed of large cubic cells of 1 km³ and a grid composed of small cubic cells of 50 m³.

**Figure S12.** A: seismicity distribution of the 1-D flat OMP routine on few cross-sections across the seismic network. Sorting if for ERH < 5 km and ERZ < 5 km (20503/24529 eartquakes). Error bars correspond to ERH and ERZ.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.