# Towards routine determination of focal mechanisms obtained from first motion P-wave arrivals

Towards routine determination of focal mechanisms obtained from first motion P-wave arrivals Summary The Bulletin of the International Seismological Centre (ISC) contains information on earthquake mechanisms collected from many different sources including national and global agencies, resulting in a satisfactory coverage over a wide magnitude range (M ∼2–9). Nevertheless, there are still a vast number of earthquakes with no reported source mechanisms especially for magnitudes up to 5. This study investigates the possibility of calculating earthquake focal mechanisms in a routine and systematic way based on P-wave first motion polarities. Any available parametric data in the ISC database is being used, as well as auto-picked polarities from waveform data up to teleseismic epicentral distances (90°) for stations that are not reported to the ISC. The determination of the earthquake mechanisms is carried out with a modified version of the HASH algorithm that is compatible with a wide range of epicentral distances and takes into account the ellipsoids defined by the ISC location errors, and the Earth’s structure uncertainties. Initially, benchmark tests for a set of ISC reviewed earthquakes (mb > 4.5) are carried out and the HASH mechanism classification scheme is used to define the mechanism quality. Focal mechanisms of quality A, B and C with an azimuthal gap up to 90° compare well to the benchmark mechanisms. Nevertheless, the majority of the obtained mechanisms fall into class D as a result of limited polarity data from stations in local/regional epicentral distances. Specifically, the computation of the minimum rotation angle between the obtained mechanisms and the benchmarks, reveals that 41  per  cent of the examined earthquakes show rotation angles up to 35°. Finally, the current technique is applied to a small set of earthquakes from the reviewed ISC bulletin where 62 earthquakes, with no previously reported source mechanisms, are successfully obtained. Body waves, Earthquake hazards, Earthquake source observations, Seismicity and tectonics 1 INTRODUCTION Earthquake focal mechanisms offer important information on the characterization of seismic source kinematics, the understanding of tectonic fault systems and the study of aftershock sequences which can be useful for hazard and stress change studies (e.g. Kilb et al. 2000; Toda & Enescu 2011). Even though a preferred representation of a seismic source is given by its moment tensor, routine waveform inversion techniques that offer this information are usually suitable for medium and large magnitude earthquakes that are detected by stations within teleseismic distances. For example, the Global Centroid Moment Tensor (GCMT) technique systematically reports moment tensor solutions for earthquakes with MW ≥ ∼5.5 (e.g. Dziewonski et al. 1981; Ekström et al. 2012). The W-phase technique was initially applied to large magnitude, global earthquakes (MW ≥ 6.5, Kanamori & Rivera 2008; Duputel et al. 2012) but there has been an effort to lower the magnitude threshold to MW ≥ 5.8 (Hayes et al. 2009), while the recently developed SCARDEC method routinely analyses earthquakes with MW ≥ 5.5 (Vallée et al. 2011). Routinely determined regional moment tensor source models are also available. Nevertheless, they cover a limited area such as the European–Mediterranean region only (e.g. Braunmiller et al. 2002; Bernardi et al. 2004; Pondrelli et al. 2011), or the broad area of Japan (National Research Institute for Earth Science and Disaster Prevention – NIED). Other local and regional moment tensor techniques using different algorithms have also been developed (e.g. Koch 1991; Stich et al. 2003; Cesca et al. 2006), but they are either not used routinely, or are applied to particular regions only. Recently Delouis (2014) proposed a method to cover a wide magnitude range (MW > 3.0), which uses strong motion data in short distances, along with broadband records. Focal mechanisms based on P-wave first motion polarities are usually preferred for the study of aftershock sequences and/or microseismicity (e.g. Schneider et al. 1987; Shearer 1998) on a local scale, where robust waveform moment tensor inversions are not always possible (De Natale et al. 1991). Moreover, polarity data can be useful in order to constrain one nodal plane or some of the parameters in moment tensor inversions (Nakanishi & Kanamori 1984; Kawakatsu & Cadena 1991). Nevertheless, the determination of focal mechanisms is not always an easy task due to numerous reasons; notably, the lack of high signal-to-noise ratios in the waveform data of small magnitude earthquakes which yield ambiguous polarity identifications, inverted polarity station channels associated with errors in the instrument response, large station azimuthal gaps, use of poor resolution and/or simplified shallow depth velocity structure in the theoretical calculations. The International Seismological Centre (ISC) is committed to collecting and making available as much of this parametric data, reported from different agencies worldwide, as possible. Nevertheless, the ISC bulletin is not always homogeneous, as the reported solutions are obtained by using different waveform data and techniques. Supporting Information Fig. S1 shows the agencies which systematically report earthquake source mechanisms to the ISC, with the greatest contributors being the HRVD/GCMT and NEIS/NEIC/USGS which provide global coverage. Other notable contributors are the Japan Meteorological Agency (JMA) and the NIED for Japan, the Pacific Northwest Seismic Network (PNSN), the Swiss Seismological Service (ZUR-RMT), and finally, the MedNet Regional Centroid Moment Tensor project (MED-RCMT) for Europe. Nevertheless, a substantial amount of earthquakes in the ISC bulletin have no reported mechanisms. The map in Fig. 1(a) shows the distribution of these earthquakes in the reviewed ISC bulletin for the time period from 1994 to 2013 (ISC database, http://www.isc.ac.uk, last accessed in January 2017, International Seismological Centre 2014). Their distribution is global with the majority of them following the tectonic plate boundaries, however, plenty of intraplate earthquakes also fall into the same category. The graph in Fig. 1(b) which shows their magnitude distribution suggests that the lack of available earthquake source mechanisms is significant up to mb ∼ 5.0 where half of the ISC relocated–reviewed earthquakes do not provide any information on the source mechanism. This gap closes gradually for larger earthquakes, primarily by the contribution of global agencies such as HRVD/GCMT and NEIS/NEIC/USGS. Fig. 2 and Supporting Information Fig. S2 summarize cases of earthquakes discussed in Fig. 1 where potentially focal mechanisms could be obtained where at least 10 stations reporting first motion polarities are available and the station azimuthal gap in the vicinity of the epicentre is kept below ∼90°. These regions are most likely to be South Europe, North India border, Japan, Indonesia, Central America and Chile. Figure 1. View largeDownload slide (a) A global map showing the ISC locations of the relocated earthquakes colour-coded by depth in the reviewed ISC bulletin from 1994 January 1 to 2013 December 31, with mb ranging between 4.0 and 5.5, having no earthquake mechanism information; (b) a histogram showing the frequency–magnitude distribution of the earthquakes presented on the map in ‘a’ in bins of 0.1 of mb. The black bars show all the relocated earthquakes by the ISC and the grey bars show those with no source mechanisms reported. The red bars represent the latter with over 10 stations that report polarities in the ISC bulletin and the green bars those where the distribution of the stations with reported polarities yields an azimuthal gap up to 90°. Note the logarithmic scale on the vertical axis. Figure 1. View largeDownload slide (a) A global map showing the ISC locations of the relocated earthquakes colour-coded by depth in the reviewed ISC bulletin from 1994 January 1 to 2013 December 31, with mb ranging between 4.0 and 5.5, having no earthquake mechanism information; (b) a histogram showing the frequency–magnitude distribution of the earthquakes presented on the map in ‘a’ in bins of 0.1 of mb. The black bars show all the relocated earthquakes by the ISC and the grey bars show those with no source mechanisms reported. The red bars represent the latter with over 10 stations that report polarities in the ISC bulletin and the green bars those where the distribution of the stations with reported polarities yields an azimuthal gap up to 90°. Note the logarithmic scale on the vertical axis. Figure 2. View largeDownload slide (a) A global map showing the locations of the ISC relocated earthquakes of Fig. 1 with depth, with over ten stations reporting first motion polarities; (b) a global map showing the earthquakes in ‘a’ with respect to the number of their associated reported first motion polarities; (c) a global map showing the locations of earthquakes in ‘a’ with respect to the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin; (d) a global map that displays the earthquakes in ‘a’ with respect to the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin. Figure 2. View largeDownload slide (a) A global map showing the locations of the ISC relocated earthquakes of Fig. 1 with depth, with over ten stations reporting first motion polarities; (b) a global map showing the earthquakes in ‘a’ with respect to the number of their associated reported first motion polarities; (c) a global map showing the locations of earthquakes in ‘a’ with respect to the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin; (d) a global map that displays the earthquakes in ‘a’ with respect to the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin. Another common occurrence is summarized in Fig. 3, where the available mechanism solutions provided by different agencies for the same earthquake show a substantial variance in the obtained nodal plane parameters. An elegant and convenient way to quantify differences in source mechanisms is the minimum rotation angle (Kagan 1990, 1991) which will be used in the remainder of this paper. Based on the concept of quaternions, the rotation angle describes the transformation of a double couple mechanism into another arbitrary mechanism through 3-D rotations. Specifically, using normalized quaternions Kagan (1990) found that the minimum rotation angle cannot exceed 120°. Thus, a preferred upper limit of ∼35° for the minimum acceptable rotation angle, that will be used in the remainder of this study, would correspond to ∼30° in combined perturbations of all nodal plane parameters (see e.g. the distance variation, blue line in fig. 1d, in Cesca et al. 2014) which is in good agreement with reported differences in strike, dip and rake among different techniques (see also Section 5). Given that the rotation angle in the distance calculation is divided by 120° (eq. 1 in Cesca et al. 2014), a 35° rotation angle limit corresponds to earthquake mechanisms similarity coefficient of 0.85 (85 per cent) in Cesca et al. (2014). Figure 3. View largeDownload slide Typical examples of earthquakes with high intraevent source mechanism variability in the ISC bulletin: (a) the 2010 May 23, $$m_{b}^{{\rm ISC}}$$ 4.8, Northern Algeria earthquake; (b) the 2012 December 7, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.2, Off east coast of Honshu earthquake; (c) the 2002 November 3, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.9, Central Alaska earthquake. The reported agencies are shown on top of the beach balls, and the values in brackets represent the rotation angle of each mechanism with respect to the HRVD/GCMT earthquake mechanism solution. Figure 3. View largeDownload slide Typical examples of earthquakes with high intraevent source mechanism variability in the ISC bulletin: (a) the 2010 May 23, $$m_{b}^{{\rm ISC}}$$ 4.8, Northern Algeria earthquake; (b) the 2012 December 7, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.2, Off east coast of Honshu earthquake; (c) the 2002 November 3, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.9, Central Alaska earthquake. The reported agencies are shown on top of the beach balls, and the values in brackets represent the rotation angle of each mechanism with respect to the HRVD/GCMT earthquake mechanism solution. Fig. 3(a) shows a typical case for strike-slip mechanisms where different solutions that might look similar at first glance, show opposite sign in rake angle that can yield a high rotation angle (41°). The rest of the examples in Fig. 3 show cases with more obvious differences (rotation angle up to ∼90°). In the case of the 2012 December 7 earthquake, $$M_{\rm W}^{{\rm GCMT}}$$ 7.2, off the east coast of Honshu (Fig. 3b) the intraevent mechanism rotation angle can be as high as 91°, and 57° for the 2002 November 3 $$M_{\rm W}^{{\rm GCMT}}$$ 7.9, Denali, Central Alaska earthquake (Fig. 3c). Later studies suggested the case of a doublet event that began with a MW 7.2, thrust earthquake and followed by an MW 7.1–7.2, normal-faulting earthquake (Harada et al. 2013; Lay et al. 2013) for the 2012 east coast of Honshu earthquake. Similarly, the 2002 Central Alaska earthquake started as a thrust event and then ruptured along the curved strike-slip Denali fault (e.g. Ozacar et al. 2003). Both earthquakes highlight the fact that the use of different data and techniques can yield different earthquake source models, which are used as input for earthquake slip distribution model determination studies (e.g. Delouis et al. 2010). In order to overcome some of the above mentioned limitations in the reviewed ISC bulletin, it is worth investigating the possibility of obtaining earthquake focal mechanisms in a routine and systematic way based on P-wave first motion polarities by using any available parametric and waveform data up to teleseismic epicentral distances (90°). A bi-product of this process could potentially be an additional set of picked polarities and direct P-wave teleseismic arrival times, that might not have been reported to the ISC. Nonetheless, it is expected that to some degree, the same differences will remain in mechanism solutions of complex sources obtained solely from first motion polarities when compared, for example, with results of waveform modelling based techniques. Moreover, the availability of stations in different epicentral distance zones can have an impact on the obtained solutions of large and complex earthquakes, as for example local/regional stations contain near-field information. In the remainder of this paper the methodology for testing this hypothesis will be presented, along with results from benchmark tests, and results from a realistic test case to a set of earthquakes ($$m_{b}^{{\rm ISC}} \ge$$ 4.5) from the reviewed ISC bulletin for the time period covering January–March 2013. 2 DATA AND METHODOLOGY 2.1 Data and resolution The purpose of this study is to examine the possibility of obtaining focal mechanisms using direct P-wave first motion polarities on a global scale. Apart from the polarity information, onset quality (impulsive or emergent), station azimuths and take-off angle calculations are also needed as input. Manually picked P-wave first motion polarities, arrival times and onsets that are reported by different agencies, are obtained from the reviewed ISC bulletin, along with their associated ISC hypocentre solutions and location errors (Bondár & Storchak 2011) for the time period from 2010 to 2013. These are usually short period pickings from stations at local and regional epicentral distances (up to ∼5°). Long period pickings from teleseismic stations (up to 90°) are also available in some cases. Furthermore, automatic picks (see details in Section 2.3) from waveform data from IRIS and EIDA data centres that are downloaded for the stations in the International Registry (IR, http://www.isc.ac.uk/registries/), for the same time period, are also included where possible. Stations beyond 100° are not allowed since the first arrival is the P-wave that is diffracted at the core–mantle boundary (Pdiff). Using first motion polarities to constrain earthquake focal mechanisms on a global scale can be challenging, since different regions are associated with different network configurations, which can have an impact on the azimuthal gap and the minimum epicentral distance where stations are picked. This is a matter of resolution since an adequate sampling of the focal sphere is a key factor in such techniques and is directly associated with the network’s geometry. Fig. 4 shows a synthetic test, where two different cases of source mechanisms recorded by the same stations (43° of azimuthal gap) are examined. In the case of a pure strike–slip earthquake, it is possible to obtain the correct mechanism even by using stations in the 30°–90° region (blue triangles in Fig. 4a). As shown in Fig. 4(b) the area surrounding the null axis is sampled adequately in this case by the current station distribution (the rotation angle with respect to the true mechanism is 4°). On the contrary, in the case of a thrust earthquake, the small azimuthal gap does not guarantee the correct solution (the rotation angle with respect to the true mechanism is 55°) as the polarities from stations in teleseismic distances tend to concentrate at the centre of the focal sphere (Fig. 4c). These polarities are associated with downgoing seismic rays. Nevertheless, the addition of stations up to regional distances (cyan triangles in Fig. 4a) may also include upgoing rays which tend to spread towards the periphery of the sphere, and hence ensure better sampling, and yield a much better constrained solution (Fig. 4d, the rotation angle with respect to the true mechanism is 14°). All the above experiments were carried out using the original HASH algorithm (Hardebeck & Shearer 2002), setting the depth uncertainty to 5 km and the presumed bad polarities to 10 per cent (more details in Sections 2.3 and 3.1), hence the small divergence from the true solution even when adding the regional data. Assuming smaller a priori errors yields better agreement between the true and the obtained nodal plane parameters. In practice, the smaller the assumed uncertainties, the smaller the fuzziness of the station locations on the focal sphere, as they are controlled from the depth uncertainty which affects the take-off angles. However, the value of the HASH method is that the obtained solution is the average of the acceptable nodal planes which fit all the possible locations of the projected stations on the focal sphere. More details will be given in the following Sections. Figure 4. View largeDownload slide (a) A global map showing two shallow (20 km) synthetic earthquakes with the same location (red circle) and their assumed source mechanisms (strike–slip [280°/86°/165°] and thrust [330°/30°/99°]), as well as the distribution of the stations used for the determination of their mechanisms. The blue triangles represent stations in teleseismic distances 30°–90° and the cyan triangles show those up to regional epicentral distances (up to 10°); (b) focal mechanism solution obtained for the strike–slip earthquake with the use of teleseismic stations only; (c) same as in ‘b’ but for the thrust earthquake instead; (d) same as in ‘c’ but the stations in regional distances have also been included in the calculations. The red circles represent compressions and the blue circles show the dilatations. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Figure 4. View largeDownload slide (a) A global map showing two shallow (20 km) synthetic earthquakes with the same location (red circle) and their assumed source mechanisms (strike–slip [280°/86°/165°] and thrust [330°/30°/99°]), as well as the distribution of the stations used for the determination of their mechanisms. The blue triangles represent stations in teleseismic distances 30°–90° and the cyan triangles show those up to regional epicentral distances (up to 10°); (b) focal mechanism solution obtained for the strike–slip earthquake with the use of teleseismic stations only; (c) same as in ‘b’ but for the thrust earthquake instead; (d) same as in ‘c’ but the stations in regional distances have also been included in the calculations. The red circles represent compressions and the blue circles show the dilatations. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. 2.2 Velocity models In conjunction with the above, robust focal mechanism computations strongly depend on the fit of the polarity observations on the focal sphere. The distribution of the polarities is controlled by the station’s azimuth and the theoretical take-off angle. The azimuth is determined by the station and epicentre locations, whereas the take-off angle depends on the choice of the velocity model. Fig. 5 compares the theoretical P-wave take-off angles obtained for three global 1-D velocity models at different source depths. In this paper a concept where the take-off angle is measured from the vertical in a manner where up = 0°, and hence, upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90° will be used throughout. The smooth 1066A model (Gilbert & Dziewonski 1975) shows the largest deviation in comparison with PREM (Dziewonski & Anderson 1981) and ak135 (Kennett et al. 1995) models for the 0°–20° epicentral distance range when the seismic source is placed within the Earth’s crust (∼≤33 km). Results obtained using PREM and ak135 are always very similar. In the ak135 model the crust consists of two uniform layers with discontinuities at 20 km and 35 km, while PREM has the crust–upper mantle boundary at 24.4 km and a transverse isotropic region from 24.4 km down to 220 km. Unrealistic positioning of layer interfaces in the 1-D velocity models, may change direct waves into refracted ones on these intra-crustal interfaces. Hardebeck & Shearer (2002) similarly highlight the importance of the vertical gradient on a 1-D seismic velocity profile. Moreover, in teleseismic distances the rays always leave the source very close to vertical (e.g. shallow depth earthquakes with take-off angles >90°) and are less prone to lateral heterogeneities. On the contrary, direct waves recorded by stations located in the vicinity of the source (up to ∼2°) are associated with near-horizontal upgoing rays that leave the source with take-off angles lower than 90°, and are affected strongly by the Earth’s 3-D structure. Thus, the use of local network data requires detailed P-wave velocity models for the Earth’s crust, where unmodelled lateral heterogeneities may introduce substantial errors in the theoretical calculations. Fig. 5 highlights this as it compares the 1-D ak135 with local/regional structure obtained from CRUST2.0 (Bassin et al. 2000) and CRUST5.1 (Mooney et al. 1998) and a local velocity model obtained from a standard model for South California (SOCAL, Shearer 1997). Figure 5. View largeDownload slide Comparison of theoretical take-off angles versus the epicentral distance (Δ) for three 1-D Earth velocity models (ak135 in blue, PREM in green, 1066a in red, left) and the 1-D ak135 (blue) global velocity model against CRUST2.0 (red), CRUST5.1 (green) and the local velocity model SOCAL (grey) for California (right). The source depth is shown on top of each subplot. Upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90°. All theoretical calculations were carried out with reference to latitude 34.2° and longitude −118.6°. Figure 5. View largeDownload slide Comparison of theoretical take-off angles versus the epicentral distance (Δ) for three 1-D Earth velocity models (ak135 in blue, PREM in green, 1066a in red, left) and the 1-D ak135 (blue) global velocity model against CRUST2.0 (red), CRUST5.1 (green) and the local velocity model SOCAL (grey) for California (right). The source depth is shown on top of each subplot. Upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90°. All theoretical calculations were carried out with reference to latitude 34.2° and longitude −118.6°. 2.3 Method The determination of focal mechanisms from P-wave first motion polarities is carried out by a modified version of the HASH algorithm originally developed by Hardebeck & Shearer (2002). HASH has become very popular (e.g. Bailey et al. 2009, 2010; Katsumata et al. 2010, 2015; Collins et al. 2012; Yang et al. 2012) as it offers an elegant way of taking into account the uncertainties in focal depth and unmodelled Earth’s structure in the theoretical take-off angle calculations. HASH is based on the widely used FPFIT algorithm (Reasenberg & Oppenheimer 1985), which is also a parameter search algorithm. The original HASH code was modified in order to connect and interact with the ISC database, from where the P-wave first motion polarities and onsets that are reported to the ISC by different agencies are obtained, along with their hypocentres and phase arrivals. The option to read from user supplied input files is also available. First motion polarity onsets and arrival times of available waveform data are auto-picked using the FilterPicker, which is an automatic broadband picker developed by Lomax et al. (2012). The picker is also implemented in the modified HASH algorithm, along with the TauP ray-tracing algorithm (Crotwell et al. 1999) in order to carry out the theoretical take-off angle calculations for various velocity models. Alternatively, the user can choose to visualize the waveforms and manually pick the first motion polarities and onsets in an interactive way. As a precautionary measure against false auto-pickings, the auto-picked arrival times are compared against reported arrival times in the ISC bulletin and theoretical arrival times based on ak135 (Kennett et al. 1995) and taking the Earth’s topography into account (Amante & Eakins 2009). Any P-wave polarities that show time residuals greater than a threshold of 4 s times the station weight associated with the ISC location, if applicable, are discarded. Given the earthquake’s location and the Earth’s velocity structure, theoretical take-off angles are calculated and the parameter search algorithm determines a group of acceptable nodal plane solutions which minimize the number of misfit impulsive polarities. The misfit function is represented as the fraction of weighted misfit polarities, and therefore, a perfect fit to the data corresponds to 0 per cent (or 0.0) and an absolute misfit to 100 per cent (or 1.0). In the presence of multiple mechanisms which fit the impulsive polarities equally well, the final set of acceptable solutions is the one that also fits the emergent polarities. In order to account for errors in the polarities, a realistic estimate expressed as the fraction of erroneous polarities is also necessary. More details are given in Section 3.1. Multiple trials are then carried out after perturbing the earthquake’s location with respect to location errors and/or the velocity models. The original HASH code takes into account the vertical error in hypocentre depths and uncertainties in 1D Earth’s structure from a set of velocity models. The algorithm then carries out take-off angle calculations for numerous trials based on random combinations of source depths and velocity models. In addition, the use of the horizontal location errors that affect both the stations’ azimuths and the theoretical take-off angles have been implemented in the calculations. This allows the algorithm to sample the location ellipsoids that are defined by the horizontal and vertical error estimations, and use the ISC locations obtained from the ISC locator (Bondár & Storchak 2011) which systematically estimates the location error. In cases where a hypocentre parameter is fixed and no error estimation is available, the median absolute difference of reported parameters is used to define the error value. As a result, new theoretical take-off angles are computed and new acceptable mechanisms are determined for each trial. Finally, the preferred solution is obtained as the average of the acceptable planes after removing any outliers. In the case where more than one clustered group of acceptable planes is defined, multiple solutions are obtained as the preferred nodal planes. Fig. 6 summarizes the main procedures of the current technique. The uncertainty estimation is expressed as the root-mean-square (RMS) difference of the acceptable nodal planes from the preferred planes and the probability that the real mechanism is close to the preferred mechanism. The probability is expressed as the fraction of acceptable solutions that are within a user-specified threshold (rotation angle) of the preferred solution. In practice, it provides an estimation of how likely it is for the correct solution to be within a chosen clustered group of nodal planes. For further details on the original HASH algorithm the reader is referred to Hardebeck & Shearer (2002, 2003). Figure 6. View largeDownload slide Workflow chart displaying the main processes of the technique presented in this study. Figure 6. View largeDownload slide Workflow chart displaying the main processes of the technique presented in this study. Fig. 7 shows a typical example of a focal mechanism solution obtained using the current technique. The algorithm carried out 50 trials with respect to errors in latitude, longitude and depth, and the errors in Earth’s structure have been taken into account by combining ak135 and PREM global velocity models. The grid angle was set to 5° for the focal mechanism search and the angle for computing the mechanisms’ probability was set to 35°. First motion polarities from the reviewed ISC bulletin and auto-picked first motion polarities for stations with available waveform data up to 90° were used in the computations. Two multiple mechanisms have been determined (Fig. 7a), with the first multiple showing a higher probability (60  per  cent) to be closer to the real mechanism as a result of lower nodal plane uncertainty. Indeed, a comparison with the available mechanisms in the ISC bulletin shows a mean rotation angle of 19°, whereas the second multiplet has a mean rotation angle of 45°, even though it fits the polarity data equally well (Fig. 7b). It would have been impossible to obtain this solution if only the reported polarities had been used as the resulting azimuthal gap would have been as high as 183°, yielding multiple solutions with rotation angles ranging from 39° to 86° when compared with the available earthquake mechanisms in the ISC bulletin. The use of auto-picked first motion polarities decreased the azimuthal gap to 60°, and even though there is a lack of stations in the vicinity of the epicentre (up to ∼2°), a fairly good solution is still obtained (ϕ1 = 347°, δ1 = 39°, λ1 = 121°, ϕ2 = 130°, δ2 = 57°, λ2 = 68°). Figure 7. View largeDownload slide (a) A typical example of the output of the current technique for the 2011 March 7, $$m_{b}^{{\rm ISC}}$$ 6.2, Bougainville–Solomon Islands earthquake. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Black dots represent compressions and white dots show dilatations. P and T principal axes are shown as red and green dots, respectively; (b) a graph showing the distribution of the misfit value (see colour scale) for the nodal plane N.P. 1 shown in ‘a’ with respect to the fault parameters (strike, dip, rake). The red and green stars represent the optimal parameters of the nodal plane N.P. 1 for the first and the second multiplet shown in ‘a’, respectively. Figure 7. View largeDownload slide (a) A typical example of the output of the current technique for the 2011 March 7, $$m_{b}^{{\rm ISC}}$$ 6.2, Bougainville–Solomon Islands earthquake. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Black dots represent compressions and white dots show dilatations. P and T principal axes are shown as red and green dots, respectively; (b) a graph showing the distribution of the misfit value (see colour scale) for the nodal plane N.P. 1 shown in ‘a’ with respect to the fault parameters (strike, dip, rake). The red and green stars represent the optimal parameters of the nodal plane N.P. 1 for the first and the second multiplet shown in ‘a’, respectively. 3 BENCHMARK TESTS In order to assess the robustness of the method, mechanism solutions of 719 earthquakes with MW > 4.5, reported by GCMT for the time period from 2010 to 2013, and relocated by the ISC have been randomly picked and used as a benchmark. Fig. 8 highlights their wide geographic distribution and their representative range in magnitude, depth and double couple mechanism types. Figure 8. View largeDownload slide (a) A global map showing the ISC locations of the earthquakes used in the benchmark tests in this study as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes shown in ‘a’, with respect to their source mechanism type. The colour scale represents the depth. Figure 8. View largeDownload slide (a) A global map showing the ISC locations of the earthquakes used in the benchmark tests in this study as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes shown in ‘a’, with respect to their source mechanism type. The colour scale represents the depth. GCMT provides an average representation of the seismic source, a centroid in space and time, whereas the methodology described in this study is associated with the nucleation point of the rupture. Moreover, both techniques follow different approaches, notably, the current technique carries out a parameter search over the seismic fault geometry based on the onset information obtained from the waveforms, while the GCMT is a waveform inversion technique solving for the centroid location and the moment tensor components. The GCMT moment tensor solutions rarely correspond to pure double couple mechanisms. The moment tensor is diagonalized and decomposed into the isotropic and deviatoric parts. Further decomposition of the deviatoric part yields a best fitting double couple and a compensated linear vector dipole. The GCMT reports the strike, dip, and rake for the best fitting double couple source. The latter will be used in the remainder of this paper and will be compared with the obtained mechanisms of the presented technique. Finally, both techniques use different Earth’s structure for the theoretical calculations and as a result of all the above an absolute agreement among their mechanism solutions cannot be expected. In spite of these fundamental differences the use of the GCMT best fitting double couple mechanisms as a benchmark is expected to provide an independent test. The mechanisms of the selected earthquakes were determined using the modified version of HASH presented in Section 2.3, based on first motion polarities for stations up to 90° as described in Section 2.1. Similarly to the example in Fig. 7, the algorithm carried out 50 trials by using a 5° grid spacing for the focal mechanism search, and the effects of polarity, location and velocity model uncertainties were tested. 3.1 Polarity errors As already mentioned in Sections 1 and 2, polarity based algorithms for focal mechanism calculations are well suited for local/regional seismic networks with a good station distribution (and small intra-station distance). As already shown the HASH algorithm performs quite satisfactorily in these cases (e.g. Hardebeck & Shearer 2002; Kilb & Hardebeck 2006; Katsumata et al. 2010, 2015). Typically HASH obtained mechanism solutions are controlled by the user supplied estimates of the location errors, the velocity model uncertainties expressed as the differences among the various velocity models being used, and an estimate of the polarity errors (different estimates can lead to slightly different results or different multiplet solutions). This latter parameter (badfrac) is supplied by the user as an estimate of the error rate for (impulsive) polarities in the data set and is used to specify the number of allowed polarity misfits for the best-fit solution. Hardebeck & Shearer (2002) studied the possible sources of errors in focal mechanism determinations and showed that even one station with the wrong polarity can be enough to substantially alter the resulting mechanism. Since the majority of the polarity data used in this paper are automatically picked and not manually reviewed, one would expect that there might be some likelihood that a few of the pickings are still erroneous, even after applying the arrival time check described in Section 2. Moreover, the goal of this study is to examine the ability of the algorithm to routinely solve earthquake mechanisms that vary widely in magnitude and for various network geometries (expressed in terms of azimuthal gap and epicentral distance) with the least possible user–algorithm interaction. Under these variable conditions a realistic estimation of the fraction of possible wrong polarities is needed, not necessarily constant, but able to adapt to various earthquake parameters. In the current study the effects of first order source parameters, such as the magnitude and the focal depth with respect to the correctness of the auto-picked polarities are investigated. The magnitude can have a direct impact in terms of signal-to-noise ratios, whereas the depth could have a somewhat secondary effect on the quality of the pick (emergent or impulsive), through different propagation effects for shallow and deeper events. Other factors such as the radiation pattern, the attenuation, the epicentral distance or the station noise level can also have a serious impact on the correctness of automatically picked polarities and phase arrival times. Numerous calculations of focal mechanisms using the benchmark set were carried out, assuming different polarity error estimates (badfrac) ranging from 1 per cent to 30 per cent. There was no clear observation of any effect from earthquakes with substantial differences in depth, possibly due to the fact that the majority of the earthquakes in the benchmark set are shallow or intermediate depth events. The results with respect to the ISC mb are summarized in Fig. 9. In practice mechanisms with dense coverage of the focal sphere and correct picked polarities that yield a low fault plane uncertainty will be more robust and somewhat independent of the user’s estimation. Indeed, earthquakes above ∼5.5 have adequate quality data and show little impact on the variation of the polarity error estimate. On the contrary, obtained mechanisms of lower magnitude earthquakes (mb 4.6–5.3) systematically show a slightly better fit when compared to the benchmark mechanisms in an error estimate (badfrac) zone ranging from 0.1 (10 per cent, mb ∼ 5.3) to 0.2 (20 per cent, mb ∼ 4.6), which can be simulated to some degree by an empirical power law:   $$badfrac = 634.23 \cdot {m_b}^{-5.22},$$ (1)where badfrac is the polarity error estimate as a function of the ISC mb. This power law is not expected to provide a full representation of all the effects that may be associated with the quality of the auto-picked polarities, but rather a bulk estimation of the expected erroneous polarities in the data set, and should only be treated as such. Preliminary tests on earthquakes with mb < 5.5, showed that the use of the error estimates determined from eq. (1) yield 10 per cent more cases showing minimum rotation angles below 35° with respect to the benchmark mechanisms, than simply using a constant error estimate of 0.1 (10 per cent), for example. In the remainder of this study eq. (1) will be used in order to carry out any focal mechanism calculations. Figure 9. View largeDownload slide Distribution of the fraction of polarities that are presumed bad (badfrac) and the $$m_{b}^{{\rm ISC}}$$. The colour scale shows the minimum rotation angle of the obtained mechanisms with respect to GCMT best fitting double couple mechanisms for the benchmark set used in this study. The grey line represents the power law described by eq. (1). Figure 9. View largeDownload slide Distribution of the fraction of polarities that are presumed bad (badfrac) and the $$m_{b}^{{\rm ISC}}$$. The colour scale shows the minimum rotation angle of the obtained mechanisms with respect to GCMT best fitting double couple mechanisms for the benchmark set used in this study. The grey line represents the power law described by eq. (1). 3.2 Earthquake location uncertainties In this section the technique described in Section 2 is applied in order to calculate focal mechanisms for the benchmark set shown in Fig. 8. Specifically, the horizontal location errors and depth errors in the ISC locations are taken into account, whereas the Earth’s structure is defined by the ak135 velocity model alone. In the simple case of a 1-D Earth model, the vertical error in event location has a direct impact on the calculated take-off angle, while horizontal errors can affect the source–station azimuth and the epicentral distance, and hence, the take-off angle. The obtained focal mechanisms are compared against the GCMT benchmark double couple mechanisms and their differences are quantified by calculating the minimum rotation angle. An acceptable match for different mechanisms is considered when the rotation angle is up to 35°, as discussed in Section 1, which is achieved for the 41  per  cent of the obtained mechanisms. Fig. 10 summarizes the results of the current experiment. Moreover, Supporting Information Fig. S3 compares the obtained mechanisms against the GCMT best fitting double couple benchmark mechanisms. It is obvious from the map in Fig. 10(a) that mid-ocean earthquakes are usually the most difficult to constrain unless they are of pure strike-slip type. The ternary plot in Fig. 10(b) shows an overall tendency for pure normal mechanisms to show systematically larger rotation angles, whereas pure strike-slip and/or oblique mechanism earthquakes are easier to constrain for the reasons described in Section 2. Mechanisms of pure thrust earthquakes are also in good agreement with the benchmark mechanisms, possibly due to the fact that these earthquakes are usually large magnitude subduction earthquakes in the benchmark set, and hence, associated with a vast distribution of high signal-to-noise ratio recordings that provide good quality polarity data. Figure 10. View largeDownload slide (a) A global map showing the locations of the earthquakes used in the experiment testing the effect of location errors as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) a ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colour scale represents the minimum rotation angle between the obtained and benchmark mechanisms. Figure 10. View largeDownload slide (a) A global map showing the locations of the earthquakes used in the experiment testing the effect of location errors as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) a ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colour scale represents the minimum rotation angle between the obtained and benchmark mechanisms. The classification scheme proposed by Hardebeck & Shearer (2002), shown in Table 1, is used in order to characterize the quality of the obtained focal mechanisms. Fig. 11 which shows their distribution, suggests that mechanisms of type A, B and C are overall well trusted solutions. Most of the cases where substantial rotation angles are observed, are cases similar to those that are presented in Fig. 3. Specifically, the case of a B mechanism showing a rotation angle of 95° in comparison with GCMT, is the case in Fig. 3(b). When the obtained solution is compared to the various NEIC mechanisms the rotation angle varies from 45° to 85°, but when it is compared to the MOS reported mechanism the rotation angle drops to 22°. The majority of D mechanisms show rotation angles above the acceptable value of 35° (63  per  cent), nevertheless, considering that 87  per  cent of the obtained mechanisms are classified as D, there are still plenty of them (228 earthquakes) which agree well with GCMT. In addition, even though the majority of the D mechanisms shall be discarded, their associated data can still be useful. Figure 11. View largeDownload slide Histograms showing the distribution of the obtained mechanisms in the experiment testing the effect of location errors, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue colour represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean (μ) and standard deviation values (σ) are also shown on every histogram subplot. Figure 11. View largeDownload slide Histograms showing the distribution of the obtained mechanisms in the experiment testing the effect of location errors, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue colour represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean (μ) and standard deviation values (σ) are also shown on every histogram subplot. Table 1. The HASH focal mechanism quality classification scheme that was applied in the current study. Quality  Misfit  RMS  Station  Mechanism      nodal plane  distribution  probability (per cent)      uncertainty (°)  ratio    A  ≤0.15  ≤25  ≥0.5  ≥80  B  ≤0.20  ≤35  ≥0.4  ≥60  C  ≤0.30  ≤45  ≥0.3  ≥50  D  all  the  rest    Quality  Misfit  RMS  Station  Mechanism      nodal plane  distribution  probability (per cent)      uncertainty (°)  ratio    A  ≤0.15  ≤25  ≥0.5  ≥80  B  ≤0.20  ≤35  ≥0.4  ≥60  C  ≤0.30  ≤45  ≥0.3  ≥50  D  all  the  rest    View Large The scatter plots in Supporting Information Fig. S4 provide information on the accuracy of the obtained mechanisms with respect to the earthquake parameters and various mechanism quality criteria. It is usually observed that robust focal mechanism determinations are possible for earthquakes with mb ≥ ∼5.0, and focal depth greater than ∼100 km given that the azimuthal gap is kept well below 120° and the take-off angle gap does not exceed 40°. This can possibly be explained as in general, the greater the magnitude, the earthquakes show a proportional increase of available polarity observations, which yields smaller azimuthal gaps, whereas deep earthquakes (>150 km) are associated with upgoing rays that offer better resolution and decrease the take-off angle gap (<30°). In terms of focal mechanism quality criteria, even though a large scatter in the obtained results is observed, it is apparent that the azimuthal gap and the misfit are probably the best parameters that define good quality mechanisms overall, whereas the nodal plane uncertainty is a good metric for the classification of the obtained mechanisms. Thus, an azimuthal gap greater than 120° or greater than 150° and a misfit value higher than 0.3 shall not be allowed. On the contrary, the station distribution ratio is always too fuzzy and scattered which makes it a less critical parameter. The station distribution ratio was originally introduced by Reasenberg & Oppenheimer (1985) in the FPFIT algorithm as a quantification parameter of the data distribution on the focal sphere. It ranges from 0 to 1, with the larger values showing more robust solutions (>0.5). The above findings are also in agreement with Kilb & Hardebeck (2006). Table 2 shows the selection criteria with respect to each quality mechanism type, where A, B and C mechanisms are only deprecated when showing excessive azimuthal and take-off angle gaps. As far as the D type mechanisms are concerned, a K-means clustering analysis was applied (Hartigan & Wong 1979) in order to distinguish between the well trusted mechanism solutions and those that should get deprecated. Since all the mechanisms that fall into type D (see Table 1) fail to match the criteria for A, B and C types, their quality parameters usually show large dispersion. The clustering analysis showed only a general trend with respect to nodal plane uncertainty, misfit, number of polarity observations and azimuthal gap. Further tweaking was needed in order to define a set of criteria (Table 2) that effectively discards D type mechanisms with large rotation angles. This allows us to select 38  per  cent of the D mechanisms, with 84  per  cent of them showing rotation angles up to 35° in comparison with the benchmark mechanisms. Detailed results for the current experiment are shown in the first two subcolumns of Table 3. Table 2. Selection criteria for HASH classified earthquake mechanisms. Quality  Misfit  RMS  Station  Mechanism  Azimuthal  Takeoff  Number of      nodal plane  distribution  probability (per cent)  gap (°)  angle gap (°)  polarities      uncertainty (°)  ratio          A          ≤150  ≤60    B          ≤150  ≤60    C          ≤120  ≤60    D  ≤0.25  ≤55  ≥0.7  ≥10  ≤90  ≤60  ≥50  D  ≤0.20  ≤55  ≥0.5  ≥10  ≤110  ≤60  ≥300  D  ≤0.20  ≤35  ≥0.6  ≥40  ≤150  ≤60  ≥50  Quality  Misfit  RMS  Station  Mechanism  Azimuthal  Takeoff  Number of      nodal plane  distribution  probability (per cent)  gap (°)  angle gap (°)  polarities      uncertainty (°)  ratio          A          ≤150  ≤60    B          ≤150  ≤60    C          ≤120  ≤60    D  ≤0.25  ≤55  ≥0.7  ≥10  ≤90  ≤60  ≥50  D  ≤0.20  ≤55  ≥0.5  ≥10  ≤110  ≤60  ≥300  D  ≤0.20  ≤35  ≥0.6  ≥40  ≤150  ≤60  ≥50  View Large Table 3. Distribution of the obtained earthquake mechanisms from benchmark tests, according to the quality classification scheme shown in Table 1, for four different scenarios testing the effects of location and velocity model uncertainties. In each case, the two subcolumns show the fraction of earthquakes that fall into each quality class and the fraction of them that show a minimum rotation angle ≤35° when compared to the benchmark earthquake mechanisms. Numbers in parentheses show the same figures after applying the selection criteria presented in Table 2.   Location and velocity model uncertainties    Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth    ak135  crust+ak135  crust+ak135/ak135  crust+ak135/ak135/PREM/1066a    Per cent  Per cent earthquakes  Per cent  Per cent earthquakes  Per cent  per cent earthquakes  Per cent  Per cent earthquakes    earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  Quality    angle ≤35°    angle ≤35°    angle ≤35°    angle ≤35°  A  2.1 (2.1)  93.3 (93.3)  2.0 (2.0)  92.9 (92.9)  2.0 (2.0)  92.9 (92.9)  1.7 (1.7)  100.0 (100.0)  B  3.9 (3.8)  78.6 (81.5)  4.2 (4.0)  76.7 (79.3)  4.0 (3.9)  79.3 (82.1)  3.2 (3.1)  87.0 (90.9)  C  7.0 (5.8)  66.0 (73.8)  6.8 (6.0)  67.4 (69.8)  5.7 (5.0)  68.3 (72.2)  5.3 (5.2)  73.7 (73.0)  D  87.1 (38.2)  36.4 (83.9)  87.1 (43.8)  35.8 (76.5)  88.3 (43.7)  36.4 (78.2)  89.9 (40.7)  36.5 (77.1)    Location and velocity model uncertainties    Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth    ak135  crust+ak135  crust+ak135/ak135  crust+ak135/ak135/PREM/1066a    Per cent  Per cent earthquakes  Per cent  Per cent earthquakes  Per cent  per cent earthquakes  Per cent  Per cent earthquakes    earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  Quality    angle ≤35°    angle ≤35°    angle ≤35°    angle ≤35°  A  2.1 (2.1)  93.3 (93.3)  2.0 (2.0)  92.9 (92.9)  2.0 (2.0)  92.9 (92.9)  1.7 (1.7)  100.0 (100.0)  B  3.9 (3.8)  78.6 (81.5)  4.2 (4.0)  76.7 (79.3)  4.0 (3.9)  79.3 (82.1)  3.2 (3.1)  87.0 (90.9)  C  7.0 (5.8)  66.0 (73.8)  6.8 (6.0)  67.4 (69.8)  5.7 (5.0)  68.3 (72.2)  5.3 (5.2)  73.7 (73.0)  D  87.1 (38.2)  36.4 (83.9)  87.1 (43.8)  35.8 (76.5)  88.3 (43.7)  36.4 (78.2)  89.9 (40.7)  36.5 (77.1)  View Large 3.3 Earthquake location and velocity model uncertainties In the current test the effects of location (same as in Section 3.2) and velocity model uncertainties are combined by calculating theoretical take-off angles for the three global velocity models discussed in Fig. 5. Moreover, detailed crustal structure is superimposed on the 1-D ak135 model, by using the CRUST5.1 velocity model (Mooney et al. 1998) for near regional stations up to an epicentral distance of 5°, and CRUST2.0 (Bassin et al. 2000) for local stations up to a distance of 2°, for each applicable source–station path (crust+ak135). The velocity model uncertainties generally have a stronger effect than the location uncertainties since the depth of the turning point of each seismic ray is controlled by the gradient of the velocity structure. As a consequence, the take-off angles of the seismic rays may vary substantially resulting in multiple focal mechanism solutions. Even though the velocity models shown in Fig. 5 do not show significant differences in the theoretical take-off angles within the 20°–90° epicentral distance range, the shorter distance, shallow depth (<35 km), source–station paths are affected the most. At the near-regional range (up to 5°), the differences between the predicted take-off angles from the global 1-D model ak135 and the local and regional crustal models (Fig. 5) are even more pronounced. As a result, a greater variability in fault parameters is observed among the preferred mechanisms. Specifically, three different experiments were carried out testing the effect of velocity structure uncertainties by calculating theoretical take-off angles when: (i) superimposing detailed crustal structure on ak135 velocity model for up to near regional epicentral distances (5°), (ii) combining the original ak135 velocity model and ak135 after superimposing detailed crustal structure on it, and (iii) all three global velocity models of Fig. 5 and ak135 after superimposing detailed crustal structure on it. The overall findings of these tests are summarized in Table 3. The results are comparable with those of the tests presented in Section 3.2. Specifically, the same 297 out of 719 earthquakes (41  per  cent) showing rotation angles up to 35° are still found, with the better determined mechanisms of types A, B and C being slightly fewer (and D mechanisms rising accordingly) at cases with high velocity structure variability. The findings for the case of superimposing more detailed crustal structure on ak135 are also summarized in Supporting Information Figs S5–S8. 4 APPLICATION TO THE REVIEWED ISC BULLETIN In this section the method discussed in Section 2.3 is applied to all the known earthquakes in the reviewed ISC bulletin with $$m_{b}^{{\rm ISC}} \ge$$ 4.5 for the time period from 2013 January 1 to March 31 (Fig. 12a). Such a search in the reviewed ISC bulletin results in 1477 earthquakes with 4.5 $$\le m_{b}^{{\rm ISC}} \le$$ 6.5 and focal depth ranging from 0.7 km to 646.3 km, where 595 earthquakes have at least one reported source mechanism, whereas 882 of them have no reported mechanisms. First motion polarities from stations up to 90° as described in Section 2.1, were used in the computations. Figure 12. View largeDownload slide (a) A global map showing the locations of the earthquakes in the reviewed ISC bulletin for the time period from 2013 January 1 to 2013 March 31 with $$m_{b}^{{\rm ISC}}$$ ≥ 4.5 to which the presented method was applied. Red circles show earthquakes with at least one reported mechanism, whereas blue circles represent those with no source mechanism information reported to the ISC bulletin; (b) a global map summarizing the results of the current application. The earthquakes in red and blue are respectively those as in ‘a’ for which the obtained mechanism solutions can be trusted. Grey circles show those that have been deprecated and green circles represent earthquakes with less than ten polarities, and hence, no solution was obtained. The size of the circles is proportional to the earthquakes’ $$m_{b}^{{\rm ISC}}$$ magnitudes. Figure 12. View largeDownload slide (a) A global map showing the locations of the earthquakes in the reviewed ISC bulletin for the time period from 2013 January 1 to 2013 March 31 with $$m_{b}^{{\rm ISC}}$$ ≥ 4.5 to which the presented method was applied. Red circles show earthquakes with at least one reported mechanism, whereas blue circles represent those with no source mechanism information reported to the ISC bulletin; (b) a global map summarizing the results of the current application. The earthquakes in red and blue are respectively those as in ‘a’ for which the obtained mechanism solutions can be trusted. Grey circles show those that have been deprecated and green circles represent earthquakes with less than ten polarities, and hence, no solution was obtained. The size of the circles is proportional to the earthquakes’ $$m_{b}^{{\rm ISC}}$$ magnitudes. In order to separate the obtained mechanisms into trusted solutions and deprecated solutions a series of automated quality tests was applied. In the first stage, the quality criteria discussed in Tables 1 and 2 were applied and any mechanisms where the azimuthal gap was greater than 150° were deprecated. In the second stage the multiplets are compared with the obtained mechanisms that have not been deprecated in the first stage. In the third stage the obtained mechanisms (including their multiplets) are compared against source mechanisms that are reported in the ISC bulletin. For those earthquakes where no reported mechanisms exist, the obtained mechanisms are compared against any reported mechanisms in the close vicinity of the hypocentre. Notably, this space is defined as ϕISC ± 2ΔϕISC, λISC ± 2ΔλISC, hISC ± 2ΔhISC, where ϕISC is the ISC latitude, ΔϕISC is the ISC latitude error, λISC is the ISC longitude, ΔλISC is the longitude error, hISC is the ISC depth and ΔhISC is the depth error. The ISC latitude and longitude errors (ΔϕISC, ΔλISC) usually range from 0.01° to 0.25°. This ensures that the maximum search distance from the epicentre is usually up to ∼55 km (e.g. 2ΔϕISC), which compares well with the typical length of some intraplate seismic faults. The depth error is limited to ± 5 km where necessary (error values are determined from the ISC locator). In all cases where mechanism solutions showed rotation angles greater than 35° when compared with the reported mechanisms, the obtained mechanisms were deprecated. Fig. 12(b) shows the geographic distribution of the final set of the obtained mechanisms after deprecating those mechanisms that did not fit the quality criteria, and/or did not compare well with the reported mechanisms in the ISC bulletin (see caption of Fig. 12 for details on the colour-scale). From the total of 1477 earthquakes, mechanism solutions were not obtained for 263 earthquakes as a result of poor data availability (number of polarities < 10). For the remaining 1214 earthquakes, well trusted mechanism solutions were obtained for 181 earthquakes (15  per  cent). Supporting Information Table S1 summarizes the obtained set of focal mechanisms using the presented technique. From the 181 earthquakes, mechanism solutions were obtained for 62 earthquakes where no previously reported mechanisms existed in the ISC bulletin. Moreover, a bi-product of this application is a set of 90108 auto-picked P-wave first motions and phase arrival times (Fig. 13) from stations with epicentral distances ranging from 0.05° to 90.0° that could potentially be used in order to enhance the existing ISC locations. Figure 13. View largeDownload slide A global map showing the locations of the earthquakes shown in Fig. 12 (red stars) and the locations of the stations (blue triangles) for which polarities and P-wave phase arrival times were automatically picked. The grey great circles highlight the associated ray-path coverage. Figure 13. View largeDownload slide A global map showing the locations of the earthquakes shown in Fig. 12 (red stars) and the locations of the stations (blue triangles) for which polarities and P-wave phase arrival times were automatically picked. The grey great circles highlight the associated ray-path coverage. Several reasons can explain the limited final number of seismic events with trusted focal mechanisms. The reported (manually picked) polarities for each event are usually limited to 10–30, depending on the magnitude and the location (mid-oceanic/continental), and many times the associated azimuthal gaps are beyond 90°. The majority of the examined earthquakes are low magnitude events. More specifically, 354 out of the 1477 earthquakes have $$m_{b}^{{\rm ISC}}$$≥ 5.0, and only 78 events have $$m_{b}^{{\rm ISC}}$$ ≥ 5.5. This substantially limits the availability of high signal-to-noise ratio recordings for the smaller magnitude earthquakes, out of which the auto-picker could pick correct first motion polarities. The automatically picked polarities usually range from a few tens to hundreds, depending on the earthquake’s location, magnitude and waveform availability. However, in order to be useful in constraining the mechanism, other factors such as the azimuthal distribution and the minimum epicentral distance have to be considered. For example, hundreds of stations in teleseismic distances having the same azimuth may not contribute as much as 20 local/regional stations with a uniform distribution on the focal sphere (small azimuthal gap). Nonetheless, an increased number of polarities is always desirable. Optimization of the auto-picker could potentially increase the performance of the technique. Further to this, almost all mid-oceanic earthquakes are shallow and lack local/regional stations (see e.g. Fig. 2c) which are essential to constrain their focal mechanisms (Fig. 12b). However, since the current technique is based on the network geometry, it is able to provide well constrained mechanism solutions in cases with good azimuthal coverage of local/regional stations even for small magnitude earthquakes where global techniques cannot be applied. Two such cases will be discussed below. Fig. 14 focuses on the north part of the New Hermides trench where the Australian plate subducts beneath the Pacific plate at a rate of 9.4 cm yr−1 (DeMets et al. 2010). The area is well known for many strong, predominantly thrust earthquakes, with 17 M7.5+ earthquakes having occurred since 1900 according to the reviewed ISC bulletin. During the studied period two substantial earthquakes occurred in the area, notably, the 2013 January 31, $$M_{\rm W}^{\rm {GCMT}}$$ 6.1, and the 2013 February 6, $$M_{\rm W}^{{\rm GCMT}}$$ 7.9, Santa Cruz Islands earthquakes, with the first located off the trench, to the northeast, showing a pure normal mechanism, whereas the second was characterized as a pure thrust earthquake which excited a tsunami with maximum run-up heights of approximately 11 m in the west part of the Nendo Island (Fritz et al. 2014). Both of these events had their mechanisms successfully determined by the presented technique and many of the other obtained mechanisms can be considered as parts of the main events’ aftershock sequences. Specifically, thrust type mechanisms, which are located in the centre of the map of Fig. 14, having focal depths ∼30 km, follow the general orientation of the New Hermides trench and compare well with those reported by other agencies. Moreover, shallow (∼30 km) pure normal type mechanisms which generally follow an N–S orientation are also in agreement with other reported mechanisms and highlight the tectonic complexity of the area, indicating N–S oriented seafloor spreading near the junction of the Solomon and the New Hermides trenches. A second cluster of pure normal mechanisms with a dominant NW–SE orientation is found to the northeast, south of the Santa Cruz Islands. Finally, strike-slip mechanisms are associated with a complex passage from the Solomon trench to the New Hermides trench (DeMets et al. 2010). Figure 14. View largeDownload slide Map (top) and cross-sections (bottom) comparing the obtained earthquake mechanisms of this study for the time period 2013 January–March and available mechanisms in the entire reviewed ISC bulletin for the north part of the New Hermides trench. White lines indicate the main subduction zone boundaries. The white arrow indicates the relative plate motion. Red contour lines on the map and red dashed lines in cross-sections show the slab’s depth as described in the Slab1.0 model (Hayes et al. 2012). The obtained mechanisms in this study with no previously reported solutions in the ISC bulletin are shown in black, and those with previously reported mechanisms are shown in grey. Their ISC locations are presented as colour-coded stars. The two large magnitude earthquakes discussed in Section 4 are plotted as the largest beach balls. The available mechanisms in the ISC bulletin for the studied area are also colour coded with respect to the earthquake depth. Cross-section traces are shown as red lines on the map. The cross-sections are made on the vertical plane and for a width of ±50 km. Figure 14. View largeDownload slide Map (top) and cross-sections (bottom) comparing the obtained earthquake mechanisms of this study for the time period 2013 January–March and available mechanisms in the entire reviewed ISC bulletin for the north part of the New Hermides trench. White lines indicate the main subduction zone boundaries. The white arrow indicates the relative plate motion. Red contour lines on the map and red dashed lines in cross-sections show the slab’s depth as described in the Slab1.0 model (Hayes et al. 2012). The obtained mechanisms in this study with no previously reported solutions in the ISC bulletin are shown in black, and those with previously reported mechanisms are shown in grey. Their ISC locations are presented as colour-coded stars. The two large magnitude earthquakes discussed in Section 4 are plotted as the largest beach balls. The available mechanisms in the ISC bulletin for the studied area are also colour coded with respect to the earthquake depth. Cross-section traces are shown as red lines on the map. The cross-sections are made on the vertical plane and for a width of ±50 km. Fig. 15 compares the obtained mechanism for the intermediate depth (164.8 km), 2013 March 9, $$m_b^{{\rm ISC}}$$ 4.6, Leeward Islands earthquake with other reported earthquake mechanisms and the regional tectonics. The main tectonic feature in the area is the Puerto Rico trench to the north-west which continues to the Lesser Antilles trench to the east, where the N. American plate subducts beneath the Caribbean plate. This transition is marked by the Anegada Passage which separates Puerto Rico from the Virgin Islands and is composed of a complex system of normal faults and strike-slip pull-apart basins (Manaker et al. 2008). The deformation across these faults is responsible for only a very small portion of the observed seismicity in the area and is mainly shallow. The seismicity to the north-west which is associated with the Puerto Rico trench is oblique to some extent, whereas the seismicity to the east is usually characterized by pure thrust mechanisms that follow the orientation of the Lesser Antilles trench (DeMets et al. 2000). No source mechanism information is available in the ISC bulletin regarding the studied earthquake, and thus, the determination of a well constrained focal mechanism is of great importance, especially due to the absence of many available earthquake mechanisms in the vicinity of the earthquake’s hypocentre. Specifically, there is only one earthquake mechanism available in the ISC bulletin (2011 January 21, Leeward Islands, $$m_b^{{\rm ISC}}$$ = 5.0, h = 163.9 km) which is located 3 km away and shows a remarkable agreement with the obtained mechanism solution of this study. Both best double couple mechanisms reveal, normal-faulting rupture with nodal planes having an NE–SW orientation. Intermediate depth (∼140 km), normal-faulting earthquakes have also been observed to the southeast (offshore Martinique Island). In these cases Laigle et al. (2013) suggested that the occurrence of these normal-fault earthquakes appears as a tear in a segmented slab. All in all, intermediate depth and deep earthquakes are special seismic events that can offer insights to deep structure source processes and deserve a more comprehensive study combining both seismological and geophysical data. Figure 15. View largeDownload slide Same as in Fig. 14 but for the northeastern Caribbean area. Figure 15. View largeDownload slide Same as in Fig. 14 but for the northeastern Caribbean area. 5 DISCUSSION This study investigated the possibility of determining focal mechanisms based on first motion P-wave polarities, in a routine and fully automated way on a global scale. The current technique is based on the commonly used HASH algorithm of Hardebeck & Shearer (2002) while some modifications were applied in order to make it more suitable for use with teleseismic distances (up to 90°). Being a parameter search algorithm which takes into account location (latitude, longitude, depth) and velocity model uncertainties, it provides a derivative free and efficient way to solve for the most possible mechanisms despite the nonlinearity of the take-off angle calculation with respect to the earthquake’s location and Earth’s structure. Focal mechanisms obtained from P-wave first motion polarities are usually determined from local network picked polarities that can achieve a good coverage and sampling of the entire focal sphere. Nevertheless, the availability of waveform data from stations in the vicinity of seismic faults is not always guaranteed, either because some networks do not share all their data, or simply because, in many cases, earthquakes occur in areas with poor station coverage. Thus, the use of stations in regional and teleseismic distances is the only option, which somewhat limits the ability to compute fault parameters for small magnitude earthquakes (<5.0). Nevertheless, first motion polarity techniques are not directly associated with the earthquake magnitude, as teleseismic waveform inversion techniques are (e.g. Dziewonski et al. 1981; Vallée et al. 2011; Duputel et al. 2012), but rather the geometry of the network (epicentral distance, azimuthal gap). As a result, techniques based on first motion polarities are preferred to study aftershock sequences and microseismicity (e.g. Schneider et al. 1987; Shearer 1998). In addition to the above, a downside of exclusively using regional and teleseismic stations is the limited sampling of the focal sphere as the seismic rays that reach these stations are downgoing and are concentrated towards the centre of the sphere. Upgoing rays that sample more adequately the focal sphere are associated with stations in short epicentral distances (<1°) for shallow earthquakes (not deeper than 15 km). Nevertheless, these rays are very sensitive to the shallow structure and their incorporation would require the use of detailed local velocity models for the Earth’s crust. On the contrary, the mechanisms of deep earthquakes are better constrained than shallow earthquakes due to the fact that, even in regional distances, the seismic rays are upgoing which yield small take-off angle gaps and offer better resolution as a consequence of more adequate sampling of the focal sphere. Moreover, the degree of lateral heterogeneities in the mantle is smaller than the lithosphere, where these rays travel almost vertically, and thus, are not too sensitive to effects from the 3-D structure. The use of simplified 1-D Earth’s structure can introduce substantial errors in the theoretical take-off angle calculations, and hence, errors in the determination of the mechanisms (Hardebeck & Shearer 2002). For the purposes of this study, 1-D global model such as ak135 which is also used in the ISC locator (Bondár & Storchak 2011), is a good compromise between accuracy and computational speed. However, the accuracy of the theoretical calculations for stations in regional distances could probably be improved by taking into account lateral heterogeneities (e.g. Simmons et al. 2012). Moreover, P and Pn phases arrive together or very close to each other in epicentral distances of ∼2°–5°, and may have opposite polarities and/or different take-off angles. However, Kilb et al. (2000) found that errors introduced due to this condition are not significant, as long as a good azimuthal coverage and good sampling of the focal sphere exists. Other possible sources of errors include wrong polarity pickings that are either a result of reversed polarity stations, or false automatic pickings. The FilterPicker that is used in this method is an automatic broadband picker, well suited for real-time, early-warning purposes (Lomax et al. 2012). After applying the residual criteria (see Section 2 for details) the automatic pickings are correct up to roughly 80 per cent when compared with manual pickings within the epicentral distance range of interest (up to 90°). The FilterPicker has a few tuning parameters which control the time window of the signal, the filtering parameters and several thresholds associated with amplitude measurements. Vassallo et al. (2012) showed that tuning of these parameters, can optimize the performance of the algorithm with respect to a specific application. Benchmark tests showed that in roughly half of the cases (41  per  cent) it is possible to constrain the source mechanism by using the minimum rotation angle (Kagan 1990, 1991) as a mean of quantification. Moreover, in many cases ( ∼ 65  per  cent) it is possible to constrain at least one nodal plane (ϕ, δ) with errors in strike and dip up to 30° (see e.g. Supporting Information Figs S3 and S7). Other studies have also shown that a similar or even higher variability is expected when comparing results obtained from different methods and data (Hjörleifsdóttir & Ekström 2010; Ferreira et al. 2011; Weston et al. 2011; Lentas et al. 2013). Even though it is not always possible to determine the entire mechanism very accurately, constraining at least one plane could prove to be a useful tool for moment tensor inversions (e.g. Nakanishi & Kanamori 1984; Kawakatsu & Cadena 1991). For example, moment tensor inversion techniques using long period surface waves suffer from a trade-off between the seismic moment and the dip angle when solving for shallow-dip, thrust, earthquakes (Kanamori & Given 1981; Tsai et al. 2011). In this case, the excitation of long period surface waves is proportional to both the seismic moment and the sine of the dip angle (Mosin(2δ)), but the two parameters cannot be resolved separately. This trade-off can have serious implications on the earthquake’s magnitude determination, thus, the robust calculation of the dip angle is of great importance. Specifically, Tsai et al. (2011) studied the 2011, MW 9.0, Tohoku earthquake and found that this trade-off can introduce an error of 0.04 in magnitude every 2° of error in dip. Lentas et al. (2014) verified this by studying several thrust earthquakes and found, on average, an error of 0.03 in moment magnitude associated with every 2° of error in dip angle. Moreover, the rake angle is always more difficult to constrain for this particular category of earthquakes, as the strike and rake are mutually correlated in a way that their difference is kept constant as a result of a geometric trade-off (Han et al. 2011, 2013). For instance, constraining the dip-slip components (Mrθ, Mrϕ) in a moment tensor inversion based on strike and dip angles determined from a P-wave polarity based technique could in theory address this issue. However, since the fault plane orientation may vary along the rupture, and given the different hypocentre location concepts that different techniques depend on (rupture nucleation point, centroid location), this procedure could potentially bias the final result for large magnitude earthquakes (M > 6.5) where the point seismic source is no longer a good approximation. As mentioned in Section 3 there are many cases of well constrained mechanisms associated with earthquakes showing excessive azimuthal gaps (>90°). Even though a good azimuthal coverage is fundamental, and always desirable, it does not necessarily guarantee that the sampling of the focal sphere is adequate. FPFIT only uses the station distribution ratio as a measure of the mechanism’s robustness (Reasenberg & Oppenheimer 1985), while the HASH method has a number of metrics to assess the quality of the determined focal mechanisms. Notably, these are the maximum azimuthal and take-off angle gaps, the station distribution ratio, the fault plane uncertainties, the average misfit and the probability of the determined mechanism to be close to the real mechanism, all implemented in a mechanism quality classification scheme (Hardebeck & Shearer 2002). Tests carried out by Kilb & Hardebeck (2006) suggested that the use of fault plane uncertainties as a mechanism quality criteria leads to a substantial increase in the fault mechanism stability. Similarly, in this study the nodal plane uncertainty and the probability were also found to be good metrics for the classification of the obtained mechanisms. Nevertheless, the azimuthal gap and the misfit are probably the best parameters that define accurate mechanism solutions overall (see Supporting Information Fig. S4). The application of the presented technique on a set of earthquakes from the reviewed ISC bulletin revealed the difficulties and limitations of attempting to determine focal mechanisms from first motion polarities on a global scale. Specifically, the current set is dominated by relatively small magnitude earthquakes, notably, 76  per  cent of them have $$m_{b}^{{\rm ISC}}$$ < 5.0. As a consequence, in most of the cases the station distribution was not ideal due to the lack of high signal-to-noise ratio recordings that would allow the auto-picker to pick a robust polarity. Only 20  per  cent of the obtained mechanisms are associated with azimuthal gaps up to 90°, and 12  per  cent of the auto-picked phases are within an epicentral distance of up to 3° which yields limited resolution. Nevertheless, it was still feasible to determine a few earthquakes with no previously reported source mechanisms. Based on this performance 100–200 source mechanisms, with no previously reported solutions, could be added to the ISC bulletin for each reviewed year. This is comparable to other ISC special bulletins, such as the IASPEI Ground Truth reference events list, for which approximately 200 events are added for each reviewed year (e.g. Bondár & McLaughlin 2009). Apart from the obvious benefit of determining earthquake mechanisms where no previously reported source information exists, the current technique can offer an alternative approach to pre-existing source models and give evidence of complex source processes that may be due to the use of different data. Moreover, the teleseismic auto-picked P-phase arrival times could be an additional set to the ISC bulletin and the polarities could also be used in order to check the validity of reported mechanisms. The latter could be a valuable addition to the ISC procedures as, for the time being, the analysis focuses only on the investigation of the reported phase arrival times and hypocentres. Any additional source information provided by the reporting agencies is not confirmed or validated. Another potential benefit from the systematic comparison of auto-picked and reported polarities with the ISC obtained and reported mechanisms could lead to the identification of possible station reversals. This evidence could be forwarded to the agencies that maintain the stations in order to further investigate any possible errors and provide solutions where necessary. Additional data, such as S/P amplitude ratios and S-wave polarities could potentially improve the accuracy of the not well constrained mechanisms, as P-wave amplitudes become larger in the vicinity of the P and T axes and smaller close to the nodal planes, while the S-wave amplitudes become larger close to the nodal planes (e.g. De Natale et al. 1991; Hardebeck & Shearer 2003). Nevertheless, their use is not always straightforward. Notably, S/P amplitude ratios need to be corrected for attenuation and path effects, and they are prone to error due to high levels of noise in the waveform data. Moreover, S-wave polarities can be affected by anisotropy (Shin et al. 1991), whereas the automatic picking of the direct S-wave phase is not always accurate due to intermixing with reflected and/or converted seismic phases. 6 CONCLUSIONS In this study, the possibility of routinely obtaining focal mechanisms from P-wave first motion polarities on a global scale was set under investigation. The main purpose is to enhance the reviewed ISC bulletin especially for cases where there is no reported information on earthquake source mechanisms. It was shown that this is a demanding task for intermediate magnitude earthquakes (mb 4–5) which show the greatest gap in source mechanism information in the bulletin, nevertheless, it may be possible in a few cases. Moreover, the obtained mechanisms could also benefit cases of vast intraevent source mechanism variability that are quite common in the ISC bulletin. Finally, the use of auto-pickings from waveform data (polarities, arrival times) that are not previously reported to the ISC (teleseismic in most cases) could also be made available and used in earthquake locations and/or tomographic studies. Acknowledgements The author wishes to thank the editor Dr Martin Schimmel, and two anonymous reviewers for their comments and suggestions which helped to improve this manuscript. The author also wants to thank Dr Bertrand Delouis for useful remarks and discussion. This research was supported by a National Science Foundation Grant (NSF Grant #1417970). The author gratefully acknowledges the availability of global seismograms from the IRIS and Orfeus data centres, as well as the availability of the HASH source code from the IRIS data centre, the TauP source code from the University of South Carolina and the FilterPicker source code from ALomax scientific. REFERENCES Amante C., Eakins B.W., 2009. ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis, NOAA Technical Memorandum NESDIS NGDC-24. Bailey I.W., Becker T.W., Ben-Zion Y., 2009. Patterns of co-seismic strain computed from southern California focal mechanisms, Geophys. J. Int. , 177, 1015– 1036. Google Scholar CrossRef Search ADS   Bailey I.W., Ben-Zion Y., Becker T.W., Holschneider M., 2010. Quantifying focal mechanism heterogeneity for fault zones in central and southern California, Geophys. J. Int. , 183, 433– 450. Google Scholar CrossRef Search ADS   Bassin C., Laske G., Masters G., 2000. The Current Limits of Resolution for Surface Wave Tomography in North America, EOS, Trans. Am. geophys. Un. , 81, F897. Bernardi F., Braunmiller J., Kradolfer U., Giardini D., 2004. Automatic regional moment tensor inversion in the European-Mediterranean region, Geophys. J. Int. , 157, 703– 716. Google Scholar CrossRef Search ADS   Bondár I., McLaughlin K.L., 2009. A new ground truth data set for seismic studies, Geophys. Res. Lett. , 80, 465– 472. Bondár I., Storchak D., 2011. Improved location procedures at the International Seismological Centre, Geophys. J. Int. , 186, 1220– 1244. Google Scholar CrossRef Search ADS   Braunmiller J., Kradolfer U., Baer M., Giardini D., 2002. Regional moment tensor determination in the European-Mediterranean area - initial results, Tectonophysics , 356, 5– 22. Google Scholar CrossRef Search ADS   Cesca S., Buforn E., Dahm T., 2006. Amplitude spectra moment tensor inversion of shallow earthquakes in Spain, Geophys. J. Int. , 166, 839– 854. Google Scholar CrossRef Search ADS   Cesca S., Sen A.T., Dahm T., 2014. Seismicity monitoring by cluster analysis of moment tensors, Geophys. J. Int. , 196, 1813– 1826. Google Scholar CrossRef Search ADS   Collins J.A., Smith D.K., McGuire J., 2012. Seismicity of the Atlantis Massif detachment fault 30°N at the Mid-Atlantic Ridge, Geochem. Geophys. Geosyst. , 13( 1), doi:10.1029/2012GC004210. Crotwell H.P., Owens T.J., Ritsema J., 1999. The TauP Toolkit: Flexible seismic travel-time and ray-path utilities, Seismol. Res. Lett. , 70, 154– 160. Google Scholar CrossRef Search ADS   Delouis B., 2014. FMNEAR: Determination of focal mechanism and first estimate of rupture directivity using near-source records and a linear distribution of point sources, Bull. seism. Soc. Am. , 104( 3), 1479– 1500. Google Scholar CrossRef Search ADS   Delouis B., Nocquet J.-M., Vallée M., 2010. Slip distribution of the February 27, 2010 MW = 8.8 Maule Earthquake, central Chile, from static and high-rate GPS, InSAR, and broadband teleseismic data, Geophys. Res. Lett. , 37, L17305, doi:10.1029/2010GL043899. Google Scholar CrossRef Search ADS   DeMets C., Jansma P.E., Mattioli G.S., Dixon T.H., Farina F., Bilham R., Calais E., Mann P., 2000. GPS geodetic constraints on Caribbean-North America plate motion, Geophys. Res. Lett. , 27, 437– 440. Google Scholar CrossRef Search ADS   DeMets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   De Natale G., Ferraro A., Virieux J., 1991. A probability method for local earthquake focal mechanisms, Geophys. Res. Lett. , 18( 4), 613– 616. Google Scholar CrossRef Search ADS   Duputel Z., Rivera L., Kanamori H., Hayes G., 2012. W phase source inversion for moderate to large earthquakes (1990–2010), Geophys. J. Int. , 189, doi: 10.1111/j.1365-246X.2012.05419.x. Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter. , 25( 4), 297– 356. Google Scholar CrossRef Search ADS   Dziewonski A.M., Chou T.A., Woodhouse J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. geophys. Res. , 86( B4), 2825– 2852. Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewonski A.M., 2012. The global CMT Project 2004–2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Int. , 200–201, 1– 9. Google Scholar CrossRef Search ADS   Ferreira A.M.G., Weston J., Funning G.J., 2011. Global compilation of interferometric synthetic aperture radar earthquake source models: 2. Effects of 3-D Earth structure. J. geophys. Res. , 116, B08409, doi:10.1029/2010JB008132. Google Scholar CrossRef Search ADS   Fritz H.M., Papantoniou A., Biukoto L., Gilly A., Wei Y., 2014. The Solomon Islands Tsunami of 6 February 2013 in the Santa Cruz Islands: Field Survey and Modeling, in EGU General Assembly , 27 April–2 May, Vienna, Austria. Gilbert F., Dziewonski A.M., 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanism from seismic spectra, Phil. Trans. R. Soc. A , 278, 187– 209. Google Scholar CrossRef Search ADS   Han S.C., Sauber J., Riva R., 2011. Contribution of satellite gravimetry to understanding seismic source processes of the 2011 Tohoku-Oki earthquake, Geophys. Res. Lett. , 38, L24312, doi:10.1029/2011GL049975. Han S.C., Riva R., Sauber J., Okal E., 2013. Source parameter inversion for recent great earthquakes from a decade-long observation of global gravity fields, J. geophys. Res. , 118, 1240– 1267. Google Scholar CrossRef Search ADS   Harada T., Murotani S., Satake K., 2013. A deep outer-rise reverse-fault earthquake immediately triggered a shallow normal-fault earthquake: the 7 December 2012 off-Sanriku earthquake (MW 7.3), Geophys. Res. Lett. , 40, 4214– 4219. Google Scholar CrossRef Search ADS   Hardebeck J.L., Shearer M.P., 2002. A new method for determining first-motion focal mechanisms, Bull. seism. Soc. Am. , 92( 6), 2264– 2276. Google Scholar CrossRef Search ADS   Hardebeck J.L., Shearer M.P., 2003. Using S/P amplitude ratios to constrain the focal mechanisms of small earthquakes, Bull. seism. Soc. Am. , 93( 6), 2434– 2444. Google Scholar CrossRef Search ADS   Hartigan J., Wong M., 1979. Algorithm AS 136: a K-means clustering algorithm, Appl. Stat. , 28( 1), 100– 108. Google Scholar CrossRef Search ADS   Hayes G.P., Rivera L., Kanamori H., 2009. Source inversion of the W-phase: Real-time implementation and extension to low magnitudes, Seismol. Res. Lett. , 80( 5), 817– 822. Google Scholar CrossRef Search ADS   Hayes G.P., Wald D.J., Johnson L., 2012. Slab1.0: A three-dimensional model of global subduction zone geometries, J. geophys. Res. , 117, B01302, doi:10.1029/2011JB008524. Hjörleifsdóttir V., Ekström G., 2010. Effects of three-dimensional Earth structure on CMT earthquake parameters, Phys. Earth planet. Inter. , 179, 178– 190. Google Scholar CrossRef Search ADS   International Seismological Centre, 2014. On-line Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent. , Thatcham, United Kingdom. Kagan Y.Y., 1990. Random stress earthquake statistics: spatial dependence, Geophys. J. Int. , 102, 573– 583. Google Scholar CrossRef Search ADS   Kagan Y.Y., 1991. 3-D rotation of double-couple earthquake sources, Geophys. J. Int. , 106, 709– 716. Google Scholar CrossRef Search ADS   Kanamori H., Given J.W., 1981. Use of long-period surface waves for rapid determination of earthquake-source parameters, Phys. Earth planet. Inter. , 27, 8– 31. Google Scholar CrossRef Search ADS   Kanamori H., Rivera L., 2008. Source inversion of W phase: speeding up seismic tsunami warning, Geophys. J. Int. , 175, 222– 238. Google Scholar CrossRef Search ADS   Katsumata K., Kosuga M., Katao H. & the Japanese University Group of the Joint Sismic Observations at, NKTZ, 2010. Focal mechanisms and stress field in the Atotsugawa fault area, central Honshu, Japan, Earth Planets Space , 62, 367– 380. Google Scholar CrossRef Search ADS   Katsumata K., Kosuga M., Katao H., Yamada T., Kato A. & the Research Group for the Joint Seismic Observations at the Nobi area, 2015. Focal mechanisms and stress field in the Nobi fault area, central Japan, Earth Planets Space , 67( 99), doi:10.1186/s40623-015-0275-2. Kawakatsu H., Cadena G.P., 1991. Focal mechanisms of the March 6, 1987 Ecuador earthquakes -CMT inversion with a first motion constraint, J. Phys. Earth. , 39, 589– 597. Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., Buland R., 1995. Constraints on seismic velocities in the Earth from traveltimes, Geophys. J. Int. , 122, 108– 124. Google Scholar CrossRef Search ADS   Kilb D., Hardebeck J.L., 2006. Fault parameter constraints using relocated earthquakes: A validation of first-motion focal-mechanism data, Bull. seism. Soc. Am. , 96( 3), 1140– 1158. Google Scholar CrossRef Search ADS   Kilb D., Gomberg J., Bodin P., 2000. Triggering of earthquake aftershocks by dynamic stresses, Nature , 408, 570– 574. Google Scholar CrossRef Search ADS PubMed  Koch K., 1991. Moment tensor inversion of local earthquake data: I. Investigation of the method and its numerical stability with model calculations, Geophys. J. Int. , 106, 305– 319. Google Scholar CrossRef Search ADS   Laigle M.et al.  , 2013. Seismic structure and activity of the north-central Lesser Antilles subduction zone from an integrated approach: Similarities with the Tohoku forearc, Tectonophysics , 603, 1– 20. Google Scholar CrossRef Search ADS   Lay T., Duputel Z., Ye L., Kanamori H., 2013. The December 7, 2012 Japan Trench intraplate doublet (MW 7.2, 7.1) and interactions between near-trench intraplate thrust and normal faulting, Phys. Earth planet. Inter. , 220, 73– 78. Google Scholar CrossRef Search ADS   Lentas K., Ferreira A.M.G., Vallée M., 2013. Assessment of SCARDEC source parameters of global large (MW ≥ 7.5) subduction earthquakes, Geophys. J. Inter. , 195, 1989– 2004. Google Scholar CrossRef Search ADS   Lentas K., Ferreira A.M.G., Clévédé E., Roch J., 2014. Source models of great earthquakes from ultra low-frequency normal mode data, Phys. Earth planet. Inter. , 233, 41– 67. Google Scholar CrossRef Search ADS   Lomax A., Satriano C., Vassallo M., 2012. Automatic picker developments and optimization: FilterPicker–a robust, broadband picker for real-time seismic monitoring and earthquake early warning, Seismol. Res. Lett. , 83( 3), 531– 540. Google Scholar CrossRef Search ADS   Manaker D.M.et al.  , 2008. Interseismic Plate coupling and strain partitioning in the Northeastern Caribbean, Geophys. J. Int. , 174, 889– 903. Google Scholar CrossRef Search ADS   Mooney W.D., Laske G., Masters G., 1998. CRUST 5.1: A global crustal model at 5°× 5°, J. geophys. Res. , 103, 727– 747. Google Scholar CrossRef Search ADS   Nakanishi I., Kanamori H., 1984. Source mechanisms of twenty-six, large shallow earthquakes (Ms ≥ 6.5) during 1980 from P-wave first motion and long-period Rayleigh wave data, Bull. seism. Soc. Am. , 74( 3), 805– 818. Ozacar A.A., Beck S.L., Christensen D.H., 2003. Source process of the 3 November 2002 Denali fault earthquake (central Alaska) from teleseismic observations, Geophys. Res. Lett. , 30( 12), doi:10.1029/2003GL017272. Pondrelli S., Salimbeni S., Morelli A., Ekström G., Postpischl L., Vannucci G., Boschi E., 2011. European-Mediterranean regional centroid moment tensor catalog: solutions for 2005–2008, Phys. Earth planet. Inter. , 185( 3–4), 74– 81. Google Scholar CrossRef Search ADS   Reasenberg P.A., Oppenheimer D., 1985. FPFIT, FPPLOT and FPPAGE: FORTRAN computer programs for calculating and displaying earthquake fault-plane solutions, US Geol. Surv. Open-File Rept., 85– 739, 109p. Schneider J.F., Pennington W.D., Meyer R.P., 1987. Microseismicity and focal mechanisms of the intermediate-depth Bucaramanga Nest, Colombia, J. gephys. Res. , 92( B13), 13 913–13 926. Shearer M.P., 1997. Improving local earthquake locations using the L1 norm and waveform cross-correlation: application to the Whittier Narrows, California, aftershock sequence, J. geophys. Res. , 102, 8269– 8283. Google Scholar CrossRef Search ADS   Shearer M.P., 1998. Evidence from a cluster of small earthquakes for a fault at 18 km depth beneath Oak Ridge, southern California, Bull. seism. Soc. Am. , 88( 6), 1327– 1336. Shin X.R., Schneider J.F., Meyer R.P., 1991. Polarities of P and S waves, and shear wave splitting observed from the Bucaramanga Nest, Colombia, J. geophys. Res. , 96( B7), 12 069–12 082. 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( B10), doi:10.1029/2012JB009525. Stich D., Ammon C.J., Morales J., 2003. Moment tensor solutions for small and moderate earthquakes in the Ibero-Maghreb region, J. geophys. Res. , 108( B3), 2148, doi:10.1029/2002JB002057. Google Scholar CrossRef Search ADS   Toda S., Enescu B., 2011. Rate/state Coulomb stress transfer model for the CSEP Japan seismicity forecast, Earth Planets Space , 63, 171– 185. Google Scholar CrossRef Search ADS   Tsai V.C., Hayes G.P., Duputel Z., 2011. Constraints on the long-period moment-dip tradeoff for the Tohoku earthquake, Geophys. Res. Lett. , 38, L00G17, doi:10.1029/2011GL049129. Google Scholar CrossRef Search ADS   Vallée M., Charléty J., Ferreira A.M.G., Delouis B., Vergoz J., 2011. SCARDEC: a new technique for the rapid determination of seismic moment magnitude, focal mechanism and source time functions for large earthquakes using body-wave deconvolution, Geophys. J. Int. , 184, 338– 358. Google Scholar CrossRef Search ADS   Vassallo M., Satriano C., Lomax A., 2012. Automatic picker developments and optimization: A strategy for improving the performances of automatic phase pickers, Seismol. Res. Lett. , 83( 3), 541– 554. Google Scholar CrossRef Search ADS   Weston J., Ferreira A.M.G., Funning G.J., 2011. Global compilation of interferometric synthetic aperture radar earthquake source models: 1. Comparisons with seismic catalogs, J. geophys. Res. , 116, B08408, doi:10.1029/2010JB008131. Google Scholar CrossRef Search ADS   Yang W., Hauksson E., Shearer P.M., 2012. Computing a large refined catalog of focal mechanisms for southern California (1981–2010): temporal stability of the style of faulting, Bull. seism. Soc. Am. , 102( 3), 1179– 1194. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Table S1. Obtained focal mechanisms for the earthquakes in the ISC reviewed bulletin, from 1/1/2013 to 31/1/2013, with $$m_b^{\rm ISC} >$$ 4.5, expressed in strike ϕ, dip δ, rake λ and P, T, B principal axes azimuth Θ and plunge Φ. The number of polarities (Nb of pol.) represents the total number including reported and auto-picked first motion polarities. In the cases where the number of available source mechanism in the ISC bulletin (Nb of mech.) is set to –, the entry in the table represents a focal mechanism solution for an earthquake with no previously reported mechanisms. Repeating origin times indicate multiple focal mechanism solutions for the same earthquake. Figure S1. Histogram showing the seismological agencies that contribute earth-quake mechanisms to the ISC Bulletin up to present. A detailed description of the agency codes can be found at http://www.isc.ac.uk/iscbulletin/agencies/. Note the logarithmic scale on the vertical axis. Figure S2. Histograms showing the distribution of ISC relocated earthquakes with over ten stations reporting first motion polarities shown in Fig. 1 in the main manuscript, with depth (a), the number of their associated reported first motion polarities (b), the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin (c), and the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin (d). Figure S3. Obtained focal mechanisms (red) versus GCMT benchmark mechanisms (black) for the experiment testing the effect of location errors. The letter on top-left represents the quality code of the mechanism solution and the value on top-right is the minimum rotation angle between the obtained and benchmark mechanisms. The mechanisms are sorted from lowest (best fit) to largest rotation angle (worst fit) from bottom to top. Figure S4. Scatter plots comparing possible tradeoffs between earthquake source parameters and focal mechanism quality metrics for the experiment testing the effect of location errors. The colourscale represents the minimum rotation angle between the obtained and GCMT benchmark mechanisms. Figure S5. (a) A global map showing the locations of the earthquakes (GCMT beachballs) used in the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The size of the beachballs is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colourscale represents the minimum rotation angle between the obtained and benchmark mechanisms. Figure S6. Histograms showing the distribution of the obtained mechanisms in the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue color represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean μ and standard deviation values σ are also shown on every histogram subplot. Figure S7. Obtained focal mechanisms (red) versus GCMT benchmark mechanisms (black) for the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The letter on top-left represents the quality code of the mechanism solution and the value on top-right is the minimum rotation angle between the obtained and benchmark mechanisms. The mechanisms are sorted from lowest (best fit) to largest rotation angle (worst fit) from bottom to top. Figure S8. Scatter plots comparing possible tradeoffs between earthquake source parameters and focal mechanism quality metrics for the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The colourscale represents the minimum rotation angle between the obtained and GCMT benchmark mechanisms. 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. © Crown copyright 2017. This article contains public sector information licensed under the Open Government Licence v3.0 (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Towards routine determination of focal mechanisms obtained from first motion P-wave arrivals

, Volume 212 (3) – Mar 1, 2018
22 pages

/lp/ou_press/towards-routine-determination-of-focal-mechanisms-obtained-from-first-LjOjnYOTlK
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx503
Publisher site
See Article on Publisher Site

### Abstract

Summary The Bulletin of the International Seismological Centre (ISC) contains information on earthquake mechanisms collected from many different sources including national and global agencies, resulting in a satisfactory coverage over a wide magnitude range (M ∼2–9). Nevertheless, there are still a vast number of earthquakes with no reported source mechanisms especially for magnitudes up to 5. This study investigates the possibility of calculating earthquake focal mechanisms in a routine and systematic way based on P-wave first motion polarities. Any available parametric data in the ISC database is being used, as well as auto-picked polarities from waveform data up to teleseismic epicentral distances (90°) for stations that are not reported to the ISC. The determination of the earthquake mechanisms is carried out with a modified version of the HASH algorithm that is compatible with a wide range of epicentral distances and takes into account the ellipsoids defined by the ISC location errors, and the Earth’s structure uncertainties. Initially, benchmark tests for a set of ISC reviewed earthquakes (mb > 4.5) are carried out and the HASH mechanism classification scheme is used to define the mechanism quality. Focal mechanisms of quality A, B and C with an azimuthal gap up to 90° compare well to the benchmark mechanisms. Nevertheless, the majority of the obtained mechanisms fall into class D as a result of limited polarity data from stations in local/regional epicentral distances. Specifically, the computation of the minimum rotation angle between the obtained mechanisms and the benchmarks, reveals that 41  per  cent of the examined earthquakes show rotation angles up to 35°. Finally, the current technique is applied to a small set of earthquakes from the reviewed ISC bulletin where 62 earthquakes, with no previously reported source mechanisms, are successfully obtained. Body waves, Earthquake hazards, Earthquake source observations, Seismicity and tectonics 1 INTRODUCTION Earthquake focal mechanisms offer important information on the characterization of seismic source kinematics, the understanding of tectonic fault systems and the study of aftershock sequences which can be useful for hazard and stress change studies (e.g. Kilb et al. 2000; Toda & Enescu 2011). Even though a preferred representation of a seismic source is given by its moment tensor, routine waveform inversion techniques that offer this information are usually suitable for medium and large magnitude earthquakes that are detected by stations within teleseismic distances. For example, the Global Centroid Moment Tensor (GCMT) technique systematically reports moment tensor solutions for earthquakes with MW ≥ ∼5.5 (e.g. Dziewonski et al. 1981; Ekström et al. 2012). The W-phase technique was initially applied to large magnitude, global earthquakes (MW ≥ 6.5, Kanamori & Rivera 2008; Duputel et al. 2012) but there has been an effort to lower the magnitude threshold to MW ≥ 5.8 (Hayes et al. 2009), while the recently developed SCARDEC method routinely analyses earthquakes with MW ≥ 5.5 (Vallée et al. 2011). Routinely determined regional moment tensor source models are also available. Nevertheless, they cover a limited area such as the European–Mediterranean region only (e.g. Braunmiller et al. 2002; Bernardi et al. 2004; Pondrelli et al. 2011), or the broad area of Japan (National Research Institute for Earth Science and Disaster Prevention – NIED). Other local and regional moment tensor techniques using different algorithms have also been developed (e.g. Koch 1991; Stich et al. 2003; Cesca et al. 2006), but they are either not used routinely, or are applied to particular regions only. Recently Delouis (2014) proposed a method to cover a wide magnitude range (MW > 3.0), which uses strong motion data in short distances, along with broadband records. Focal mechanisms based on P-wave first motion polarities are usually preferred for the study of aftershock sequences and/or microseismicity (e.g. Schneider et al. 1987; Shearer 1998) on a local scale, where robust waveform moment tensor inversions are not always possible (De Natale et al. 1991). Moreover, polarity data can be useful in order to constrain one nodal plane or some of the parameters in moment tensor inversions (Nakanishi & Kanamori 1984; Kawakatsu & Cadena 1991). Nevertheless, the determination of focal mechanisms is not always an easy task due to numerous reasons; notably, the lack of high signal-to-noise ratios in the waveform data of small magnitude earthquakes which yield ambiguous polarity identifications, inverted polarity station channels associated with errors in the instrument response, large station azimuthal gaps, use of poor resolution and/or simplified shallow depth velocity structure in the theoretical calculations. The International Seismological Centre (ISC) is committed to collecting and making available as much of this parametric data, reported from different agencies worldwide, as possible. Nevertheless, the ISC bulletin is not always homogeneous, as the reported solutions are obtained by using different waveform data and techniques. Supporting Information Fig. S1 shows the agencies which systematically report earthquake source mechanisms to the ISC, with the greatest contributors being the HRVD/GCMT and NEIS/NEIC/USGS which provide global coverage. Other notable contributors are the Japan Meteorological Agency (JMA) and the NIED for Japan, the Pacific Northwest Seismic Network (PNSN), the Swiss Seismological Service (ZUR-RMT), and finally, the MedNet Regional Centroid Moment Tensor project (MED-RCMT) for Europe. Nevertheless, a substantial amount of earthquakes in the ISC bulletin have no reported mechanisms. The map in Fig. 1(a) shows the distribution of these earthquakes in the reviewed ISC bulletin for the time period from 1994 to 2013 (ISC database, http://www.isc.ac.uk, last accessed in January 2017, International Seismological Centre 2014). Their distribution is global with the majority of them following the tectonic plate boundaries, however, plenty of intraplate earthquakes also fall into the same category. The graph in Fig. 1(b) which shows their magnitude distribution suggests that the lack of available earthquake source mechanisms is significant up to mb ∼ 5.0 where half of the ISC relocated–reviewed earthquakes do not provide any information on the source mechanism. This gap closes gradually for larger earthquakes, primarily by the contribution of global agencies such as HRVD/GCMT and NEIS/NEIC/USGS. Fig. 2 and Supporting Information Fig. S2 summarize cases of earthquakes discussed in Fig. 1 where potentially focal mechanisms could be obtained where at least 10 stations reporting first motion polarities are available and the station azimuthal gap in the vicinity of the epicentre is kept below ∼90°. These regions are most likely to be South Europe, North India border, Japan, Indonesia, Central America and Chile. Figure 1. View largeDownload slide (a) A global map showing the ISC locations of the relocated earthquakes colour-coded by depth in the reviewed ISC bulletin from 1994 January 1 to 2013 December 31, with mb ranging between 4.0 and 5.5, having no earthquake mechanism information; (b) a histogram showing the frequency–magnitude distribution of the earthquakes presented on the map in ‘a’ in bins of 0.1 of mb. The black bars show all the relocated earthquakes by the ISC and the grey bars show those with no source mechanisms reported. The red bars represent the latter with over 10 stations that report polarities in the ISC bulletin and the green bars those where the distribution of the stations with reported polarities yields an azimuthal gap up to 90°. Note the logarithmic scale on the vertical axis. Figure 1. View largeDownload slide (a) A global map showing the ISC locations of the relocated earthquakes colour-coded by depth in the reviewed ISC bulletin from 1994 January 1 to 2013 December 31, with mb ranging between 4.0 and 5.5, having no earthquake mechanism information; (b) a histogram showing the frequency–magnitude distribution of the earthquakes presented on the map in ‘a’ in bins of 0.1 of mb. The black bars show all the relocated earthquakes by the ISC and the grey bars show those with no source mechanisms reported. The red bars represent the latter with over 10 stations that report polarities in the ISC bulletin and the green bars those where the distribution of the stations with reported polarities yields an azimuthal gap up to 90°. Note the logarithmic scale on the vertical axis. Figure 2. View largeDownload slide (a) A global map showing the locations of the ISC relocated earthquakes of Fig. 1 with depth, with over ten stations reporting first motion polarities; (b) a global map showing the earthquakes in ‘a’ with respect to the number of their associated reported first motion polarities; (c) a global map showing the locations of earthquakes in ‘a’ with respect to the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin; (d) a global map that displays the earthquakes in ‘a’ with respect to the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin. Figure 2. View largeDownload slide (a) A global map showing the locations of the ISC relocated earthquakes of Fig. 1 with depth, with over ten stations reporting first motion polarities; (b) a global map showing the earthquakes in ‘a’ with respect to the number of their associated reported first motion polarities; (c) a global map showing the locations of earthquakes in ‘a’ with respect to the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin; (d) a global map that displays the earthquakes in ‘a’ with respect to the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin. Another common occurrence is summarized in Fig. 3, where the available mechanism solutions provided by different agencies for the same earthquake show a substantial variance in the obtained nodal plane parameters. An elegant and convenient way to quantify differences in source mechanisms is the minimum rotation angle (Kagan 1990, 1991) which will be used in the remainder of this paper. Based on the concept of quaternions, the rotation angle describes the transformation of a double couple mechanism into another arbitrary mechanism through 3-D rotations. Specifically, using normalized quaternions Kagan (1990) found that the minimum rotation angle cannot exceed 120°. Thus, a preferred upper limit of ∼35° for the minimum acceptable rotation angle, that will be used in the remainder of this study, would correspond to ∼30° in combined perturbations of all nodal plane parameters (see e.g. the distance variation, blue line in fig. 1d, in Cesca et al. 2014) which is in good agreement with reported differences in strike, dip and rake among different techniques (see also Section 5). Given that the rotation angle in the distance calculation is divided by 120° (eq. 1 in Cesca et al. 2014), a 35° rotation angle limit corresponds to earthquake mechanisms similarity coefficient of 0.85 (85 per cent) in Cesca et al. (2014). Figure 3. View largeDownload slide Typical examples of earthquakes with high intraevent source mechanism variability in the ISC bulletin: (a) the 2010 May 23, $$m_{b}^{{\rm ISC}}$$ 4.8, Northern Algeria earthquake; (b) the 2012 December 7, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.2, Off east coast of Honshu earthquake; (c) the 2002 November 3, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.9, Central Alaska earthquake. The reported agencies are shown on top of the beach balls, and the values in brackets represent the rotation angle of each mechanism with respect to the HRVD/GCMT earthquake mechanism solution. Figure 3. View largeDownload slide Typical examples of earthquakes with high intraevent source mechanism variability in the ISC bulletin: (a) the 2010 May 23, $$m_{b}^{{\rm ISC}}$$ 4.8, Northern Algeria earthquake; (b) the 2012 December 7, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.2, Off east coast of Honshu earthquake; (c) the 2002 November 3, $$M_{{\rm W}}^{{\rm GCMT}}$$ 7.9, Central Alaska earthquake. The reported agencies are shown on top of the beach balls, and the values in brackets represent the rotation angle of each mechanism with respect to the HRVD/GCMT earthquake mechanism solution. Fig. 3(a) shows a typical case for strike-slip mechanisms where different solutions that might look similar at first glance, show opposite sign in rake angle that can yield a high rotation angle (41°). The rest of the examples in Fig. 3 show cases with more obvious differences (rotation angle up to ∼90°). In the case of the 2012 December 7 earthquake, $$M_{\rm W}^{{\rm GCMT}}$$ 7.2, off the east coast of Honshu (Fig. 3b) the intraevent mechanism rotation angle can be as high as 91°, and 57° for the 2002 November 3 $$M_{\rm W}^{{\rm GCMT}}$$ 7.9, Denali, Central Alaska earthquake (Fig. 3c). Later studies suggested the case of a doublet event that began with a MW 7.2, thrust earthquake and followed by an MW 7.1–7.2, normal-faulting earthquake (Harada et al. 2013; Lay et al. 2013) for the 2012 east coast of Honshu earthquake. Similarly, the 2002 Central Alaska earthquake started as a thrust event and then ruptured along the curved strike-slip Denali fault (e.g. Ozacar et al. 2003). Both earthquakes highlight the fact that the use of different data and techniques can yield different earthquake source models, which are used as input for earthquake slip distribution model determination studies (e.g. Delouis et al. 2010). In order to overcome some of the above mentioned limitations in the reviewed ISC bulletin, it is worth investigating the possibility of obtaining earthquake focal mechanisms in a routine and systematic way based on P-wave first motion polarities by using any available parametric and waveform data up to teleseismic epicentral distances (90°). A bi-product of this process could potentially be an additional set of picked polarities and direct P-wave teleseismic arrival times, that might not have been reported to the ISC. Nonetheless, it is expected that to some degree, the same differences will remain in mechanism solutions of complex sources obtained solely from first motion polarities when compared, for example, with results of waveform modelling based techniques. Moreover, the availability of stations in different epicentral distance zones can have an impact on the obtained solutions of large and complex earthquakes, as for example local/regional stations contain near-field information. In the remainder of this paper the methodology for testing this hypothesis will be presented, along with results from benchmark tests, and results from a realistic test case to a set of earthquakes ($$m_{b}^{{\rm ISC}} \ge$$ 4.5) from the reviewed ISC bulletin for the time period covering January–March 2013. 2 DATA AND METHODOLOGY 2.1 Data and resolution The purpose of this study is to examine the possibility of obtaining focal mechanisms using direct P-wave first motion polarities on a global scale. Apart from the polarity information, onset quality (impulsive or emergent), station azimuths and take-off angle calculations are also needed as input. Manually picked P-wave first motion polarities, arrival times and onsets that are reported by different agencies, are obtained from the reviewed ISC bulletin, along with their associated ISC hypocentre solutions and location errors (Bondár & Storchak 2011) for the time period from 2010 to 2013. These are usually short period pickings from stations at local and regional epicentral distances (up to ∼5°). Long period pickings from teleseismic stations (up to 90°) are also available in some cases. Furthermore, automatic picks (see details in Section 2.3) from waveform data from IRIS and EIDA data centres that are downloaded for the stations in the International Registry (IR, http://www.isc.ac.uk/registries/), for the same time period, are also included where possible. Stations beyond 100° are not allowed since the first arrival is the P-wave that is diffracted at the core–mantle boundary (Pdiff). Using first motion polarities to constrain earthquake focal mechanisms on a global scale can be challenging, since different regions are associated with different network configurations, which can have an impact on the azimuthal gap and the minimum epicentral distance where stations are picked. This is a matter of resolution since an adequate sampling of the focal sphere is a key factor in such techniques and is directly associated with the network’s geometry. Fig. 4 shows a synthetic test, where two different cases of source mechanisms recorded by the same stations (43° of azimuthal gap) are examined. In the case of a pure strike–slip earthquake, it is possible to obtain the correct mechanism even by using stations in the 30°–90° region (blue triangles in Fig. 4a). As shown in Fig. 4(b) the area surrounding the null axis is sampled adequately in this case by the current station distribution (the rotation angle with respect to the true mechanism is 4°). On the contrary, in the case of a thrust earthquake, the small azimuthal gap does not guarantee the correct solution (the rotation angle with respect to the true mechanism is 55°) as the polarities from stations in teleseismic distances tend to concentrate at the centre of the focal sphere (Fig. 4c). These polarities are associated with downgoing seismic rays. Nevertheless, the addition of stations up to regional distances (cyan triangles in Fig. 4a) may also include upgoing rays which tend to spread towards the periphery of the sphere, and hence ensure better sampling, and yield a much better constrained solution (Fig. 4d, the rotation angle with respect to the true mechanism is 14°). All the above experiments were carried out using the original HASH algorithm (Hardebeck & Shearer 2002), setting the depth uncertainty to 5 km and the presumed bad polarities to 10 per cent (more details in Sections 2.3 and 3.1), hence the small divergence from the true solution even when adding the regional data. Assuming smaller a priori errors yields better agreement between the true and the obtained nodal plane parameters. In practice, the smaller the assumed uncertainties, the smaller the fuzziness of the station locations on the focal sphere, as they are controlled from the depth uncertainty which affects the take-off angles. However, the value of the HASH method is that the obtained solution is the average of the acceptable nodal planes which fit all the possible locations of the projected stations on the focal sphere. More details will be given in the following Sections. Figure 4. View largeDownload slide (a) A global map showing two shallow (20 km) synthetic earthquakes with the same location (red circle) and their assumed source mechanisms (strike–slip [280°/86°/165°] and thrust [330°/30°/99°]), as well as the distribution of the stations used for the determination of their mechanisms. The blue triangles represent stations in teleseismic distances 30°–90° and the cyan triangles show those up to regional epicentral distances (up to 10°); (b) focal mechanism solution obtained for the strike–slip earthquake with the use of teleseismic stations only; (c) same as in ‘b’ but for the thrust earthquake instead; (d) same as in ‘c’ but the stations in regional distances have also been included in the calculations. The red circles represent compressions and the blue circles show the dilatations. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Figure 4. View largeDownload slide (a) A global map showing two shallow (20 km) synthetic earthquakes with the same location (red circle) and their assumed source mechanisms (strike–slip [280°/86°/165°] and thrust [330°/30°/99°]), as well as the distribution of the stations used for the determination of their mechanisms. The blue triangles represent stations in teleseismic distances 30°–90° and the cyan triangles show those up to regional epicentral distances (up to 10°); (b) focal mechanism solution obtained for the strike–slip earthquake with the use of teleseismic stations only; (c) same as in ‘b’ but for the thrust earthquake instead; (d) same as in ‘c’ but the stations in regional distances have also been included in the calculations. The red circles represent compressions and the blue circles show the dilatations. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. 2.2 Velocity models In conjunction with the above, robust focal mechanism computations strongly depend on the fit of the polarity observations on the focal sphere. The distribution of the polarities is controlled by the station’s azimuth and the theoretical take-off angle. The azimuth is determined by the station and epicentre locations, whereas the take-off angle depends on the choice of the velocity model. Fig. 5 compares the theoretical P-wave take-off angles obtained for three global 1-D velocity models at different source depths. In this paper a concept where the take-off angle is measured from the vertical in a manner where up = 0°, and hence, upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90° will be used throughout. The smooth 1066A model (Gilbert & Dziewonski 1975) shows the largest deviation in comparison with PREM (Dziewonski & Anderson 1981) and ak135 (Kennett et al. 1995) models for the 0°–20° epicentral distance range when the seismic source is placed within the Earth’s crust (∼≤33 km). Results obtained using PREM and ak135 are always very similar. In the ak135 model the crust consists of two uniform layers with discontinuities at 20 km and 35 km, while PREM has the crust–upper mantle boundary at 24.4 km and a transverse isotropic region from 24.4 km down to 220 km. Unrealistic positioning of layer interfaces in the 1-D velocity models, may change direct waves into refracted ones on these intra-crustal interfaces. Hardebeck & Shearer (2002) similarly highlight the importance of the vertical gradient on a 1-D seismic velocity profile. Moreover, in teleseismic distances the rays always leave the source very close to vertical (e.g. shallow depth earthquakes with take-off angles >90°) and are less prone to lateral heterogeneities. On the contrary, direct waves recorded by stations located in the vicinity of the source (up to ∼2°) are associated with near-horizontal upgoing rays that leave the source with take-off angles lower than 90°, and are affected strongly by the Earth’s 3-D structure. Thus, the use of local network data requires detailed P-wave velocity models for the Earth’s crust, where unmodelled lateral heterogeneities may introduce substantial errors in the theoretical calculations. Fig. 5 highlights this as it compares the 1-D ak135 with local/regional structure obtained from CRUST2.0 (Bassin et al. 2000) and CRUST5.1 (Mooney et al. 1998) and a local velocity model obtained from a standard model for South California (SOCAL, Shearer 1997). Figure 5. View largeDownload slide Comparison of theoretical take-off angles versus the epicentral distance (Δ) for three 1-D Earth velocity models (ak135 in blue, PREM in green, 1066a in red, left) and the 1-D ak135 (blue) global velocity model against CRUST2.0 (red), CRUST5.1 (green) and the local velocity model SOCAL (grey) for California (right). The source depth is shown on top of each subplot. Upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90°. All theoretical calculations were carried out with reference to latitude 34.2° and longitude −118.6°. Figure 5. View largeDownload slide Comparison of theoretical take-off angles versus the epicentral distance (Δ) for three 1-D Earth velocity models (ak135 in blue, PREM in green, 1066a in red, left) and the 1-D ak135 (blue) global velocity model against CRUST2.0 (red), CRUST5.1 (green) and the local velocity model SOCAL (grey) for California (right). The source depth is shown on top of each subplot. Upgoing rays have take-off angles lower than 90°, whereas downgoing rays show take-off angles greater than 90°. All theoretical calculations were carried out with reference to latitude 34.2° and longitude −118.6°. 2.3 Method The determination of focal mechanisms from P-wave first motion polarities is carried out by a modified version of the HASH algorithm originally developed by Hardebeck & Shearer (2002). HASH has become very popular (e.g. Bailey et al. 2009, 2010; Katsumata et al. 2010, 2015; Collins et al. 2012; Yang et al. 2012) as it offers an elegant way of taking into account the uncertainties in focal depth and unmodelled Earth’s structure in the theoretical take-off angle calculations. HASH is based on the widely used FPFIT algorithm (Reasenberg & Oppenheimer 1985), which is also a parameter search algorithm. The original HASH code was modified in order to connect and interact with the ISC database, from where the P-wave first motion polarities and onsets that are reported to the ISC by different agencies are obtained, along with their hypocentres and phase arrivals. The option to read from user supplied input files is also available. First motion polarity onsets and arrival times of available waveform data are auto-picked using the FilterPicker, which is an automatic broadband picker developed by Lomax et al. (2012). The picker is also implemented in the modified HASH algorithm, along with the TauP ray-tracing algorithm (Crotwell et al. 1999) in order to carry out the theoretical take-off angle calculations for various velocity models. Alternatively, the user can choose to visualize the waveforms and manually pick the first motion polarities and onsets in an interactive way. As a precautionary measure against false auto-pickings, the auto-picked arrival times are compared against reported arrival times in the ISC bulletin and theoretical arrival times based on ak135 (Kennett et al. 1995) and taking the Earth’s topography into account (Amante & Eakins 2009). Any P-wave polarities that show time residuals greater than a threshold of 4 s times the station weight associated with the ISC location, if applicable, are discarded. Given the earthquake’s location and the Earth’s velocity structure, theoretical take-off angles are calculated and the parameter search algorithm determines a group of acceptable nodal plane solutions which minimize the number of misfit impulsive polarities. The misfit function is represented as the fraction of weighted misfit polarities, and therefore, a perfect fit to the data corresponds to 0 per cent (or 0.0) and an absolute misfit to 100 per cent (or 1.0). In the presence of multiple mechanisms which fit the impulsive polarities equally well, the final set of acceptable solutions is the one that also fits the emergent polarities. In order to account for errors in the polarities, a realistic estimate expressed as the fraction of erroneous polarities is also necessary. More details are given in Section 3.1. Multiple trials are then carried out after perturbing the earthquake’s location with respect to location errors and/or the velocity models. The original HASH code takes into account the vertical error in hypocentre depths and uncertainties in 1D Earth’s structure from a set of velocity models. The algorithm then carries out take-off angle calculations for numerous trials based on random combinations of source depths and velocity models. In addition, the use of the horizontal location errors that affect both the stations’ azimuths and the theoretical take-off angles have been implemented in the calculations. This allows the algorithm to sample the location ellipsoids that are defined by the horizontal and vertical error estimations, and use the ISC locations obtained from the ISC locator (Bondár & Storchak 2011) which systematically estimates the location error. In cases where a hypocentre parameter is fixed and no error estimation is available, the median absolute difference of reported parameters is used to define the error value. As a result, new theoretical take-off angles are computed and new acceptable mechanisms are determined for each trial. Finally, the preferred solution is obtained as the average of the acceptable planes after removing any outliers. In the case where more than one clustered group of acceptable planes is defined, multiple solutions are obtained as the preferred nodal planes. Fig. 6 summarizes the main procedures of the current technique. The uncertainty estimation is expressed as the root-mean-square (RMS) difference of the acceptable nodal planes from the preferred planes and the probability that the real mechanism is close to the preferred mechanism. The probability is expressed as the fraction of acceptable solutions that are within a user-specified threshold (rotation angle) of the preferred solution. In practice, it provides an estimation of how likely it is for the correct solution to be within a chosen clustered group of nodal planes. For further details on the original HASH algorithm the reader is referred to Hardebeck & Shearer (2002, 2003). Figure 6. View largeDownload slide Workflow chart displaying the main processes of the technique presented in this study. Figure 6. View largeDownload slide Workflow chart displaying the main processes of the technique presented in this study. Fig. 7 shows a typical example of a focal mechanism solution obtained using the current technique. The algorithm carried out 50 trials with respect to errors in latitude, longitude and depth, and the errors in Earth’s structure have been taken into account by combining ak135 and PREM global velocity models. The grid angle was set to 5° for the focal mechanism search and the angle for computing the mechanisms’ probability was set to 35°. First motion polarities from the reviewed ISC bulletin and auto-picked first motion polarities for stations with available waveform data up to 90° were used in the computations. Two multiple mechanisms have been determined (Fig. 7a), with the first multiple showing a higher probability (60  per  cent) to be closer to the real mechanism as a result of lower nodal plane uncertainty. Indeed, a comparison with the available mechanisms in the ISC bulletin shows a mean rotation angle of 19°, whereas the second multiplet has a mean rotation angle of 45°, even though it fits the polarity data equally well (Fig. 7b). It would have been impossible to obtain this solution if only the reported polarities had been used as the resulting azimuthal gap would have been as high as 183°, yielding multiple solutions with rotation angles ranging from 39° to 86° when compared with the available earthquake mechanisms in the ISC bulletin. The use of auto-picked first motion polarities decreased the azimuthal gap to 60°, and even though there is a lack of stations in the vicinity of the epicentre (up to ∼2°), a fairly good solution is still obtained (ϕ1 = 347°, δ1 = 39°, λ1 = 121°, ϕ2 = 130°, δ2 = 57°, λ2 = 68°). Figure 7. View largeDownload slide (a) A typical example of the output of the current technique for the 2011 March 7, $$m_{b}^{{\rm ISC}}$$ 6.2, Bougainville–Solomon Islands earthquake. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Black dots represent compressions and white dots show dilatations. P and T principal axes are shown as red and green dots, respectively; (b) a graph showing the distribution of the misfit value (see colour scale) for the nodal plane N.P. 1 shown in ‘a’ with respect to the fault parameters (strike, dip, rake). The red and green stars represent the optimal parameters of the nodal plane N.P. 1 for the first and the second multiplet shown in ‘a’, respectively. Figure 7. View largeDownload slide (a) A typical example of the output of the current technique for the 2011 March 7, $$m_{b}^{{\rm ISC}}$$ 6.2, Bougainville–Solomon Islands earthquake. The grey nodal planes are the acceptable mechanisms, whereas the thick black nodal planes represent the preferred solutions. Black dots represent compressions and white dots show dilatations. P and T principal axes are shown as red and green dots, respectively; (b) a graph showing the distribution of the misfit value (see colour scale) for the nodal plane N.P. 1 shown in ‘a’ with respect to the fault parameters (strike, dip, rake). The red and green stars represent the optimal parameters of the nodal plane N.P. 1 for the first and the second multiplet shown in ‘a’, respectively. 3 BENCHMARK TESTS In order to assess the robustness of the method, mechanism solutions of 719 earthquakes with MW > 4.5, reported by GCMT for the time period from 2010 to 2013, and relocated by the ISC have been randomly picked and used as a benchmark. Fig. 8 highlights their wide geographic distribution and their representative range in magnitude, depth and double couple mechanism types. Figure 8. View largeDownload slide (a) A global map showing the ISC locations of the earthquakes used in the benchmark tests in this study as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes shown in ‘a’, with respect to their source mechanism type. The colour scale represents the depth. Figure 8. View largeDownload slide (a) A global map showing the ISC locations of the earthquakes used in the benchmark tests in this study as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes shown in ‘a’, with respect to their source mechanism type. The colour scale represents the depth. GCMT provides an average representation of the seismic source, a centroid in space and time, whereas the methodology described in this study is associated with the nucleation point of the rupture. Moreover, both techniques follow different approaches, notably, the current technique carries out a parameter search over the seismic fault geometry based on the onset information obtained from the waveforms, while the GCMT is a waveform inversion technique solving for the centroid location and the moment tensor components. The GCMT moment tensor solutions rarely correspond to pure double couple mechanisms. The moment tensor is diagonalized and decomposed into the isotropic and deviatoric parts. Further decomposition of the deviatoric part yields a best fitting double couple and a compensated linear vector dipole. The GCMT reports the strike, dip, and rake for the best fitting double couple source. The latter will be used in the remainder of this paper and will be compared with the obtained mechanisms of the presented technique. Finally, both techniques use different Earth’s structure for the theoretical calculations and as a result of all the above an absolute agreement among their mechanism solutions cannot be expected. In spite of these fundamental differences the use of the GCMT best fitting double couple mechanisms as a benchmark is expected to provide an independent test. The mechanisms of the selected earthquakes were determined using the modified version of HASH presented in Section 2.3, based on first motion polarities for stations up to 90° as described in Section 2.1. Similarly to the example in Fig. 7, the algorithm carried out 50 trials by using a 5° grid spacing for the focal mechanism search, and the effects of polarity, location and velocity model uncertainties were tested. 3.1 Polarity errors As already mentioned in Sections 1 and 2, polarity based algorithms for focal mechanism calculations are well suited for local/regional seismic networks with a good station distribution (and small intra-station distance). As already shown the HASH algorithm performs quite satisfactorily in these cases (e.g. Hardebeck & Shearer 2002; Kilb & Hardebeck 2006; Katsumata et al. 2010, 2015). Typically HASH obtained mechanism solutions are controlled by the user supplied estimates of the location errors, the velocity model uncertainties expressed as the differences among the various velocity models being used, and an estimate of the polarity errors (different estimates can lead to slightly different results or different multiplet solutions). This latter parameter (badfrac) is supplied by the user as an estimate of the error rate for (impulsive) polarities in the data set and is used to specify the number of allowed polarity misfits for the best-fit solution. Hardebeck & Shearer (2002) studied the possible sources of errors in focal mechanism determinations and showed that even one station with the wrong polarity can be enough to substantially alter the resulting mechanism. Since the majority of the polarity data used in this paper are automatically picked and not manually reviewed, one would expect that there might be some likelihood that a few of the pickings are still erroneous, even after applying the arrival time check described in Section 2. Moreover, the goal of this study is to examine the ability of the algorithm to routinely solve earthquake mechanisms that vary widely in magnitude and for various network geometries (expressed in terms of azimuthal gap and epicentral distance) with the least possible user–algorithm interaction. Under these variable conditions a realistic estimation of the fraction of possible wrong polarities is needed, not necessarily constant, but able to adapt to various earthquake parameters. In the current study the effects of first order source parameters, such as the magnitude and the focal depth with respect to the correctness of the auto-picked polarities are investigated. The magnitude can have a direct impact in terms of signal-to-noise ratios, whereas the depth could have a somewhat secondary effect on the quality of the pick (emergent or impulsive), through different propagation effects for shallow and deeper events. Other factors such as the radiation pattern, the attenuation, the epicentral distance or the station noise level can also have a serious impact on the correctness of automatically picked polarities and phase arrival times. Numerous calculations of focal mechanisms using the benchmark set were carried out, assuming different polarity error estimates (badfrac) ranging from 1 per cent to 30 per cent. There was no clear observation of any effect from earthquakes with substantial differences in depth, possibly due to the fact that the majority of the earthquakes in the benchmark set are shallow or intermediate depth events. The results with respect to the ISC mb are summarized in Fig. 9. In practice mechanisms with dense coverage of the focal sphere and correct picked polarities that yield a low fault plane uncertainty will be more robust and somewhat independent of the user’s estimation. Indeed, earthquakes above ∼5.5 have adequate quality data and show little impact on the variation of the polarity error estimate. On the contrary, obtained mechanisms of lower magnitude earthquakes (mb 4.6–5.3) systematically show a slightly better fit when compared to the benchmark mechanisms in an error estimate (badfrac) zone ranging from 0.1 (10 per cent, mb ∼ 5.3) to 0.2 (20 per cent, mb ∼ 4.6), which can be simulated to some degree by an empirical power law:   $$badfrac = 634.23 \cdot {m_b}^{-5.22},$$ (1)where badfrac is the polarity error estimate as a function of the ISC mb. This power law is not expected to provide a full representation of all the effects that may be associated with the quality of the auto-picked polarities, but rather a bulk estimation of the expected erroneous polarities in the data set, and should only be treated as such. Preliminary tests on earthquakes with mb < 5.5, showed that the use of the error estimates determined from eq. (1) yield 10 per cent more cases showing minimum rotation angles below 35° with respect to the benchmark mechanisms, than simply using a constant error estimate of 0.1 (10 per cent), for example. In the remainder of this study eq. (1) will be used in order to carry out any focal mechanism calculations. Figure 9. View largeDownload slide Distribution of the fraction of polarities that are presumed bad (badfrac) and the $$m_{b}^{{\rm ISC}}$$. The colour scale shows the minimum rotation angle of the obtained mechanisms with respect to GCMT best fitting double couple mechanisms for the benchmark set used in this study. The grey line represents the power law described by eq. (1). Figure 9. View largeDownload slide Distribution of the fraction of polarities that are presumed bad (badfrac) and the $$m_{b}^{{\rm ISC}}$$. The colour scale shows the minimum rotation angle of the obtained mechanisms with respect to GCMT best fitting double couple mechanisms for the benchmark set used in this study. The grey line represents the power law described by eq. (1). 3.2 Earthquake location uncertainties In this section the technique described in Section 2 is applied in order to calculate focal mechanisms for the benchmark set shown in Fig. 8. Specifically, the horizontal location errors and depth errors in the ISC locations are taken into account, whereas the Earth’s structure is defined by the ak135 velocity model alone. In the simple case of a 1-D Earth model, the vertical error in event location has a direct impact on the calculated take-off angle, while horizontal errors can affect the source–station azimuth and the epicentral distance, and hence, the take-off angle. The obtained focal mechanisms are compared against the GCMT benchmark double couple mechanisms and their differences are quantified by calculating the minimum rotation angle. An acceptable match for different mechanisms is considered when the rotation angle is up to 35°, as discussed in Section 1, which is achieved for the 41  per  cent of the obtained mechanisms. Fig. 10 summarizes the results of the current experiment. Moreover, Supporting Information Fig. S3 compares the obtained mechanisms against the GCMT best fitting double couple benchmark mechanisms. It is obvious from the map in Fig. 10(a) that mid-ocean earthquakes are usually the most difficult to constrain unless they are of pure strike-slip type. The ternary plot in Fig. 10(b) shows an overall tendency for pure normal mechanisms to show systematically larger rotation angles, whereas pure strike-slip and/or oblique mechanism earthquakes are easier to constrain for the reasons described in Section 2. Mechanisms of pure thrust earthquakes are also in good agreement with the benchmark mechanisms, possibly due to the fact that these earthquakes are usually large magnitude subduction earthquakes in the benchmark set, and hence, associated with a vast distribution of high signal-to-noise ratio recordings that provide good quality polarity data. Figure 10. View largeDownload slide (a) A global map showing the locations of the earthquakes used in the experiment testing the effect of location errors as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) a ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colour scale represents the minimum rotation angle between the obtained and benchmark mechanisms. Figure 10. View largeDownload slide (a) A global map showing the locations of the earthquakes used in the experiment testing the effect of location errors as GCMT beach balls. The size of the beach balls is proportional to their MW magnitudes; (b) a ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colour scale represents the minimum rotation angle between the obtained and benchmark mechanisms. The classification scheme proposed by Hardebeck & Shearer (2002), shown in Table 1, is used in order to characterize the quality of the obtained focal mechanisms. Fig. 11 which shows their distribution, suggests that mechanisms of type A, B and C are overall well trusted solutions. Most of the cases where substantial rotation angles are observed, are cases similar to those that are presented in Fig. 3. Specifically, the case of a B mechanism showing a rotation angle of 95° in comparison with GCMT, is the case in Fig. 3(b). When the obtained solution is compared to the various NEIC mechanisms the rotation angle varies from 45° to 85°, but when it is compared to the MOS reported mechanism the rotation angle drops to 22°. The majority of D mechanisms show rotation angles above the acceptable value of 35° (63  per  cent), nevertheless, considering that 87  per  cent of the obtained mechanisms are classified as D, there are still plenty of them (228 earthquakes) which agree well with GCMT. In addition, even though the majority of the D mechanisms shall be discarded, their associated data can still be useful. Figure 11. View largeDownload slide Histograms showing the distribution of the obtained mechanisms in the experiment testing the effect of location errors, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue colour represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean (μ) and standard deviation values (σ) are also shown on every histogram subplot. Figure 11. View largeDownload slide Histograms showing the distribution of the obtained mechanisms in the experiment testing the effect of location errors, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue colour represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean (μ) and standard deviation values (σ) are also shown on every histogram subplot. Table 1. The HASH focal mechanism quality classification scheme that was applied in the current study. Quality  Misfit  RMS  Station  Mechanism      nodal plane  distribution  probability (per cent)      uncertainty (°)  ratio    A  ≤0.15  ≤25  ≥0.5  ≥80  B  ≤0.20  ≤35  ≥0.4  ≥60  C  ≤0.30  ≤45  ≥0.3  ≥50  D  all  the  rest    Quality  Misfit  RMS  Station  Mechanism      nodal plane  distribution  probability (per cent)      uncertainty (°)  ratio    A  ≤0.15  ≤25  ≥0.5  ≥80  B  ≤0.20  ≤35  ≥0.4  ≥60  C  ≤0.30  ≤45  ≥0.3  ≥50  D  all  the  rest    View Large The scatter plots in Supporting Information Fig. S4 provide information on the accuracy of the obtained mechanisms with respect to the earthquake parameters and various mechanism quality criteria. It is usually observed that robust focal mechanism determinations are possible for earthquakes with mb ≥ ∼5.0, and focal depth greater than ∼100 km given that the azimuthal gap is kept well below 120° and the take-off angle gap does not exceed 40°. This can possibly be explained as in general, the greater the magnitude, the earthquakes show a proportional increase of available polarity observations, which yields smaller azimuthal gaps, whereas deep earthquakes (>150 km) are associated with upgoing rays that offer better resolution and decrease the take-off angle gap (<30°). In terms of focal mechanism quality criteria, even though a large scatter in the obtained results is observed, it is apparent that the azimuthal gap and the misfit are probably the best parameters that define good quality mechanisms overall, whereas the nodal plane uncertainty is a good metric for the classification of the obtained mechanisms. Thus, an azimuthal gap greater than 120° or greater than 150° and a misfit value higher than 0.3 shall not be allowed. On the contrary, the station distribution ratio is always too fuzzy and scattered which makes it a less critical parameter. The station distribution ratio was originally introduced by Reasenberg & Oppenheimer (1985) in the FPFIT algorithm as a quantification parameter of the data distribution on the focal sphere. It ranges from 0 to 1, with the larger values showing more robust solutions (>0.5). The above findings are also in agreement with Kilb & Hardebeck (2006). Table 2 shows the selection criteria with respect to each quality mechanism type, where A, B and C mechanisms are only deprecated when showing excessive azimuthal and take-off angle gaps. As far as the D type mechanisms are concerned, a K-means clustering analysis was applied (Hartigan & Wong 1979) in order to distinguish between the well trusted mechanism solutions and those that should get deprecated. Since all the mechanisms that fall into type D (see Table 1) fail to match the criteria for A, B and C types, their quality parameters usually show large dispersion. The clustering analysis showed only a general trend with respect to nodal plane uncertainty, misfit, number of polarity observations and azimuthal gap. Further tweaking was needed in order to define a set of criteria (Table 2) that effectively discards D type mechanisms with large rotation angles. This allows us to select 38  per  cent of the D mechanisms, with 84  per  cent of them showing rotation angles up to 35° in comparison with the benchmark mechanisms. Detailed results for the current experiment are shown in the first two subcolumns of Table 3. Table 2. Selection criteria for HASH classified earthquake mechanisms. Quality  Misfit  RMS  Station  Mechanism  Azimuthal  Takeoff  Number of      nodal plane  distribution  probability (per cent)  gap (°)  angle gap (°)  polarities      uncertainty (°)  ratio          A          ≤150  ≤60    B          ≤150  ≤60    C          ≤120  ≤60    D  ≤0.25  ≤55  ≥0.7  ≥10  ≤90  ≤60  ≥50  D  ≤0.20  ≤55  ≥0.5  ≥10  ≤110  ≤60  ≥300  D  ≤0.20  ≤35  ≥0.6  ≥40  ≤150  ≤60  ≥50  Quality  Misfit  RMS  Station  Mechanism  Azimuthal  Takeoff  Number of      nodal plane  distribution  probability (per cent)  gap (°)  angle gap (°)  polarities      uncertainty (°)  ratio          A          ≤150  ≤60    B          ≤150  ≤60    C          ≤120  ≤60    D  ≤0.25  ≤55  ≥0.7  ≥10  ≤90  ≤60  ≥50  D  ≤0.20  ≤55  ≥0.5  ≥10  ≤110  ≤60  ≥300  D  ≤0.20  ≤35  ≥0.6  ≥40  ≤150  ≤60  ≥50  View Large Table 3. Distribution of the obtained earthquake mechanisms from benchmark tests, according to the quality classification scheme shown in Table 1, for four different scenarios testing the effects of location and velocity model uncertainties. In each case, the two subcolumns show the fraction of earthquakes that fall into each quality class and the fraction of them that show a minimum rotation angle ≤35° when compared to the benchmark earthquake mechanisms. Numbers in parentheses show the same figures after applying the selection criteria presented in Table 2.   Location and velocity model uncertainties    Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth    ak135  crust+ak135  crust+ak135/ak135  crust+ak135/ak135/PREM/1066a    Per cent  Per cent earthquakes  Per cent  Per cent earthquakes  Per cent  per cent earthquakes  Per cent  Per cent earthquakes    earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  Quality    angle ≤35°    angle ≤35°    angle ≤35°    angle ≤35°  A  2.1 (2.1)  93.3 (93.3)  2.0 (2.0)  92.9 (92.9)  2.0 (2.0)  92.9 (92.9)  1.7 (1.7)  100.0 (100.0)  B  3.9 (3.8)  78.6 (81.5)  4.2 (4.0)  76.7 (79.3)  4.0 (3.9)  79.3 (82.1)  3.2 (3.1)  87.0 (90.9)  C  7.0 (5.8)  66.0 (73.8)  6.8 (6.0)  67.4 (69.8)  5.7 (5.0)  68.3 (72.2)  5.3 (5.2)  73.7 (73.0)  D  87.1 (38.2)  36.4 (83.9)  87.1 (43.8)  35.8 (76.5)  88.3 (43.7)  36.4 (78.2)  89.9 (40.7)  36.5 (77.1)    Location and velocity model uncertainties    Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth  Lat./Lon./Depth    ak135  crust+ak135  crust+ak135/ak135  crust+ak135/ak135/PREM/1066a    Per cent  Per cent earthquakes  Per cent  Per cent earthquakes  Per cent  per cent earthquakes  Per cent  Per cent earthquakes    earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  earthquakes  with rot.  Quality    angle ≤35°    angle ≤35°    angle ≤35°    angle ≤35°  A  2.1 (2.1)  93.3 (93.3)  2.0 (2.0)  92.9 (92.9)  2.0 (2.0)  92.9 (92.9)  1.7 (1.7)  100.0 (100.0)  B  3.9 (3.8)  78.6 (81.5)  4.2 (4.0)  76.7 (79.3)  4.0 (3.9)  79.3 (82.1)  3.2 (3.1)  87.0 (90.9)  C  7.0 (5.8)  66.0 (73.8)  6.8 (6.0)  67.4 (69.8)  5.7 (5.0)  68.3 (72.2)  5.3 (5.2)  73.7 (73.0)  D  87.1 (38.2)  36.4 (83.9)  87.1 (43.8)  35.8 (76.5)  88.3 (43.7)  36.4 (78.2)  89.9 (40.7)  36.5 (77.1)  View Large 3.3 Earthquake location and velocity model uncertainties In the current test the effects of location (same as in Section 3.2) and velocity model uncertainties are combined by calculating theoretical take-off angles for the three global velocity models discussed in Fig. 5. Moreover, detailed crustal structure is superimposed on the 1-D ak135 model, by using the CRUST5.1 velocity model (Mooney et al. 1998) for near regional stations up to an epicentral distance of 5°, and CRUST2.0 (Bassin et al. 2000) for local stations up to a distance of 2°, for each applicable source–station path (crust+ak135). The velocity model uncertainties generally have a stronger effect than the location uncertainties since the depth of the turning point of each seismic ray is controlled by the gradient of the velocity structure. As a consequence, the take-off angles of the seismic rays may vary substantially resulting in multiple focal mechanism solutions. Even though the velocity models shown in Fig. 5 do not show significant differences in the theoretical take-off angles within the 20°–90° epicentral distance range, the shorter distance, shallow depth (<35 km), source–station paths are affected the most. At the near-regional range (up to 5°), the differences between the predicted take-off angles from the global 1-D model ak135 and the local and regional crustal models (Fig. 5) are even more pronounced. As a result, a greater variability in fault parameters is observed among the preferred mechanisms. Specifically, three different experiments were carried out testing the effect of velocity structure uncertainties by calculating theoretical take-off angles when: (i) superimposing detailed crustal structure on ak135 velocity model for up to near regional epicentral distances (5°), (ii) combining the original ak135 velocity model and ak135 after superimposing detailed crustal structure on it, and (iii) all three global velocity models of Fig. 5 and ak135 after superimposing detailed crustal structure on it. The overall findings of these tests are summarized in Table 3. The results are comparable with those of the tests presented in Section 3.2. Specifically, the same 297 out of 719 earthquakes (41  per  cent) showing rotation angles up to 35° are still found, with the better determined mechanisms of types A, B and C being slightly fewer (and D mechanisms rising accordingly) at cases with high velocity structure variability. The findings for the case of superimposing more detailed crustal structure on ak135 are also summarized in Supporting Information Figs S5–S8. 4 APPLICATION TO THE REVIEWED ISC BULLETIN In this section the method discussed in Section 2.3 is applied to all the known earthquakes in the reviewed ISC bulletin with $$m_{b}^{{\rm ISC}} \ge$$ 4.5 for the time period from 2013 January 1 to March 31 (Fig. 12a). Such a search in the reviewed ISC bulletin results in 1477 earthquakes with 4.5 $$\le m_{b}^{{\rm ISC}} \le$$ 6.5 and focal depth ranging from 0.7 km to 646.3 km, where 595 earthquakes have at least one reported source mechanism, whereas 882 of them have no reported mechanisms. First motion polarities from stations up to 90° as described in Section 2.1, were used in the computations. Figure 12. View largeDownload slide (a) A global map showing the locations of the earthquakes in the reviewed ISC bulletin for the time period from 2013 January 1 to 2013 March 31 with $$m_{b}^{{\rm ISC}}$$ ≥ 4.5 to which the presented method was applied. Red circles show earthquakes with at least one reported mechanism, whereas blue circles represent those with no source mechanism information reported to the ISC bulletin; (b) a global map summarizing the results of the current application. The earthquakes in red and blue are respectively those as in ‘a’ for which the obtained mechanism solutions can be trusted. Grey circles show those that have been deprecated and green circles represent earthquakes with less than ten polarities, and hence, no solution was obtained. The size of the circles is proportional to the earthquakes’ $$m_{b}^{{\rm ISC}}$$ magnitudes. Figure 12. View largeDownload slide (a) A global map showing the locations of the earthquakes in the reviewed ISC bulletin for the time period from 2013 January 1 to 2013 March 31 with $$m_{b}^{{\rm ISC}}$$ ≥ 4.5 to which the presented method was applied. Red circles show earthquakes with at least one reported mechanism, whereas blue circles represent those with no source mechanism information reported to the ISC bulletin; (b) a global map summarizing the results of the current application. The earthquakes in red and blue are respectively those as in ‘a’ for which the obtained mechanism solutions can be trusted. Grey circles show those that have been deprecated and green circles represent earthquakes with less than ten polarities, and hence, no solution was obtained. The size of the circles is proportional to the earthquakes’ $$m_{b}^{{\rm ISC}}$$ magnitudes. In order to separate the obtained mechanisms into trusted solutions and deprecated solutions a series of automated quality tests was applied. In the first stage, the quality criteria discussed in Tables 1 and 2 were applied and any mechanisms where the azimuthal gap was greater than 150° were deprecated. In the second stage the multiplets are compared with the obtained mechanisms that have not been deprecated in the first stage. In the third stage the obtained mechanisms (including their multiplets) are compared against source mechanisms that are reported in the ISC bulletin. For those earthquakes where no reported mechanisms exist, the obtained mechanisms are compared against any reported mechanisms in the close vicinity of the hypocentre. Notably, this space is defined as ϕISC ± 2ΔϕISC, λISC ± 2ΔλISC, hISC ± 2ΔhISC, where ϕISC is the ISC latitude, ΔϕISC is the ISC latitude error, λISC is the ISC longitude, ΔλISC is the longitude error, hISC is the ISC depth and ΔhISC is the depth error. The ISC latitude and longitude errors (ΔϕISC, ΔλISC) usually range from 0.01° to 0.25°. This ensures that the maximum search distance from the epicentre is usually up to ∼55 km (e.g. 2ΔϕISC), which compares well with the typical length of some intraplate seismic faults. The depth error is limited to ± 5 km where necessary (error values are determined from the ISC locator). In all cases where mechanism solutions showed rotation angles greater than 35° when compared with the reported mechanisms, the obtained mechanisms were deprecated. Fig. 12(b) shows the geographic distribution of the final set of the obtained mechanisms after deprecating those mechanisms that did not fit the quality criteria, and/or did not compare well with the reported mechanisms in the ISC bulletin (see caption of Fig. 12 for details on the colour-scale). From the total of 1477 earthquakes, mechanism solutions were not obtained for 263 earthquakes as a result of poor data availability (number of polarities < 10). For the remaining 1214 earthquakes, well trusted mechanism solutions were obtained for 181 earthquakes (15  per  cent). Supporting Information Table S1 summarizes the obtained set of focal mechanisms using the presented technique. From the 181 earthquakes, mechanism solutions were obtained for 62 earthquakes where no previously reported mechanisms existed in the ISC bulletin. Moreover, a bi-product of this application is a set of 90108 auto-picked P-wave first motions and phase arrival times (Fig. 13) from stations with epicentral distances ranging from 0.05° to 90.0° that could potentially be used in order to enhance the existing ISC locations. Figure 13. View largeDownload slide A global map showing the locations of the earthquakes shown in Fig. 12 (red stars) and the locations of the stations (blue triangles) for which polarities and P-wave phase arrival times were automatically picked. The grey great circles highlight the associated ray-path coverage. Figure 13. View largeDownload slide A global map showing the locations of the earthquakes shown in Fig. 12 (red stars) and the locations of the stations (blue triangles) for which polarities and P-wave phase arrival times were automatically picked. The grey great circles highlight the associated ray-path coverage. Several reasons can explain the limited final number of seismic events with trusted focal mechanisms. The reported (manually picked) polarities for each event are usually limited to 10–30, depending on the magnitude and the location (mid-oceanic/continental), and many times the associated azimuthal gaps are beyond 90°. The majority of the examined earthquakes are low magnitude events. More specifically, 354 out of the 1477 earthquakes have $$m_{b}^{{\rm ISC}}$$≥ 5.0, and only 78 events have $$m_{b}^{{\rm ISC}}$$ ≥ 5.5. This substantially limits the availability of high signal-to-noise ratio recordings for the smaller magnitude earthquakes, out of which the auto-picker could pick correct first motion polarities. The automatically picked polarities usually range from a few tens to hundreds, depending on the earthquake’s location, magnitude and waveform availability. However, in order to be useful in constraining the mechanism, other factors such as the azimuthal distribution and the minimum epicentral distance have to be considered. For example, hundreds of stations in teleseismic distances having the same azimuth may not contribute as much as 20 local/regional stations with a uniform distribution on the focal sphere (small azimuthal gap). Nonetheless, an increased number of polarities is always desirable. Optimization of the auto-picker could potentially increase the performance of the technique. Further to this, almost all mid-oceanic earthquakes are shallow and lack local/regional stations (see e.g. Fig. 2c) which are essential to constrain their focal mechanisms (Fig. 12b). However, since the current technique is based on the network geometry, it is able to provide well constrained mechanism solutions in cases with good azimuthal coverage of local/regional stations even for small magnitude earthquakes where global techniques cannot be applied. Two such cases will be discussed below. Fig. 14 focuses on the north part of the New Hermides trench where the Australian plate subducts beneath the Pacific plate at a rate of 9.4 cm yr−1 (DeMets et al. 2010). The area is well known for many strong, predominantly thrust earthquakes, with 17 M7.5+ earthquakes having occurred since 1900 according to the reviewed ISC bulletin. During the studied period two substantial earthquakes occurred in the area, notably, the 2013 January 31, $$M_{\rm W}^{\rm {GCMT}}$$ 6.1, and the 2013 February 6, $$M_{\rm W}^{{\rm GCMT}}$$ 7.9, Santa Cruz Islands earthquakes, with the first located off the trench, to the northeast, showing a pure normal mechanism, whereas the second was characterized as a pure thrust earthquake which excited a tsunami with maximum run-up heights of approximately 11 m in the west part of the Nendo Island (Fritz et al. 2014). Both of these events had their mechanisms successfully determined by the presented technique and many of the other obtained mechanisms can be considered as parts of the main events’ aftershock sequences. Specifically, thrust type mechanisms, which are located in the centre of the map of Fig. 14, having focal depths ∼30 km, follow the general orientation of the New Hermides trench and compare well with those reported by other agencies. Moreover, shallow (∼30 km) pure normal type mechanisms which generally follow an N–S orientation are also in agreement with other reported mechanisms and highlight the tectonic complexity of the area, indicating N–S oriented seafloor spreading near the junction of the Solomon and the New Hermides trenches. A second cluster of pure normal mechanisms with a dominant NW–SE orientation is found to the northeast, south of the Santa Cruz Islands. Finally, strike-slip mechanisms are associated with a complex passage from the Solomon trench to the New Hermides trench (DeMets et al. 2010). Figure 14. View largeDownload slide Map (top) and cross-sections (bottom) comparing the obtained earthquake mechanisms of this study for the time period 2013 January–March and available mechanisms in the entire reviewed ISC bulletin for the north part of the New Hermides trench. White lines indicate the main subduction zone boundaries. The white arrow indicates the relative plate motion. Red contour lines on the map and red dashed lines in cross-sections show the slab’s depth as described in the Slab1.0 model (Hayes et al. 2012). The obtained mechanisms in this study with no previously reported solutions in the ISC bulletin are shown in black, and those with previously reported mechanisms are shown in grey. Their ISC locations are presented as colour-coded stars. The two large magnitude earthquakes discussed in Section 4 are plotted as the largest beach balls. The available mechanisms in the ISC bulletin for the studied area are also colour coded with respect to the earthquake depth. Cross-section traces are shown as red lines on the map. The cross-sections are made on the vertical plane and for a width of ±50 km. Figure 14. View largeDownload slide Map (top) and cross-sections (bottom) comparing the obtained earthquake mechanisms of this study for the time period 2013 January–March and available mechanisms in the entire reviewed ISC bulletin for the north part of the New Hermides trench. White lines indicate the main subduction zone boundaries. The white arrow indicates the relative plate motion. Red contour lines on the map and red dashed lines in cross-sections show the slab’s depth as described in the Slab1.0 model (Hayes et al. 2012). The obtained mechanisms in this study with no previously reported solutions in the ISC bulletin are shown in black, and those with previously reported mechanisms are shown in grey. Their ISC locations are presented as colour-coded stars. The two large magnitude earthquakes discussed in Section 4 are plotted as the largest beach balls. The available mechanisms in the ISC bulletin for the studied area are also colour coded with respect to the earthquake depth. Cross-section traces are shown as red lines on the map. The cross-sections are made on the vertical plane and for a width of ±50 km. Fig. 15 compares the obtained mechanism for the intermediate depth (164.8 km), 2013 March 9, $$m_b^{{\rm ISC}}$$ 4.6, Leeward Islands earthquake with other reported earthquake mechanisms and the regional tectonics. The main tectonic feature in the area is the Puerto Rico trench to the north-west which continues to the Lesser Antilles trench to the east, where the N. American plate subducts beneath the Caribbean plate. This transition is marked by the Anegada Passage which separates Puerto Rico from the Virgin Islands and is composed of a complex system of normal faults and strike-slip pull-apart basins (Manaker et al. 2008). The deformation across these faults is responsible for only a very small portion of the observed seismicity in the area and is mainly shallow. The seismicity to the north-west which is associated with the Puerto Rico trench is oblique to some extent, whereas the seismicity to the east is usually characterized by pure thrust mechanisms that follow the orientation of the Lesser Antilles trench (DeMets et al. 2000). No source mechanism information is available in the ISC bulletin regarding the studied earthquake, and thus, the determination of a well constrained focal mechanism is of great importance, especially due to the absence of many available earthquake mechanisms in the vicinity of the earthquake’s hypocentre. Specifically, there is only one earthquake mechanism available in the ISC bulletin (2011 January 21, Leeward Islands, $$m_b^{{\rm ISC}}$$ = 5.0, h = 163.9 km) which is located 3 km away and shows a remarkable agreement with the obtained mechanism solution of this study. Both best double couple mechanisms reveal, normal-faulting rupture with nodal planes having an NE–SW orientation. Intermediate depth (∼140 km), normal-faulting earthquakes have also been observed to the southeast (offshore Martinique Island). In these cases Laigle et al. (2013) suggested that the occurrence of these normal-fault earthquakes appears as a tear in a segmented slab. All in all, intermediate depth and deep earthquakes are special seismic events that can offer insights to deep structure source processes and deserve a more comprehensive study combining both seismological and geophysical data. Figure 15. View largeDownload slide Same as in Fig. 14 but for the northeastern Caribbean area. Figure 15. View largeDownload slide Same as in Fig. 14 but for the northeastern Caribbean area. 5 DISCUSSION This study investigated the possibility of determining focal mechanisms based on first motion P-wave polarities, in a routine and fully automated way on a global scale. The current technique is based on the commonly used HASH algorithm of Hardebeck & Shearer (2002) while some modifications were applied in order to make it more suitable for use with teleseismic distances (up to 90°). Being a parameter search algorithm which takes into account location (latitude, longitude, depth) and velocity model uncertainties, it provides a derivative free and efficient way to solve for the most possible mechanisms despite the nonlinearity of the take-off angle calculation with respect to the earthquake’s location and Earth’s structure. Focal mechanisms obtained from P-wave first motion polarities are usually determined from local network picked polarities that can achieve a good coverage and sampling of the entire focal sphere. Nevertheless, the availability of waveform data from stations in the vicinity of seismic faults is not always guaranteed, either because some networks do not share all their data, or simply because, in many cases, earthquakes occur in areas with poor station coverage. Thus, the use of stations in regional and teleseismic distances is the only option, which somewhat limits the ability to compute fault parameters for small magnitude earthquakes (<5.0). Nevertheless, first motion polarity techniques are not directly associated with the earthquake magnitude, as teleseismic waveform inversion techniques are (e.g. Dziewonski et al. 1981; Vallée et al. 2011; Duputel et al. 2012), but rather the geometry of the network (epicentral distance, azimuthal gap). As a result, techniques based on first motion polarities are preferred to study aftershock sequences and microseismicity (e.g. Schneider et al. 1987; Shearer 1998). In addition to the above, a downside of exclusively using regional and teleseismic stations is the limited sampling of the focal sphere as the seismic rays that reach these stations are downgoing and are concentrated towards the centre of the sphere. Upgoing rays that sample more adequately the focal sphere are associated with stations in short epicentral distances (<1°) for shallow earthquakes (not deeper than 15 km). Nevertheless, these rays are very sensitive to the shallow structure and their incorporation would require the use of detailed local velocity models for the Earth’s crust. On the contrary, the mechanisms of deep earthquakes are better constrained than shallow earthquakes due to the fact that, even in regional distances, the seismic rays are upgoing which yield small take-off angle gaps and offer better resolution as a consequence of more adequate sampling of the focal sphere. Moreover, the degree of lateral heterogeneities in the mantle is smaller than the lithosphere, where these rays travel almost vertically, and thus, are not too sensitive to effects from the 3-D structure. The use of simplified 1-D Earth’s structure can introduce substantial errors in the theoretical take-off angle calculations, and hence, errors in the determination of the mechanisms (Hardebeck & Shearer 2002). For the purposes of this study, 1-D global model such as ak135 which is also used in the ISC locator (Bondár & Storchak 2011), is a good compromise between accuracy and computational speed. However, the accuracy of the theoretical calculations for stations in regional distances could probably be improved by taking into account lateral heterogeneities (e.g. Simmons et al. 2012). Moreover, P and Pn phases arrive together or very close to each other in epicentral distances of ∼2°–5°, and may have opposite polarities and/or different take-off angles. However, Kilb et al. (2000) found that errors introduced due to this condition are not significant, as long as a good azimuthal coverage and good sampling of the focal sphere exists. Other possible sources of errors include wrong polarity pickings that are either a result of reversed polarity stations, or false automatic pickings. The FilterPicker that is used in this method is an automatic broadband picker, well suited for real-time, early-warning purposes (Lomax et al. 2012). After applying the residual criteria (see Section 2 for details) the automatic pickings are correct up to roughly 80 per cent when compared with manual pickings within the epicentral distance range of interest (up to 90°). The FilterPicker has a few tuning parameters which control the time window of the signal, the filtering parameters and several thresholds associated with amplitude measurements. Vassallo et al. (2012) showed that tuning of these parameters, can optimize the performance of the algorithm with respect to a specific application. Benchmark tests showed that in roughly half of the cases (41  per  cent) it is possible to constrain the source mechanism by using the minimum rotation angle (Kagan 1990, 1991) as a mean of quantification. Moreover, in many cases ( ∼ 65  per  cent) it is possible to constrain at least one nodal plane (ϕ, δ) with errors in strike and dip up to 30° (see e.g. Supporting Information Figs S3 and S7). Other studies have also shown that a similar or even higher variability is expected when comparing results obtained from different methods and data (Hjörleifsdóttir & Ekström 2010; Ferreira et al. 2011; Weston et al. 2011; Lentas et al. 2013). Even though it is not always possible to determine the entire mechanism very accurately, constraining at least one plane could prove to be a useful tool for moment tensor inversions (e.g. Nakanishi & Kanamori 1984; Kawakatsu & Cadena 1991). For example, moment tensor inversion techniques using long period surface waves suffer from a trade-off between the seismic moment and the dip angle when solving for shallow-dip, thrust, earthquakes (Kanamori & Given 1981; Tsai et al. 2011). In this case, the excitation of long period surface waves is proportional to both the seismic moment and the sine of the dip angle (Mosin(2δ)), but the two parameters cannot be resolved separately. This trade-off can have serious implications on the earthquake’s magnitude determination, thus, the robust calculation of the dip angle is of great importance. Specifically, Tsai et al. (2011) studied the 2011, MW 9.0, Tohoku earthquake and found that this trade-off can introduce an error of 0.04 in magnitude every 2° of error in dip. Lentas et al. (2014) verified this by studying several thrust earthquakes and found, on average, an error of 0.03 in moment magnitude associated with every 2° of error in dip angle. Moreover, the rake angle is always more difficult to constrain for this particular category of earthquakes, as the strike and rake are mutually correlated in a way that their difference is kept constant as a result of a geometric trade-off (Han et al. 2011, 2013). For instance, constraining the dip-slip components (Mrθ, Mrϕ) in a moment tensor inversion based on strike and dip angles determined from a P-wave polarity based technique could in theory address this issue. However, since the fault plane orientation may vary along the rupture, and given the different hypocentre location concepts that different techniques depend on (rupture nucleation point, centroid location), this procedure could potentially bias the final result for large magnitude earthquakes (M > 6.5) where the point seismic source is no longer a good approximation. As mentioned in Section 3 there are many cases of well constrained mechanisms associated with earthquakes showing excessive azimuthal gaps (>90°). Even though a good azimuthal coverage is fundamental, and always desirable, it does not necessarily guarantee that the sampling of the focal sphere is adequate. FPFIT only uses the station distribution ratio as a measure of the mechanism’s robustness (Reasenberg & Oppenheimer 1985), while the HASH method has a number of metrics to assess the quality of the determined focal mechanisms. Notably, these are the maximum azimuthal and take-off angle gaps, the station distribution ratio, the fault plane uncertainties, the average misfit and the probability of the determined mechanism to be close to the real mechanism, all implemented in a mechanism quality classification scheme (Hardebeck & Shearer 2002). Tests carried out by Kilb & Hardebeck (2006) suggested that the use of fault plane uncertainties as a mechanism quality criteria leads to a substantial increase in the fault mechanism stability. Similarly, in this study the nodal plane uncertainty and the probability were also found to be good metrics for the classification of the obtained mechanisms. Nevertheless, the azimuthal gap and the misfit are probably the best parameters that define accurate mechanism solutions overall (see Supporting Information Fig. S4). The application of the presented technique on a set of earthquakes from the reviewed ISC bulletin revealed the difficulties and limitations of attempting to determine focal mechanisms from first motion polarities on a global scale. Specifically, the current set is dominated by relatively small magnitude earthquakes, notably, 76  per  cent of them have $$m_{b}^{{\rm ISC}}$$ < 5.0. As a consequence, in most of the cases the station distribution was not ideal due to the lack of high signal-to-noise ratio recordings that would allow the auto-picker to pick a robust polarity. Only 20  per  cent of the obtained mechanisms are associated with azimuthal gaps up to 90°, and 12  per  cent of the auto-picked phases are within an epicentral distance of up to 3° which yields limited resolution. Nevertheless, it was still feasible to determine a few earthquakes with no previously reported source mechanisms. Based on this performance 100–200 source mechanisms, with no previously reported solutions, could be added to the ISC bulletin for each reviewed year. This is comparable to other ISC special bulletins, such as the IASPEI Ground Truth reference events list, for which approximately 200 events are added for each reviewed year (e.g. Bondár & McLaughlin 2009). Apart from the obvious benefit of determining earthquake mechanisms where no previously reported source information exists, the current technique can offer an alternative approach to pre-existing source models and give evidence of complex source processes that may be due to the use of different data. Moreover, the teleseismic auto-picked P-phase arrival times could be an additional set to the ISC bulletin and the polarities could also be used in order to check the validity of reported mechanisms. The latter could be a valuable addition to the ISC procedures as, for the time being, the analysis focuses only on the investigation of the reported phase arrival times and hypocentres. Any additional source information provided by the reporting agencies is not confirmed or validated. Another potential benefit from the systematic comparison of auto-picked and reported polarities with the ISC obtained and reported mechanisms could lead to the identification of possible station reversals. This evidence could be forwarded to the agencies that maintain the stations in order to further investigate any possible errors and provide solutions where necessary. Additional data, such as S/P amplitude ratios and S-wave polarities could potentially improve the accuracy of the not well constrained mechanisms, as P-wave amplitudes become larger in the vicinity of the P and T axes and smaller close to the nodal planes, while the S-wave amplitudes become larger close to the nodal planes (e.g. De Natale et al. 1991; Hardebeck & Shearer 2003). Nevertheless, their use is not always straightforward. Notably, S/P amplitude ratios need to be corrected for attenuation and path effects, and they are prone to error due to high levels of noise in the waveform data. Moreover, S-wave polarities can be affected by anisotropy (Shin et al. 1991), whereas the automatic picking of the direct S-wave phase is not always accurate due to intermixing with reflected and/or converted seismic phases. 6 CONCLUSIONS In this study, the possibility of routinely obtaining focal mechanisms from P-wave first motion polarities on a global scale was set under investigation. The main purpose is to enhance the reviewed ISC bulletin especially for cases where there is no reported information on earthquake source mechanisms. It was shown that this is a demanding task for intermediate magnitude earthquakes (mb 4–5) which show the greatest gap in source mechanism information in the bulletin, nevertheless, it may be possible in a few cases. Moreover, the obtained mechanisms could also benefit cases of vast intraevent source mechanism variability that are quite common in the ISC bulletin. Finally, the use of auto-pickings from waveform data (polarities, arrival times) that are not previously reported to the ISC (teleseismic in most cases) could also be made available and used in earthquake locations and/or tomographic studies. Acknowledgements The author wishes to thank the editor Dr Martin Schimmel, and two anonymous reviewers for their comments and suggestions which helped to improve this manuscript. The author also wants to thank Dr Bertrand Delouis for useful remarks and discussion. This research was supported by a National Science Foundation Grant (NSF Grant #1417970). The author gratefully acknowledges the availability of global seismograms from the IRIS and Orfeus data centres, as well as the availability of the HASH source code from the IRIS data centre, the TauP source code from the University of South Carolina and the FilterPicker source code from ALomax scientific. REFERENCES Amante C., Eakins B.W., 2009. ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis, NOAA Technical Memorandum NESDIS NGDC-24. Bailey I.W., Becker T.W., Ben-Zion Y., 2009. Patterns of co-seismic strain computed from southern California focal mechanisms, Geophys. J. Int. , 177, 1015– 1036. Google Scholar CrossRef Search ADS   Bailey I.W., Ben-Zion Y., Becker T.W., Holschneider M., 2010. Quantifying focal mechanism heterogeneity for fault zones in central and southern California, Geophys. J. Int. , 183, 433– 450. Google Scholar CrossRef Search ADS   Bassin C., Laske G., Masters G., 2000. The Current Limits of Resolution for Surface Wave Tomography in North America, EOS, Trans. Am. geophys. Un. , 81, F897. Bernardi F., Braunmiller J., Kradolfer U., Giardini D., 2004. Automatic regional moment tensor inversion in the European-Mediterranean region, Geophys. J. Int. , 157, 703– 716. Google Scholar CrossRef Search ADS   Bondár I., McLaughlin K.L., 2009. A new ground truth data set for seismic studies, Geophys. Res. Lett. , 80, 465– 472. Bondár I., Storchak D., 2011. Improved location procedures at the International Seismological Centre, Geophys. J. Int. , 186, 1220– 1244. Google Scholar CrossRef Search ADS   Braunmiller J., Kradolfer U., Baer M., Giardini D., 2002. Regional moment tensor determination in the European-Mediterranean area - initial results, Tectonophysics , 356, 5– 22. Google Scholar CrossRef Search ADS   Cesca S., Buforn E., Dahm T., 2006. Amplitude spectra moment tensor inversion of shallow earthquakes in Spain, Geophys. J. Int. , 166, 839– 854. Google Scholar CrossRef Search ADS   Cesca S., Sen A.T., Dahm T., 2014. Seismicity monitoring by cluster analysis of moment tensors, Geophys. J. Int. , 196, 1813– 1826. Google Scholar CrossRef Search ADS   Collins J.A., Smith D.K., McGuire J., 2012. Seismicity of the Atlantis Massif detachment fault 30°N at the Mid-Atlantic Ridge, Geochem. Geophys. Geosyst. , 13( 1), doi:10.1029/2012GC004210. Crotwell H.P., Owens T.J., Ritsema J., 1999. The TauP Toolkit: Flexible seismic travel-time and ray-path utilities, Seismol. Res. Lett. , 70, 154– 160. Google Scholar CrossRef Search ADS   Delouis B., 2014. FMNEAR: Determination of focal mechanism and first estimate of rupture directivity using near-source records and a linear distribution of point sources, Bull. seism. Soc. Am. , 104( 3), 1479– 1500. Google Scholar CrossRef Search ADS   Delouis B., Nocquet J.-M., Vallée M., 2010. Slip distribution of the February 27, 2010 MW = 8.8 Maule Earthquake, central Chile, from static and high-rate GPS, InSAR, and broadband teleseismic data, Geophys. Res. Lett. , 37, L17305, doi:10.1029/2010GL043899. Google Scholar CrossRef Search ADS   DeMets C., Jansma P.E., Mattioli G.S., Dixon T.H., Farina F., Bilham R., Calais E., Mann P., 2000. GPS geodetic constraints on Caribbean-North America plate motion, Geophys. Res. Lett. , 27, 437– 440. Google Scholar CrossRef Search ADS   DeMets C., Gordon R.G., Argus D.F., 2010. Geologically current plate motions, Geophys. J. Int. , 181, 1– 80. Google Scholar CrossRef Search ADS   De Natale G., Ferraro A., Virieux J., 1991. A probability method for local earthquake focal mechanisms, Geophys. Res. Lett. , 18( 4), 613– 616. Google Scholar CrossRef Search ADS   Duputel Z., Rivera L., Kanamori H., Hayes G., 2012. W phase source inversion for moderate to large earthquakes (1990–2010), Geophys. J. Int. , 189, doi: 10.1111/j.1365-246X.2012.05419.x. Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth planet. Inter. , 25( 4), 297– 356. Google Scholar CrossRef Search ADS   Dziewonski A.M., Chou T.A., Woodhouse J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, J. geophys. Res. , 86( B4), 2825– 2852. Google Scholar CrossRef Search ADS   Ekström G., Nettles M., Dziewonski A.M., 2012. The global CMT Project 2004–2010: Centroid-moment tensors for 13,017 earthquakes, Phys. Earth planet. Int. , 200–201, 1– 9. Google Scholar CrossRef Search ADS   Ferreira A.M.G., Weston J., Funning G.J., 2011. Global compilation of interferometric synthetic aperture radar earthquake source models: 2. Effects of 3-D Earth structure. J. geophys. Res. , 116, B08409, doi:10.1029/2010JB008132. Google Scholar CrossRef Search ADS   Fritz H.M., Papantoniou A., Biukoto L., Gilly A., Wei Y., 2014. The Solomon Islands Tsunami of 6 February 2013 in the Santa Cruz Islands: Field Survey and Modeling, in EGU General Assembly , 27 April–2 May, Vienna, Austria. Gilbert F., Dziewonski A.M., 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanism from seismic spectra, Phil. Trans. R. Soc. A , 278, 187– 209. Google Scholar CrossRef Search ADS   Han S.C., Sauber J., Riva R., 2011. Contribution of satellite gravimetry to understanding seismic source processes of the 2011 Tohoku-Oki earthquake, Geophys. Res. Lett. , 38, L24312, doi:10.1029/2011GL049975. Han S.C., Riva R., Sauber J., Okal E., 2013. Source parameter inversion for recent great earthquakes from a decade-long observation of global gravity fields, J. geophys. Res. , 118, 1240– 1267. Google Scholar CrossRef Search ADS   Harada T., Murotani S., Satake K., 2013. A deep outer-rise reverse-fault earthquake immediately triggered a shallow normal-fault earthquake: the 7 December 2012 off-Sanriku earthquake (MW 7.3), Geophys. Res. Lett. , 40, 4214– 4219. Google Scholar CrossRef Search ADS   Hardebeck J.L., Shearer M.P., 2002. A new method for determining first-motion focal mechanisms, Bull. seism. Soc. Am. , 92( 6), 2264– 2276. Google Scholar CrossRef Search ADS   Hardebeck J.L., Shearer M.P., 2003. Using S/P amplitude ratios to constrain the focal mechanisms of small earthquakes, Bull. seism. Soc. Am. , 93( 6), 2434– 2444. Google Scholar CrossRef Search ADS   Hartigan J., Wong M., 1979. Algorithm AS 136: a K-means clustering algorithm, Appl. Stat. , 28( 1), 100– 108. Google Scholar CrossRef Search ADS   Hayes G.P., Rivera L., Kanamori H., 2009. Source inversion of the W-phase: Real-time implementation and extension to low magnitudes, Seismol. Res. Lett. , 80( 5), 817– 822. Google Scholar CrossRef Search ADS   Hayes G.P., Wald D.J., Johnson L., 2012. Slab1.0: A three-dimensional model of global subduction zone geometries, J. geophys. Res. , 117, B01302, doi:10.1029/2011JB008524. Hjörleifsdóttir V., Ekström G., 2010. Effects of three-dimensional Earth structure on CMT earthquake parameters, Phys. Earth planet. Inter. , 179, 178– 190. Google Scholar CrossRef Search ADS   International Seismological Centre, 2014. On-line Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent. , Thatcham, United Kingdom. Kagan Y.Y., 1990. Random stress earthquake statistics: spatial dependence, Geophys. J. Int. , 102, 573– 583. Google Scholar CrossRef Search ADS   Kagan Y.Y., 1991. 3-D rotation of double-couple earthquake sources, Geophys. J. Int. , 106, 709– 716. Google Scholar CrossRef Search ADS   Kanamori H., Given J.W., 1981. Use of long-period surface waves for rapid determination of earthquake-source parameters, Phys. Earth planet. Inter. , 27, 8– 31. Google Scholar CrossRef Search ADS   Kanamori H., Rivera L., 2008. Source inversion of W phase: speeding up seismic tsunami warning, Geophys. J. Int. , 175, 222– 238. Google Scholar CrossRef Search ADS   Katsumata K., Kosuga M., Katao H. & the Japanese University Group of the Joint Sismic Observations at, NKTZ, 2010. Focal mechanisms and stress field in the Atotsugawa fault area, central Honshu, Japan, Earth Planets Space , 62, 367– 380. Google Scholar CrossRef Search ADS   Katsumata K., Kosuga M., Katao H., Yamada T., Kato A. & the Research Group for the Joint Seismic Observations at the Nobi area, 2015. Focal mechanisms and stress field in the Nobi fault area, central Japan, Earth Planets Space , 67( 99), doi:10.1186/s40623-015-0275-2. Kawakatsu H., Cadena G.P., 1991. Focal mechanisms of the March 6, 1987 Ecuador earthquakes -CMT inversion with a first motion constraint, J. Phys. Earth. , 39, 589– 597. Google Scholar CrossRef Search ADS   Kennett B.L.N., Engdahl E.R., Buland R., 1995. Constraints on seismic velocities in the Earth from traveltimes, Geophys. J. Int. , 122, 108– 124. Google Scholar CrossRef Search ADS   Kilb D., Hardebeck J.L., 2006. Fault parameter constraints using relocated earthquakes: A validation of first-motion focal-mechanism data, Bull. seism. Soc. Am. , 96( 3), 1140– 1158. Google Scholar CrossRef Search ADS   Kilb D., Gomberg J., Bodin P., 2000. Triggering of earthquake aftershocks by dynamic stresses, Nature , 408, 570– 574. Google Scholar CrossRef Search ADS PubMed  Koch K., 1991. Moment tensor inversion of local earthquake data: I. Investigation of the method and its numerical stability with model calculations, Geophys. J. Int. , 106, 305– 319. Google Scholar CrossRef Search ADS   Laigle M.et al.  , 2013. Seismic structure and activity of the north-central Lesser Antilles subduction zone from an integrated approach: Similarities with the Tohoku forearc, Tectonophysics , 603, 1– 20. Google Scholar CrossRef Search ADS   Lay T., Duputel Z., Ye L., Kanamori H., 2013. The December 7, 2012 Japan Trench intraplate doublet (MW 7.2, 7.1) and interactions between near-trench intraplate thrust and normal faulting, Phys. Earth planet. Inter. , 220, 73– 78. Google Scholar CrossRef Search ADS   Lentas K., Ferreira A.M.G., Vallée M., 2013. Assessment of SCARDEC source parameters of global large (MW ≥ 7.5) subduction earthquakes, Geophys. J. Inter. , 195, 1989– 2004. Google Scholar CrossRef Search ADS   Lentas K., Ferreira A.M.G., Clévédé E., Roch J., 2014. Source models of great earthquakes from ultra low-frequency normal mode data, Phys. Earth planet. Inter. , 233, 41– 67. Google Scholar CrossRef Search ADS   Lomax A., Satriano C., Vassallo M., 2012. Automatic picker developments and optimization: FilterPicker–a robust, broadband picker for real-time seismic monitoring and earthquake early warning, Seismol. Res. Lett. , 83( 3), 531– 540. Google Scholar CrossRef Search ADS   Manaker D.M.et al.  , 2008. Interseismic Plate coupling and strain partitioning in the Northeastern Caribbean, Geophys. J. Int. , 174, 889– 903. Google Scholar CrossRef Search ADS   Mooney W.D., Laske G., Masters G., 1998. CRUST 5.1: A global crustal model at 5°× 5°, J. geophys. Res. , 103, 727– 747. Google Scholar CrossRef Search ADS   Nakanishi I., Kanamori H., 1984. Source mechanisms of twenty-six, large shallow earthquakes (Ms ≥ 6.5) during 1980 from P-wave first motion and long-period Rayleigh wave data, Bull. seism. Soc. Am. , 74( 3), 805– 818. Ozacar A.A., Beck S.L., Christensen D.H., 2003. Source process of the 3 November 2002 Denali fault earthquake (central Alaska) from teleseismic observations, Geophys. Res. Lett. , 30( 12), doi:10.1029/2003GL017272. Pondrelli S., Salimbeni S., Morelli A., Ekström G., Postpischl L., Vannucci G., Boschi E., 2011. European-Mediterranean regional centroid moment tensor catalog: solutions for 2005–2008, Phys. Earth planet. Inter. , 185( 3–4), 74– 81. Google Scholar CrossRef Search ADS   Reasenberg P.A., Oppenheimer D., 1985. FPFIT, FPPLOT and FPPAGE: FORTRAN computer programs for calculating and displaying earthquake fault-plane solutions, US Geol. Surv. Open-File Rept., 85– 739, 109p. Schneider J.F., Pennington W.D., Meyer R.P., 1987. Microseismicity and focal mechanisms of the intermediate-depth Bucaramanga Nest, Colombia, J. gephys. Res. , 92( B13), 13 913–13 926. Shearer M.P., 1997. Improving local earthquake locations using the L1 norm and waveform cross-correlation: application to the Whittier Narrows, California, aftershock sequence, J. geophys. Res. , 102, 8269– 8283. Google Scholar CrossRef Search ADS   Shearer M.P., 1998. Evidence from a cluster of small earthquakes for a fault at 18 km depth beneath Oak Ridge, southern California, Bull. seism. Soc. Am. , 88( 6), 1327– 1336. Shin X.R., Schneider J.F., Meyer R.P., 1991. Polarities of P and S waves, and shear wave splitting observed from the Bucaramanga Nest, Colombia, J. geophys. Res. , 96( B7), 12 069–12 082. 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( B10), doi:10.1029/2012JB009525. Stich D., Ammon C.J., Morales J., 2003. Moment tensor solutions for small and moderate earthquakes in the Ibero-Maghreb region, J. geophys. Res. , 108( B3), 2148, doi:10.1029/2002JB002057. Google Scholar CrossRef Search ADS   Toda S., Enescu B., 2011. Rate/state Coulomb stress transfer model for the CSEP Japan seismicity forecast, Earth Planets Space , 63, 171– 185. Google Scholar CrossRef Search ADS   Tsai V.C., Hayes G.P., Duputel Z., 2011. Constraints on the long-period moment-dip tradeoff for the Tohoku earthquake, Geophys. Res. Lett. , 38, L00G17, doi:10.1029/2011GL049129. Google Scholar CrossRef Search ADS   Vallée M., Charléty J., Ferreira A.M.G., Delouis B., Vergoz J., 2011. SCARDEC: a new technique for the rapid determination of seismic moment magnitude, focal mechanism and source time functions for large earthquakes using body-wave deconvolution, Geophys. J. Int. , 184, 338– 358. Google Scholar CrossRef Search ADS   Vassallo M., Satriano C., Lomax A., 2012. Automatic picker developments and optimization: A strategy for improving the performances of automatic phase pickers, Seismol. Res. Lett. , 83( 3), 541– 554. Google Scholar CrossRef Search ADS   Weston J., Ferreira A.M.G., Funning G.J., 2011. Global compilation of interferometric synthetic aperture radar earthquake source models: 1. Comparisons with seismic catalogs, J. geophys. Res. , 116, B08408, doi:10.1029/2010JB008131. Google Scholar CrossRef Search ADS   Yang W., Hauksson E., Shearer P.M., 2012. Computing a large refined catalog of focal mechanisms for southern California (1981–2010): temporal stability of the style of faulting, Bull. seism. Soc. Am. , 102( 3), 1179– 1194. Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Table S1. Obtained focal mechanisms for the earthquakes in the ISC reviewed bulletin, from 1/1/2013 to 31/1/2013, with $$m_b^{\rm ISC} >$$ 4.5, expressed in strike ϕ, dip δ, rake λ and P, T, B principal axes azimuth Θ and plunge Φ. The number of polarities (Nb of pol.) represents the total number including reported and auto-picked first motion polarities. In the cases where the number of available source mechanism in the ISC bulletin (Nb of mech.) is set to –, the entry in the table represents a focal mechanism solution for an earthquake with no previously reported mechanisms. Repeating origin times indicate multiple focal mechanism solutions for the same earthquake. Figure S1. Histogram showing the seismological agencies that contribute earth-quake mechanisms to the ISC Bulletin up to present. A detailed description of the agency codes can be found at http://www.isc.ac.uk/iscbulletin/agencies/. Note the logarithmic scale on the vertical axis. Figure S2. Histograms showing the distribution of ISC relocated earthquakes with over ten stations reporting first motion polarities shown in Fig. 1 in the main manuscript, with depth (a), the number of their associated reported first motion polarities (b), the minimum epicentral distance of the associated stations that report first motion polarities in the ISC bulletin (c), and the associated station azimuthal gaps calculated by taking into account the stations that report first motion polarities in the ISC bulletin (d). Figure S3. Obtained focal mechanisms (red) versus GCMT benchmark mechanisms (black) for the experiment testing the effect of location errors. The letter on top-left represents the quality code of the mechanism solution and the value on top-right is the minimum rotation angle between the obtained and benchmark mechanisms. The mechanisms are sorted from lowest (best fit) to largest rotation angle (worst fit) from bottom to top. Figure S4. Scatter plots comparing possible tradeoffs between earthquake source parameters and focal mechanism quality metrics for the experiment testing the effect of location errors. The colourscale represents the minimum rotation angle between the obtained and GCMT benchmark mechanisms. Figure S5. (a) A global map showing the locations of the earthquakes (GCMT beachballs) used in the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The size of the beachballs is proportional to their MW magnitudes; (b) A ternary plot showing the distribution of the earthquakes used in the same experiment with respect to their source mechanism type. The colourscale represents the minimum rotation angle between the obtained and benchmark mechanisms. Figure S6. Histograms showing the distribution of the obtained mechanisms in the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model, with respect to their quality classification codes (A, B, C, D), versus different quality criteria. Quality codes are shown on the top left corner of each subplot. Blue color represents all the obtained mechanisms, whereas red highlights the mechanisms which have rotation angles up to 35° when compared to the GCMT benchmark mechanisms. Mean μ and standard deviation values σ are also shown on every histogram subplot. Figure S7. Obtained focal mechanisms (red) versus GCMT benchmark mechanisms (black) for the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The letter on top-left represents the quality code of the mechanism solution and the value on top-right is the minimum rotation angle between the obtained and benchmark mechanisms. The mechanisms are sorted from lowest (best fit) to largest rotation angle (worst fit) from bottom to top. Figure S8. Scatter plots comparing possible tradeoffs between earthquake source parameters and focal mechanism quality metrics for the experiment testing the effects of location errors and detailed crustal structure (using CRUST2.0 and CRUST5.1) superimposed on the 1D ak135 velocity model. The colourscale represents the minimum rotation angle between the obtained and GCMT benchmark mechanisms. 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. © Crown copyright 2017. This article contains public sector information licensed under the Open Government Licence v3.0 (http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/).

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial