On the source of the 8 May 1939 Azores earthquake – tsunami observations and numerical modelling

ABSTRACT On 8 May 1939, an earthquake (Ms7.1) occurred near the Azores archipelago, with an epicentre located close to the western end of the Gloria fault. Previous studies present different epicentre locations spreading over a large area, and two different types of focal mechanisms. Given these uncertainties, the interpretation of the seismological information in a complex tectonic environment between the Gloria Fault and the Terceira Ridge is a matter of debate. The event caused a small tsunami recorded in the Azores Islands. In this study, we use the tsunami observations and tsunami numerical modelling to select the earthquake fault rupture that best fits the tsunami observations. We consider the different focal mechanism solutions, perform tsunami numerical modelling, and compute synthetic tsunami waveforms at the tide gauge locations. We find that an earthquake caused by a low-angle dipping fault with dominant strike–slip movement generates a tsunami that reproduces well the record at Ponta Delgada tide gauge. Finally, in areas where earthquakes are rare, the study of ancient earthquakes must use all information available, namely tsunami observations and mareograph data.


Introduction
The Azores Islands constitute an archipelago of volcanic origin located near the triple junction between Eurasian, African, and North American plates. The area has recurrent seismic activity manifested in a significant number of micro-earthquakes and a few moderate and big events (Andrade et al. 2006;Bezzeghoud et al. 2008). The Islands chain location, along the Mid-Atlantic Ridge, presents a suite of inter-related morphological, seismic, and volcanic features that favour tsunami generation and impact (Andrade et al. 2006). Very little is known about the tsunamigenic history of the region (Andrade et al. 2006). The historical data on earthquakes and tsunamis are available in the form of written accounts and geologic evidence. Since the settlement of the archipelago in the fifteenth century, several tsunamis struck Azorean coastal zones, and most of them were generated by earthquakes (Andrade et al. 2006).
The location of the archipelago and the frequency of the seismic events justified the implementation of both seismic and tide gauge networks at the Azores archipelago. Since the installation of the tide gauge network in 1924 in the Azores archipelago (GLOSS 2016), the sensors in the Azores CONTACT C. Reis claudiavdreis@sapo.pt Some documents also include tsunami observations. According to Moreira (1968), a small tsunami followed the 1939 earthquake. The tide stations of Ponta Delgada and Angra do Heroismo, respectively, in São Miguel and Terceira Islands (see Figure 2), recorded the tsunami event.
At Ponta Delgada, the maximum peak-to-peak amplitude was 16 cm. The first movement in the tide gauge record was a downward movement at 02:05 am (UTC). The record of Angra do Heroismo, in Terceira Island, is unreadable. The newspaper 'Correio dos Açores' quoting information from the Azores Meteorological Service reports three tsunami waves of small amplitude, recorded in Angra do Heroismo. Baptista and Miranda (2009) compiled the Portuguese catalogue of tsunamis  Udias et al. (1976) and by Buforn et al. (1988) for the 1939 event. For both focal mechanisms, letters 'A' and 'B' correspond to nodal planes A and B, respectively. Stars are the epicentral locations proposed by several authors (complementary information in Table 1), the best constrained locations are numbered (4, 5 and 7). and reported that the 1939 earthquake triggered a tsunami recorded by Ponta Delgada and Angra do Heroismo tide gauges.
In this study, we investigate which of the proposed fault planes is more compatible with the tsunami observations. To do this, we digitized the record of Ponta Delgada and de-tided it to compare with the synthetic waveform (see Figure 3). We use the focal mechanisms published for the 1939 event to compute sea-bottom deformations using the half-space elastic deformation theory (Mansinha & Smiley 1971;Okada 1985) implemented in Mirone suite (Luis 2007). Assuming that the initial sea-surface displacement mimics the sea-bottom deformation, we use them as initial conditions to run the tsunami simulations. We employ a validated nonlinear shallow water code À NSWING (Miranda et al. 2014). We computed synthetic tsunami waveforms at grid points close to Ponta Delgada and Angra do Heroismo. The grid points for waveform computation are the points closest to the locations of the real mareographs whose water depth is not less than 10 m. We compared the synthetic tsunami waveforms with the observed data. Based on this comparison we select the fault model that best matches the tsunami records.

Geological context and tsunamigenic sources
The Azores Islands are located in the North Atlantic Ocean, between 37 N and 40 N latitude and 25 W and 31 W longitude, close to the triple junction of the North America, Eurasia, and Africa plates ( Figure 1). The Azores archipelago is subdivided into three island groups. The western group (the Flores and Corvo Islands) is inserted on the North American plate. The central group (Terceira, Graciosa, São Jorge, Pico and Faial Islands) and eastern group (São Miguel and Santa Maria Islands and the Formigas Islets) are both on the EurasiaÀAfrican plate junction (Figure 2(a)).
The central and eastern groups stand in a region that absorbs the deformation induced by the eastward differential displacement of the Eurasian and African plates, east of the Azores triple junction, morphologically expressed by the Azores Plateau (Matias et al. 2007, and references therein).
The Azores Plateau is a shallow (less than 2000 m) triangular topographic zone. It is located where the Terceira Rift, a NWÀSE to WNWÀESE fault system with a normal dextral component, intersects the Mid-Atlantic Ridge, with an approximate NÀS direction and a general triangular shape, as illustrated in Figure 1 (Matias et al. 2007. and references therein). The Azores Plateau is surrounded by the nearby abyssal plains, with depths over 3500 m.
The Azores-Gibraltar Fracture Zone is less prominent, with an extension from Azores to the Strait of Gibraltar, delimiting the EurasianÀAfrican plate boundary. The plate boundary is supposed to follow a large linear feature, denominated Gloria Fault (GF) between 24 W and 17 W, characterized by a normal sinistral displacement in the NNWÀSSE trend (Matias et al. 2007, and references therein). The GF is represented in Figure 1.
The volcanic nature of the Islands and the active tectonic structures in the area contribute to the archipelago seismicity. The present seismic activity recorded by the seismic network is characterized by a large number of small-magnitude events (magnitude below 3). Moderate-and strong-magnitude earthquakes are frequent in the area (Bezzeghoud et al. 2014). Recently, the seismic network recorded several events of magnitude close to 6 (e.g. the most recent are the 27 June 1997, 9 July 1998, 5 April 2007, and 7 April 2007. Since the installation of the seismic network, a few important events struck the Azores region with magnitude 7 or larger: the 20 May 1931, the 8 May 1939, the 25 November 1941, and the 1 January 1980. The tsunami vulnerability of the area results from (1) its location in the Atlantic, which makes the islands vulnerable to far-field events; and (2) their geotectonic setting that includes different sources of tsunami triggers, namely earthquakes, volcanoes, submarine landslides, and aerial slope mass movements (Latter 1981;Beg et 2000;De Lange et al. 2001;Freundt 2003;Dawson et al. 2004;Andrade et al. 2006).
According to the available documents, some tsunami impacted the Azorean coasts since the settling of the islands in the fifteenth century. Most of these events are of tectonic origin. About half of them are distal source events (Baptista &Miranda 2009).

Epicentre and focal mechanism
The 8 May 1939 earthquake, Ms7.1, occurred close to the western end of the GF. Many epicentre locations have been proposed for this event. These locations are spread over a vast area that precludes a simple tectonic interpretation (Figure 2(b), Table 1). In the following discussion, we will give preference to the epicentre locations that have been recently revised, EHB, ISC-GEM, and IPMA (as in Table 1). Udias et al. (1976) presented the first focal mechanism for the 1939 earthquake (see Table 2). These authors used 26 P-wave polarities and manually fit two nodal planes so that only one polarity was wrong. They also took into consideration that 4 readings (out of 26) were evaluated as doubtful. The focal mechanism proposed consists of a nearly pure strikeÀslip solution, with one of the planes following the GF trend.
Later on, Buforn et al. (1988) revised the event (see Table 2). These authors added one polarity reading and revised a few of the waveforms. The P-wave polarity distribution remained very similar but the solution proposed was considerably different. The reason for that was the use of an automatic algorithm to compute the focal mechanism using the maximum-likelihood method proposed by Brillinger et al. (1980). This method compares the polarities observed with the P-wave expected amplitudes and computes them from a likelihood function. The method will tend to locate the nodal planes as far as possible from the observations, even at the price of increasing the number of wrong polarities. The automatic method has the additional advantage of estimating errors for the best-fit solution, assuming a linear error model. Plane A of the Buforn et al.'s (1988) focal mechanism (Table 2) consists of a low-angle dipping fault with a dominant right lateral strikeÀslip movement, while plane B represents a high-angle dipping fault with a dominant dipÀslip normal movement. The overall polarity score is lower than the one given by Udias et al. (1976) with four incorrect polarities for the best-fit solution.
Since Udias et al. (1976) did not publish the values of the rake angles for their focal mechanism, Moreira (1985) presented the same solution, with plane A and plane P angles recomputed from the P and T axis parameters, which were published before. The catalogue of Udias et al. (1989) also reproduces the results of Udias et al. (1976) with the rake angle computed from the two nodal plane strike and slip angles. All these parameters are reported in Table 2.
In summary, with very similar sets of P-wave polarities, we have two very different focal mechanisms proposed. One mechanism results from the best fit of polarities and the other mechanism is obtained by an automatic algorithm using the maximum-likelihood method applied to polarities and amplitudes. (1) Rake angles computed from nodal plane angles.
(2) Rake angles computed from P and T axes by the author.
(3) Rake angles computed from nodal plane angles by the authors.
In the following sections, we will evaluate the contribution of tsunami observations and modelling to differentiate the fault planes proposed for the 1939 earthquake. Both nodal planes from the two focal mechanisms (SSÀSS and SSÀNF in Table 2) will be considered. These two possible focal mechanisms for the 1939 earthquake have been published in earthquake source compilations, like Udias et al. (1989) or Vannucci and Gasperini (2003).

Tsunami observations
Several documents present brief descriptions of the 8 May 1939 tsunami. According to Moreira (1968), a small tsunami followed the 8 May 1939 earthquake recorded at Ponta Delgada, in São Miguel Island (see Figure 2(b) for location) with a 'double-amplitude wave' (i.e. peak-to-peak wave amplitude) of 16 cm and the first descendant movement at 02:05 am (UTC). According to Moreira (1968), the tide record of Angra do Heroismo in Terceira Island was unreadable. The local newspaper, 'Correio dos Açores', quoting information from the Azores Meteorological Service, refers to three waves recorded by the Angra do Heroismo tide gauge at 02:25, 02:29, and 02:45 am (UTC) following the earthquake felt at 01:47 am(UTC). Baptista and Miranda (2009) compiled the Portuguese catalogue of tsunamis in which they mentioned that the 1939 earthquake was responsible for triggering the tsunami recorded at Ponta Delgada and Angra do Heroismo tide gauges.

Methodology and data
We performed tsunami numerical modelling using as initial condition the sea-surface elevation caused by the co-seismic displacement of the earthquake. To do that, we computed the co-seismic displacement of the ocean bottom and transferred it to the water column, assuming that the water is incompressible. The methodology includes (1) compilation of the bathymetric data-set for tsunami simulations; (2) use of scaling laws to infer the fault dimensions (length and width) and the average co-seismic slip; (3) computation of the initial sea-surface elevation; (4) modelling tsunami propagation and extraction of synthetic tsunami waveforms at the grid node closest to the position of the old tide gauge of Ponta Delgada; and (5) comparison of the synthetic waveforms with the observations.

Data: bathymetric data and tide records
The bathymetric model used in this study is generated from the GEBCO (General Bathymetric Chart of the Oceans) 30-arcsec grid and from the 50 m horizontal resolution digital elevation model available at http://w3.ualg.pt/»jluis/mirone/data-links.html. The 30-arcsec grid extends from 36 N to 41 N and from 22 W to 33 W and covers the earthquake source area and the Azores archipelago western, central, and eastern Islands groups. To guarantee a good representation of the bathymetric features in the study area, close to the tide stations, we improve the GEBCO bathymetric model by integrating the 50 m resolution grid data. We built a 200 m resolution final bathymetric model to simulate the tsunami from the source to the coast. This final model was generated from an extrapolation of the 50 m resolution data that encompasses the area close to São Miguel Island and an interpolation of the GEBCO 30 arcsec grid data around the rest of the Azores Islands, including Terceira.
The record of Ponta Delgada was digitized, linear interpolated, and de-tided to isolate the tsunami signal. To remove the tide signal, we used a polynomial fitting of order 9 (Figure 3).

Fault models
To initiate the tsunami simulations we need the initial sea-surface elevation caused by the sea-bottom deformation due to the earthquake. To compute the sea-bottom deformation we need the data presented in Table 2 plus the length, width, and slip along the fault. These parameters are obtained through the scaling laws that we propose for two fault models: (1) strikeÀslip earthquake at the GF and (2) normal faults in the Azores.

StrikeÀslip model at the Gloria Fault
The fault model for a strikeÀslip earthquake at the GF is based on the study of two large earthquakes (M > 7) placed in the fault, the 25 November 1941, Ms8.4 (Udias et al. 1976;Lynnes & Ruff 1985, and re-evaluated by Baptista et al. 2011) and the 26 May 1975, Ms7.9 (Buforn et al. 1988Argus et al. 1989).
The GF is in many ways similar to an Oceanic Fracture Zone (OFZ) since it represents an old and actual plate boundary where the Eurasia and African plates have slid without major compressional or extensional features. However, this is not a common transform fault since it lacks the extensional ridges at both ends. GF also departs from all known OFZ since it offsets very old lithosphere from 50 Ma in the west to 120 Ma in the east, close to the Tore-Madeira Rise (M€ uller et al. 2008). Current plate kinematic models show that the GF trend follows very closely the relative velocity between Nubia and Eurasia that has an estimated value of 5 mm/year (Fernandes et al. 2003).
The seismogenic thickness of the lithosphere in the GF area is inferred from the comparison of plate cooling models and the well-constrained focal depths of oceanic earthquakes (e.g. McKenzie et al. 2005). Considering the ages of the lithosphere involved, we use an average value of 40 km for the maximum rupture thickness in the western end of the GF.
In the re-evaluation of the 25 November 1941 earthquake, by Baptista et al. (2011), the location and focal mechanism were recomputed. The source parameters used for the hydrodynamic modelling were a shear modulus of 3 £ 10 10 Pa, a fault length of 200 km, a fault width of 45 km, and an average slip of 12 m. The inferred moment magnitude is 8.3. We remark that increasing the shear modulus to 5 £ 10 10 Pa would increase the moment magnitude Mw to 8.4, without affecting the tsunami modelling. The stress drop will also increase from 5.6 MPa (m D 3 £ 10 10 Pa) to 9.6 MPa (m D 5 £ 10 10 Pa).
A summary of the main seismotectonic parameters related to the two largest events that occurred during instrumental times in the GF area is presented in Table 3.
Considering two global fault scaling compilations, the one by Wells and Coppersmith (1994) and the one by Stirling et al. (2002), and using a shear modulus of m D 5 £ 10 10 Pa, more appropriate for the oceanic lithosphere, we conclude that both compilations fail to reproduce the characteristics of large earthquakes in the GF: the aspect ratio is larger and variable in the compilations (e.g. 6.0 or 5.0 for Mw7.1) and the stress drop is much smaller (e.g. 2.0 or 2.25 MPa for Mw7.1). These deviant characteristics of large earthquakes can be explained by the fact that these global compilations are dominated by strikeÀslip earthquakes on continental fracture zones.
For the GF strikeÀslip earthquakes, we propose a local fault model based on the semi-empirical scaling law of Manighetti et al. (2007) (equation (1)): Here D max is the maximum displacement on the fault, L is the fault length, and W Sat is the saturation width, which we assimilate to the maximum rupture thickness discussed above. We will also assume that D max D 2D as in Manighetti et al. (2007). The a parameter is proportional to a constant static stress drop. Considering the data of large earthquakes in the GF and other oceanic intra-plate strikeÀslip events, we propose a scaling law between maximum displacement and fault length using the equation of Manighetti et al. (2007), depicted above, with a D 50 £1 0 ¡5 and W Sat D 40 km. We further consider that the aspect ratio of GF earthquakes is 4 for events that do not rupture the whole brittle lithosphere. Using these considerations we infer the fault model in the GF given in Table 4 for Mw7.0 and Mw7.27 earthquakes. The fault static stress drop is computed using the relationships proposed by Kanamori and Anderson (1975) (equation (2)), for strikeÀslip faults:

Normal faults model in the Azores
As before, considering a rupture of oceanic lithosphere, we assume that the shear modulus (m) is 5.0 £ 10 10 Pa. Contrary to the GF area, we have no well-constrained normal fault earthquakes in the Azores archipelago to allow a definition of a local semi-empirical scaling law. For this reason, to obtain the scaling relationship between magnitude and the fault dimensions (and average slip), we will use as a reference the database of great earthquakes with a known extended source compiled by Stirling et al. (2002). From this compilation, we extract the faults with pure strikeÀslip, pure normal, or a mixture of the two regimes, that we consider representative for the Azores Plateau. For this data-set, which includes incomplete information on subsurface rupture length, we extract an average relationship between surface and subsurface lengths. The scaling factor found is 1.65, L Subsurface D 1:65L Superficial . This allows us to estimate the length of the fault as the subsurface rupture for all the earthquakes where the width and average displacement on the fault are known. Using this data, we infer the following empirical relationships between moment magnitude and fault length (L) (equation (3)) and between moment magnitude and fault width (W) (equation (4)): log.L/ D 0:698Mw ¡ 3:049 (3) log.W/ D 0:157Mw C 0:044 These empirical relationships are given in Table 4 for two magnitudes, Mw D 7.0 and Mw D 7.25.

Tsunami numerical modelling
Tsunami numerical simulations, presented here, include the generation and propagation, and the waveform computation closest to the tide gauge location where the recorded signal is available (i.e. Ponta Delgada tide gauge). We assume that the earthquake rupture is instantaneous, and we compute the seabed deformation using the half-space elastic theory first introduced by Mansinha and Smiley (1971) and later extended by Okada (1985). Assuming that the water is incompressible, the vertical sea-bottom displacement is then transferred to the free ocean surface. The tsunami propagation is modelled using the in-house developed and benchmarked NSWING (Non-Linear Shallow Water model With Nested Grids) code (Miranda et al. 2014). NSWING is a finite difference tsunami code that solves the nonlinear shallow water equations, using the second- Table 4. Scaling laws for the length (L), width (W), and average slip (u), assuming m D 5 £ 10 10 Pa, for (1) the strikeÀslip models in the Gloria Fault and (2) for the SS and NF mechanism codes considered in Table 2 order discretization and a staggered explicit leap-frog scheme, in Cartesian or spherical coordinates, like in COMCOT code (Liu et al. 1998). The code incorporates an open radiative boundary condition, Coriolis acceleration, bottom friction, a moving boundary scheme to model run-up, and the possibility to use multiple different levels of nested grids. An open radiative boundary condition is used on the outward limit of the external grid, whenever it does not correspond to dry land. All tsunami simulations were computed for mean sea level condition. We used one single grid (200 m resolution, see Section 4.1) to simulate the tsunami from the source to the target points. We ran the numerical code for a total time of 4 h with a one-second time step to meet the CourantÀFriedrichs-Lewy numerical stability criterion.

Results
For each nodal plane of the two representative focal mechanisms described in Section 3.1 and presented in Table 2, we define a fault model based on the scaling laws introduced in Section 4.2 that is used to generate a tsunami that is propagated to the target points in the Azores Islands. In the following sections, we present the results in terms of tsunami initial deformations, tsunami maximum wave height distributions, and synthetic waveforms close to Ponta Delgada. We then compare the simulated tsunami waveforms with the tsunami observations. A summary of the parameters used in the tsunami modelling is given in Table 5. All simulations are done considering that the fault ruptures up to the sea bottom.

SSÀSS focal mechanism
This focal mechanism, proposed by Udias et al. (1976) with rake computed by Moreira (1985), comprises two fault planes, both nearly pure strikeÀslip (Table 2 and Figure 4(a) and 4(b), representing plane A and plane B, respectively).  Figure 4. Initial sea-surface deformation using Udias et al. (1976) focal mechanism and the fault parameters listed in Table 5, associated to (a) the nodal plane A and (b) the nodal plane B.

Plane A À SS
This nodal plane follows the trend of the GF on its western boundary (close to EÀW). Figures 4(a) and 5(a) and 5(b) depict the tsunami generation, the maximum wave amplitude distribution, and the comparison between the observed and the synthetic waveforms, respectively.
The analysis of these results shows that a tsunami wave of about 0.46 m was generated at the source area (Figure 4(a)) and propagated towards the Azores surrounding Islands (Figure 5(a)). Most tsunami energy is steered towards the São Miguel and Santa Maria Islands (Figure 5(a)). The tsunami reaches all the islands in the Azores with amplitudes between 0.02 m at Flores and 0.17 m at São Miguel ( Figure 5(a)). The comparison of the recorded and simulated tsunami at Ponta Delgada tide gauge ( Figure 5(b)) shows a good agreement between both signals in terms of arrival time and the first wave amplitude. However, in terms of signal frequency and number of waves, the signals are very different. The synthetic tsunami suffers a significant attenuation after the arrival of the first wave.  (Figure 4(a)). Respective maximum tsunami wave amplitude simulated in each island (in metres); (b) comparison between the recorded and simulated tsunami waveforms at Ponta Delgada tide gauge.

Plane B À SS
This nodal plane is also a nearly pure strikeÀslip with an azimuth close to NÀS. Figures 4(b) and 6(a) and 6(b) show for this fault plane the tsunami generation, the maximum wave amplitudes distribution, and a comparison between the recorded and simulated tsunami signals at Ponta Delgada. Figure 4(b) shows that a wave of 0.55 m (maximum) is generated close to the fault trace. The maximum wave amplitude distribution map (Figure 6(a)) indicates that most tsunami energy is steered eastwards from the source area. São Miguel and Santa Maria are the islands most affected by the tsunami with waves of about 0.17 and 0.21 m in amplitudes, respectively (Figure 6(a)).
The results obtained for planes A and B for the first arrival time and first wave amplitudes are similar, but observations of the number of waves and the signal period are different (see Figure 6(b)).  (Figure 4(b)). Respective maximum tsunami wave amplitude simulated in each island (in metres); (b) comparison between the recorded and simulated tsunami waveforms at Ponta Delgada tide gauge.

SSÀNF focal mechanism
This is the focal mechanism proposed by Buforn et al. (1988) applying an automatic procedure to a set of P-wave polarities similar to the ones of Udias et al. (1976). In Figure 7(a) and 7(b), plane A and plane B, are depicted, respectively, together with the tsunami generation.

Plane A À SS
This nodal plane represents a low-angle dipping fault with a predominant strikeÀslip motion. The tsunami modelling results including the tsunami generation, maximum wave amplitudes, and comparison between simulated and recorded tsunami signals at Ponta Delgada tide gauge, are presented in Figures 7(a) and 8(a) and 8(b). Figure 7(a) shows that the earthquake rupture causes a maximum vertical sea-surface displacement of 0.15 m. Figure 8(a) depicts the maximum wave amplitudes distribution with values ranging from 0.03 m at Flores to about 0.35 m at Santa Maria. The comparison of the recorded and simulated tsunami at Ponta Delgada tide gauge (Figure 8(b)) shows good agreement in terms of arrival time and signal 'shape'. The simulated signal fits relatively well with the registered one except for the first wave that is slightly higher than the recorded one (Figure 8(b)). Also, the attenuation of the simulated tsunami signal compares well with the observed one.

Plane B À NF
This nodal plane represents a high-angle dipping fault with a predominant normal dipÀslip motion. Figure 7(b) shows that a wave of maximum 0.91 m in amplitude was generated at the source. Maximum wave amplitudes distribution (Figure 9(a)) clearly shows that the largest tsunami impact occurred at Santa Maria and São Miguel Islands with waves of about 0.9 and 0.5 m in amplitude, respectively. In contrast to the models for the SSÀSS focal mechanism (Section 5.1), the tsunami signal generated by this normal fault shows no significant attenuation after the arrival of the first wave (Figure 9(b)). The comparison of the synthetic and observed de-tided records shows a good agreement in the tsunami arrival time and also in the 'shape' of the signals (Figure 9(b)). However, the simulated signal has relatively higher amplitudes than the recorded one (Figure 9(b)).

Discussion and proposed fault model
The tsunami modelling results for the SSÀSS focal mechanism (Udias et al. 1976) showed similar results for both nodal planes, namely a relatively good agreement with recorded signal in terms of arrival time, polarity, and first wave amplitude. On the other hand, the synthetic waveforms show differences regarding the number of significant waves. The waves that follow the first pulse on the synthetic waveform are considerably attenuated, with one peak for plane A and two peaks for plane Figure 7. Initial sea-surface deformation using Buforn et al. (1988) focal mechanism and fault parameters listed in Table 5, associated to (a) the nodal plane A and (b) the nodal plane B.
B, against the seven peaks presented in the Ponta Delgada mareogram. We suppose that this disagreement can be associated with (1) the surrounding sea-bottom morphology (the source area is relatively shallow when compared with the remaining area, as represented in Figure 2), or to (2) the focal mechanism. To understand the possible influence of each one, we simulate an SS-type rupture occurring in a shallower bathymetry. The obtained results show that if the fault associated with the SS solution is positioned in a shallow bathymetry, the attenuation effect is not noticeable. Therefore, the bathymetry is the main factor influencing the attenuation (number of important wave peaks). However, we have to consider that, despite the complex tectonics of this area (intersection of the GF with the Terceira Rift), pure strikeÀslip mechanisms are not found in the shallow area of the Azores Plateau (e.g. Bezzeghoud et al. 2014). Furthermore, the needed displacement of the fault rupture renders the fault planes incompatible with the revised epicentres for the 1939 earthquake.  (Figure 7(a)). Respective maximum tsunami wave amplitude simulated in each island (in metres); (b) comparison between the recorded and simulated tsunami waveforms at Ponta Delgada tide gauge.
Regarding the SSÀNF focal mechanism (Buforn et al. 1988), the arrival time is well estimated when we consider the plane A (low-angle strikeÀslip fault). Also, the polarity of synthetic waveforms is compatible with the observations. However, the simulated amplitudes exceed approximately twice those observed by the tide gauge, and the frequency of the synthetic waveform is higher than the recorded one. This mismatch in frequencies between the simulated and recorded signals may be due to the relatively small fault width obtained, as inferred from the scaling laws, as well as to the low-dip angle value from the proposed focal mechanism.
The comparison between the observed waveform at Ponta Delgada and the synthetic tsunami waveforms for all nodal planes of the two focal mechanisms as proposed by Udias et al. (1976) and Buforn et al. (1988), with the corresponding scaling laws, show that none of these solutions reproduce the original mareogram.  (Figure 7(b)). Respective maximum tsunami wave amplitude simulated in each island (in metres); (b) comparison between the recorded and simulated tsunami waveforms at Ponta Delgada tide gauge.
In summary, two types of discrepancies were found: (1) the attenuation observed for the nodal planes of the SSÀSS focal mechanism, and (2) the high frequency and amplitude for the nodal planes of the SSÀNF focal mechanism. To overcome the first discrepancy, it is required to displace the fault to the shallow bathymetry area to the SW area, which is incompatible with the epicentre area and less likely in the regional geodynamic framework. The second discrepancy can be solved by decreasing the dip of the fault and increasing the width, as discussed in the following section.
The analysis of the tsunami results and comparison with the tsunami signal at Ponta Delgada tide gauge and Terceira enabled us to select the plane A of the focal mechanism, proposed by Buforn et al. (1988) as the possible rupture solution that best approximates the observed tsunami. This conclusion was based on five main criteria. They were (1) tsunami amplitude of the waveform; (2) tsunami signal polarity of the first wave; (3) tsunami arrival time; (4) number of significant waves: the number of initial waves with similar peak to peak values; and (5) geological evidence.
Taking the SSÀNF mechanism, nodal plane A and the earthquake magnitude as initial conditions for this analysis, we adjusted the fault model (length, width, and slip) to better simulate the recorded signal. We have seen above that the simulated tsunami waveform at Ponta Delgada shows a signal with relatively shorter periods when compared to the recorded one (Figure 8(b)). By trial and error, we concluded that a larger fault width with a slightly lower dip angle (25 ) overcomes the difference in periods between the simulated and recorded signals. We ran multiple tsunami simulations using different values of fault length, width, and co-seismic slip. We found that a fault model of dimensions (L £ W) 40.5 km £ 40 km and with a co-seismic slip of 1.5 m can generate a tsunami that best fits the recorded tsunami signal. The fault dimensions and parameters correspond to a seismic moment of 1.2 £ 10 20 Nm, which implies a magnitude Mw7.3. The fault model parameters used for the best fit are provided in Table 5.
In Figures 10 and 11, we present the results of the numerical modelling for the selected and adjusted fault model. The simulation of the initial sea-surface displacement (tsunami generation) ( Figure 10) indicates that the maximum wave height in the source area reaches 0.25 m. Among the Azores Islands, Santa Maria and São Miguel are most impacted by the resulting tsunami (Figure 11(a)).  The Angra do Heroismo tide gauge, although considered unreadable, was also virtually simulated and compared with the newspaper reports.
The tsunami simulations show a 0.2 m wave amplitude at Angra do Heroismoa (Figure 12), which is in agreement with the report in 'Correio dos Açores'.

Conclusions
In this study, we used the tsunami modelling and observations to constrain the earthquake rupture. This is particularly important for old earthquakes where only P-wave polarities are observed. The use of P-wave data only does not allow for selection of the fault plane. The additional information provided by tsunami observations and numerical simulation proves to be a good tool to solve this ambiguity.
The 8 May 1939 earthquake in the Azores area is a good example where this methodology can be useful. In this case, very similar polarity data-sets resulted in two different focal mechanisms proposed, when different methods are applied, which adds to the fault rupture ambiguity.
The major seismotectonic questions addressed were as follows: (1) is the best focal mechanism a strikeÀslip or normal fault? (2) Can we distinguish between the two nodal planes and select one as the fault plane? (3) Is the preferred fault plane coherent with the revised epicentral parameters?
The main findings of this study are as follows: The comparison between the tide gauge record and the synthetic waveforms shows that (1) the SSÀSS focal mechanism (two nodal planes nearly pure strikeÀslip as proposed by Udias et al. (1976)) shows good agreement in terms of polarity, arrival time, and first wave amplitude, but attenuation after the first wave is incompatible with the Ponta Delgada recorded waveform, and (2) the SSÀNF focal mechanism (proposed by Buforn et al. (1988)) shows agreement in terms of polarity (for both planes A and B) and arrival time (for plane A, a low-angle nearly pure strikeÀslip), but presents higher frequency of the wave due to the values of the fault geometry and dip, suggested by the scaling laws.
Multiple simulation tests considering all four nodal planes of the SSÀSS and SSÀNF focal mechanisms, fault locations, and geometries highlight that the attenuation observed in the synthetic waveforms was caused by the bathymetric features separating the source and tide gauge location.
We found that a fault model based on the Buforn et al. (1988) nodal plane A with dimensions (L £ W) 40.5 km £ 40 km and a co-seismic slip of 1.5 m may generate a tsunami that best fits the recorded tsunami signal. The 8 May 1939 tsunami is best explained from an earthquake generated at a low-angle dipping fault with dominant strikeÀslip displacement. The fault is located in the area of shallow bathymetry, but it is compatible with the best constrained epicentres published (Table 1). The focal mechanism of the preferred fault model is only slightly different from the solution proposed by Buforn et al. (1988), plane A. As mentioned above, Buforn et al. (1988) proposed a focal mechanism using an automatic procedure based on the maximum-likelihood method of Brillinger et al. (1980). In this method, the likelihood function is evaluated using a 'precision parameter r' that can be assimilated to a signal/noise ratio. The choice of this parameter affects the results obtained by the method. These nodal planes that are obtained by the Brillinger et al. (1980) method using the Buforn et al. (1988) polarities, with the parameter r ranging between 1 and 10. If we compare these solutions with the best-fit one proposed in this paper, we see that Buforn et al. (1988) solution corresponds to r D 7, while the best-fit solution falls inside the set of possible solutions for variable r.
The geometry and slip movement of the preferred solution inferred by tsunami modelling (Figure 10) is unusual in the Azores Archipelago setting. However, we must bear in mind that the 8 May 1939 earthquake occurred at the intersection of the GF with the Terceira Rift which creates a complex and poorly known tectonic framework, as already outlined by Udias et al. (1976).
Finally, it is important to note that large earthquakes are rare in the Azores and that the original seismograms for events occurred before the WWSSN are hard to find. Thus, the use of all supplementary data, namely tsunami observations and tide-gauge records, constitute valuable data to shed light on the rupture mechanism.