Very Fast Characterization of Focal Mechanism Parameters Through W-Phase Centroid Inversion in the Context of Tsunami Warning

Most of the tsunami potential seismic sources in the NEAM region are in a magnitude range of 6.5≤Mw≤7.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6.5 \le M_{w} \le 7.5$$\end{document} (e.g. the tsunami triggered by the Boumerdes earthquake of 2003 with Mw=6.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{w}=6.9$$\end{document}). The CENtre d’ALerte aux Tsunamis (CENALT), in operation since 2012 as the French National Tsunami Warning Centre (NTWC) and Candidate Tsunami Service Provider (CTSP), has to issue warning messages within 15 min of earthquake origin time. These warnings are based on the seismic source parameters (Mw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{w}$$\end{document} magnitude, focal depth and type of fault), which are computed by focal mechanisms and centroid inversion methods. The W-phase method, developed by Kanamori and Rivera, allows quick computation of seismic source parameters due to the early arrival time between P-waves and surface waves, and is therefore particularly useful for monitoring. We assess the W-phase method with 29 events of magnitude Mw≥\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_w \ge$$\end{document} 5.8 for the period 2010–2015 in the NEAM region. Results with 10 min of signal length are in good agreement compared to the Global Centroid Moment Tensor (GCMT) catalog.


Introduction
Since the occurrence of the Sumatra mega earthquake in December 2004 (M w ¼ 9:1) that had over 200,000 casualties, UNESCO has coordinated the implementation of tsunami warning systems in all oceanic basins that could be affected by tsunami waves. In the past, large tsunamis have been generated in the Mediterranean sea and the Atlantic ocean-for example the 1755 Lisbon earthquake ROGER et al. (2010) and the 1908 Messina earthquake TINTI et al. (1999) both generated tsunamis that claimed over 10,000 victims. More recently, the 2003 Boumerdes earthquake of Magnitude 6.9 generated a tsunami that hit the Balearic coast 40 min after origin time with waves reaching 2 m in amplitude (ALASSET et al. 2006;SAHAL et al. 2009). Although no victims from the tsunami were reported, the tsunami has caused widespread infrastructural damage to harbours, roads and housing.
In past decades the Mediterranean and Atlantic coasts have seen rapid demographic and infrastructural growth due to tourism and economic development, making these regions increasingly vulnerable to tsunami-related infrastructural damage and causualities. The CENALT (ROUDIL et al. 2013;SCHINDELÉ et al. 2015) is the French National Tsunami Warning Centre (NTWC) and CTSP for the North East Atlantic, Mediterranean and adjacent sea (NEAM) region. The CENALT is responsible for issuing alerts in cases of tsunami risk to the French authorities, international warning centers, and tsunami focal points in the NEAM region within 15 min of the earthquake origin time. A tsunami warning system includes two phases : the first phase consists in the determination of the sources parameters (location, magnitude, focal mechanism). The second phase, rely on the tsunami wave propagation forecast GAILLER et al. (2013). The warning level is currently based on a decision matrix depending on the magnitude and the location of the hypocenter. Most of tsunamigenic seismic sources in the NEAM region are in the magnitude range between 6.0 and 7.5 and the response time is very short as a tsunami can reach the balearic islands within 40 min of the earthquake origin time. In this context, it is important to have a fast and precise method to analyze the seismic source to assess the tsunami risk.
The W-phase algorithm (KANAMORI and RIVERA 2008;HAYES et al. 2009;DUPUTEL et al. 2012) has been used at the CENALT since 2012. Although the method was initially designed to characterize large earthquakes M w ! 7:0, many tsunami warning centres world-wide use the W-phase algorithm for smaller magnitudes. The USGS has applied the W-phase method to invert seismic sources as small as M w ¼ 5:8, and results have shown good agreement with the GCMT solutions HAYES et al. (2009).
In this paper, we applied the W-phase algorithm on 29 seismic events in the NEAM region of M w ! 5:8, selected as part of the European ASTARTE project, for the period 2010-2015. To assess the accuracy of the W-phase method, seismic traces are processed in time windows of 10 and 20 min from the origin time. The W-phase algorithm is used in a regional mode which means that only stations between 2 and 30 are taken into account and a 1D regional model is used for the Green's functions calculation. In addition to the magnitude and the location used for the decision matrix, we also compute the focal mechanism which can be useful to distinguish normal and thrust faults from strike-slip faults as the latter ones are less prone to tsunami generation. Results are then compared to the GCMT catalog (DZIEWONSKI and ANDERSON 1981;EKSTRÖ M et al. 2005).

Method
The W-phase is a very long period [100; 1000] s and high group velocity [4.5; 9] km s -1 phase that has been first observed by KANAMORI (1993) on the vertical displacement record of the 1992 Nicaragua earthquake (M w ¼ 7:7). It represents the total elastic field (near-field and far-field terms). From a ray theory point of view, the W-phase can be interpreted as the superposition of different phases such as P, PP, S, SS, SP, etc. From a normal mode theory point of view, the W-phase can be interpreted as the superposition of the fundamental, 1st, 2nd, and 3rd overtones of the spheroidal modes.
After deconvolution in displacement and bandpass filtering, the W-phase is observable on the 3 components of large-band records CUMMINS (1997). Figure 1 shows the W-phase observed and the corresponding synthetic for the 2015/02/13 Atlantic ridge M = 7.1 earthquake. To cope with a low signal-tonoise ratio (SNR) at smaller magnitudes, a moving frequency band relying on the preliminary magnitude of the event is applied (HAYES et al. 2009;DUPUTEL et al. 2012). The frequency band is shifted towards shorter periods as the magnitude decreases (Fig. 2).
The W-phase algorithm is a centroid moment tensor inversion method based on the so called Wphase. At CENALT, the W-phase algorithm is integrated to the software tools SeisComP3 WEBER et al. (2007), which gives the preliminary determination of epicenters (PDE). To improve the response time, the W-phase algorithm has been configured to operate on a regional scale which involves two modifications compared to a teleseismic configuration: (1) the Green's functions database is computed for a distance range of 0 D 30 with the discrete wave number method (DWNM) HERRMANN (2013) along with a European 1D model, provided by the CENALT, and (2) the time window is set to take into account stations as close as 2 . The W-phase time window starts at the theoretical P-wave arrival and ends before the surface waves. For small distances, however, the Wphase signal is too short to be used. In such cases the time window of the W-phase is widened between 0 D 12 following the relation: Figure 1 Comparison of observed W-phase (black lines) and the corresponding synthetics (red lines) of the BHN and BHE components at the station LOR from the RD (CEA/DASE) network for the 2015-2-13 Mw = 7.1 Atlantic earthquake. The two red dots correspond to the time window over which the W-phase is inverted Figure 3 shows the difference between the Wphase time window uses in a regional and global scale. Stations within 0 D 21 are available at t 0 þ 10 min and stations within 0 D 46 are available at t 0 þ 20 min. The main drawback of widened the W-phase time window between 2 and 12 is the risk to include a saturated signals in the time window, especially at very short distances. Beyond 12 , the end time of the W-phase time window is set in order to avoid the fundamental branch of the surface waves which are more sensitive to the crust and so to the shallow heterogeneities.
The moment tensor inversion is performed assuming no net volume change (M rr þ M hh þ M // ¼ 0), the inverse problem is then a linear problem solved by the least square method: where m corresponds to the 5 elements of the moment tensor, the 3 spatial coordinates, and a centroid time. u corresponds to the concatenated observed W-phase at every station used in the process and U are the kernel functions resulting from the convolution of the Green functions and a triangular source time function S(t) characterizes by an half duration (h d ) and a time shift (s s ). The half duration is estimated by a scaling law from the preliminary seismic moment tensor M 0 EK-STRÖ M and ENGDAHL (1989) : with M 0 in dyn cm and h d in s. In real-time, the preliminary magnitude, source location and origin time are determined by body wave arrivals from SeisComP3 between 2 and 9 min of the earthquake origin time, depending on the location and the stations coverage (t 3 min for a Mediterranean earthquake and t 9 min for an Atlantic earthquake). A poor estimate of the preliminary magnitude can lead to a poor half duration estimation. To circumvent this issue, the W-phase inversion is executed twice. During the first execution, only the grid search in time is done, which gives us a new seismic source half-duration. The PDE, updated with the half-duration, is then used as an input for the second inversion during which a gridsearch in time and space is done.

Data and Analysis
29 recent earthquakes of magnitude M w ! 5:8, with different types of focal mechanisms, inside the NEAM region, occurred from 2010 to 2015 were The corner frequencies used for bandpass filtering (butterworth, 4th order, causal) depends on the preliminary magnitude obtained by SeisComP3 Figure 3 W-phase time window for regional inversion. The bottom red line corresponds to the first P waves arrivals (t p ). The top red line corresponds to the end of the W-phase time window in a regional mode. The red dash line corresponds to the end of the W-phase time window used for a global scale inversion. The vertical black dash linecorreponds to the maximum epicentral distance that can reach with 10 min of signal Vol. 173, (2016) Very fast characterization of focal mechanism parameters 3883 selected by the ASTARTE partners based on the GCMT catalog. A list of the 29 earthquakes and their respective characteristics are reported in Table 1. Figure 4 shows the geographical location of the 29 earthquakes. They are distributed in 4 geographical areas : Atlantic (green), Mediterranean (pink), North (blue), East (brown). Figure 5 shows the seismic station network used for the W-phase source inversion for which a good azimuthal coverage and a high signal to noise ratio (SNR) were the main criteria. In total, 172 seismic stations from 25 seismic networks have been taken into account. All stations in Fig. 5 have been used at least once. The same seismic station network have been requested for every earthquake but due to the change of the global seismic network some stations were not yet available in 2010 or other were not available anymore in 2015. Three components broad band data have been collected from IRIS and GFZ web-service. Other data from the KOERI seismic network, and the French seismic network have been added. High sampling rate data (i.e. between 20 and 120 Hz) are used by SeisComP3 to detect and pick P waves arrivals. The records are then downsampled to 1 Hz to be compatible with the W-phase algorithm.
The decision matrix to assess the tsunami hazard in the context of a tsunami warning center is currently based on the earthquake location (lat, lon, depth) and the magnitude. On top of the magnitude and the centroid location, it is very important to know the focal mechanism in order to assess the tsunami potential of an earthquake. A strike-slip fault will barely generate a tsunami, compared to an event of the same magnitude with reverse or normal faulting OKAL (1988). A ternary diagram FROHLICH and AP-PERSON (1992) is used to discriminate the type of rupture (strike-slip, normal, thrust). The focal mechanism is represented on a triangle diagram where the three vertices correspond to pure strike-slip (top), reverse (right), normal (left) fault. The moment tensor is expressed in terms of its principal axis which are defined by an eigenvalue, a plunge angle, and an azimuth angle. The location of each points on the diagram depends only on the plunge angle of the 3 principal axes with sin 2 d T þ sin 2 d N þ sin 2 d P ¼ 1 where the T (Tensional) axis correspond to the direction of maximum dilatation, the P (Compressional) axis correspond to the direction of the maximum compression and the N axis is the Null axis. A focal mechanism is considered as a strike-slip fault when d N ! 60 , a normal fault when d P ! 60 and a thrust fault when d T ! 50 .
The ternary diagram representation described above doesn't give a similarity between two solutions of focal mechanism. A method developped by RIVERA and KANAMORI (2014) evaluate the similarity/difference between two focal mechanism regardless the scalar moment is used to compare the focal geometry retrieve by the W-phase inversion and the one given by the GCMT. This method is similar to the Kagan angle KAGAN (1991) which evaluate the minimum rotation angle to match two focal mechanism. However, the Kagan angle can be counterintuitive to interpret and can be bias when the double couple percentage is low. The Focal Mechanism Correlation (hereafter, abbreviated to FMC) takes into account the 2 normalized moment tensors, where ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi M i : M i p is the Frobenius Norm. We define the difference between the 2 moment tensors, Then the FMC is defined by : The FMC varies between 0 (opposite) and 1 (identical). A good correlation between 2 focal geometry when FMC C0.75.

Results
At CENALT, the SeisComP3 software is used for the real-time data acquisition, detection, location, and magnitude estimation based on the P waves phases. In the NEAM region, a PDE solution is given between 2 and 9 min after the origin time of an earthquake. In such a short time, the preliminary magnitude and location can suffer from a lack of accuracy. For the purpose of this study, the PDE 3884 J. Roch et al. Pure Appl. Geophys. Table 1 GCMT solutions and W-phase solutions at t 0 þ 20 and t 0 þ 10 min The SeisComP3 magnitude is shown next to the GCMT magnitude The color code is the same as defined in the Fig. 3 The column S/C/c correspond respectively to the number of stations (S), the number of channels (C) and the azimuthal gap (c) in degrees Vol. 173, (2016) Very fast characterization of focal mechanism parameters 3885 depth has been set to the GCMT depth. The GCMT and SeisComP3 magnitudes are shown in Table 1. The comparison of the two magnitudes can be useful to evaluate the results of the inversion as the W-phase pass-band filter depends on the preliminary magnitude. The W-phase magnitude and focal mechanisms are compared to the GCMT solution for every earthquake. To do so, all events have been processed by the W-phase source inversion within t 0 þ 20 and t 0 þ 10 min of the origin time. The goal of this study is two-fold : first, assessing the W-phase inversion results regarding the GCMT and second, compared the W-phase results within 10 and 20 min of the earthquake origin time. All earthquakes have been processed successfully, meaning a magnitude and a focal mechanism have been computed. However, results obtained with less than 5 stations and/or with azimuthal gap greater than 270 are poorly constrained and should be considered with cautious, if not, to be discard. Table 1 shows the W-phase solutions next to the GCMT solutions.

Magnitude Comparison
In this section, we compare the W-phase magnitude, obtained with 10 and 20 min of signal length with the GCMT magnitude. Figure 6 shows the distribution of magnitude difference between the Wphase magnitude M wp  earthquakes origin time, 86 % of the events have a magnitude within AE0:2 U and 69 % of the events have a magnitude within AE0:1 of the GCMT magnitude. When considering solution within 10 min of the earthquakes origin time, 83 % of the events have a magnitude within AE0:2 U and 66 % of the events have a magnitude within AE0:1 of the GCMT magnitude. Only 3 % difference between magnitudes obtained at t 0 þ 20 and t 0 þ 10 min are observed which comfort us in the use of the W-phase algorithm within 10 min of the earthquake origin time. Earthquakes with a jDM w j ! 0:2 U are located on the NEAM region boundaries where the stations coverage is sparse and so too few stations were used during the inversion (N stat 5). Taking into account 10 min of signals length allow the W-phase algorithm to include stations as far as 21 from the epicentral distance. Atlantic earthquakes occur on the Atlantic ridge which suffer from a lack of seismic stations within 21 which explains the poor quality of the result for these events. 3 out of 6 W-phase magnitudes for the Atlantic region are improved at t 0 þ 20 min (Fig. 6, green markers and Table 1 which include the number of stations taken into account). Some earthquakes in the East region suffer also from a poor seismic stations coverage (Fig. 6, brown markers and Table 1) which make the W-phase solution less accurate. Except for earthquakes occurring in the NEAM region boundaries (Atlantic and East), W-phase magnitudes are quite accurate (within AE0:2 U of the GCMT magnitude) even when the difference between PDE magnitude obtained from SeisComP3 and the GCMT magnitude is significant (Table 1). One events (2013-06-16-Crete) shows a difference greater than 0.2 U although the number of stations/components and azimuthal gap taken into account are good enough (N stations ! 15 and c 146). Moreover, the focal mechanism is also in a good agreement with the GCMT mechanism (see Sect. 4.2). This event shows also a large difference in magnitude between the GCMT magnitude and the preliminary magnitude with jDMj ¼ 0:3 (Table 1). Changing the preliminary magnitude to match the GCMT magnitude didn't improve the final solution. However, the 2013-06-16-Crete earthquake exhibit a half duration higher than expected by the scaling law for such magnitude. This event is a slow earthquake as shown by the SCARDEC inversion VALLÉE et al. (2011) with a source time function length of 18 s and slightly asymetric which is consistent with the half duration obtain by the Wphase algorithm (h d ¼ 6 s at 20 min and h d ¼ 7 s at 10 min).

Focal Mechanism Comparison
In this section, we first compare the strike and dip angles between the W-phase solutions and the GCMT solutions and then we determine the percentage of strike-slip, normal and thrust component represented by a ternary diagram. The strike and dip angles are key parameters in terms of tsunami waves modelisation. The strike and dip angles are linked to the directivity of the tsunami and the wave amplitude. Figure 7 shows the distribution of strike and dip angles difference between the W-phase nodal plans and the GCMT ones. The distribution of strike angle difference has been determined by minimizing the angle between the 2 nodal plans of the W-phase solutions and the 2 nodal plans of the GCMT solutions. The distribution of dip angle difference is associated to the minimum strike difference of every earthquakes.
Earthquakes with a difference in strike and/or dip greater than 20 correspond to solutions obtained with less than 5 stations or with a large azimuthal gap. In average 76 % of the W-phase strikes are within AE10 in 20 and 10 min. The rate increase to 90 and 79 % within AE20 in 20 and 10 min respectively. Normal and thrust faults are the most dangerous case in terms of tsunami generation and hazard. Table 2 shows the percentage of normal and thrust earthquakes (15 out of 29) within 10 and 20 of the GCMT angles for the two time length series.
Comparing the nodal plan angles is useful in the context of tsunami wave propagation but it can be biais when the non double couple component of the rupture process is important. To assess the similarity of the whole source geometry between the W-phase and the GCMT, we compute the Focal Mechanism Correlation (FMC) for the two time series length as defined by RIVERA and KANAMORI (2014). As the moment tensors are normalized, no assumption is made over the magnitude so only the focal geometry is taking into account. A good correlation between focal mechanism is assumed for a FMC parameter greater that 0.75. The Fig. 8 shows the FMC of the set of events within 20 and 10 min of the earthquake origin time.
Focal mechanisms are better constrained, although the difference is small, within 20 min rather than within 10 min of the earthquake origin time. All events with a FMC \0:75 occured on the NEAM region boundary (Atlantic ocean or East region) where the station coverage is poor. Regardless the length of the time series, 7 out of 10 events with a FMC less than 0.75 have a small number of stations.
Another important aspect in the context of tsunami warning center is to evaluate the percentage of strike-slip, normal, thrust component of the rupture. As discuss before, a strike-slip event generate a small tsunami, if not, compare to a normal or thrust event with the same seismic moment. A ternary  . Triangles correspond to the difference of strike/dip obtained with t 0 þ 10 min of signals. Squares correspond to the difference of strike/dip obtained with t 0 þ 20 min of signals. The color of each marker is related to the location of the event described before. The red dash lines represent a difference of AE20 and the red dot lines represent a difference of AE10 . Red dates correspond to normal or thrust faulting Vol. 173, (2016) Very fast characterization of focal mechanism parameters 3889 diagram have been proposed by FROHLICH and APPERSON (1992) to determine the fraction of normal, strike-slip, and thrust fault components of an earthquake. The ternary diagram provides a useful representation in a real-time context to quickly determine the type of rupture involved. We follow the same convention of FROHLICH and APPERSON (1992), a rupture is defined as a strike-slip if dN ! 60 , a normal fault if dP ! 60 and a thrust fault if dT ! 50 . This information is important in real-time to evaluate the tsunami risk. Figure 9 shows the proportions of normal, strike-slip, and thrust faulting for the set of events for the GCMT (a) and the W-phase algorithm at t 0 þ 20 min (b-d), and t 0 þ 10 min (e-g). The GCMT solution have been separated in three groups delimited by the black dash lines. For a sake of clarity, the W-phase solution has been represented for each group. Figure 9b, e represent the percentage of thrust component according to the W-phase solution for the 4 events (green) considered as mostly thrust fault for the GCMT solution. Figure   stations. In Fig. 9d, g, the events (3,9,22,27) are all Atlantic earthquakes which are also poorly constraints (less than 5 stations). It is important to note that no event considered as normal or thrust fault according to the GCMT has been identified as a strike-slip fault by the W-phase algorithm which will be the worst case scenario as it will underestimate the tsunami potential of the earthquake.

Conclusions
The W-phase method has been established as a useful addition to the CENALT software tools and tsunami monitoring procedures. The inversion algorithm has been successfully applied to 29 earthquakes of magnitude as small as M w ¼ 5:8 between 2010 and 2015 with 10 and 20 min of signal length. The goal of (a) (b) (c) (d) (e) (f) (g) Figure 9 Ternary diagrams representing the percentage of strike-slip, normal and thrust components of the 29 earthquakes under study for the GCMT (a) and the W-phase algorithm at t 0 þ 20 min (b-d), and t 0 þ 10 min (e-g). The GCMT solution have been separated in three groups delimited by the black dash lines. Black beachballs around the diagram represent an example of focal mechanism. The number on the top of each dots correspond to the earthquake index of occurrence in our list Vol. 173, (2016) Very fast characterization of focal mechanism parameters 3891 this study was to assess the W-phase algorithm when dealing with only 10 min of signal length after origin time. The regional configuration (2 D 30 ) of the W-phase method alongside with the use of regional Green's functions produces rapid and reliable magnitudes (within AE0:2 U of the GCMT magnitude). Both magnitude and focal mechanism retrieved within 10 and 20 min of the earthquake origin time are in good agreement with the GCMT catalog. The discrepancy between the results obtained with 10 and 20 min are very small and most of the main differences can be explain by the seismic stations distribution, especially near the NEAM region border. For the purpose of this study, we use all stations available and let the W-phase algorithm select the signals. Improved results might be obtained by selecting a primary seismic station network with good SNR and a more homogeneous coverage.