New study on the 1941 Gloria Fault earthquake and tsunami

The M ∼ 8.3–8.4 25 November 1941 was one of the largest submarine strike-slip earthquakes ever recorded in the Northeast (NE) Atlantic basin. This event occurred along the Eurasia–Nubia plate boundary between the Azores and the Strait of Gibraltar. After the earthquake, the tide stations in the NE Atlantic recorded a small tsunami with maximum amplitudes of 40 cm peak to through in the Azores and Madeira islands. In this study, we present a re-evaluation of the earthquake epicentre location using seismological data not included in previous studies. We invert the tsunami travel times to obtain a preliminary tsunami source location using the backward ray tracing (BRT) technique. We invert the tsunami waveforms to infer the initial sea surface displacement using empirical Green’s functions, without prior assumptions about the geometry of the source. The results of the BRT simulation locate the tsunami source quite close to the new epicentre. This fact suggests that the co-seismic deformation of the earthquake induced the tsunami. The waveform inversion of tsunami data favours the conclusion that the earthquake ruptured an approximately 160 km segment of the plate boundary, in the eastern section of the Gloria Fault between −20.249 and −18.630 E. The results presented here contribute to the evaluation of tsunami hazard in the Northeast Atlantic basin.


Introduction
The western segment of the Eurasia-Nubia plate boundary extends from the Mid-Atlantic Ridge in the Azores towards the Strait of Gibraltar.Between −24 and −20 • E, the plate boundary is supposed to follow a prominent morphological feature, the Gloria Fault, first mapped by Laughton et al. (1972) (see Fig. 1 for location).The development of the Gloria Fault as a segment of the Eurasia-Nubia plate boundary occurred ∼ 27 Myr ago (Luis and Miranda, 2008).The relative interplate motion is slow, at ∼ 5 mm yr −1 (Fernandes et al., 2003), with unusual seismic activity.In spite of this, a few large earthquakes have been located close to this plate boundary.The predominant focal mechanism is strike-slip (Buforn et al., 1998(Buforn et al., , 2004;;Argus et al., 1989).Among these events are those that occurred on 25 November 1941 (Udias et al., 1976;Lynnes and Ruff, 1985), 9 June 1969 (Argus et al., 1989), 26 May 1975 (Buforn et al., 1998;Argus et al., 1989;Kaabouben et al., 2008) and 17 October 1983 (Argus et al., 1989).Some of these earthquakes generated tsunamis, described in historical documents and recorded by tide stations, namely on 31 March 1761, 25 November 1941and 26 May 1975(Baptista et al., 2006;Kaabouben et al., 2008;Baptista and Miranda, 2009).
In this study, we present a reanalysis of the seismic data of the 25 November 1941 event and the first comprehensive analysis of the associated tsunami.We use a newly acquired set of historic seismograms to re-evaluate the position of the epicentre and the computation of the seismic moment.We digitized, de-tided and filtered all tide records available in the Northeast Atlantic basin.The use of seismic and tsunami data provides a better understanding of this earthquake and tsunami.
In this study, we present a relocation of the earthquake source using the phases published by the ISC bulletin complemented with readings from nearby stations AVE in Morocco (−7.4133To compute the hypocentre, we used "HYPOCENTRE", a damped least square algorithm for earthquake location (Lienert et al., 1986), running under a SEISAN environment, a seismic analysis package containing a complete set of programs and a simple database for analysing earthquake data (Ottemöller et al., 2011).
Additionally, we tested two velocity models to compute the epicentre.The IASPEI91 velocity model (Kennett and Engdahl, 1991) locates the epicentre at −18.965 • E, 37.344 • N. The regional velocity model by Carrilho et al. (2004) for the Portuguese area locates the epicentre at −19.038 • E, 37.405 • N, which is our preferred solution.The focal depths given by these models are 12 km (root mean square -rms = 2.5 s) and 5.3 km (rms = 2.2 s) respectively.
We computed the scalar seismic moment using the Dineva et al. (2002) approach.Original analogue seismograms were digitized, and the seismic moment is computed independently from the spectra of body waves' ground motion for each component.Twenty-six seismograms from 14 seismic stations were digitized.Using P waves from 13 stations, we obtain a seismic moment of 3.96 1021 Nm, while using S waves from seven stations, the seismic moment gives 2.96 × 10 21 Nm.These seismic moment values correspond to a moment magnitude of 8.3 ± 0.2 and 8.2 ± 0.3 respectively.To re-evaluate the focal mechanism, we used the polarities from 72 seismograms.The parameters of the focal mechanism are strike 76 • , dip 88 • and rake angle −161 • .The focal mechanism is depicted in Fig. 1, and an example of a newly digitized seismogram is shown in Fig. 2.

Tsunami data analysis
After the earthquake, the tide stations in Portugal mainland, Morocco, Madeira and Azores islands (see Fig. 1 for locations) recorded a small tsunami.In the United Kingdom, Newlyn (−5.5 • E, 50.1 • N) tide station (not shown in Fig. 1) recorded the tsunami.Ponta Delgada tide station, in the Azores, recorded a maximum amplitude of 0.2 m.The Portuguese press reported tsunami observations close to Lisbon and the overtopping of shallow areas close to Oporto (Leixões) (see Baptista and Miranda (2009) for details).Haslett and Bryant (2008) described tsunami observations along the English Channel at Newlyn (Cornwall, United Kingdom) and Le Havre (France).Debrach (1946) and Rothé (1951) reported damage to the submarine cables Brest-Casablanca and Brest-Dakar.In Spain, the Bulletin of the Seismic Observatory in Almeria reported tsunami wave observations in the Canary Islands, but the tide record of Tenerife is unreadable.In this study, we selected the records of seven tide stations presented in Table 1.We digitized the original paper records except those of Casablanca and Essaouira (Morocco).Debrach (1946) presents drawings of these tide records: Casablanca-Petite Darse, Casablanca-Jetee Transversale and Essaouira.We could not find the original records of these stations.Instead, we had to rely on the drawings reproduced in Debrach (1946).We discarded the station Casablanca-Petite-Darse because the position of the tide gauge is uncertain.In this study, we use Casablanca-Jetee Transversale, here named Casablanca.We did not use the record of Newlyn because of the quality of the record.
All digitized records were linearly interpolated into a 1 min time step.For de-tiding, we used a least squares algorithm to fit the interpolated data into polynomials whose degrees guaranteed an adjusted R-square parameter above 0.995.To further ensure the goodness of our de-tiding we applied a bandpass filter.The parameters of the filter are as follows: end of the first stop band at 91 min, beginning of first passband at 89 min, end of the passband at 5 min and beginning of the second stop band at 2 min.The filter's magnitude response is of 30 dB attenuation in the first stopband, a passband ripple of 1 dB and a second stopband attenuation of 40 dB.We used Morlet's continuous wavelet analysis.Wavelet analysis has the advantage of providing information when a particular wave component is present, even if it appears only once in the time series.Hence, it can be better suited to the analysis of transient signals like tsunami waves.In this study, we present the results of wavelet analysis (Eq.1), using the Morlet mother function (Eq.2): The Morlet mother function is given by a plane wave modulated by a Gaussian function: where ω 0 is the non-dimensional frequency, here taken to be 6 to satisfy the admissibility condition.Figure 3 depicts the tsunami signals and the results of wavelet analysis.Madeira station presents the smallest dominant period close to 6 min.Ponta Delgada, Cascais and Lagos present a dominant period close to 12 min.Leixões presents a dominant period close to 20 min but this period is present in the spectrum well before the arrival of the tsunami -in this station the signal-to-noise ratio is low, when compared with the other records.Casablanca presents a 22.5 min dominant period.In the case of Essaouira, the dominant period is 45 min.

Preliminary location of the tsunami source
We used the travel time of the first wave arrival at each station to compute a preliminary location of the tsunami source.
To do this, we used the backward ray tracing (BRT) technique (Gjevik et al., 1997).We defined a set of points of in-terest (POIs), one per tide station.The position of each POI is the grid node closest to the actual location of the tide station, at a depth not less than 10 m to avoid strong non-linear effects.Table 1 shows the location and depth of POIs used in this study.We back-propagated the wave fronts from each POI using the algorithm implemented in the Mirone suite (Luis, 2007).The location of the tsunami source is the minimum of the averaged travel time square errors (in the least squares sense).The spatial distribution of the misfit is given by where t o i and t p i are the observed and predicted travel time for each POI, and n is the total number of POIs.
To compute the BRT solution, we used the tsunami travel times shown in Table 1 and a bathymetric grid with a horizontal resolution of 0.1 • .In Fig. 4, we present the minimum of the averaged travel time error, ε(x, y) = 7.2 min, at −19.0 • E, 37.35 • N.This location is 10 km to the SE of the seismic epicentre, −19.038 • E, 37.405 • N, presented in Sect. 2 (see Fig. 4).
The BRT method has limitations inherent to the linear shallow water approximation, implying an overestimation of the speed in shallow waters.Nevertheless, in the case of the 1941 tsunami, the good azimuthal coverage ensures that these limitations are partially averaged out.

Inversion methodology
To estimate the initial sea surface displacement (tsunami source) we need to invert the tsunami waveforms.Moreover, for tsunamis generated by earthquakes, we assume that the initial sea surface displacement mimics the elastic deformation of the seafloor, thus providing information on the earthquake focal mechanism.
The problem of tsunami waveform inversion without assuming a fault model was firstly addressed by Aida (1972in Satake, 1987).Since then, many studies have focused on the use of waveform inversion methods to estimate the tsunami source.These studies can be broadly divided into two categories: with a priori assumptions on the fault model (e.g.Satake, 1987Satake, , 1993;;Hirata et al., 2003;Titov et al., 2005), and without a priori assumptions on the source, namely Baba et al. (2005), Satake et al. (2005), Tsushima et al. (2009), Wu and Ho (2011) and Yasuda and Mase (2012).
Here, we use the method proposed by Tsushima et al. (2009), Koike (2011) and Miranda et al. (2014) based on the use of empirical Green's functions to efficiently perform the linear shallow water forward problem, where the synthesis of the tsunami waveform η m (t t ) at the m point of interest along the coast for the time span t t is given by a superposition principle: where (nl × nc) are the set of unit sources and h ij is the amplitude of the initial water displacement attributed to the ij unit source.G m ij (t t ) is the response to the unitary source of the ith element of the sea surface at the j th location, and m is the total number of unitary sources.The waveform inversion is simply obtained by solving We incorporate the suggestion by Tsushima et al. (2009) of defining an influence area around the minimum of the BRT.
The initial water displacement is forced to become null as the distance to the minimum increases, therefore allowing the computation of a spatially compact solution.To do this, to Eq. ( 5), we added additional equations of the form 0 0 0 0. ..q. ..0 where q is a positive number given by and r > r max is the distance between each parameter (unit source) and the minimum of BRT.r max is a cut-off distance after which the initial water displacement converges to zero.Equations ( 5) and ( 6) form a mixed-determined linear inversion problem (Menke, 2012) which can be solved in the sense of the least squares weighted minimum norm, by In Eq. ( 8), D is the first-order derivative operator and W e is a diagonal matrix that weights the likelihood that we attribute to the different observations.Each observation corresponds to a single measurement made by a tide gauge.Different values of p correspond to various levels of smoothness of the final solution.

Application to the 25 November 1941 event
In this work, the study area is shown in Figs. 1 and 3, of which the limits are −26.28 to −5.020 • E and 31.08 to 42.65 • N. The Green's functions database covers the area defined by −22.5 to −17.4 • E and 35.5 to 38.6 • N.This area has a total of 51 × 31 unit prisms.The dimensions of the unitary prisms are approximately 0.12 • × 0.12 • .The time span for the calculation was fixed as 180 min, and the time step as 1 min.
Experience shows that giving the same weight to all data points is misleading as waves that arrive later carry information resulting from reflections and interference, which are very hard to simulate with the inversion procedure.Here, we used a simple strategy of attributing a larger weight to the first incoming wave.We used only two different weights (Eq.8) for the data set.Data between the origin time of the earthquake and the first tsunami wave weighted 100 times more than the remaining data values, except for the Essaouira data.
Using the results of the BRT, we selected the area within the 10 min contour (corresponding to an average travel time error of less than 10 min) as the influence area (see Fig. 4).We performed several tests to assess the best choice of the parameters (p, r max ) that describe the relative quality of the data, the smoothness in the final solution and the size of the influence area, respectively.Parameter p is non-dimensional, while r is measured in grid units, and corrected for the latitude.
In Fig. 5, we show a "chessboard" synthetic source, with alternating +1/−1 m displacements, to study the resolving power of the inversion procedure.The test reproduces the geometry of the inversion procedure (same location of the POIs, and identical values of p, and r max ).The results presented in Fig. 4 show that the inversion algorithm recovers most of the chessboard information inside the influence area.Different values of r max do not change the results for the inner area qualitatively.Nevertheless, the value of p must be tuned by "trial and error" to get the best trade-off between resolution and smoothness of the solution.The geometry of the tide gauge network, with no stations SW of the source area, decreases the resolution of the solution in that direction and increases the sensitivity towards changes in the smoothness parameter p.
Figure 6 shows the solutions for the same set of parameters used in the chessboard test.In all cases the maxima of the initial displacement field correspond to segment BC (see Fig. 4), matching the 10 min error contour of the BRT and the earthquake epicentre.The lateral extent of the source is roughly 160 km.If we look at the solutions presented in Fig. 5 more closely, there is an apparent rotation when compared with the strike of segment BC.This is not a result of the azimuthal distribution of the tide gauge stations as the chessboard test shows (Fig. 4).It can be the result of the apparent contradiction between the records from Cascais (see lines 230-231) and Lagos tide stations, but we have no independent assessment of this problem to weight the station data differently.
Figure 7 depicts the comparison between the forward computation using the initial displacement field depicted in Fig. 6 Figure 5. Chessboard test of the inversion approach.The source has alternating +1/−1 m displacements, and the synthetic waveforms are computed with a shallow water code for the locations that correspond to our tide gauge network.Inversion results for different values of r max (Eq.4) and p (Eq. 5) are depicted.
for the preferred solution (p = 0.3, r = 8) and the observed data.

Discussion and conclusions
The seismological data of the 25 November 1941 earthquake lead to an epicentre location at −19.038 • E, 37.4050 • N, and a moment magnitude of 8.3.The new epicentre lies west of the epicentre presented by Guttenberg and Richter (1949).However, it is important to note that both determinations lie along segment BC (see Fig. 4) of the plate boundary, whereas Udias et al. (1976) and Antunes (1944) locate the epicentre further north.
The analysis of the tsunami data shows a clear signal of small amplitude in most of the stations, with an apparent change in the frequency content except for Leixões.The spectral analysis of the tsunami records shows different frequency contents.These differences are mainly caused by differences in local morphological conditions and to the relative position of the station and the seismic source.Madeira station located approximately across-strike shows the smallest dominant period of 6 min.Leixões, located along-strike, shows the largest dominant period of 20 min.Ponta Delgada, Cascais and Lagos show dominant periods close to 11 min (see Figs. 1 and 3).These results are consistent with a tsunami source located along the Nubia-Eurasia plate boundary.
The BRT simulation locates the minimum of the averaged travel time (7.2 min), at −19.000 • E, 37.350 • N, quite close to the epicentre location.This result is consistent with the fact that the earthquake's co-seismic deformation induced the tsunami.It also shows the coherence of the tsunami data set.
Because of the limitations inherent to the use of the linear shallow water approximation and the location of the tide gauges inside the harbours, we opted to give different weight to the data points corresponding to up the first incoming wave.The set of solutions shown in Fig. 6, corresponding to different values of p (smoothness) and r max (influence area), show a similar deformation pattern.All solutions show a smooth profile, with a maximum upward displacement of +0.6 m compatible with the initial sea surface displacement computed with Okada (1985) formulae for the co-seismic deformation due to a fault in an elastic half-space (see Fig. 6a).Moreover, all inverse-problem solutions favour the conclusion that the earthquake ruptured approximately 160 km of the plate boundary between −20.249 and −18.630 • E.
The results presented in Fig. 7 show that the inversion technique used here can fairly reproduce the first incoming wave for most of the observations presented in Table 1.In Cascais, the synthetic wave arrives 8 min earlier than the observed wave.The location of the tide station inside the old harbour in a very shallow area (less than 4 m) might be responsible for this discrepancy.In addition, the smoothness coefficient used to compute the inverse solution may contribute to a "spreading" of the lateral extent of the source, thus producing earlier arrivals at some points on the coast.The observation in Essaouira is incompatible with all solutions presented in Fig. 6.The lack of the original tide records of Morocco and the poor quality of their reproduction in Debrach (1946) does not allow for further investiga-tion.However, we cannot exclude the hypothesis of a local second tsunami source close to the coast of Morocco.Debrach (1946) and Rothé (1951) report damage to the submarine cables Brest-Casablanca and Brest-Dakar after the earthquake.This fact may suggest the occurrence of a submarine landslide close to the coast of Morocco.
In spite of the limitations inherent to the use of old instrumental records, some of them with a low-amplitude signalto-noise ratio, the inversion of the tsunami waveforms provides an independent estimate of the extension of the tsunami source and indirectly, the extent of the seismic source.
The seismic source solution (strike 76 • , rake 161 • ) also favours the conclusion that the earthquake ruptured segment BC (see Fig. 4) of the plate boundary.This position coincides with a segment of the plate boundary that is marked by a change in strike when compared to the Gloria Fault between −23.960 and −20.249 • E (AB in Fig. 4), and to the eastern segment CD (see Fig. 4) between −18.142 and −15.972 • E. While the relative motion between Eurasia and Nubia generates pure strike-slip events along AB and CD, there is a significant thrust component associated with the seismic rupture within segment BC.This is consistent with both (seismic and tsunami) studies.Segments AB and CD follow the small circle of the Eurasia-Nubia rotation pole, while segment BC makes an angle of approximately 12 • with the small circle of Nuvel 1A.
We cannot discard the fact that a local landslide occurred because the break of the submarine cables is documented.We do not have information either on the location of the submarine cable nor where it broke.This lack of information makes the modelling of the landslide very difficult.We think that the occurrence of a landslide close to the coast of Africa (Morocco and Senegal) might influence the signal of the Morocco tide records.However, our strongest tsunami observation is in the Azores (Ponta Delgada); this station is too far away from the submarine cable (and to the possible landslide location), therefore we conclude that the tsunami observed there is mostly due to the earthquake's co-seismic deformation.
Due to the small Eurasia-Nubia relative motion in this segment of the plate boundary (5 mm yr −1 ), the recurrence of similar events does not happen often.However, other parts of the same plate boundary can generate similar events, with noticeable effects in the North Atlantic coasts, and the tsunami-genic potential of the Gloria Fault is relevant for the quantification of tsunami hazard in the Northeast Atlantic basin.

Data availability
Copies of the paper tide records can be obtained from the authors.

Figure 1 .
Figure 1.General overview of the Northeast Atlantic at the latitude of Iberia and focal mechanism of the earthquake.The yellow star represents the epicentre of the earthquake.Black dots show the location of the tide stations.

Figure 2 .
Figure 2. Top panel: digitized seismograms corresponding to the N-S and E-W components of the 1941 Atlantic earthquake as recorded by the Wiechert seismograph at Uppsala (UPP).Bottom left panel: ground motion spectra of the P waves, as marked at the top figure (blue and red, respectively), recorded at UPP, and the corresponding fit to the Brune model.Bottom right panel: ground motion spectra of the S waves, as marked at the top figure (blue and red, respectively), recorded at UPP, and the corresponding fit to the Brune model.

Figure 3 .
Figure 3.For each station: tsunami signals (de-tided and filtered record) (left panels), wavelet amplitude spectrum (middle panels) and global wavelet spectrum (GWS) (right panels).The initial time, 18:04 UTC, is the time of the earthquake.

Figure 4 .
Figure 4. Tsunami source location given by BRT.The red star represents the minimum of the backward ray tracing misfit (∼ 7 min) at −19 • E, 37.35 • N, less than 10 km away, to the SE of the seismic epicentre −19.038 • E, 37.405 • N (yellow circle).Black dots represent the tide stations.Segments AB and CD follow the small circle of the Eurasia-Nubia rotation pole; the angle between segment AB and BC is 12 • with the small circle of Nuvel 1A.

Figure 6 .
Figure 6.(a) Co-seismic deformation using the focal mechanism parameters proposed in this study: L is 170 km, W is 45 km; strike is 255 • , rake is 161 • , dip is 88 • , slip is 8 m, yellow triangle depicts epicentre; (b)-(j) initial displacement field of the 1941 Gloria Fault tsunami, obtained from waveform inversion for different values of r and p parameters; dashed line sketches the approximate location of the plate boundary.

Figure 7 .
Figure 7.Comparison between observed and synthetic waveforms.The solution shown corresponds to r = 8, p = 03.Red Line depicts the synthetic waveform.Black line depicts the observed waveform.Time zero is the time of the earthquake.

Table 1 .
POIs (points of interest) of the 1941 tsunami used in this study.Longitude, latitude and depth refer to the POI.POIs are named after the closest tide station.