The intraplate Maranhão earthquake of 2017 January 3, northern Brazil: evidence for uniform regional stresses along the Brazilian equatorial margin

The intraplate Maranhão earthquake of 2017 January 3, northern Brazil: evidence for uniform... Summary Lithospheric stresses in intraplate regions can be characterized by many different wavelengths. In some areas, stresses vary over short distances of less than ∼100 km, but in other regions uniform stresses can be recognized for more than ∼1000 km or so. However, not all intraplate regions are well sampled with stress measurements to allow a good characterization of the lithospheric stresses. On 2017 January 3, a magnitude mb 4 earthquake occurred near the equatorial coast of the Maranhão State, an aseismic area of northern Brazil. Despite the few permanent stations in northern Brazil, a well-constrained strike-slip mechanism was obtained from regional moment-tensor inversion. A detailed analysis of the backazimuths of aftershocks recorded by the closest station (∼40 km away) allowed the identification of the fault plane to be the NNW–SSE trending nodal plane. An estimate of the rupture length, about 2 km, was also possible. The strike-slip mechanism has coast-parallel P axis and coast-perpendicular T axis, in agreement with most of the focal mechanisms found further to the east. The coast parallel P axis is also similar to the SHmax orientations from breakouts measurements further along the coast. The Maranhão earthquake fills an important gap of stress indicators in northern Brazil and suggests that the intraplate stress field is uniform along the 2000 km long northern coast. South America, Waveform inversion, Earthquake source observations, Dynamics: seismotectonics INTRODUCTION The observed intraplate stress field is often used as a key test for numerical models of mantle flow (e.g. Lithgow-Bertelloni & Guynn 2004; Naliboff et al. 2012). However, defining regional patterns of the intraplate stress field can be a challenge due to the lack of suitable earthquakes with well-constrained focal mechanisms. Assumpção et al. (2016) showed that, to a first order, a large region in mid-plate South America (from the Guyana shield to Central Brazil, about ∼2000 km long) seems to be characterized by NW–SE oriented maximum horizontal stresses (SHmax). Large regions of uniform stress field have been observed in other parts of the world, called ‘first-order patterns’, usually explained by predominance of far-field plate boundary forces (Heidbach et al. 2010). In South America more data are needed to confirm this uniform NW–SE compressional field, which has not been adequately explained yet by numerical models of global lithospheric stresses. Very few focal mechanisms in northern Brazil have been published due to the sparsity of seismic stations and relatively low seismicity. Here, we join recent efforts (Carvalho et al. 2016) to complement the database of stress indicators in this part of South America. On 2017 January 3, an earthquake with magnitude mb∼4 occurred in an aseismic area of the Parnaíba Basin, near the northern coast of Brazil (Fig. 1a), and was recorded by several stations of the Brazilian Seismographic Network—RSBR, which is composed of broadband stations (50 Hz to 120 s velocity response), with relatively low noise (Bianchi et al. 2015). Despite the large distances between most RSBR stations, a well-constrained strike-slip focal mechanism was obtained with coast-parallel P axis and coast-perpendicular T axis, as will be described in detail below. This kind of focal mechanism is similar to most mechanisms further to the east, and is consistent with SHmax orientations from breakout measurements further to the West along the equatorial coast (Fig. 1b). Figure 1. View largeDownload slide Seismicity and stresses in northern Brazil. (a) Epicentres of the Brazilian catalogue (Assumpção et al. 2014) with magnitudes >3.5 mb. Size of red circles indicates mb magnitude. Yellow star is the epicentre of the Maranhão 2017 event. Geological provinces are cratons (pink), foldbelts (grey) and intracratonic basins (yellow). GS, Guyana shield; CBS, Central Brazil shield (both part of the Amazon Craton); SFC, São Francisco Craton; BB, Borborema; TP, Tocantins provinces; PnB, Parnaíba Basin; PcB, Parecis Basin (from CPRM, Brazilian Geological Survey). (b) Maximum horizontal stress (SHmax) orientations from focal mechanisms (bars with white circle), in situ (IS, open bar) and breakout data (BO, solid bar). Bar colour denotes stress regime: blue for inverse faulting (T), green for strike-slip (SS) and red for normal faulting (N). Bar length denotes data quality, (C1, C2 or D) as used by the World Stress Map (http://www.world-stress-map.org/) and Assumpção et al. (2016). The strike-slip result of the Maranhão 2017 event is already included in the figure. Figure 1. View largeDownload slide Seismicity and stresses in northern Brazil. (a) Epicentres of the Brazilian catalogue (Assumpção et al. 2014) with magnitudes >3.5 mb. Size of red circles indicates mb magnitude. Yellow star is the epicentre of the Maranhão 2017 event. Geological provinces are cratons (pink), foldbelts (grey) and intracratonic basins (yellow). GS, Guyana shield; CBS, Central Brazil shield (both part of the Amazon Craton); SFC, São Francisco Craton; BB, Borborema; TP, Tocantins provinces; PnB, Parnaíba Basin; PcB, Parecis Basin (from CPRM, Brazilian Geological Survey). (b) Maximum horizontal stress (SHmax) orientations from focal mechanisms (bars with white circle), in situ (IS, open bar) and breakout data (BO, solid bar). Bar colour denotes stress regime: blue for inverse faulting (T), green for strike-slip (SS) and red for normal faulting (N). Bar length denotes data quality, (C1, C2 or D) as used by the World Stress Map (http://www.world-stress-map.org/) and Assumpção et al. (2016). The strike-slip result of the Maranhão 2017 event is already included in the figure. In addition, several aftershocks were well recorded by the nearest station (ROSB, ∼40 km away). A detailed analysis of the backazimuths and S–P differences allowed the determination of the rupture size of the main shock to be around 2 km. In this paper, we describe the determination of the focal mechanism by regional moment-tensor inversion. Data from a few local stations installed for aftershock studies are also consistent with the focal mechanism. The pattern of coast-parallel SHmax, along most of the equatorial margin, is then discussed. FOCAL MECHANISM The Maranhão earthquake had a teleseismic magnitude mb 4.0 and Ms 3.3 (International Data Centre). The Brazilian regional magnitude, equivalent to the teleseismic mb (Assumpção 1983; Assumpção et al. 2014), was mR 4.6 measured by 10 broadband stations from RSBR in the range of 200–1000 km from the epicentre. P-wave first motions were clear at only a few regional stations. For the main shock, 11 polarities were read from unfiltered seismograms of both vertical and radial components. We also included two polarities from a small aftershock of February 2 recorded by local stations installed a few weeks after the main event (stations VGM4 and VGM6, Fig. 2a). No reliable P-wave polarity was recorded at teleseismic stations due to the small magnitude and the strike-slip focal mechanism (as we show later). The P-wave polarities alone (Fig. 2a) were not enough to constrain the mechanism, although a strike-slip component is suggested. Figure 2. View largeDownload slide Focal mechanism solution. (a) P-wave polarity used to qualify the moment-tensor solutions. Push and pull are indicated by crosses and circles, respectively. Large and small symbols denote more and less reliable first-motion. Data from one aftershock recorded at two local stations (VGM4 and VGM6, large ‘X’ symbols) were also included. (b) Map of stations used in the moment-tensor inversion. (c) Frequency range test for the moment-tensor inversion: Plot of variance reduction (VR) for all four stations versus frequency band. The horizontal bar denotes the frequency band used with the DC solution in the middle. The beachball colours denote the fraction of polarity fit shown in the scale on the right. Note the stability of the fault-plane solutions fitting ∼80 per cent of the polarities Figure 2. View largeDownload slide Focal mechanism solution. (a) P-wave polarity used to qualify the moment-tensor solutions. Push and pull are indicated by crosses and circles, respectively. Large and small symbols denote more and less reliable first-motion. Data from one aftershock recorded at two local stations (VGM4 and VGM6, large ‘X’ symbols) were also included. (b) Map of stations used in the moment-tensor inversion. (c) Frequency range test for the moment-tensor inversion: Plot of variance reduction (VR) for all four stations versus frequency band. The horizontal bar denotes the frequency band used with the DC solution in the middle. The beachball colours denote the fraction of polarity fit shown in the scale on the right. Note the stability of the fault-plane solutions fitting ∼80 per cent of the polarities A regional moment-tensor inversion was carried out with the ISOLA code (Sokos & Zahradnik 2008; 2013). The deviatoric moment tensor is retrieved by least-squares fitting of the waveforms for a chosen frequency range. The centroid depth and origin time are found by grid-search. The fit of the solution is expressed by the variance reduction (VR) where VR = 1 means a perfect match between synthetic and observed traces. The reliability of the result can be quantified by the condition number (CN): low values, less than 5, imply that the moment-tensor is mathematically well resolved; larger values indicate an ill-posed problem, that is, the final solution may have no physical meaning. We used the four closest stations of the Brazilian Seismic Network (Fig. 2b). The nearest station, ROSB, is 40 km away. The other three distant stations (NBCL, PRPB and SMTB) are between 600 and 750 km away. To overcome the difficulty in using distant stations in regional moment-tensor inversion of small events, path-specific velocity models were derived from Love- and Rayleigh-wave group velocity dispersion following the method of Dias et al. (2016). Rayleigh and Love dispersion curves were retrieved from vertical and transverse components, respectively from the range of ∼ 6–30 s for the distant stations. They were inverted into shear-wave velocity models through the linearized inversion technique of Julià et al. (2000). During the inversion, we fixed the Moho depth according to the data compilation from Assumpção et al. (2013). For the ROSB path, we used the velocity model constrained by the location of one aftershock recorded at five local stations (see Supporting Information Table S2). These local data constrain the average crustal velocities, with the crustal thickness of 40 km taken from receiver functions at station ROSB (Assumpção et al. 2015). Crustal thickness at ROSB is not too critical as there is no mantle (Pn or Sn) phase recorded at 40 km distance. Using the path specific velocity models the Green's functions for the ISOLA code were calculated for a fixed epicentre and depths ranging from 2.5 to 19.5 km, allowing a centroid time shift of ±5 s for the waveform fitting. Fig. 2(c) shows the moment-tensor inversion for all four stations for many different frequency bands between 0.01 and 0.15 Hz with a width of at least 0.03 Hz. For each frequency band, we find the centroid time and depth that produces the best waveform fit. The P-wave polarities are not used directly in the moment-tensor inversion but they are only used as a quality check of the possible solutions. In Fig. 2(c), the beachballs are color-coded according to their polarity fit (PF). For example, PF = 0.85 means that 85 per cent of the polarities are satisfied by the mechanism. All polarities are shown in Fig. 2(a). Stations with sharp P-wave arrivals are represented as larger symbols. Note that a PF = 1 is impossible for the polarities in the dataset of Fig 2(a). Fig. 2(c) also shows that a satisfactory inversion (i.e. with a VR ≥ 0.5) can be obtained for frequencies up to 0.15 Hz, and the best fit is obtained for the frequency range 0.07–0.10 Hz, in which just surface waves are modelled. The main point of Fig. 2(c) is that the focal mechanism solution is very stable, that is, all solutions from different frequency ranges exhibit the same strike-slip mechanism. However, as pointed out by Dias et al. (2016), stability of the solution is a necessary condition but is not sufficient to guarantee that the focal mechanism is correct. Another feature from Fig. 2(c) is the presence of two ‘flipped’ solutions (red beachballs, i.e. solution with poor polarity fit). By flipped, we mean that those solutions have the compressional and dilatational quadrants swapped compared to the green solutions, that is, the green and red solutions have opposite rake angles. Flipped solutions appear when the frequency band used in the inversion produces seismograms close to sinusoidal, and so a time shift of half-period would also match the seismogram by changing the rake by 180° (like a cycle skip). One way to verify whether the solution is flipped or not is to check the fit to the first motion polarities. Also, checking whether the centroid time was shifted by more than 2–3 s can be an indication of cycle skipping. In conclusion, the main reason for the reliability of the strike-slip mechanisms (green solutions in Fig. 2c) is the fact that those solutions are able to explain simultaneously both the waveforms at four stations as well as the polarities, in addition to having a small condition number (CN = 1.8). The proximity and the high signal/noise ratio of station ROSB (40 km) allowed us to go to frequencies higher than 0.15 Hz aiming to fit both the P and S waves. Since we have clear body waves in the seismogram it is important to check whether our focal mechanism is able to explain them. Fig. 3 exhibits our preferred solution, a strike-slip mechanism with ESE-oriented P axis and NNE-oriented T axis obtained using the 0.07–0.10 Hz band for the distant stations (NBCL, PRPB and SMTB, fitting mostly surface waves), and the 0.07–0.4 Hz band for ROSB (fitting P and S waves). The effects of choosing such a wider frequency range are shown later in the text. Figure 3. View largeDownload slide Focal mechanism with ISOLA regional moment-tensor inversion. The black traces are the observed displacement waveforms (Z, vertical; N, north; E, east), and the red traces are the synthetics. Frequency band: 0.07–0.40 Hz for stations ROSB (39 km), and 0.07–0.10 Hz for NBCL (630 km), PRPB (740 km) and SMTB (755 km). The P- and S-wave arrival times, origin time, station azimuth and distance, as well as the amplitude and timescale for each station are indicated. The grey area in the beachball shows the range of possible solutions within the threshold of 90 per cent. Legend gives the information about focal mechanism: decomposition into DC, CLVD and VOL components, moment magnitude Mw, centroid depth (km), centroid time shift with respect to origin time (CT), nodal planes (NP1 and NP2), condition number (CN), variance reduction (VR), maximum Kagan angle among the solutions. The numbers in parentheses are the range used to search for possible solutions. For the Kagan angle, the number in parentheses is the median of the set. Figure 3. View largeDownload slide Focal mechanism with ISOLA regional moment-tensor inversion. The black traces are the observed displacement waveforms (Z, vertical; N, north; E, east), and the red traces are the synthetics. Frequency band: 0.07–0.40 Hz for stations ROSB (39 km), and 0.07–0.10 Hz for NBCL (630 km), PRPB (740 km) and SMTB (755 km). The P- and S-wave arrival times, origin time, station azimuth and distance, as well as the amplitude and timescale for each station are indicated. The grey area in the beachball shows the range of possible solutions within the threshold of 90 per cent. Legend gives the information about focal mechanism: decomposition into DC, CLVD and VOL components, moment magnitude Mw, centroid depth (km), centroid time shift with respect to origin time (CT), nodal planes (NP1 and NP2), condition number (CN), variance reduction (VR), maximum Kagan angle among the solutions. The numbers in parentheses are the range used to search for possible solutions. For the Kagan angle, the number in parentheses is the median of the set. This preferred solution and its uncertainty were determined as follows: for each solution from the time-depth grid search (Supporting Information Fig. S1), we calculated its polarity agreement (PF). The best solution was the one with the largest product VR*PF. VR*PF = 1 means a perfect fit between data and synthetic and a perfect match for all polarities, while an index of zero indicates that VR or/and PF are equal to zero. The best solution (highlighted as the largest beachball in Supporting Information Fig. S1) had VR = 0.76 and PF = 0.79, thus VR*PF = 0.60. The waveform fits and the source parameters of this preferred solution are shown in Fig. 3. The best mechanism is a pure strike-slip fault (strike/dip/rake = 340°/83°/−2°). The small CN (1.8) means that it is a well-conditioned problem. The Centroid Time (CT) was only −0.2 s, close to the estimated origin the time of the event. To estimate the uncertainty of this preferred solution, we chose a threshold of 90 per cent of the best VR*PF index. That is, all solutions with VR*PF > 0.9 of maximum VR*PF in the time-depth grid were taken as possible solutions for the Maranhão focal mechanism, that is, all solutions with VR*PF > 0.54 are accepted as possible solutions. These solutions are shown in Supporting Information Fig. S1 as larger beachballs. We use the Kagan angle (Kagan 1991) to measure the difference between these possible solutions and the best (preferred) solution. Zahradnik & Custódio (2012) consider focal mechanisms as equivalents for K-angles less than 30° and considerably different for >40°. The maximum Kagan angle of the suite of possible solutions is only 9°. The depth of the preferred solution is 7.5 km (largest beachball in Supporting Information Fig. S1, VR*PF = 0.60), but it is not well constrained since the depths of the acceptable (VR*PF > 0.54) solutions range from 5.0 to 13.5 km (larger beachballs in Supporting Information Fig. S1). As will be seen later, one small aftershock recorded by five local temporary stations had a depth of about 12 km, consistent with the possible depth range of the main shock. A large non-DC component (CLVD) of 59 per cent was found for the best solution but is not well constrained and is probably due to low signal-to-noise ratio of the distant stations. Solutions with 100 per cent DC are also present in the acceptable range of solutions (VR*PF > 0.54) as can be seen in Supporting Information Fig. S2 (plot VR*PF × DC per cent). Most solutions have high DC per cent. An Mw 4.3 was obtained with the moment-tensor inversion. This is consistent with the mb 4.0 and mR 4.6 values. As noted by Marco et al. (2010), strike-slip events tend to have smaller teleseismic mb magnitude compared with the regional mR value because the P-wave take-off angles are close to the vertically oriented null axis. All the results above were retrieved for the inversions up to 0.4 Hz at ROSB station. To test the effect of this wider frequency range, we compared how depth, centroid time shift, double-couple percentage and Mw change when inverting all stations for the 0.07–0.10 Hz range. This is shown in Supporting Information Fig. S2. VR*PF index decreases when using the higher frequency range. This is due to the fact that by using a wider frequency range we are fitting more data, and it is expected that the final solution has a lower VR. The distribution of possible depths is the same: the optimum depth is 7.5 km (VR*PF = 0.65 for the 0.07–0.10 Hz range and VR*PF = 0.60 for the 0.07–0.40 Hz range) and the depths satisfying the acceptable threshold of 0.90 are mainly between 5.0 and 12.5 km for both cases (Supporting Information Fig. S2). The optimum centroid time shift (CT) for the narrower frequency range is −0.2 s, close to −0.1 s suggested by the wider frequency band. In both cases the CT uncertainty remains less than ±0.5 s. A denser concentration of higher DC per cent was found for the lower frequencies (0.07–0.10 Hz) with the highest values around 80 per cent, whilst for higher frequencies (0.07–0.4 Hz) DC per cent exhibits highest values around 60 per cent. However, for both cases, solutions with DC > 90 per cent are also found. Results for the Mw do not change much: they are centred at 4.34 for the narrower frequency band and 4.32 for the wider frequency band. These tests show that the regional moment tensor solutions are stable and do not depend very much on the frequency band used for ROSB. EVENT DEPTH The regional moment-tensor inversion (Figs 2 and 3, and Supporting Information Fig. S1) indicated good solutions in the depth range of 5 to 13 km. The hypocentral determination of an aftershock with five local temporary stations gave a depth of 12 km (Supporting Information Figs S3 and S4 and Supporting Information Table S2). We tried to confirm the 12 km depth looking for pP phases at teleseismic distances. However, due to the small teleseismic magnitude (4.0 mb) and the strike-slip mechanism, a search for the depth phases pP and sP at teleseismic stations was not very successful. Only one station, the highly sensitive QSPA at the South Pole, showed signals that could be interpreted as the pP + sP phases. We calculated synthetic teleseismic P-waves at QSPA (87° distance) with the hudson96 code from Herrmann's (2013) package. The code is designed for teleseismic waveform modelling and it is based on the stationary phase approximation of Hudson (1969). Examples of the use of this code can be found in Ford et al. (2012) and Bobrov et al. (2016). One advantage of this code lies in the possibility of using different velocity models for the source, path and receiver. We used a crustal model at the source having a 300m-thick sedimentary layer, as expected for this area of the Parnaíba Basin (Otaviano C.P. Neto, Petrobras, personal communication, 2017). This is consistent with P to S and S to P converted phases beneath the local stations (see Supporting Information Fig. S4). The AK135 model (Kennett et al. 1995) was used for the upper mantle and receiver structure. The preferred fault plane solution (Figs 2 and 3) was used as a reference. We performed a grid-search for focal depths between 1 and 19 km, as done for the moment tensor inversion, and varied the focal mechanism allowing maximum Kagan angle = 9°. According to Cormier (1982), the attenuation factor (t*) for teleseismic P waves is around t* = 1 s. We allowed t* to vary from 0.1 up to 1.5. Hosseini (2016), in his fig. 2.9, showed a compilation of the half-durations for M > 4.5 earthquakes from the GCMT catalogue and obtained the relationship: half-duration = 0.00272 × exp(1.134 × M) for events between 2005 and 2014. The magnitude of the Maranhão earthquake, according to our regional modelling, is Mw 4.3, which corresponds to an expected duration of 0.7 s. In our search, we allowed duration from 0.1 up to 1.4 s. The result of the vertical component of QSPA inversion is shown in Fig. 4. The three histograms are the depth, t* attenuation and source duration, respectively. Although the signal-to-noise ratio at QSPA is far from ideal, the best estimate from this inversion is strike/dip/rake = 336°/84°/−10° (Kagan angle of 8° with respect to the solution of Figs 2 and 3). The best correlation coefficient (0.71) was found for a hypocentral depth of 12 or 13 km, generally consistent with the regional moment tensor inversion results (5–14 km) and the depth of the aftershock located with the local stations (12 km). The last two histograms show that the strong trade-off between t* and duration precludes any resolution to solve for t* and duration parameters independently. Figure 4. View largeDownload slide Teleseismic P-wave modelling at station QSPA (87° distance). All seismograms were bandpass filtered 1–3 Hz. (a) All tested focal mechanism solutions (grey lines) based on Fig. 3, allowing max Kagan angle of 9°. The black mechanism provided the best correlation (0.71) with the waveform. (b) Waveform fit (vertical component) for all ∼15 000 synthetic traces (grey lines) which had correlation coefficient better than 0.64 with the observed trace (black line), together with the depth phases pP and sP. The histograms (grey coded by the correlation coefficient) show how the depth (c), attenuation t* (d) and duration (e) vary in the grid-search. Despite the poor signal-to-noise ratio, it can be seen that a depth of about 12–13 km (histogram c) is well constrained. We do not have resolution for the t* and duration parameters. Figure 4. View largeDownload slide Teleseismic P-wave modelling at station QSPA (87° distance). All seismograms were bandpass filtered 1–3 Hz. (a) All tested focal mechanism solutions (grey lines) based on Fig. 3, allowing max Kagan angle of 9°. The black mechanism provided the best correlation (0.71) with the waveform. (b) Waveform fit (vertical component) for all ∼15 000 synthetic traces (grey lines) which had correlation coefficient better than 0.64 with the observed trace (black line), together with the depth phases pP and sP. The histograms (grey coded by the correlation coefficient) show how the depth (c), attenuation t* (d) and duration (e) vary in the grid-search. Despite the poor signal-to-noise ratio, it can be seen that a depth of about 12–13 km (histogram c) is well constrained. We do not have resolution for the t* and duration parameters. FAULT PLANE ORIENTATION FROM AFTERSHOCKS AT ROSB The Maranhão earthquake of 2017 January 3 was followed by tens of aftershocks well recorded by station ROSB, 40 km away. The epicentres of 24 events (the main shock and 23 aftershocks) were located by measuring the backazimuth (abbreviated as BAz) with the P-wave particle motion at ROSB. (Supporting Information Table S1 lists all aftershocks used in this work). Correlation of the P- and S-wave signals, after filtering all events with the same 2 to 8 Hz passband, provided accurate relative S–P traveltime differences. Fig. 5 shows two examples of P-wave particle motions. Although the particle motions do not show a clear difference on a visual inspection, the small difference in BAz (131.6° and 134.1°) is significant, as will be discussed below. Supporting Information Figs S5 and S6 show the correlated P- and S-wave signals. Figure 5. View largeDownload slide Two examples of P-wave particle motion used to get backazimuths. Panels (a) and (b) show north and east components for events 2 and 5, respectively: raw signal and 2–8 Hz bandpass filtered. Vertical marks ‘P’ denote the P-wave arrivals after cross-correlation; ‘A’ and ‘F’ denote the time windows (0.05 s before and 0.35 s after the P arrival) used for the particle motions shown in panels (c) and (d), and the backazimuth calculations. The magnitudes of events 2 and 5 were 3.1 and 2.2, respectively. Figure 5. View largeDownload slide Two examples of P-wave particle motion used to get backazimuths. Panels (a) and (b) show north and east components for events 2 and 5, respectively: raw signal and 2–8 Hz bandpass filtered. Vertical marks ‘P’ denote the P-wave arrivals after cross-correlation; ‘A’ and ‘F’ denote the time windows (0.05 s before and 0.35 s after the P arrival) used for the particle motions shown in panels (c) and (d), and the backazimuth calculations. The magnitudes of events 2 and 5 were 3.1 and 2.2, respectively. The BAz at ROSB were measured by finding the eigenvectors of the matrix of cross-correlations of the three components, using the tools of the Seismic Handler system (Stammler 1993). The epicentral distances were calculated from the S–P difference assuming an average Vp = 6.1 km s−1 and a Vp/Vs ratio of 1.73 (best value from a Wadati diagram with the local stations and ROSB). We assume that the events occurred in the upper crust (say, 12 km depth, according to the aftershock location using the local stations and the teleseismic modelling at QSPA) so that the depths of the events have little effect on the S–P times. We are more interested in the relative locations of the aftershocks, so that errors in the average Vp between the epicentral area and ROSB will affect the absolute locations of all events but not their positions relative to one another. The BAz measured for any event can vary by a few degrees depending on the filter passband and the chosen time-window. However, measuring the BAz of all events with the same bandpass filter and using exactly the same time window after correlating the highly similar P-wave signals, ensured that the results have enough resolution to allow a good determination of their relative epicentres, as shown in Fig. 6. Figure 6. View large Download slide (a) Backazimuth estimation for the main shock (event 1) and 23 aftershocks, using nine different filter bands and time windows. Large/small circles indicate magnitudes larger/smaller than 1.0, as indicated in the rightmost column. Event numbers in the left side correspond to Supporting Information Table S1. (b) Epicentres for all events with magnitude >1 determined with the BAz and S–P derived distance at ROSB station (all events assumed at the same fixed depth). Location is the medium value of nine different filter/time-window combinations and error bars are standard deviations from (a). The grey star is the 4.6 mR main shock; square is event no. 4 used as a reference in the cross-correlations. Except for three events (probably off the main fault), the tight cluster of 14 events with roughly N–S orientation indicates the probable rupture length. Figure 6. View large Download slide (a) Backazimuth estimation for the main shock (event 1) and 23 aftershocks, using nine different filter bands and time windows. Large/small circles indicate magnitudes larger/smaller than 1.0, as indicated in the rightmost column. Event numbers in the left side correspond to Supporting Information Table S1. (b) Epicentres for all events with magnitude >1 determined with the BAz and S–P derived distance at ROSB station (all events assumed at the same fixed depth). Location is the medium value of nine different filter/time-window combinations and error bars are standard deviations from (a). The grey star is the 4.6 mR main shock; square is event no. 4 used as a reference in the cross-correlations. Except for three events (probably off the main fault), the tight cluster of 14 events with roughly N–S orientation indicates the probable rupture length. To estimate the uncertainties of the relative locations using the BAz at ROSB, we repeated the same analysis for different frequency bands and time windows. We used six different frequency bands for the BAz determinations: 2–8, 2–10, 3–10, 2–14, and 5–16 Hz. Different time windows were used, ranging from 0.2 to 0.4 s around the maximum energy of the P waves. A total of nine measurements were carried out. For each measurement, the signal within the time-window was cross-correlated with that of the reference event (no. 4, magnitude 2.2, Supporting Information Table S1). A similar procedure (different frequencies and time-windows) was employed to align the S arrivals. Fig. 6(a) shows the nine BAz results for each of the 24 events. The medium values of the BAz and S–P times were used to locate the events. The standard deviations of the BAz were less than 1.7° for all events with magnitudes >1. For example, the average BAz for events 2 and 5 (with particle motion shown in Fig. 5) are 130.1 ± 1.4 and 134.9 ± 1.6, respectively. Fig. 6(b) shows the relative epicentres calculated with the median BAz (similar results are obtained with the average BAz). The error bars were given by the standard deviations of the BAz and S–P times. The delay times between the correlated P-wave signals from different events are usually less than 0.01 s. The same happens with the S arrivals. This means that the largest uncertainty in the relative location comes from the BAz (Fig. 6). The tight alignment of epicentres, with an approximate N–S orientation (Fig. 6b), seems well constrained in view of the relatively small epicentral uncertainties. The same pattern is seen with the epicentres calculated with our preferred frequency band (2–8 Hz, with the best signal to noise ratio) as shown in Supporting Information Fig. S7. The consistency of our epicentres, using different filter bands and time windows, strongly indicates that the N–S alignment is the orientation of the fault plane. Fig. 6 also shows that two events are located outside the main fault, perhaps in secondary nearby faults with similar focal mechanisms. Together with the P-wave polarity data and regional moment-tensor inversion, a left-lateral strike-slip mechanism on an NNW–SSE oriented fault is confirmed. It is interesting to observe that the mapped geological faults in this part of the Parnaíba basin do not show an NNW–SSE predominance, as indicated in Supporting Information Fig. S3. The P axis is oriented in the WNW–ESE direction (parallel to the equatorial margin) and the T axis NNE–SSW, perpendicular to the coast. This kind of mechanism is the most common in the states of Ceará and Rio Grande do Norte (further to the east), and suggests that the same stress field (coast parallel compression and coast-perpendicular extension) prevails along most of the Brazilian equatorial margin. RUPTURE LENGTH The tight cluster of epicentres in Fig. 6 most likely indicates that the rupture length of the main shock reached about 2 km. Nuttli's (1983) relation for intraplate earthquakes predicts a magnitude mb 4.4 for a 2 km long fault, which is not inconsistent with our average magnitude of 4.3 (mean of mb and mR, averaging out effects of radiation, as discussed earlier). On the other hand, the relation of Wells & Coppersmith (1994) suggests that a 2 km length would correspond to a magnitude Mw = 4.8 ± 0.2, higher than our Mw 4.3 obtained with the regional moment-tensor inversion. Using the magnitude Mw 4.3 and a rupture radius of 1 km (2 km length), we obtain a stress drop of Δσ = (7/16) Mo/r3 = 1.5 MPa. This is near the lower limit but still within the expected range of 1–40 MPa for intraplate events (Allmann & Shearer 2009). The Maranhão event seems to have a lower stress drop than expected for the average intraplate earthquake. DISCUSSION AND CONCLUSIONS This earthquake in Maranhão state provided a rare chance to sample the stress field in an aseismic area near the equatorial coast of Brazil. Despite the large distances of the regional stations, a moment-tensor inversion using path-specific models (Dias et al. 2016) was able to constrain the focal mechanism as a strike-slip event with a coast-parallel P axis and a coast-perpendicular T axis. The highly correlated P-wave signals from tens of aftershocks allowed the estimate of backazimuths with enough resolution to provide reliable relative epicentres showing the rupture to be about 2 km long with a roughly N–S orientation, in agreement with one of the nodal planes of the focal mechanism solution (Fig. 1b). Strike-slip mechanisms with coast-parallel P axes are the most frequent type of faulting in NE Brazil (further to the east, Fig. 1b) as shown by Ferreira et al. (1998) and Assumpção et al. (2016). The lateral density variations across the continent/oceanic lithosphere cause spreading stresses in the continent, that is, coast-perpendicular extension (as modelled by Stein et al. 1989), and as observed in Brazil (Assumpção 1998; Oliveira et al. 2015; 2016) and many other continental margins around the world (Tingay et al. 2006). SHmax measurements from borehole breakouts in other parts of the northern coast (Lima et al. 1997) also show the same ESE–WNW, coast-parallel trend (Fig. 1b). This seems to indicate that the whole Brazilian equatorial margin is part of a large, uniform stress province. Worldwide compilation of stress measurements (Heidbach et al. 2010) has shown that first-order stress provinces (up to 2000 km long) are found in several intraplate areas (e.g. in North America and NE Asia). These have been explained by the prevalence of far-field plate-boundary forces (e.g. Tingay et al. 2006; Heidbach et al. 2010). Shorter-wavelength provinces (the so-called second-order patterns with ∼100 km dimensions) can be caused by local forces due to large structural lithospheric contrasts (Tingay et al. 2006), and flexural effects. Near the Amazon river fan, for example, Watts et al. (2009) showed that the load of the sediments in the continental shelf should cause high flexural stresses with a coast-perpendicular extension onshore. Sediment load and continent/ocean transition are usually considered to have a ‘secondary’ effect on the stress field, compared to plate-boundary forces. However, the long extent of the NW–SE oriented Brazilian equatorial margin may cause an apparent first-order, ∼2000 km long stress province, which must be taken into account when modelling the stresses resulting from far-field plate boundary forces. The Maranhão earthquake was related to stresses commonly found in passive continental margins, far from present plate boundary processes. For this reason, it is a typical Stable Continental Region (SCR) earthquake, as defined by Johnston (1989). Because plate-boundary forces (such as Mid-Atlantic ridge-push and collisional forces at the convergence between Nazca and South American plates) can cause stresses across an entire plate (e.g. Heidbach et al. 2010), SCR earthquakes have been thought to follow the same strain loading/release cycles as plate-boundary events, the only difference being the much lower strain rates and longer recurrence intervals. However, an alternative hypothesis has been proposed (e.g. Stein & Liu 2009; Liu et al. 2010; Calais et al. 2016): intraplate stresses are near critical and roughly constant, with earthquakes occurring by localized small transient variations of stress or fault strength. There would be no significant slow strain accumulation with a long-term seismic cycle. Earthquakes along the entire Brazilian northern coast are subject to typical continental margin stresses, mainly due to lateral density contrasts across continental-oceanic crust, which are ‘non-renewable’ and do not follow seismic cycles. If these stresses predominate over long-range stresses from plate boundary forces, then the Maranhão earthquake and the seismicity along the northern Brazilian coast would favour the new hypothesis of a constant pre-stressed lithosphere. ACKNOWLEDGEMENTS We thank Petrobras for supporting the installation of the Brazilian Seismic Network (RSBR), and CPRM (Brazilian Geological Survey) for the present operation and maintenance of the network. Adriano Botelho and Marcelo Rocha helped with the installation of the local stations for the aftershock study. This work was carried out with the support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) grant 30.6547/2013-9 and CPRM grant 23809. RSBR data are freely available at www.sismo.iag.usp.br, (last accessed 2018 January 14). We thank Stephen Hicks and an anonymous reviewer for the many helpful comments and suggestions to improve the paper. REFERENCES Allmann B.P., Shearer P.M., 2009. Global variations of stress drop for moderate to large earthquakes, J. geophys. Res.  114 B01310, doi:10.1029/2008JB005821. https://doi.org/10.1029/2008JB005821 Google Scholar CrossRef Search ADS   Assumpção M., 1983. A regional magnitude scale for Brazil, Bull. seism. Soc. Am.  73 237– 246. Assumpção M., 1998. Seismicity and stresses in the Brazilian passive margin, Bull. seism. Soc. Am.  88( 1), 160– 169. Assumpção M., Feng M., Tassara A., Julià J., 2013. Models of crustal thickness for South America from seismic refraction, receiver functions and surface wave tomography, Tectonophysics  609 82– 96. https://doi.org/10.1016/j.tecto.2012.11.014 Google Scholar CrossRef Search ADS   Assumpção M., Ferreira J., Barros L., Bezerra H., França G., Barbosa J., Dourado J., 2014. Intraplate seismicity in Brazil, in Intraplate Earthquakes  . 50– 71, ed. Talwani P., Cambridge Univ. Press, Cambridge, doi:10.1017/CBO9781139628921.004. Google Scholar CrossRef Search ADS   Assumpção M., Bianchi M., Albuquerque D.F., França G.S.L., Barros L.V., 2015. Crustal thickness map in South America: an updated version, in 1st Brazilian Symposium of Seismology  Brasília, Dezembro 2015, Expanded Abstract. Assumpção M., Dias F.L., Zevallos I., Naliboff J.B., 2016. Intraplate stress field in South America from earthquake focal mechanisms, J. South Am. Earth Sci.  71 278– 295. https://doi.org/10.1016/j.jsames.2016.07.005 Google Scholar CrossRef Search ADS   Barros L.V., Assumpção M., 2011. Basement depths in the Parecis Basin (Amazon) with receiver functions from small local earthquakes in the Porto dos Gaúchos Seismic Zone, J. South Am. Earth Sci.  32 142– 151. Google Scholar CrossRef Search ADS   Bianchi M. et al.   2015. The Brazilian Seismographic Network: Historical overview and current status, Summary. Bull. Int. Seismol. Cent.  49( 1–6), 70– 90, http://www.isc.ac.uk/iscbulletin/summary/. Bobrov D.I., Kitov I.O., Rozhkov M.V., Friberg P., 2016. Towards global seismic monitoring of underground nuclear explosions using waveform cross correlation. Part I: Grand master events, Seism. Instrum.  52( 1), 43– 59. https://doi.org/10.3103/S0747923916010035 Google Scholar CrossRef Search ADS   Calais E., Camelbeeck T., Stein S., Liu M., Craig T.J., 2016. A new paradigm for large earthquakes in stable continental plate interiors, Geophys. Res. Lett.  43 10 621–10 637. Google Scholar CrossRef Search ADS   Carvalho J., Barros L.V., Zahradnik J. 2016. Focal mechanisms and moment magnitudes of micro-earthquakes in central Brazil by waveform inversion with quality assessment and inference of the local stress field, J. South Am. Earth Sci.  71 333– 343. https://doi.org/10.1016/j.jsames.2015.07.020 Google Scholar CrossRef Search ADS   Cormier V.F., 1982. The effect of attenuation on seismic body waves, Bull. seism. Soc. Am.  72( 6B), S169– S200. Dias F.L., Zahradnik J., Assumpção M., 2016. Path-specific, dispersion-based velocity models and moment tensors of moderate events recorded at few distant stations: examples from Brazil and Greece, J. South Am. Earth Sci.  71 344– 358. https://doi.org/10.1016/j.jsames.2016.07.004 Google Scholar CrossRef Search ADS   Ferreira J.M., Oliveira R.T., Takeya M.K., Assumpção M., 1998. Superposition of local and regional stresses in northeast Brazil: evidence from focal mechanisms around the Potiguar marginal basin, Geophys. J. Int.  134 2, 341– 355. https://doi.org/10.1046/j.1365-246x.1998.00563.x Google Scholar CrossRef Search ADS   Ford S.R., Walter W.R., Dreger D.S., 2012. Event discrimination using regional moment tensors with teleseismic-P constraints, Bull. seism. Soc. Am.  102( 2), 867– 872. https://doi.org/10.1785/0120110227 Google Scholar CrossRef Search ADS   Heidbach O., Tingay M., Barth A., Reinecker J., Kurfes D., Müller B., 2010. Global crustal stress pattern based on the World Stress Map database release 2008, Tectonophysics ., 482( 1–4), 3– 15. https://doi.org/10.1016/j.tecto.2009.07.023 Google Scholar CrossRef Search ADS   Herrmann R.B., 2013. Computer Programs in Seismology: An evolving tool for instruction and research, Seismol. Res. Lett.  84( 6), 1081– 1088. https://doi.org/10.1785/0220110096 Google Scholar CrossRef Search ADS   Hosseini K., 2016. Global multiple-frequency seismic tomography using teleseismic and core-diffracted body waves, Doctoral thesis , Ludwig Maximilians University, http://nbn-resolving.de/urn:nbn:de:bvb:19-195973. Hudson J.A. ( 1969). A quantitative evaluation of seismic signals at teleseismic distances—II: Body waves and surface waves from an extended source, Geophys. J. Int.  18( 4), 353– 370. https://doi.org/10.1111/j.1365-246X.1969.tb03574.x Google Scholar CrossRef Search ADS   Johnston A.C., 1989. The seismicity of ‘Stable Continental Interiors’, in Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound  . 299– 327, eds Gregersen S., Basham P.W., Springer. Google Scholar CrossRef Search ADS   Juliá J., Ammon C.J., Herrmann R.B., Correig A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, Geophys. J. Int.  143 99– 112. Google Scholar CrossRef Search ADS   Kagan Y., 1991. 3-D rotation of double-couple earthquake sources, Geophys. Int ., 106( 3), 709– 716. https://doi.org/10.1111/j.1365-246X.1991.tb06343.x 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( 1), 108– 124. https://doi.org/10.1111/j.1365-246X.1995.tb03540.x Google Scholar CrossRef Search ADS   Lima C., Nascimento E., Assumpção M., 1997. Stress orientations in Brazilian sedimentary basins from breakout analysis: implications for force models in the South American plate, Geophys. J. Int.  130( 1), 112– 124. https://doi.org/10.1111/j.1365-246X.1997.tb00991.x Google Scholar CrossRef Search ADS   Lithgow-Bertelloni C., Guynn J.H., 2004. Origin of the lithospheric stress field, J. geophys. Res.  109( B1), B01408, doi:10.1029/2003JB002467. https://doi.org/10.1029/2003JB002467 Google Scholar CrossRef Search ADS   Liu M., Stein S., Wang H., 2011. 2000 years of migrating earthquakes in North China: how earthquakes in midcontinents differ from those at plate boundaries, Lithosphere  3( 2), 128– 132. https://doi.org/10.1130/L129.1 Google Scholar CrossRef Search ADS   Marco E., Barcelos A., Assumpção M., 2010. Magnitude determination of regional earthquakes in Brazil: a comparison with the new IASPEI magnitude formulas, in AGU Meeting of the Americas  Iguazu Falls, Brazil, Abstract. Naliboff J.B., Lithgow-Bertelonni C., Ruff L.J., Koker N., 2012. The effects of lithospheric thickness and density structure on Earth's stress field, Geophys. J. Int.  188 1, 1– 17. https://doi.org/10.1111/j.1365-246X.2011.05248.x Google Scholar CrossRef Search ADS   Nuttli O.W., 1983. Average seismic source-parameter relations for mid-plate earthquakes, Bull. seism. Soc. Am.  73( 2), 519– 535. Oliveira P.H.S., Ferreira J.M., Bezerra F.H.R., Assumpção M., Nascimento A.F., Sousa M.O.L., Menezes E.A.S., 2015. Influence of the continental margin on the stress field and seismicity in the intraplate Acaraú seismic zone, NE Brazil, Geophys. J. Int.  202( 3), 1453– 1462. https://doi.org/10.1093/gji/ggv211 Google Scholar CrossRef Search ADS   Sokos E.N., Zahradnik J., 2008. ISOLA a Fortran code and a Matlab GUI to perform multiple-point source inversion of seismic data, Comput. Geosci.  34( 8), 967– 977. https://doi.org/10.1016/j.cageo.2007.07.005 Google Scholar CrossRef Search ADS   Sokos E.N., Zahradnik J., 2013. Evaluating centroid-moment-tensor uncertainty in the new version of ISOLA software, Seismol. Res. Lett.  84( 4), 656– 665. https://doi.org/10.1785/0220130002 Google Scholar CrossRef Search ADS   Stammler K. 1993. Seismichandler—programmable multichannel data handler for interactive and automatic processing of seismological analyses, Comput. Geosci.  19 2, 135– 140. https://doi.org/10.1016/0098-3004(93)90110-Q Google Scholar CrossRef Search ADS   Stein S., Liu M., 2009. Long aftershock sequences within continents and implications for earthquake hazard assessment, Nat. Geosci.  462( 7269), 87– 89. Stein S., Cloetingh S., Sleep N.H., Wortel R., ( 1989). Passive margin earthquakes, stresses and rheology, in Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound , pp. 231– 259, eds Gregersen S., Basham P.W., Kluwer Academic. Google Scholar CrossRef Search ADS   Tingay M.R.P, Müller B., Reinecker J., Heidbach O., 2006. State and origin of the present-day stress field in sedimentary basins: new results from the World Stress Map project, paper presented at the 41st U.S. Symposium on Rock Mechanics (USRMS), 06–1049, Golden, CO. Watts A.B., Rodger M., Peirce C., Greenroyd C.J., Hobbs R.W., 2009. Seismic structure, gravity anomalies, and flexure of the Amazon continental margin, NE Brazil, J. geophys. Res.  114 B07103, doi:10.1029/2008JB006259. https://doi.org/10.1029/2008JB006259 Google Scholar CrossRef Search ADS   Wells D.L., Coppersmith K.J., 1994. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, Bull. seism. Soc. Am.  84 974– 1002. Zahradnik J., Custódio S., 2012. Moment tensor resolvability: application to southwest Iberia, Bull. seism. Soc. Am.  102 1235– 1254. https://doi.org/10.1785/0120110216 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Grid search for centroid time shift and depth (‘source position’). The depth varies from 2.5 to 17 km and time shift from −5.0 and 5.0 s. The contours indicate the VR*PF index. The preferred solution is the mechanism with highest PF*VR value, highlighted with the largest beachball in the grid. For the uncertainty, we use a threshold of 0.9 for VR and PF: all solutions with VR*PF ≥ 0.9*(VR*PF)MAX (subscript MAX indicates the highest PF*VR) were considered and are shown as larger beachballs in the middle of the grid. Figure S2. Sensitivity to filter bands. Comparison between the moment tensor parameters from different frequency ranges: 0.07–0.10 (grey dots) and 0.07–0.40 Hz (black dots). Each dot represents an acceptable solution, that is, VR*PF ≥ 0.9*(VR*PF)MAX. Top-left histogram is depth, top-right represents the centroid time shift (CT), bottom-left shows double-couple percentage (DC) and bottom-right moment magnitude (Mw). Figure S3. Map of the temporary local stations (solid triangles, with numbers depicting stations VGM1 to VGM6) and the permanent station ROSB used to locate the aftershock of 2017 February 2 (yellow star, event 24 in Table S1). Solid lines are faults from the CPRM (Brazilian Geological Survey) database. Solid lines in the inset show the major geological provinces as in Fig. 1 of the main text. Figure S4. Seismograms of the 2017 February 2 aftershock (magnitude 0.6) at the three local stations VGM2, 4, 6. Sp and Ps are converted phases at the base of a ∼0.5 km thick sedimentary layer. Each station shows the vertical (Z), radial (R) and transverse (T) components. The Ps phase is best seen in the radial component; the Sp precursor in the vertical component. The time differences are roughly the same: Ps-P = S-Sp. See Barros & Assumpção (2011) for similar observations. Figure S5. Alignment of P arrival (left: 0.3 s before to 1.0 s after the P) and S arrival (right: 0.4 s before to 0.6 after), both at the vertical component for the first 15 aftershocks in Table S1. Bandpass filtered 2–8 Hz. Marks denoted by ‘S’ are preliminary attempts to identify the absolute arrivals (less accurate) picked before the correlation. The correlated arrivals for the S waves are denoted by ‘T1’. Figure S6. Vertical component traces aligned by the correlated S arrival (mark ‘T1’), showing significant differences in the S–P times. Events are 1 (main shock) to 15 from top to bottom. The two more distant events (larger S–P times) are nos. 8 and 15. No. 14 is the closest event. Figure S7. Backazimuths, dip angles, S–P and epicentres, calculated with 2–8 Hz bandpass filter and a time window of −0.05 to 0.35 s around the P onset. Small circles denote magnitudes <1 mR. The grey star is the 4.3 Mw main shock. Top left: S–P × backazimuth; top right: S–P × dip angle; bottom: epicentres determined with the BAz and the S–P derived distance (all events assumed at the same fixed depth). Except for three events (probably off the main fault), the main shock and 11 aftershocks show good correlation between S–P, BAz and Dip angle. The tight cluster of 12 events with roughly N–S orientation indicates the probable rupture length. Table S1. Main shock and aftershocks used in the epicentre determination. Two backazimuths (Baz) are shown: BAz-2–8 were measured with the 2–8 Hz bandpass filter and window −0.05 s before to 0.35 s after the P arrival; Baz-M is the medium value of nine different filter/window combinations. Magnitudes mR of the aftershocks were determined by comparing the P-wave rms amplitudes with that of the main shock. Table S2. Result of hypocentral determination of the 2017-Feb-02 07:15 aftershock (No. 24 in Table S1) using a best-fitting local model with Vp/Vs = 1.73 (from a Wadati diagram). The sedimentary layer (∼0.5 km thick, Vp = 3.0 km s−1, Vp/Vs = 2.0), similar to values of the Parecis Basin (Barros & Assumpção 2011), was replaced by station corrections (column ‘Delay’). We used the HYPO71 code modified for millisecond accuracy. Stations VGM1,5 had problems with GPS clock and only S–P times were used. 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. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

The intraplate Maranhão earthquake of 2017 January 3, northern Brazil: evidence for uniform regional stresses along the Brazilian equatorial margin

Loading next page...
 
/lp/ou_press/the-intraplate-maranh-o-earthquake-of-2017-january-3-northern-brazil-win5Pu6T8D
Publisher
Oxford University Press
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx560
Publisher site
See Article on Publisher Site

Abstract

Summary Lithospheric stresses in intraplate regions can be characterized by many different wavelengths. In some areas, stresses vary over short distances of less than ∼100 km, but in other regions uniform stresses can be recognized for more than ∼1000 km or so. However, not all intraplate regions are well sampled with stress measurements to allow a good characterization of the lithospheric stresses. On 2017 January 3, a magnitude mb 4 earthquake occurred near the equatorial coast of the Maranhão State, an aseismic area of northern Brazil. Despite the few permanent stations in northern Brazil, a well-constrained strike-slip mechanism was obtained from regional moment-tensor inversion. A detailed analysis of the backazimuths of aftershocks recorded by the closest station (∼40 km away) allowed the identification of the fault plane to be the NNW–SSE trending nodal plane. An estimate of the rupture length, about 2 km, was also possible. The strike-slip mechanism has coast-parallel P axis and coast-perpendicular T axis, in agreement with most of the focal mechanisms found further to the east. The coast parallel P axis is also similar to the SHmax orientations from breakouts measurements further along the coast. The Maranhão earthquake fills an important gap of stress indicators in northern Brazil and suggests that the intraplate stress field is uniform along the 2000 km long northern coast. South America, Waveform inversion, Earthquake source observations, Dynamics: seismotectonics INTRODUCTION The observed intraplate stress field is often used as a key test for numerical models of mantle flow (e.g. Lithgow-Bertelloni & Guynn 2004; Naliboff et al. 2012). However, defining regional patterns of the intraplate stress field can be a challenge due to the lack of suitable earthquakes with well-constrained focal mechanisms. Assumpção et al. (2016) showed that, to a first order, a large region in mid-plate South America (from the Guyana shield to Central Brazil, about ∼2000 km long) seems to be characterized by NW–SE oriented maximum horizontal stresses (SHmax). Large regions of uniform stress field have been observed in other parts of the world, called ‘first-order patterns’, usually explained by predominance of far-field plate boundary forces (Heidbach et al. 2010). In South America more data are needed to confirm this uniform NW–SE compressional field, which has not been adequately explained yet by numerical models of global lithospheric stresses. Very few focal mechanisms in northern Brazil have been published due to the sparsity of seismic stations and relatively low seismicity. Here, we join recent efforts (Carvalho et al. 2016) to complement the database of stress indicators in this part of South America. On 2017 January 3, an earthquake with magnitude mb∼4 occurred in an aseismic area of the Parnaíba Basin, near the northern coast of Brazil (Fig. 1a), and was recorded by several stations of the Brazilian Seismographic Network—RSBR, which is composed of broadband stations (50 Hz to 120 s velocity response), with relatively low noise (Bianchi et al. 2015). Despite the large distances between most RSBR stations, a well-constrained strike-slip focal mechanism was obtained with coast-parallel P axis and coast-perpendicular T axis, as will be described in detail below. This kind of focal mechanism is similar to most mechanisms further to the east, and is consistent with SHmax orientations from breakout measurements further to the West along the equatorial coast (Fig. 1b). Figure 1. View largeDownload slide Seismicity and stresses in northern Brazil. (a) Epicentres of the Brazilian catalogue (Assumpção et al. 2014) with magnitudes >3.5 mb. Size of red circles indicates mb magnitude. Yellow star is the epicentre of the Maranhão 2017 event. Geological provinces are cratons (pink), foldbelts (grey) and intracratonic basins (yellow). GS, Guyana shield; CBS, Central Brazil shield (both part of the Amazon Craton); SFC, São Francisco Craton; BB, Borborema; TP, Tocantins provinces; PnB, Parnaíba Basin; PcB, Parecis Basin (from CPRM, Brazilian Geological Survey). (b) Maximum horizontal stress (SHmax) orientations from focal mechanisms (bars with white circle), in situ (IS, open bar) and breakout data (BO, solid bar). Bar colour denotes stress regime: blue for inverse faulting (T), green for strike-slip (SS) and red for normal faulting (N). Bar length denotes data quality, (C1, C2 or D) as used by the World Stress Map (http://www.world-stress-map.org/) and Assumpção et al. (2016). The strike-slip result of the Maranhão 2017 event is already included in the figure. Figure 1. View largeDownload slide Seismicity and stresses in northern Brazil. (a) Epicentres of the Brazilian catalogue (Assumpção et al. 2014) with magnitudes >3.5 mb. Size of red circles indicates mb magnitude. Yellow star is the epicentre of the Maranhão 2017 event. Geological provinces are cratons (pink), foldbelts (grey) and intracratonic basins (yellow). GS, Guyana shield; CBS, Central Brazil shield (both part of the Amazon Craton); SFC, São Francisco Craton; BB, Borborema; TP, Tocantins provinces; PnB, Parnaíba Basin; PcB, Parecis Basin (from CPRM, Brazilian Geological Survey). (b) Maximum horizontal stress (SHmax) orientations from focal mechanisms (bars with white circle), in situ (IS, open bar) and breakout data (BO, solid bar). Bar colour denotes stress regime: blue for inverse faulting (T), green for strike-slip (SS) and red for normal faulting (N). Bar length denotes data quality, (C1, C2 or D) as used by the World Stress Map (http://www.world-stress-map.org/) and Assumpção et al. (2016). The strike-slip result of the Maranhão 2017 event is already included in the figure. In addition, several aftershocks were well recorded by the nearest station (ROSB, ∼40 km away). A detailed analysis of the backazimuths and S–P differences allowed the determination of the rupture size of the main shock to be around 2 km. In this paper, we describe the determination of the focal mechanism by regional moment-tensor inversion. Data from a few local stations installed for aftershock studies are also consistent with the focal mechanism. The pattern of coast-parallel SHmax, along most of the equatorial margin, is then discussed. FOCAL MECHANISM The Maranhão earthquake had a teleseismic magnitude mb 4.0 and Ms 3.3 (International Data Centre). The Brazilian regional magnitude, equivalent to the teleseismic mb (Assumpção 1983; Assumpção et al. 2014), was mR 4.6 measured by 10 broadband stations from RSBR in the range of 200–1000 km from the epicentre. P-wave first motions were clear at only a few regional stations. For the main shock, 11 polarities were read from unfiltered seismograms of both vertical and radial components. We also included two polarities from a small aftershock of February 2 recorded by local stations installed a few weeks after the main event (stations VGM4 and VGM6, Fig. 2a). No reliable P-wave polarity was recorded at teleseismic stations due to the small magnitude and the strike-slip focal mechanism (as we show later). The P-wave polarities alone (Fig. 2a) were not enough to constrain the mechanism, although a strike-slip component is suggested. Figure 2. View largeDownload slide Focal mechanism solution. (a) P-wave polarity used to qualify the moment-tensor solutions. Push and pull are indicated by crosses and circles, respectively. Large and small symbols denote more and less reliable first-motion. Data from one aftershock recorded at two local stations (VGM4 and VGM6, large ‘X’ symbols) were also included. (b) Map of stations used in the moment-tensor inversion. (c) Frequency range test for the moment-tensor inversion: Plot of variance reduction (VR) for all four stations versus frequency band. The horizontal bar denotes the frequency band used with the DC solution in the middle. The beachball colours denote the fraction of polarity fit shown in the scale on the right. Note the stability of the fault-plane solutions fitting ∼80 per cent of the polarities Figure 2. View largeDownload slide Focal mechanism solution. (a) P-wave polarity used to qualify the moment-tensor solutions. Push and pull are indicated by crosses and circles, respectively. Large and small symbols denote more and less reliable first-motion. Data from one aftershock recorded at two local stations (VGM4 and VGM6, large ‘X’ symbols) were also included. (b) Map of stations used in the moment-tensor inversion. (c) Frequency range test for the moment-tensor inversion: Plot of variance reduction (VR) for all four stations versus frequency band. The horizontal bar denotes the frequency band used with the DC solution in the middle. The beachball colours denote the fraction of polarity fit shown in the scale on the right. Note the stability of the fault-plane solutions fitting ∼80 per cent of the polarities A regional moment-tensor inversion was carried out with the ISOLA code (Sokos & Zahradnik 2008; 2013). The deviatoric moment tensor is retrieved by least-squares fitting of the waveforms for a chosen frequency range. The centroid depth and origin time are found by grid-search. The fit of the solution is expressed by the variance reduction (VR) where VR = 1 means a perfect match between synthetic and observed traces. The reliability of the result can be quantified by the condition number (CN): low values, less than 5, imply that the moment-tensor is mathematically well resolved; larger values indicate an ill-posed problem, that is, the final solution may have no physical meaning. We used the four closest stations of the Brazilian Seismic Network (Fig. 2b). The nearest station, ROSB, is 40 km away. The other three distant stations (NBCL, PRPB and SMTB) are between 600 and 750 km away. To overcome the difficulty in using distant stations in regional moment-tensor inversion of small events, path-specific velocity models were derived from Love- and Rayleigh-wave group velocity dispersion following the method of Dias et al. (2016). Rayleigh and Love dispersion curves were retrieved from vertical and transverse components, respectively from the range of ∼ 6–30 s for the distant stations. They were inverted into shear-wave velocity models through the linearized inversion technique of Julià et al. (2000). During the inversion, we fixed the Moho depth according to the data compilation from Assumpção et al. (2013). For the ROSB path, we used the velocity model constrained by the location of one aftershock recorded at five local stations (see Supporting Information Table S2). These local data constrain the average crustal velocities, with the crustal thickness of 40 km taken from receiver functions at station ROSB (Assumpção et al. 2015). Crustal thickness at ROSB is not too critical as there is no mantle (Pn or Sn) phase recorded at 40 km distance. Using the path specific velocity models the Green's functions for the ISOLA code were calculated for a fixed epicentre and depths ranging from 2.5 to 19.5 km, allowing a centroid time shift of ±5 s for the waveform fitting. Fig. 2(c) shows the moment-tensor inversion for all four stations for many different frequency bands between 0.01 and 0.15 Hz with a width of at least 0.03 Hz. For each frequency band, we find the centroid time and depth that produces the best waveform fit. The P-wave polarities are not used directly in the moment-tensor inversion but they are only used as a quality check of the possible solutions. In Fig. 2(c), the beachballs are color-coded according to their polarity fit (PF). For example, PF = 0.85 means that 85 per cent of the polarities are satisfied by the mechanism. All polarities are shown in Fig. 2(a). Stations with sharp P-wave arrivals are represented as larger symbols. Note that a PF = 1 is impossible for the polarities in the dataset of Fig 2(a). Fig. 2(c) also shows that a satisfactory inversion (i.e. with a VR ≥ 0.5) can be obtained for frequencies up to 0.15 Hz, and the best fit is obtained for the frequency range 0.07–0.10 Hz, in which just surface waves are modelled. The main point of Fig. 2(c) is that the focal mechanism solution is very stable, that is, all solutions from different frequency ranges exhibit the same strike-slip mechanism. However, as pointed out by Dias et al. (2016), stability of the solution is a necessary condition but is not sufficient to guarantee that the focal mechanism is correct. Another feature from Fig. 2(c) is the presence of two ‘flipped’ solutions (red beachballs, i.e. solution with poor polarity fit). By flipped, we mean that those solutions have the compressional and dilatational quadrants swapped compared to the green solutions, that is, the green and red solutions have opposite rake angles. Flipped solutions appear when the frequency band used in the inversion produces seismograms close to sinusoidal, and so a time shift of half-period would also match the seismogram by changing the rake by 180° (like a cycle skip). One way to verify whether the solution is flipped or not is to check the fit to the first motion polarities. Also, checking whether the centroid time was shifted by more than 2–3 s can be an indication of cycle skipping. In conclusion, the main reason for the reliability of the strike-slip mechanisms (green solutions in Fig. 2c) is the fact that those solutions are able to explain simultaneously both the waveforms at four stations as well as the polarities, in addition to having a small condition number (CN = 1.8). The proximity and the high signal/noise ratio of station ROSB (40 km) allowed us to go to frequencies higher than 0.15 Hz aiming to fit both the P and S waves. Since we have clear body waves in the seismogram it is important to check whether our focal mechanism is able to explain them. Fig. 3 exhibits our preferred solution, a strike-slip mechanism with ESE-oriented P axis and NNE-oriented T axis obtained using the 0.07–0.10 Hz band for the distant stations (NBCL, PRPB and SMTB, fitting mostly surface waves), and the 0.07–0.4 Hz band for ROSB (fitting P and S waves). The effects of choosing such a wider frequency range are shown later in the text. Figure 3. View largeDownload slide Focal mechanism with ISOLA regional moment-tensor inversion. The black traces are the observed displacement waveforms (Z, vertical; N, north; E, east), and the red traces are the synthetics. Frequency band: 0.07–0.40 Hz for stations ROSB (39 km), and 0.07–0.10 Hz for NBCL (630 km), PRPB (740 km) and SMTB (755 km). The P- and S-wave arrival times, origin time, station azimuth and distance, as well as the amplitude and timescale for each station are indicated. The grey area in the beachball shows the range of possible solutions within the threshold of 90 per cent. Legend gives the information about focal mechanism: decomposition into DC, CLVD and VOL components, moment magnitude Mw, centroid depth (km), centroid time shift with respect to origin time (CT), nodal planes (NP1 and NP2), condition number (CN), variance reduction (VR), maximum Kagan angle among the solutions. The numbers in parentheses are the range used to search for possible solutions. For the Kagan angle, the number in parentheses is the median of the set. Figure 3. View largeDownload slide Focal mechanism with ISOLA regional moment-tensor inversion. The black traces are the observed displacement waveforms (Z, vertical; N, north; E, east), and the red traces are the synthetics. Frequency band: 0.07–0.40 Hz for stations ROSB (39 km), and 0.07–0.10 Hz for NBCL (630 km), PRPB (740 km) and SMTB (755 km). The P- and S-wave arrival times, origin time, station azimuth and distance, as well as the amplitude and timescale for each station are indicated. The grey area in the beachball shows the range of possible solutions within the threshold of 90 per cent. Legend gives the information about focal mechanism: decomposition into DC, CLVD and VOL components, moment magnitude Mw, centroid depth (km), centroid time shift with respect to origin time (CT), nodal planes (NP1 and NP2), condition number (CN), variance reduction (VR), maximum Kagan angle among the solutions. The numbers in parentheses are the range used to search for possible solutions. For the Kagan angle, the number in parentheses is the median of the set. This preferred solution and its uncertainty were determined as follows: for each solution from the time-depth grid search (Supporting Information Fig. S1), we calculated its polarity agreement (PF). The best solution was the one with the largest product VR*PF. VR*PF = 1 means a perfect fit between data and synthetic and a perfect match for all polarities, while an index of zero indicates that VR or/and PF are equal to zero. The best solution (highlighted as the largest beachball in Supporting Information Fig. S1) had VR = 0.76 and PF = 0.79, thus VR*PF = 0.60. The waveform fits and the source parameters of this preferred solution are shown in Fig. 3. The best mechanism is a pure strike-slip fault (strike/dip/rake = 340°/83°/−2°). The small CN (1.8) means that it is a well-conditioned problem. The Centroid Time (CT) was only −0.2 s, close to the estimated origin the time of the event. To estimate the uncertainty of this preferred solution, we chose a threshold of 90 per cent of the best VR*PF index. That is, all solutions with VR*PF > 0.9 of maximum VR*PF in the time-depth grid were taken as possible solutions for the Maranhão focal mechanism, that is, all solutions with VR*PF > 0.54 are accepted as possible solutions. These solutions are shown in Supporting Information Fig. S1 as larger beachballs. We use the Kagan angle (Kagan 1991) to measure the difference between these possible solutions and the best (preferred) solution. Zahradnik & Custódio (2012) consider focal mechanisms as equivalents for K-angles less than 30° and considerably different for >40°. The maximum Kagan angle of the suite of possible solutions is only 9°. The depth of the preferred solution is 7.5 km (largest beachball in Supporting Information Fig. S1, VR*PF = 0.60), but it is not well constrained since the depths of the acceptable (VR*PF > 0.54) solutions range from 5.0 to 13.5 km (larger beachballs in Supporting Information Fig. S1). As will be seen later, one small aftershock recorded by five local temporary stations had a depth of about 12 km, consistent with the possible depth range of the main shock. A large non-DC component (CLVD) of 59 per cent was found for the best solution but is not well constrained and is probably due to low signal-to-noise ratio of the distant stations. Solutions with 100 per cent DC are also present in the acceptable range of solutions (VR*PF > 0.54) as can be seen in Supporting Information Fig. S2 (plot VR*PF × DC per cent). Most solutions have high DC per cent. An Mw 4.3 was obtained with the moment-tensor inversion. This is consistent with the mb 4.0 and mR 4.6 values. As noted by Marco et al. (2010), strike-slip events tend to have smaller teleseismic mb magnitude compared with the regional mR value because the P-wave take-off angles are close to the vertically oriented null axis. All the results above were retrieved for the inversions up to 0.4 Hz at ROSB station. To test the effect of this wider frequency range, we compared how depth, centroid time shift, double-couple percentage and Mw change when inverting all stations for the 0.07–0.10 Hz range. This is shown in Supporting Information Fig. S2. VR*PF index decreases when using the higher frequency range. This is due to the fact that by using a wider frequency range we are fitting more data, and it is expected that the final solution has a lower VR. The distribution of possible depths is the same: the optimum depth is 7.5 km (VR*PF = 0.65 for the 0.07–0.10 Hz range and VR*PF = 0.60 for the 0.07–0.40 Hz range) and the depths satisfying the acceptable threshold of 0.90 are mainly between 5.0 and 12.5 km for both cases (Supporting Information Fig. S2). The optimum centroid time shift (CT) for the narrower frequency range is −0.2 s, close to −0.1 s suggested by the wider frequency band. In both cases the CT uncertainty remains less than ±0.5 s. A denser concentration of higher DC per cent was found for the lower frequencies (0.07–0.10 Hz) with the highest values around 80 per cent, whilst for higher frequencies (0.07–0.4 Hz) DC per cent exhibits highest values around 60 per cent. However, for both cases, solutions with DC > 90 per cent are also found. Results for the Mw do not change much: they are centred at 4.34 for the narrower frequency band and 4.32 for the wider frequency band. These tests show that the regional moment tensor solutions are stable and do not depend very much on the frequency band used for ROSB. EVENT DEPTH The regional moment-tensor inversion (Figs 2 and 3, and Supporting Information Fig. S1) indicated good solutions in the depth range of 5 to 13 km. The hypocentral determination of an aftershock with five local temporary stations gave a depth of 12 km (Supporting Information Figs S3 and S4 and Supporting Information Table S2). We tried to confirm the 12 km depth looking for pP phases at teleseismic distances. However, due to the small teleseismic magnitude (4.0 mb) and the strike-slip mechanism, a search for the depth phases pP and sP at teleseismic stations was not very successful. Only one station, the highly sensitive QSPA at the South Pole, showed signals that could be interpreted as the pP + sP phases. We calculated synthetic teleseismic P-waves at QSPA (87° distance) with the hudson96 code from Herrmann's (2013) package. The code is designed for teleseismic waveform modelling and it is based on the stationary phase approximation of Hudson (1969). Examples of the use of this code can be found in Ford et al. (2012) and Bobrov et al. (2016). One advantage of this code lies in the possibility of using different velocity models for the source, path and receiver. We used a crustal model at the source having a 300m-thick sedimentary layer, as expected for this area of the Parnaíba Basin (Otaviano C.P. Neto, Petrobras, personal communication, 2017). This is consistent with P to S and S to P converted phases beneath the local stations (see Supporting Information Fig. S4). The AK135 model (Kennett et al. 1995) was used for the upper mantle and receiver structure. The preferred fault plane solution (Figs 2 and 3) was used as a reference. We performed a grid-search for focal depths between 1 and 19 km, as done for the moment tensor inversion, and varied the focal mechanism allowing maximum Kagan angle = 9°. According to Cormier (1982), the attenuation factor (t*) for teleseismic P waves is around t* = 1 s. We allowed t* to vary from 0.1 up to 1.5. Hosseini (2016), in his fig. 2.9, showed a compilation of the half-durations for M > 4.5 earthquakes from the GCMT catalogue and obtained the relationship: half-duration = 0.00272 × exp(1.134 × M) for events between 2005 and 2014. The magnitude of the Maranhão earthquake, according to our regional modelling, is Mw 4.3, which corresponds to an expected duration of 0.7 s. In our search, we allowed duration from 0.1 up to 1.4 s. The result of the vertical component of QSPA inversion is shown in Fig. 4. The three histograms are the depth, t* attenuation and source duration, respectively. Although the signal-to-noise ratio at QSPA is far from ideal, the best estimate from this inversion is strike/dip/rake = 336°/84°/−10° (Kagan angle of 8° with respect to the solution of Figs 2 and 3). The best correlation coefficient (0.71) was found for a hypocentral depth of 12 or 13 km, generally consistent with the regional moment tensor inversion results (5–14 km) and the depth of the aftershock located with the local stations (12 km). The last two histograms show that the strong trade-off between t* and duration precludes any resolution to solve for t* and duration parameters independently. Figure 4. View largeDownload slide Teleseismic P-wave modelling at station QSPA (87° distance). All seismograms were bandpass filtered 1–3 Hz. (a) All tested focal mechanism solutions (grey lines) based on Fig. 3, allowing max Kagan angle of 9°. The black mechanism provided the best correlation (0.71) with the waveform. (b) Waveform fit (vertical component) for all ∼15 000 synthetic traces (grey lines) which had correlation coefficient better than 0.64 with the observed trace (black line), together with the depth phases pP and sP. The histograms (grey coded by the correlation coefficient) show how the depth (c), attenuation t* (d) and duration (e) vary in the grid-search. Despite the poor signal-to-noise ratio, it can be seen that a depth of about 12–13 km (histogram c) is well constrained. We do not have resolution for the t* and duration parameters. Figure 4. View largeDownload slide Teleseismic P-wave modelling at station QSPA (87° distance). All seismograms were bandpass filtered 1–3 Hz. (a) All tested focal mechanism solutions (grey lines) based on Fig. 3, allowing max Kagan angle of 9°. The black mechanism provided the best correlation (0.71) with the waveform. (b) Waveform fit (vertical component) for all ∼15 000 synthetic traces (grey lines) which had correlation coefficient better than 0.64 with the observed trace (black line), together with the depth phases pP and sP. The histograms (grey coded by the correlation coefficient) show how the depth (c), attenuation t* (d) and duration (e) vary in the grid-search. Despite the poor signal-to-noise ratio, it can be seen that a depth of about 12–13 km (histogram c) is well constrained. We do not have resolution for the t* and duration parameters. FAULT PLANE ORIENTATION FROM AFTERSHOCKS AT ROSB The Maranhão earthquake of 2017 January 3 was followed by tens of aftershocks well recorded by station ROSB, 40 km away. The epicentres of 24 events (the main shock and 23 aftershocks) were located by measuring the backazimuth (abbreviated as BAz) with the P-wave particle motion at ROSB. (Supporting Information Table S1 lists all aftershocks used in this work). Correlation of the P- and S-wave signals, after filtering all events with the same 2 to 8 Hz passband, provided accurate relative S–P traveltime differences. Fig. 5 shows two examples of P-wave particle motions. Although the particle motions do not show a clear difference on a visual inspection, the small difference in BAz (131.6° and 134.1°) is significant, as will be discussed below. Supporting Information Figs S5 and S6 show the correlated P- and S-wave signals. Figure 5. View largeDownload slide Two examples of P-wave particle motion used to get backazimuths. Panels (a) and (b) show north and east components for events 2 and 5, respectively: raw signal and 2–8 Hz bandpass filtered. Vertical marks ‘P’ denote the P-wave arrivals after cross-correlation; ‘A’ and ‘F’ denote the time windows (0.05 s before and 0.35 s after the P arrival) used for the particle motions shown in panels (c) and (d), and the backazimuth calculations. The magnitudes of events 2 and 5 were 3.1 and 2.2, respectively. Figure 5. View largeDownload slide Two examples of P-wave particle motion used to get backazimuths. Panels (a) and (b) show north and east components for events 2 and 5, respectively: raw signal and 2–8 Hz bandpass filtered. Vertical marks ‘P’ denote the P-wave arrivals after cross-correlation; ‘A’ and ‘F’ denote the time windows (0.05 s before and 0.35 s after the P arrival) used for the particle motions shown in panels (c) and (d), and the backazimuth calculations. The magnitudes of events 2 and 5 were 3.1 and 2.2, respectively. The BAz at ROSB were measured by finding the eigenvectors of the matrix of cross-correlations of the three components, using the tools of the Seismic Handler system (Stammler 1993). The epicentral distances were calculated from the S–P difference assuming an average Vp = 6.1 km s−1 and a Vp/Vs ratio of 1.73 (best value from a Wadati diagram with the local stations and ROSB). We assume that the events occurred in the upper crust (say, 12 km depth, according to the aftershock location using the local stations and the teleseismic modelling at QSPA) so that the depths of the events have little effect on the S–P times. We are more interested in the relative locations of the aftershocks, so that errors in the average Vp between the epicentral area and ROSB will affect the absolute locations of all events but not their positions relative to one another. The BAz measured for any event can vary by a few degrees depending on the filter passband and the chosen time-window. However, measuring the BAz of all events with the same bandpass filter and using exactly the same time window after correlating the highly similar P-wave signals, ensured that the results have enough resolution to allow a good determination of their relative epicentres, as shown in Fig. 6. Figure 6. View large Download slide (a) Backazimuth estimation for the main shock (event 1) and 23 aftershocks, using nine different filter bands and time windows. Large/small circles indicate magnitudes larger/smaller than 1.0, as indicated in the rightmost column. Event numbers in the left side correspond to Supporting Information Table S1. (b) Epicentres for all events with magnitude >1 determined with the BAz and S–P derived distance at ROSB station (all events assumed at the same fixed depth). Location is the medium value of nine different filter/time-window combinations and error bars are standard deviations from (a). The grey star is the 4.6 mR main shock; square is event no. 4 used as a reference in the cross-correlations. Except for three events (probably off the main fault), the tight cluster of 14 events with roughly N–S orientation indicates the probable rupture length. Figure 6. View large Download slide (a) Backazimuth estimation for the main shock (event 1) and 23 aftershocks, using nine different filter bands and time windows. Large/small circles indicate magnitudes larger/smaller than 1.0, as indicated in the rightmost column. Event numbers in the left side correspond to Supporting Information Table S1. (b) Epicentres for all events with magnitude >1 determined with the BAz and S–P derived distance at ROSB station (all events assumed at the same fixed depth). Location is the medium value of nine different filter/time-window combinations and error bars are standard deviations from (a). The grey star is the 4.6 mR main shock; square is event no. 4 used as a reference in the cross-correlations. Except for three events (probably off the main fault), the tight cluster of 14 events with roughly N–S orientation indicates the probable rupture length. To estimate the uncertainties of the relative locations using the BAz at ROSB, we repeated the same analysis for different frequency bands and time windows. We used six different frequency bands for the BAz determinations: 2–8, 2–10, 3–10, 2–14, and 5–16 Hz. Different time windows were used, ranging from 0.2 to 0.4 s around the maximum energy of the P waves. A total of nine measurements were carried out. For each measurement, the signal within the time-window was cross-correlated with that of the reference event (no. 4, magnitude 2.2, Supporting Information Table S1). A similar procedure (different frequencies and time-windows) was employed to align the S arrivals. Fig. 6(a) shows the nine BAz results for each of the 24 events. The medium values of the BAz and S–P times were used to locate the events. The standard deviations of the BAz were less than 1.7° for all events with magnitudes >1. For example, the average BAz for events 2 and 5 (with particle motion shown in Fig. 5) are 130.1 ± 1.4 and 134.9 ± 1.6, respectively. Fig. 6(b) shows the relative epicentres calculated with the median BAz (similar results are obtained with the average BAz). The error bars were given by the standard deviations of the BAz and S–P times. The delay times between the correlated P-wave signals from different events are usually less than 0.01 s. The same happens with the S arrivals. This means that the largest uncertainty in the relative location comes from the BAz (Fig. 6). The tight alignment of epicentres, with an approximate N–S orientation (Fig. 6b), seems well constrained in view of the relatively small epicentral uncertainties. The same pattern is seen with the epicentres calculated with our preferred frequency band (2–8 Hz, with the best signal to noise ratio) as shown in Supporting Information Fig. S7. The consistency of our epicentres, using different filter bands and time windows, strongly indicates that the N–S alignment is the orientation of the fault plane. Fig. 6 also shows that two events are located outside the main fault, perhaps in secondary nearby faults with similar focal mechanisms. Together with the P-wave polarity data and regional moment-tensor inversion, a left-lateral strike-slip mechanism on an NNW–SSE oriented fault is confirmed. It is interesting to observe that the mapped geological faults in this part of the Parnaíba basin do not show an NNW–SSE predominance, as indicated in Supporting Information Fig. S3. The P axis is oriented in the WNW–ESE direction (parallel to the equatorial margin) and the T axis NNE–SSW, perpendicular to the coast. This kind of mechanism is the most common in the states of Ceará and Rio Grande do Norte (further to the east), and suggests that the same stress field (coast parallel compression and coast-perpendicular extension) prevails along most of the Brazilian equatorial margin. RUPTURE LENGTH The tight cluster of epicentres in Fig. 6 most likely indicates that the rupture length of the main shock reached about 2 km. Nuttli's (1983) relation for intraplate earthquakes predicts a magnitude mb 4.4 for a 2 km long fault, which is not inconsistent with our average magnitude of 4.3 (mean of mb and mR, averaging out effects of radiation, as discussed earlier). On the other hand, the relation of Wells & Coppersmith (1994) suggests that a 2 km length would correspond to a magnitude Mw = 4.8 ± 0.2, higher than our Mw 4.3 obtained with the regional moment-tensor inversion. Using the magnitude Mw 4.3 and a rupture radius of 1 km (2 km length), we obtain a stress drop of Δσ = (7/16) Mo/r3 = 1.5 MPa. This is near the lower limit but still within the expected range of 1–40 MPa for intraplate events (Allmann & Shearer 2009). The Maranhão event seems to have a lower stress drop than expected for the average intraplate earthquake. DISCUSSION AND CONCLUSIONS This earthquake in Maranhão state provided a rare chance to sample the stress field in an aseismic area near the equatorial coast of Brazil. Despite the large distances of the regional stations, a moment-tensor inversion using path-specific models (Dias et al. 2016) was able to constrain the focal mechanism as a strike-slip event with a coast-parallel P axis and a coast-perpendicular T axis. The highly correlated P-wave signals from tens of aftershocks allowed the estimate of backazimuths with enough resolution to provide reliable relative epicentres showing the rupture to be about 2 km long with a roughly N–S orientation, in agreement with one of the nodal planes of the focal mechanism solution (Fig. 1b). Strike-slip mechanisms with coast-parallel P axes are the most frequent type of faulting in NE Brazil (further to the east, Fig. 1b) as shown by Ferreira et al. (1998) and Assumpção et al. (2016). The lateral density variations across the continent/oceanic lithosphere cause spreading stresses in the continent, that is, coast-perpendicular extension (as modelled by Stein et al. 1989), and as observed in Brazil (Assumpção 1998; Oliveira et al. 2015; 2016) and many other continental margins around the world (Tingay et al. 2006). SHmax measurements from borehole breakouts in other parts of the northern coast (Lima et al. 1997) also show the same ESE–WNW, coast-parallel trend (Fig. 1b). This seems to indicate that the whole Brazilian equatorial margin is part of a large, uniform stress province. Worldwide compilation of stress measurements (Heidbach et al. 2010) has shown that first-order stress provinces (up to 2000 km long) are found in several intraplate areas (e.g. in North America and NE Asia). These have been explained by the prevalence of far-field plate-boundary forces (e.g. Tingay et al. 2006; Heidbach et al. 2010). Shorter-wavelength provinces (the so-called second-order patterns with ∼100 km dimensions) can be caused by local forces due to large structural lithospheric contrasts (Tingay et al. 2006), and flexural effects. Near the Amazon river fan, for example, Watts et al. (2009) showed that the load of the sediments in the continental shelf should cause high flexural stresses with a coast-perpendicular extension onshore. Sediment load and continent/ocean transition are usually considered to have a ‘secondary’ effect on the stress field, compared to plate-boundary forces. However, the long extent of the NW–SE oriented Brazilian equatorial margin may cause an apparent first-order, ∼2000 km long stress province, which must be taken into account when modelling the stresses resulting from far-field plate boundary forces. The Maranhão earthquake was related to stresses commonly found in passive continental margins, far from present plate boundary processes. For this reason, it is a typical Stable Continental Region (SCR) earthquake, as defined by Johnston (1989). Because plate-boundary forces (such as Mid-Atlantic ridge-push and collisional forces at the convergence between Nazca and South American plates) can cause stresses across an entire plate (e.g. Heidbach et al. 2010), SCR earthquakes have been thought to follow the same strain loading/release cycles as plate-boundary events, the only difference being the much lower strain rates and longer recurrence intervals. However, an alternative hypothesis has been proposed (e.g. Stein & Liu 2009; Liu et al. 2010; Calais et al. 2016): intraplate stresses are near critical and roughly constant, with earthquakes occurring by localized small transient variations of stress or fault strength. There would be no significant slow strain accumulation with a long-term seismic cycle. Earthquakes along the entire Brazilian northern coast are subject to typical continental margin stresses, mainly due to lateral density contrasts across continental-oceanic crust, which are ‘non-renewable’ and do not follow seismic cycles. If these stresses predominate over long-range stresses from plate boundary forces, then the Maranhão earthquake and the seismicity along the northern Brazilian coast would favour the new hypothesis of a constant pre-stressed lithosphere. ACKNOWLEDGEMENTS We thank Petrobras for supporting the installation of the Brazilian Seismic Network (RSBR), and CPRM (Brazilian Geological Survey) for the present operation and maintenance of the network. Adriano Botelho and Marcelo Rocha helped with the installation of the local stations for the aftershock study. This work was carried out with the support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) grant 30.6547/2013-9 and CPRM grant 23809. RSBR data are freely available at www.sismo.iag.usp.br, (last accessed 2018 January 14). We thank Stephen Hicks and an anonymous reviewer for the many helpful comments and suggestions to improve the paper. REFERENCES Allmann B.P., Shearer P.M., 2009. Global variations of stress drop for moderate to large earthquakes, J. geophys. Res.  114 B01310, doi:10.1029/2008JB005821. https://doi.org/10.1029/2008JB005821 Google Scholar CrossRef Search ADS   Assumpção M., 1983. A regional magnitude scale for Brazil, Bull. seism. Soc. Am.  73 237– 246. Assumpção M., 1998. Seismicity and stresses in the Brazilian passive margin, Bull. seism. Soc. Am.  88( 1), 160– 169. Assumpção M., Feng M., Tassara A., Julià J., 2013. Models of crustal thickness for South America from seismic refraction, receiver functions and surface wave tomography, Tectonophysics  609 82– 96. https://doi.org/10.1016/j.tecto.2012.11.014 Google Scholar CrossRef Search ADS   Assumpção M., Ferreira J., Barros L., Bezerra H., França G., Barbosa J., Dourado J., 2014. Intraplate seismicity in Brazil, in Intraplate Earthquakes  . 50– 71, ed. Talwani P., Cambridge Univ. Press, Cambridge, doi:10.1017/CBO9781139628921.004. Google Scholar CrossRef Search ADS   Assumpção M., Bianchi M., Albuquerque D.F., França G.S.L., Barros L.V., 2015. Crustal thickness map in South America: an updated version, in 1st Brazilian Symposium of Seismology  Brasília, Dezembro 2015, Expanded Abstract. Assumpção M., Dias F.L., Zevallos I., Naliboff J.B., 2016. Intraplate stress field in South America from earthquake focal mechanisms, J. South Am. Earth Sci.  71 278– 295. https://doi.org/10.1016/j.jsames.2016.07.005 Google Scholar CrossRef Search ADS   Barros L.V., Assumpção M., 2011. Basement depths in the Parecis Basin (Amazon) with receiver functions from small local earthquakes in the Porto dos Gaúchos Seismic Zone, J. South Am. Earth Sci.  32 142– 151. Google Scholar CrossRef Search ADS   Bianchi M. et al.   2015. The Brazilian Seismographic Network: Historical overview and current status, Summary. Bull. Int. Seismol. Cent.  49( 1–6), 70– 90, http://www.isc.ac.uk/iscbulletin/summary/. Bobrov D.I., Kitov I.O., Rozhkov M.V., Friberg P., 2016. Towards global seismic monitoring of underground nuclear explosions using waveform cross correlation. Part I: Grand master events, Seism. Instrum.  52( 1), 43– 59. https://doi.org/10.3103/S0747923916010035 Google Scholar CrossRef Search ADS   Calais E., Camelbeeck T., Stein S., Liu M., Craig T.J., 2016. A new paradigm for large earthquakes in stable continental plate interiors, Geophys. Res. Lett.  43 10 621–10 637. Google Scholar CrossRef Search ADS   Carvalho J., Barros L.V., Zahradnik J. 2016. Focal mechanisms and moment magnitudes of micro-earthquakes in central Brazil by waveform inversion with quality assessment and inference of the local stress field, J. South Am. Earth Sci.  71 333– 343. https://doi.org/10.1016/j.jsames.2015.07.020 Google Scholar CrossRef Search ADS   Cormier V.F., 1982. The effect of attenuation on seismic body waves, Bull. seism. Soc. Am.  72( 6B), S169– S200. Dias F.L., Zahradnik J., Assumpção M., 2016. Path-specific, dispersion-based velocity models and moment tensors of moderate events recorded at few distant stations: examples from Brazil and Greece, J. South Am. Earth Sci.  71 344– 358. https://doi.org/10.1016/j.jsames.2016.07.004 Google Scholar CrossRef Search ADS   Ferreira J.M., Oliveira R.T., Takeya M.K., Assumpção M., 1998. Superposition of local and regional stresses in northeast Brazil: evidence from focal mechanisms around the Potiguar marginal basin, Geophys. J. Int.  134 2, 341– 355. https://doi.org/10.1046/j.1365-246x.1998.00563.x Google Scholar CrossRef Search ADS   Ford S.R., Walter W.R., Dreger D.S., 2012. Event discrimination using regional moment tensors with teleseismic-P constraints, Bull. seism. Soc. Am.  102( 2), 867– 872. https://doi.org/10.1785/0120110227 Google Scholar CrossRef Search ADS   Heidbach O., Tingay M., Barth A., Reinecker J., Kurfes D., Müller B., 2010. Global crustal stress pattern based on the World Stress Map database release 2008, Tectonophysics ., 482( 1–4), 3– 15. https://doi.org/10.1016/j.tecto.2009.07.023 Google Scholar CrossRef Search ADS   Herrmann R.B., 2013. Computer Programs in Seismology: An evolving tool for instruction and research, Seismol. Res. Lett.  84( 6), 1081– 1088. https://doi.org/10.1785/0220110096 Google Scholar CrossRef Search ADS   Hosseini K., 2016. Global multiple-frequency seismic tomography using teleseismic and core-diffracted body waves, Doctoral thesis , Ludwig Maximilians University, http://nbn-resolving.de/urn:nbn:de:bvb:19-195973. Hudson J.A. ( 1969). A quantitative evaluation of seismic signals at teleseismic distances—II: Body waves and surface waves from an extended source, Geophys. J. Int.  18( 4), 353– 370. https://doi.org/10.1111/j.1365-246X.1969.tb03574.x Google Scholar CrossRef Search ADS   Johnston A.C., 1989. The seismicity of ‘Stable Continental Interiors’, in Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound  . 299– 327, eds Gregersen S., Basham P.W., Springer. Google Scholar CrossRef Search ADS   Juliá J., Ammon C.J., Herrmann R.B., Correig A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, Geophys. J. Int.  143 99– 112. Google Scholar CrossRef Search ADS   Kagan Y., 1991. 3-D rotation of double-couple earthquake sources, Geophys. Int ., 106( 3), 709– 716. https://doi.org/10.1111/j.1365-246X.1991.tb06343.x 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( 1), 108– 124. https://doi.org/10.1111/j.1365-246X.1995.tb03540.x Google Scholar CrossRef Search ADS   Lima C., Nascimento E., Assumpção M., 1997. Stress orientations in Brazilian sedimentary basins from breakout analysis: implications for force models in the South American plate, Geophys. J. Int.  130( 1), 112– 124. https://doi.org/10.1111/j.1365-246X.1997.tb00991.x Google Scholar CrossRef Search ADS   Lithgow-Bertelloni C., Guynn J.H., 2004. Origin of the lithospheric stress field, J. geophys. Res.  109( B1), B01408, doi:10.1029/2003JB002467. https://doi.org/10.1029/2003JB002467 Google Scholar CrossRef Search ADS   Liu M., Stein S., Wang H., 2011. 2000 years of migrating earthquakes in North China: how earthquakes in midcontinents differ from those at plate boundaries, Lithosphere  3( 2), 128– 132. https://doi.org/10.1130/L129.1 Google Scholar CrossRef Search ADS   Marco E., Barcelos A., Assumpção M., 2010. Magnitude determination of regional earthquakes in Brazil: a comparison with the new IASPEI magnitude formulas, in AGU Meeting of the Americas  Iguazu Falls, Brazil, Abstract. Naliboff J.B., Lithgow-Bertelonni C., Ruff L.J., Koker N., 2012. The effects of lithospheric thickness and density structure on Earth's stress field, Geophys. J. Int.  188 1, 1– 17. https://doi.org/10.1111/j.1365-246X.2011.05248.x Google Scholar CrossRef Search ADS   Nuttli O.W., 1983. Average seismic source-parameter relations for mid-plate earthquakes, Bull. seism. Soc. Am.  73( 2), 519– 535. Oliveira P.H.S., Ferreira J.M., Bezerra F.H.R., Assumpção M., Nascimento A.F., Sousa M.O.L., Menezes E.A.S., 2015. Influence of the continental margin on the stress field and seismicity in the intraplate Acaraú seismic zone, NE Brazil, Geophys. J. Int.  202( 3), 1453– 1462. https://doi.org/10.1093/gji/ggv211 Google Scholar CrossRef Search ADS   Sokos E.N., Zahradnik J., 2008. ISOLA a Fortran code and a Matlab GUI to perform multiple-point source inversion of seismic data, Comput. Geosci.  34( 8), 967– 977. https://doi.org/10.1016/j.cageo.2007.07.005 Google Scholar CrossRef Search ADS   Sokos E.N., Zahradnik J., 2013. Evaluating centroid-moment-tensor uncertainty in the new version of ISOLA software, Seismol. Res. Lett.  84( 4), 656– 665. https://doi.org/10.1785/0220130002 Google Scholar CrossRef Search ADS   Stammler K. 1993. Seismichandler—programmable multichannel data handler for interactive and automatic processing of seismological analyses, Comput. Geosci.  19 2, 135– 140. https://doi.org/10.1016/0098-3004(93)90110-Q Google Scholar CrossRef Search ADS   Stein S., Liu M., 2009. Long aftershock sequences within continents and implications for earthquake hazard assessment, Nat. Geosci.  462( 7269), 87– 89. Stein S., Cloetingh S., Sleep N.H., Wortel R., ( 1989). Passive margin earthquakes, stresses and rheology, in Earthquakes at North-Atlantic Passive Margins: Neotectonics and Postglacial Rebound , pp. 231– 259, eds Gregersen S., Basham P.W., Kluwer Academic. Google Scholar CrossRef Search ADS   Tingay M.R.P, Müller B., Reinecker J., Heidbach O., 2006. State and origin of the present-day stress field in sedimentary basins: new results from the World Stress Map project, paper presented at the 41st U.S. Symposium on Rock Mechanics (USRMS), 06–1049, Golden, CO. Watts A.B., Rodger M., Peirce C., Greenroyd C.J., Hobbs R.W., 2009. Seismic structure, gravity anomalies, and flexure of the Amazon continental margin, NE Brazil, J. geophys. Res.  114 B07103, doi:10.1029/2008JB006259. https://doi.org/10.1029/2008JB006259 Google Scholar CrossRef Search ADS   Wells D.L., Coppersmith K.J., 1994. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, Bull. seism. Soc. Am.  84 974– 1002. Zahradnik J., Custódio S., 2012. Moment tensor resolvability: application to southwest Iberia, Bull. seism. Soc. Am.  102 1235– 1254. https://doi.org/10.1785/0120110216 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Grid search for centroid time shift and depth (‘source position’). The depth varies from 2.5 to 17 km and time shift from −5.0 and 5.0 s. The contours indicate the VR*PF index. The preferred solution is the mechanism with highest PF*VR value, highlighted with the largest beachball in the grid. For the uncertainty, we use a threshold of 0.9 for VR and PF: all solutions with VR*PF ≥ 0.9*(VR*PF)MAX (subscript MAX indicates the highest PF*VR) were considered and are shown as larger beachballs in the middle of the grid. Figure S2. Sensitivity to filter bands. Comparison between the moment tensor parameters from different frequency ranges: 0.07–0.10 (grey dots) and 0.07–0.40 Hz (black dots). Each dot represents an acceptable solution, that is, VR*PF ≥ 0.9*(VR*PF)MAX. Top-left histogram is depth, top-right represents the centroid time shift (CT), bottom-left shows double-couple percentage (DC) and bottom-right moment magnitude (Mw). Figure S3. Map of the temporary local stations (solid triangles, with numbers depicting stations VGM1 to VGM6) and the permanent station ROSB used to locate the aftershock of 2017 February 2 (yellow star, event 24 in Table S1). Solid lines are faults from the CPRM (Brazilian Geological Survey) database. Solid lines in the inset show the major geological provinces as in Fig. 1 of the main text. Figure S4. Seismograms of the 2017 February 2 aftershock (magnitude 0.6) at the three local stations VGM2, 4, 6. Sp and Ps are converted phases at the base of a ∼0.5 km thick sedimentary layer. Each station shows the vertical (Z), radial (R) and transverse (T) components. The Ps phase is best seen in the radial component; the Sp precursor in the vertical component. The time differences are roughly the same: Ps-P = S-Sp. See Barros & Assumpção (2011) for similar observations. Figure S5. Alignment of P arrival (left: 0.3 s before to 1.0 s after the P) and S arrival (right: 0.4 s before to 0.6 after), both at the vertical component for the first 15 aftershocks in Table S1. Bandpass filtered 2–8 Hz. Marks denoted by ‘S’ are preliminary attempts to identify the absolute arrivals (less accurate) picked before the correlation. The correlated arrivals for the S waves are denoted by ‘T1’. Figure S6. Vertical component traces aligned by the correlated S arrival (mark ‘T1’), showing significant differences in the S–P times. Events are 1 (main shock) to 15 from top to bottom. The two more distant events (larger S–P times) are nos. 8 and 15. No. 14 is the closest event. Figure S7. Backazimuths, dip angles, S–P and epicentres, calculated with 2–8 Hz bandpass filter and a time window of −0.05 to 0.35 s around the P onset. Small circles denote magnitudes <1 mR. The grey star is the 4.3 Mw main shock. Top left: S–P × backazimuth; top right: S–P × dip angle; bottom: epicentres determined with the BAz and the S–P derived distance (all events assumed at the same fixed depth). Except for three events (probably off the main fault), the main shock and 11 aftershocks show good correlation between S–P, BAz and Dip angle. The tight cluster of 12 events with roughly N–S orientation indicates the probable rupture length. Table S1. Main shock and aftershocks used in the epicentre determination. Two backazimuths (Baz) are shown: BAz-2–8 were measured with the 2–8 Hz bandpass filter and window −0.05 s before to 0.35 s after the P arrival; Baz-M is the medium value of nine different filter/window combinations. Magnitudes mR of the aftershocks were determined by comparing the P-wave rms amplitudes with that of the main shock. Table S2. Result of hypocentral determination of the 2017-Feb-02 07:15 aftershock (No. 24 in Table S1) using a best-fitting local model with Vp/Vs = 1.73 (from a Wadati diagram). The sedimentary layer (∼0.5 km thick, Vp = 3.0 km s−1, Vp/Vs = 2.0), similar to values of the Parecis Basin (Barros & Assumpção 2011), was replaced by station corrections (column ‘Delay’). We used the HYPO71 code modified for millisecond accuracy. Stations VGM1,5 had problems with GPS clock and only S–P times were used. 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. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.

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
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

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.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off