The Horseshoe Abyssal plain Thrust could be the source of the 1755 Lisbon earthquake and tsunami

The southwest Iberia margin is widely believed to have hosted the 1755 Great Lisbon earthquake and ensuing tsunami, one of the most destructive natural events in European history. Here we combine geophysical data and numerical tsunami modelling to investigate the source and mechanism responsible for this event. We find that an intra-plate, lithospheric¬-scale thrust fault located at the Horseshoe Abyssal Plain coincides with the location and focal mechanisms of the largest regional earthquakes and is likely to have suitable dimensions and fault-rock properties to account for the magnitude of the 1755 event. We present tsunami simulations with this fault as the source, and find that they reproduce reported tsunami energy propagation patterns, arrival-times and run up heights more successfully than other modelled sources. We propose that a reverse dip-slip mechanism on the northwest verging Horseshoe Abyssal plain Thrust, combined with the two-state mechanical behaviour of serpentinite, is the most likely candidate for the source of the 1755 Great Lisbon earthquake and for other recent large regional earthquakes. The 1755 Great Lisbon earthquake potentially had its source on the Horseshoe Abyssal plain Thrust, southwest of Iberia, according to a comparison of geophysical observations, historical reports and numerical modelling of the tsunami that was generated

T he present-day Eurasia-Africa plate boundary off southwest Iberia is dominated by a northwest (NW)-southeast (SE) trending convergence (3.8-5.6 mm yr −1 ) 1 (Fig. 1). The plate boundary is diffuse 2 and deformation is accommodated along a 600 km-long deformation band 2 by a series of thrust and strike-slip faults 3,4 . It is one of the most seismically active regions in Western Europe, showing a constant background of seismic activity with low to moderate magnitude events of up to magnitude (M w )~6.0 [5][6][7][8] , mostly concentrating at upper mantle depths [6][7][8] . Besides the background seismicity, this area has also experienced large magnitude events in historical and instrumental times, such as the 1755 Lisbon (estimated M w 8.5+) [9][10][11][12][13][14] and the 1969 Horseshoe (M w 7.9) 15 ones.
The 1755 Great Lisbon earthquake has the largest documented felt area of any shallow earthquake 9,16 . The ground motion was felt as far away as Hamburg, Great Britain, the Azores and Cape Verde Islands 9 , producing a remarkable agitation of lakes, rivers and harbours in a rhythmic way as far as Great Britain, Switzerland, Finland and Sweden 17,18 , and it was one of the most destructive known earthquake in European history 19 . The associated tsunami, with run-up heights reported to reach 5-15 m 15 , devastated the coasts of SW Iberia and NW Morocco. The event is thought to have caused near 20,000 victims 10 combining the effects of ground-shaking, the related fires and ensuring tsunami, which were reported to hit as far as England, Newfoundland, the Caribbean and Brazil 19 . Given such evidence, the scientific community widely supports an estimated M w ≥ 8.5 [9][10][11][12][13][14] for the Lisbon earthquake, although a recent work has also suggested a lower magnitude 20 . The source and the generating mechanism of the Lisbon earthquake have been the subject of numerous studies and the matter of a continuous debate during the last decades  . A number of candidate sources have been proposed 11,[22][23][24][25][26][27][28] , although none of them explains satisfactorily the inferred seismic moment of the event, as well as the arrival times and, particularly, the run-up heights of the tsunami referred to in historical reports.
The southwest Iberian margin has undergone a complex geodynamic history spanning from the Mesozoic extension to the Neogene-to-present-day convergence. Martínez-Loriente et al. 30 used P-wave seismic tomography models obtained along two regional controlled-source wide-angle seismic profiles (profiles P1 and P2 in Fig. 1), to show the presence of three different oceanic domains in the area originated in three different opening systems (Fig. 1a). First, the Seine Abyssal Plain (SAP) domain, which is made of oceanic crust generated during the first stages of seafloor spreading in the NE Central Atlantic (Early Jurassic). Second, the Gulf of Cadiz (GC) domain, formed by Jurassic-age oceanic crust generated at the Alpine-Tethys oblique-spreading system between Iberia and Africa (i.e., the Western Tethys) 31 . Third, the Gorringe Bank (GB) domain, composed of unroofed mantle rocks with little synchronous magmatism. These rocks, which constitute the basement of the GB, southern Tagus and northern Horseshoe Abyssal plains (HAP), were exhumed during the first stages of North Atlantic opening (Lower Cretaceous) 32 . The present-day plate configuration was established~20 Ma ago, when the plate boundary between Africa and Eurasia shifted from north to south Iberia 33 (Fig. 1b). The Neogene NW-SE convergence resulted on the emplacement of large allochthonous masses during the  Table 1). White and grey circles show epicentral locations of regional earthquakes with M w > 3 for the period 1960-2020 57  Tortonian 34 and caused the reactivation of major Mesozoic extensional faults 3 . These data also allowed identifying the presence of the Horseshoe Abyssal plain Thrust (HAT), a deep thrust fault located in the HAP, away from any well-established subduction zone.
The deformation resulting from the present-day NW-SE convergence (3.8-5.6 mm/yr) −1 is mainly accommodated by two types of active faults. The first includes regional NE-SW trending thrusts, such as the HAT 30 , the GB Fault (GBF) 9,32 or the Marquês de Pombal Fault (MPF) 24,25 , whereas the second comprises large west northwest-east southeast trending dextral strike-slip faults (i.e., the South West Iberian Margin (SWIM) faults) 28 , such as the so-called Lineaments North and South 3,35,36 (LS) (Fig. 1). According to Martínez-Loriente et al. 30 , the SAP and GC oceanic basement domains are separated by the LS strike-slip fault, whereas the GC and GB domains are limited by the HAT. This means that the HAT separates exhumed mantle rocks (NW) from oceanic lithosphere (SE) (Figs. 1 and 2).
In this work, we show that the HAT has the dimensions, geometry, and fault-rock bulk properties required to explain the deep regional seismicity, and that it has the potential to generate earthquakes as large as the 1755 Lisbon one. In addition, we have performed tsunami simulations to show that the HAT source can generate transatlantic tsunamis with characteristics compatible with those referred in historical records of the 1755 event. Finally, we discuss the regional and global implications linked with the presence of this deep thrust and the bulk-rock properties in terms of seismic and tsunami hazard assessment.

Results and discussion
Origin of deep seismicity offshore SW Iberia. The seismic activity in the study area mostly concentrates north of the LS, between the Gorringe and Guadalquivir banks (Fig. 1). Regional micro-seismicity recorded by temporal networks of Ocean Bottom Seismometers (OBS) 6-8 has allowed identifying three main seismicity clusters located in the GB, HAP, and Sao Vicente Canyon (SVC) 8 displaying a combination of thrust and strike-slip focal mechanisms. These events, which concentrate at upper mantle depths (i.e., between 30 and 60 km depth), appear to be associated with inverted rift structures identified at the transition zones between the different rheological domains that localize regional stresses and hence seismic strain.
Previous works on the 1755 earthquake have proposed a number of candidate sources for this event, such as the GB 9 , MPF 11,22,24,25 , the subduction beneath Gibraltar 26 , the SWIM strike-slip faults 28 , Horseshoe Fault 24 , São Vicente Fault 24 or some composite systems 22,24 . However, due to either their limited surface areas, their locations and geometry, or both, they do not appear to be capable of accounting for the large seismic moment, the overall reach and energy pattern of the tsunami, as well as the historical recorded amplitudes and arrival times. Among all the proposed candidates, the only tsunami simulation considering far-field observations 27 Table 1). Serpentinization degree is indicated as gradation of colour.
Supplementary Table 1). Both earthquakes share the focal mechanism, location and a deep hypocentre, suggesting that they might have been generated by the same fault system. The M w 7.9 1969 earthquake generated a small tsunami with centimetric wave heights on the SW Portuguese coasts 39 . For comparison, several metre wave heights were reported in the same area after the 1755 earthquake, indicating that the energy released by this event was much larger than that of the 1969 one (between 15 13 and 30 19 times larger according to previous works). The analysis of peak ground accelerations estimated from reported intensities of the 1755 earthquake supports a NE-SW oriented reverse dip-slip mechanism 37 . This would be similar to that of the 17 December 2009 earthquake (M w = 5.5; Supplementary Table 1) 37 associated with the MPF, located at the mouth of the SVC where the HAT changes its orientation in depth and would be linked to the MPF (Figs. 1 and 2, and Supplementary Fig. 1).
As stated above, recent work shows evidence for the presence of a deep structure in the HAP 30 between the SAP and GB domains with the location and geometry required to be the source of these large earthquakes. This structure was interpreted as a NE-SW trending, NW-verging blind-deep thrust located in the middle of the HAP (Figs. 1 and 2). Although the HAT was initially interpreted to extend to the SVC, the re-analysing of available geophysical data, in combination with recent seismological studies [5][6][7][8]37,38 , suggest that it rather spans from the MPF to the LS (~200 km). Its location and geometry makes it a likely source candidate for part of the deep seismicity recorded by OBS networks in the area (Figs. 1 and 2).
Earthquakes result from stick-slip frictional instability and its frictional behaviour is commonly described by a rate-and statevariable constitutive law 40 . At steady-state, the frictional response depends on a slip velocity-dependent friction parameter. Depending on the sign of this parameter, the material is either velocity strengthening (i.e., earthquakes cannot nucleate and any earthquake propagating there will tend to stop) or velocity weakening (i.e., earthquakes can nucleate and propagate there). This parameter depends on rock lithology, temperature and pressure. For example, granite is velocity weakening below 300°C and velocity strengthening above, due to the onset of plasticity of quartz, the most ductile mineral in granite. In oceanic basalt, frictional stability is controlled by feldspar transition, which occurs at~450°C 32 . This temperature is typically assumed to mark the downdip limit of inter-plate seismicity in subduction zones.
According to the interpretation of the local controlled-source seismic tomography models 30 , the HAT separates serpentinized peridotite basement at the NW (GB domain) from oceanic lithosphere at the SE (GC domain), so the deep frictional stability, which concentrates at upper mantle depths, should be controlled by the rheological transition of olivine, which occurs at~600°C 40,41 . Oceanic plate cooling models indicate that this isotherm is located at 50-55 km depth in a~140 Myr-old lithosphere 42 , in agreement with the location of the deepest earthquakes in the area 6 (Fig. 2). This pre-eminence of deep seismicity strongly supports that it is controlled by peridotite/ serpentinite rheology. The high peridotite rigidity (~70 GPa) and the potential large stress drop, combined with the large downdip extension of the seismogenic layer (>60 km for a 45°thrust), makes that modest-length thrust faults such as the HAT could potentially account for large earthquakes 5 . Therefore, if we assume a focal depth of 50-60 km on peridotite rock (which has a rigidity of~70 GPa) and an upward propagation over lower rigidity rocks such as serpentinite (for an average rigidity of~50 GPa), a dip angle of 45°and a realistic slip-to-length ratio of 8.5 × 10 −5 , a M w 8.5+ earthquake like that of the 1755 Lisbon one requires a fault length of just 160-220 km (Fig. 3).
In this rupture scenario, seismic rupture at shallower levels should be controlled by serpentinite, which is more ductile than peridotite. Experimental results show that, at crustal pressure, temperature conditions, serpentinite is velocity strengthening at plate tectonic rates but it is unstable at higher velocities. This means that although serpentinized regions are not a likely site of earthquake nucleation, seismic rupture of deeper earthquakes could propagate through them 43 . This two-states mechanical behaviour of serpentinite could explain two apparently contradictory observations: the lack of shallow earthquakes in the region 6-8 and the occurrence of large, occasionally devastating tsunamis requiring shallow slip and substantial seafloor deformation 44 .
Simulation of the tsunami caused by a HAT-like source. Although multichannel seismic data and TOPAS sub-bottom profiles across the HAT (Fig. 1) do not provide clear evidence of the thrust reaching the seafloor, the OBS records 30 show a laterally coherent deep seismic phase consistent with a reflection at a NW-verging structure reaching at least 25-30 km depth, which is interpreted to extend to at least 60 km depth based on seismological data [5][6][7][8]36 . Such a structure can totalize~16,700 km 2 of surface area (Figs. 1 and 2, and Supplementary Fig. 1). The controlled-source seismic tomography model shows that the dip of the blind thrust changes with depth, from 20°to 30°from the surface to 25-30 km depth, to 40°-45°deeper than this level, as indicated by focal mechanism solutions of deep seismicity on this area [5][6][7][8]36 (Fig. 2 and Supplementary Fig. 1).
We have used this fault plane as a source for tsunami simulations (Methods), in order to verify if such a feature is compatible with the reported characteristics of the largest tsunamigenic earthquakes in the area, and in particular with those of the 1755 Lisbon one (Fig. 4). To test this, the simulated travel times, amplitudes and tsunami propagation resulting from the HAT source are compared with the available historical records of the 1755 event [10][11][12][13][14][16][17][18][19]27,[45][46][47][48] and with previous works 11,26,27 (Figs. 5  and 6, and Supplementary Tables 2 and 3).
The modelled tsunami propagation fits better than any other proposed candidate with the information included in the available historical chronicles (Figs. 4, 5 and 6, and Supplementary Tables 2  and 3). The HAT source scenario reproduces wave arrivals to all locations reported to be hit by the 1755 tsunami including the distant ones such as the Ireland, Great Britain, Newfoundland, Barbados or Brazilian coasts. In addition, it correctly predicts the zones protected from the tsunami impact, such as Vigo in West Iberia or cities like Boston, Baltimore, and New York in the eastern US coast. The bathymetric roughness has an important effect on tsunami propagation, dispersing, amplifying or diminishing the tsunami waves 49 . For instance, and according to our simulation, Vigo was protected from the tsunami by the rías (i.e., sea-inlet or estuary), southern Florida by the Bahamas Banks, and the Madeira Tore-Rise acted as scatter for tsunami energy protecting the US East coast. In contrast, the Madeira island scattered the energy towards the Caribbean, whereas the Great Meteor seamount scattered it towards Brazil.
Comparing arrival times and tsunami amplitude in particular locations is problematic due to the inaccuracy and incompatibility of some arrival times in historical reports. For instance, they indicate a 50 min difference in the arrival time between Lisbon and Porto Novo, which are just 50 km apart, 35 min difference between Penzance and Newlyn, which are 2 km apart, or 30 min difference between Porto Santo and Madeira Islands, which are 50 km apart. Likewise, they report huge wave heights reaching~15 m in Lisbon, Cadiz, Tangiers and Azores, which are hard to believe 26,27 . Historical observations contain information on estimated run-up heights (i.e., maximum height of the inundation), but tsunami simulations provide maximum wave height (H m ) before reaching the coast. Therefore, we are comparing different features. To overcome this issue, we have normalized our results to 1 m (Ĥ m ) below seafloor according to the Green's Law 47 (Methods). Almost all previous works 13,15,21,22 assumed a run-up of 15 m in Cadiz, although the French scientist Louis Godin, who was in Cadiz when the tsunami hit the city, reported a wave height of 5 m 50 . Shallow water estuaries are highly affected by non-linear propagation effects, which can induce notorious modelling errors. Therefore, we included in our comparison Oerias and Vila Real de San Antonio instead of Lisbon and Huelva (respectively), which are affected by estuaries.
Overall, our results fit better than any previous modelled source of the reported observations (Figs. 5 and 6, and Supplementary Tables 2 and 3). In the localities with information available for comparison in the four simulations, the  improvement is remarkable (~80%) for the maximum wave heights, whereas the arrival times, which suffer from evident inconsistencies, are comparable or slightly better on average. The HAT modelling results show lower deviations in the estimated wave heights from historical reports in all the reported comparable cases (i.e., 24 localities) except for Porto, Cadiz and Tangier. In the comparison between the only two sources capable of generating transatlantic tsunamis, the HAT always obtains a better fit than Source 8, both in the far-and near-fields (except in Porto, where the difference is minimal). Regarding arrival times, the HAT simulation obtains good results in the far-field, with average differences of <20%, and acceptable results in the nearfield, where we obtain values within or close to the reported time range, except in localities with incompatible values for the differential arrival times (i.e., Lisbon, Safi, Madeira and Azores) as also pointed out in previous works 26,27 . It is noteworthy that only one of the previous works proposed a source capable of generating a transatlantic tsunami and considered near-and far-field historical observations. Their results favour a tsunami generated by a west northwest-east southeast trending fault located in the centre of the HAP 27 , so that it is at the same location as the HAT but with a different strike direction. However, the proposed orientation was chosen ad hoc for the purposes of the modelling, but it does not correspond to any known structure in this area and is nearly parallel to the plate convergence direction.
Regional and global implications. In summary, the HAT can explain the location and focal mechanism of the largest instrumental events having occurred in the area, and has the dimensions and the inferred fault-rock bulk properties at focal depths to account for such a giant earthquake. In addition, the tsunami simulation fits better than any other proposed source of the historical reports from the 1755 Lisbon event, particularly the observed wave heights.
Our work has relevant implications at regional level in the fields of active tectonics, petrophysics, seismology and for the seismic and tsunami hazard assessment. We show that the deep regional seismicity is controlled by peridotite/serpentinite rheology, allowing modest-length thrust faults to generate anomalously large earthquakes (Fig. 3). This suggests that there is a number of active thrusts in the area 3,4,25 similar to the HAT that might also have the potential of generating greater earthquakes than previously thought. One example is the GBF, a 155 km-long, NW-verging thrust fault located in between the Tagus and the HAP that thrust two exhumed mantle blocs 32 (Fig. 1). Therefore, the earthquake and tsunami hazard associated to these structures must be re-evaluated in light of the new findings. Globally, the M w 8.5+ Lisbon earthquake is a case study evidencing that great destructive tsunamigenic earthquakes can also occur in intra-plate settings (i.e., away from a well-defined plate boundary), challenging the widely extended idea in the scientific community that these events are restricted to subduction zones. There are some recent examples worldwide of great earthquakes in intra-plate settings. The most relevant ones are two large earthquakes of M w 8.6 and M w 8.2 associated with modest-length strike-slip faults occurred in the Wharton Basin (Indian Ocean) in April 2012 51 . This has been the largest intra-plate earthquake sequence that has been instrumentally recorded 52 . The events initiated ∼400 km to the SW of northern Sumatra, Indonesia, on the oceanic side of the Sunda megathrust, within the Wharton Basin and west of the Investigator Fracture Zone 51 . According to these authors, the WNW-trending faults experienced large slip and stress drops. The slip was deep and extended into the oceanic mantle (at least up to 60 km depth). The event had such a large magnitude, because these faults experienced very large slip in highrigidity material 51 . Therefore, the occurrence of great earthquakes in intra-plate settings and the possible associated tsunami should be considered in geohazard assessments worldwide.

Methods
Tsunami modelling. To evaluate the seismic and tsunamigenic potential of the HAT and estimate if this structure would generate great earthquakes (M w ≥ 8.5) and transatlantic tsunamis, we first specify fault parameters and physical properties, considering that: log 10 ðM 0 Þ À 6:07 ð1Þ where M w is seismic magnitude, M 0 is seismic moment, μ is rigidity, A is fault area and Ḋ is slip rate. M 0 is 7.16 × 10 21 Nm for M w = 8.5. A is 16.730 km 2 (Supplementary Fig. 1) and μ is 50 GPa. This is an average value assuming that the earthquake nucleated in fresh peridotite at~50 km depth (~70 GPa) and propagated along the fault plane up to shallower levels where serpentinite and basalts predominate (~30 GPa) 5 . Therefore, the slip value derived from the scaling relations that is required to generate a large earthquake (M w = 8.5+) is 8.5 m. We carried out a succession of tsunami simulations using the HAT fault plane (Supplementary Fig. 1) with a rake of 90°, the parameters described above and varying the slip from 7 to 20 m. Our results show that a slip of 15 m is required to generate a transatlantic tsunami that affects all localities that reported the 1755 event. This value is similar to that assumed in previous works 9, 26,27 (Fig. 3b). Tsunami simulations, for the equivalent of 12 h of propagation, were performed using the finite-volume Tsunami-HySEA software 53 . Tsunami-HySEA is a Graphics Processing Unit-optimized code that solves non-linear shallow water equations in spherical or Cartesian coordinate system 53 . The computational grid has a spatial resolution of 2 arcmin and covers all the Atlantic from Newfoundland to Brazil, considering open boundary conditions. The initial seawater elevation is considered equal to the vertical displacement of the sea bottom, computed as a Volterra's formulation of the elastic dislocation theory applied to a source buried in an elastic half-space, using the Okada's fault deformation model 54 . The rupture is considered instantaneous and the initial velocity field is zero. Results of tsunami modelling were processed to obtain regional tsunami wave propagation patterns and maximum wave height (H m ) profiles at selected coastal points. H m were extracted at~50 m depth. They are subsequently extrapolated at the 1 m depth using the Green's law approximation for one-dimensional shoaling as the fourth root of the 50 m, i.e., by applying a multiplicative amplification factor of~2.66 (Supplementary Table 2). Arrival times at the coast were calculated with the TTT software (http://www.geoware-online.com/tsunami.html), using the 2 arcmin digital bathymetry (SRTM15, https://topex.ucsd.edu/WWW_html/srtm15_plus. html), 64 Huygens nodes, and a set of source points on a rectangle obtained from the projection of the fault on the Earth's surface (Supplementary Table 3).
The leading wave of a tsunami is considered to be a shallow water gravity wave, where the ocean depth is negligible compared to the wavelength. In deep waters, the propagation velocity is approximately proportional to ffiffiffiffiffiffiffiffiffi g h p (3), where g is the acceleration of gravity and h is the water depth in metres. In deep waters, the tsunami wave travels at a velocity of~800 km/h, whereas at a water depth of 1000 m the velocity drops to 360 km/h. In shallow continental shelf, the wave propagation speed drops dramatically. Under these conditions, an excellent bathymetric grid is required for accurately modelling amplitudes and tsunami arrival times. Due to the enormous dimensions of our tsunami scenario and the unavailability of high-resolution bathymetry for all coastal areas, we used a grid size of 2 arcmin for the simulations.
When comparing modelled results with historical information, it should be considered that some reported observations are incompatible in many cases for any single source, so they are not credible, such as Lisbon, Porto Novo, Penzance or Madeira. In addition, some delays can be due to poorly constrained shallow shelf bathymetry and to the potential effects of landslides (i.e., multiple tsunami sources). However, the turbidite related to the Lisbon earthquake with a thin volume (1.3 km 3 ) is unlikely to have contributed importantly to the 1755 tsunami 29 . In addition, an increase of 2-3 m should be considered due to the high tide at the time of the tsunami 27,45,55,56 . In order to facilitate the comparison of results, we have calculated the relative variation of the maximum wave height (ΔH m ) and the arrival time (Δt a ) between the historical values reported (H obs m and t obs a ) and every simulation (H calc m and t calc a ) at each location of Tables 2 and 3 ( Figs. 5b and 6b).
and Δt a ¼ jt obs a À t calc a j t obs a ð4Þ It is noteworthy that to enable the comparison of results in Fig. 5b, the wave heights have been approached to zero in the localities where the impact of the 1755 tsunami was not reported (identified with "NR" in Table 2 and in Fig. 5).

Data availability
The bathymetric data used for tsunami simulations are available from GEBCO (https:// www.gebco.net/data_and_products/gridded_bathymetry_data/). Detailed bathymetry used to generate Figs. 1 and 2 is published in Zitellini et al 28 . The epicentral location of the seismicity shown in Fig. 1 is from the IGN Seismic Catalogue (https://www.ign.es/ web/ign/portal/sis-catalogo-terremotos). The HAT fault plane needed to reproduce the tsunami simulation is available in the figshare repository: https://doi.org/10.6084/m9. figshare.14589315.