How Alpine seismicity relates to lithospheric strength

Despite the amount of research focussed on the Alpine orogen, different hypotheses still exist regarding varying spatial seismicity distribution patterns throughout the region. Previous measurement-constrained regional 3D models of lithospheric density distribution and thermal field facilitate the generation of a data-based rheological model of the region. In this study, we compute the long-term lithospheric strength and compare its spatial variation to observed seismicity patterns. We demonstrate how strength maxima within the crust (~ 1 GPa) and upper mantle (> 2 GPa) occur at temperatures characteristic of the onset of crystal plasticity in those rocks (crust: 200–400 °C; mantle: ~ 600 °C), with almost all seismicity occurring in these regions. Correlation in the northern and southern forelands between crustal and lithospheric strengths and seismicity show different patterns of event distribution, reflecting their different tectonic settings. Seismicity in the plate boundary setting of the southern foreland corresponds to the integrated lithospheric strength, occurring mainly in the weaker domains surrounding the strong Adriatic plate. In the intraplate setting of the northern foreland, seismicity correlates to modelled crustal strength, and it mainly occurs in the weaker and warmer crust beneath the Upper Rhine Graben. We, therefore, suggest that seismicity in the upper crust is linked to weak crustal domains, which are more prone to localise deformation promoting failure and, depending on the local properties of the fault, earthquakes at relatively lower levels of accumulated stress than their neighbouring stronger counterparts. Upper mantle seismicity at depths greater than modelled brittle conditions, can be either explained by embrittlement of the mantle due to grain-size sensitive deformation within domains of active or recent slab cooling, or by dissipative weakening mechanisms, such as thermal runaway from shear heating and/or dehydration reactions within an overly ductile mantle. Results generated in this study are available for open access use to further discussions on the region.


Introduction
The present-day seismicity distribution is still poorly understood and presently debated in the Alpine orogen and its northern and southern forelands (Deichmann 1992;Bonjer 1997;Cattaneo et al. 1999;Singer et al. 2014;Eva et al. 2015). Therefore, any additional constraints on the controlling factors affecting event localisation are of interest. Local models of seismic focal mechanisms have been used with some success to assess the stress regime on individual faults thereby offering important insights into the local dynamics driving the observed localisation of seismicity (e.g. Bonjer 1997). However, despite its relevance in the current ongoing scientific debate, the correlation between monitored seismicity and its localisation in space with respect to the long-term stress state of the whole lithosphere is still lacking.
Previous works in the region have attempted a quantification of the long-term mechanical state of the lithosphere mainly relying on 2D sections across the Alpine chain such as the EGT (Okaya et al. 1996) and Transalp (Willingshofer and Cloetingh 2003) profiles. Relatively few models exist that have attempted to compute lithospheric strength variations across the entire orogen (e.g. Tesauro et al. 2011;Marotta and Splendore 2014). In this contribution we propose a revision of the current understanding of the long-term lithospheric strength of the lithosphere in the Alpine area in light of higher resolution, region specific 3D geophysical models. To this end, we make use of recently published results derived from a gravity and seismically constrained structural and density model of the region (Spooner et al. 2019a) that has been verified by secondary sources (Magrin and Rossi 2020) along with a wellbore measurement-constrained thermal field (Spooner et al. 2020a) to arrive at an updated model of the rheological configuration of the study area.
After a detailed summary of the methodology used to compute lithospheric strengths, we dedicate the second part of the manuscript to a critical discussion of the correlation between modelled lithospheric strength variations and a comprehensive dataset of regional seismicity to investigate the potential role that lithospheric heterogeneity plays in shaping the observed localisation of seismicity throughout the Alps and their forelands. In depth analysis leads us to partially revise the main conclusions derived from a recent numerical work that attributes the seismicity depth distribution across the region solely to the relatively low plate convergence rate, with negligible influence from tectonic inheritance (Dal Zilio et al. 2018).
We also quantify the role of viscosity contrasts caused by the presence of a laterally heterogeneous lithospheric configuration in limiting the maximum depth of seismicity via dissipative creep mechanisms. As a result, the 3D distribution of effective viscosities calculated within this work, are able to complement the ongoing debate on the relative impact of glacial isostatic adjustment versus tectonic and mantle dynamic processes as causes of the observed presentday kinematic state of the Alpine region (e.g. Norton and Hampel 2010; Chéry et al., 2016;Mey et al. 2017;Sternai et al. 2019).

Geological and seismic history
Large scale crustal differentiation within the northern foreland of the Alps (the European plate) primarily results from the Carboniferous age Variscan orogeny (Franke 2000), such as the juxtaposition of the Moldanubian and Saxothuringian terrains (Babuška and Plomerová 1992;Freymark et al. 2017) and the assemblage of crystalline basement presently exposed in the Vosges, Black Forest and Bohemian massifs. Heterogeneity within the Alpine orogen, largely stems from the ongoing collision of the Adriatic plate with the European plate from the Cretaceous to present-day (Handy et al. 2010). Fault systems such as the Periadriatic Lineament, that runs east to west through the spine of the orogen, and the Giudicarie Fault that offsets it 100 km to the north-east, formed to accommodate this convergence. Traditionally, the Alpine crust is split up according to its margin of origin prior to orogenesis, such as the Helvetic Alps, deriving from the proximal domain of the European margin and the Southern and Eastern Alps that represent proximal (Friuli-Ionian) to distal (Canavese) parts of the early Mesozoic Adriatic continental margin. (Schmid et al. 1989. The terms proximal and distal are used here with respect to Alpine Tethys. The foreland basins related to the orogen are the Po Basin and Veneto-Friuli Plain of the southern foreland and the Molasse Basin of the northern foreland. Also in the northern foreland, the Upper Rhine Graben formed as part of the European Cenozoic Rift System in the Eocene (Dèzes at el. 2004). The locations of all relevant tectonic features within the region can be found in Fig. 1.
Seismicity across the Alpine region largely occurs within the shallow upper crust along local faults (e.g. Thouvenot et al. 2016). Deeper events in the forelands are observed, however no consensus yet exists for their occurrence, with explanations ranging from thickened continental crust (Schmid and Kissling 2000) to slab steepening during the late collisional stage of orogenesis (Singer et al. 2014), suggesting the importance of both lithospheric structure and regional dynamics in the localisation of seismicity. Different styles of deformation are apparent along the European-Adriatic plate boundary. Ongoing indentation occurs in the Eastern Alps (Scharf et al. 2013;Restivo et al. 2016) where, geodetic measurements indicate that the rigid (i.e. mechanically stiff) Adriatic plate is currently rotating counterclockwise and indenting the Eastern Alps at ~ 2 mm/year (Vrabec and Fodor 2006;D'Agostino et al. 2008;Serpelloni et al. 2016). In the Central and Western Alps, where the convergence rate is lower (≤ 1 mm/year, Nocquet and Calais 2004; Tesauro et al. 2005), post-collisional extension and isostatic uplift occur (Sue and Tricart 2003). The relative impact of glacial isostatic adjustment versus tectonic and mantle dynamic processes, such as tearing of the Alpine slab beneath the Central and Western Alps (Fox et al. 2015(Fox et al. , 2016, as causes of the observed present-day kinematic state of the Alpine region is still debated (e.g. Norton and Hampel 2010;Chéry et al. 2016;Mey et al. 2017;Sternai et al. 2019). However, some indications from the present-day uplift rates suggest that the seismicity in the eastern Alpine domain could be instead related to still ongoing viscous relaxation from the waning of the last ice cap.
As such, whilst some seismicity within the Alps corresponds to the convergence of the European and Adriatic plates in North-East Italy (Chiarabba et al. 2005;Restivo et al. 2016), seismicity across the Alpine region is largely characterised by rather diffuse, shallow seismicity that does not necessarily follow major faults or plate boundaries (e.g. Deichmann 1992). Numerous destructive seismic events (Stucchi et al. 2012;Grünthal et al. 2013), such as the Magnitude 6.6 Basel earthquake in 1356 AD (Ustaszewski and Schmid 2007), the 1348 Villach earthquake (Reinecker and Lenhardt 1999) and events above the crustal root of the 1 3 Central Alps in southern Switzerland (Fäh et al. 2011) lie in areas of the orogen and forelands with low amounts of horizontal surface strain (Sánchez et al. 2018). These observations suggest that factors other than only convergence rate, such as the strength and rheology of the crust, affect earthquake activity in the Alpine region.

Structural model and thermal field
Lithospheric structure and density were taken from the model of Spooner et al. (2019a), where they were obtained by integrating previous geoscientific datasets in the region. Topography and bathymetry ( Fig. 1) came from ETOPO1 (Amante and Eakins 2009), and the seismically derived Lithosphere-Asthenosphere Boundary (referred to as LAB hereon) originates from Geissler et al (2010). Other subsurface lithospheric depths were constrained from the use of numerous published deep seismic surveys (e.g. Brückl et al. 2007;Hetényi et al. 2018), existing structural models of smaller subsets of the study area (e.g. Ebbing 2002;Przybycin et al. 2014;Freymark et al. 2017) and European plate crustal models (Tesauro et al. 2008;Molinari and Morelli 2011). Densities were constrained using forward 3D gravity modelling in IGMAS + (Schmidt et al. 2010(Schmidt et al. , 2020 and the global satellite gravity model EIGEN-6C4 (Förste et al. 2014), with the lithospheric layers split into domains of different density to account for lateral heterogeneity.
The resulting structural and density model (Spooner et al. 2019a), comprises an area of 660 km × 620 km ( Fig. 1) that is focussed on the central Alpine region, with a horizontal grid resolution of 20 km × 20 km. The vast majority of the Alps and their forelands are included, with the Central and Eastern Alps and the northern foreland being the best-covered regions, however the most extreme western and eastern portions of the Alpine chain are not covered as input data of a requisite quality were lacking (Spooner et al. 2019a). The model represents the highest spatial resolution among available models and conforms to both seismic and gravity-based observations. Five model layers that represent key lithospheric structural and density contrasts were used for the rheological calculations: (1) unconsolidated sediments (mostly Quaternary strata); (2) consolidated sediments (mostly Mesozoic strata); (3) upper crystalline crust (quartzite); (4) lower crystalline crust (gabbro); and (5) lithospheric mantle. The thicknesses of these layers and the domains of different density that they are comprised of are shown in Figs. 2 and 3c. Slabs and subduction interfaces are not considered within this work as no consistent model of their geometry or polarity throughout the whole study area is available at present, as debate is still ongoing in the Eastern Alps (e.g. Handy et al. 2021).
The temperature distribution input to the rheological calculations was obtained from a thermal model of the Alps and their forelands (Spooner et al. 2020a) generated from the same structural model utilised here (Spooner et al. 2019a). The 3D conductive steady-state thermal field was computed using the numerical simulator GOLEM, that can calculate coupled thermal-hydraulic-mechanical processes (Cacace and Jacquey 2017). Therefore, steady-state conductive heat transport was assumed as the main mechanism and specific thermal parameters were assigned to domains of the structural model, to fit a compiled dataset of measured sub-surface temperatures (Przybycin et al. 2015 and references therein;Freymark et al. 2017 and references;Trumpy and Manzella 2017). Depth maps of the calculated 200 °C, 400 °C, 600 °C, and 800 °C steady-state isotherms are plotted in Fig. 4.
The Alpine specific thermal field generated by Spooner et al. (2020a) was judged as preferential to European plate scale thermal models (e.g. Limberger et al. 2018) due to its higher degree of resolution in the resolved geological configuration of the plate. With the LAB exerting a firstorder control on temperature gradients calculated within these models, the definition of its depth and temperature is crucial to the resultant thermal field. The LAB utilised in Limberger et al. (2018) was derived from a European plate model by Tesauro et al. (2009), which has been smoothed, thus resulting in largely homogenous depths of 125 km for the forelands and 150 km for the orogen, whilst using a constant temperature of 1200 °C. The LAB varies more significantly in Spooner et al. (2020a) with depths taken from Geissler et al. (2010) ranging from 70 km beneath the URG to 140 km beneath the Adriatic Sea, and temperatures ranging from 1250 to 1400 °C converted from the shear wave velocity dataset of Schaeffer and Lebedev (2013). As such whilst lower temperature isotherms from Limberger et al. (2018) in the Alps (200 °C at 5 km) and Po Basin (200 °C at 10 km) are at a similar depth in Spooner et al. (2020a: Alps, 200 °C at 5 km;Po Basin, 200 °C at 9 km) higher temperature isotherms vary significantly (Limberger et al. 2018

Seismicity catalogue
Recorded seismic events from the ISC catalogue (International Seismological Centre 2020) were compared to the regional catalogues of Switzerland and Italy, regions with the highest seismicity within the study area, to assess the most suitable catalogue or combination of catalogues for use in this study. Figure 5 shows the spatial distribution of events from the ISC catalogue along with events from the SED catalogue in Switzerland (Fäh et al. 2011) and the INGV catalogue in Italy (Rovida et al. 2020). Events from each catalogue larger than magnitude 2 that occurred between November 2010 and November 2017 have been plotted, showing that the spatial coverage of all three catalogues is very similar across the whole region. Three distinct sub-sections were defined to better ascertain differences between the catalogues at depth. The first section compares the ISC and SED catalogues within Switzerland, the second section the ISC, SED and INGV catalogues along the Swiss-Italian border and the third section the ISC and INGV catalogues in Italy. The generated histograms, with a bin size of 500 m, show that in Switzerland and Italy, the ISC catalogue replicates the seismicity trends at depth of the respective regional catalogues, and portrays a greater number of events. In the border region, the two regional catalogues exhibit different trends however the ISC catalogue represents both trends, again with a larger number of events. As such, to maximise the number of seismic events available for interpretation, the ISC catalogue was used rather than regional catalogues.
Additionally, the depth distribution of events also indicates that depth fixing is an issue with the ISC catalogue. As the aim of this work is to compare strength at depth with the occurrence of seismicity, the ISC catalogue was filtered to remove events where fixed depths were assigned, where depth error estimates were absent or where the depth error was in excess of 20% of the event's hypocentre depth (allowing errors of < 3 km at a depth of 15 km), to both maximise the accuracy and number of useable events. The catalogue was further filtered to remove events smaller than magnitude 2 in an effort to remove the effects of smaller events 'clustering' around seismic stations observed in similar seismic catalogues (e.g. González 2016), whilst maintaining coverage across the entire study area. The threshold of magnitude 2 was selected after carrying out catalogue completeness analysis shown in Fig. 6, that shows the variability of completeness through time every 5 years starting from 1975. The magnitude of completeness for each period was estimated by two methods: the maximum curvature approach and goodness-of-fit test (please refer to Wiemer and Wyss (2000) for both), and the average value of two estimates was taken for each period. The estimated magnitude of completeness for the whole period from 1975 is around 1.5, while the average value for 5 years bins is 1.8. Therefore, the value of 2 represents a reasonable threshold that ensures the completeness from a statistical point of view. The resulting filtered seismicity catalogue contained 4,405 seismic events (shown in Figs. 7,9,11 and 12) from March 1964 to November 2017.

Method
The yield strength of the lithosphere (maximum differential stress achievable prior to failure (Goetze and Evans 1979)  Fig. 1a caption). The ISC catalogue is compared to the SED catalogue in a polygon cover-ing Switzerland in panel 1, it is compared to the INGV catalogue in a polygon covering Italy in panel 3, and it is compared to both the SED and INGV catalogues in a polygon covering the Swiss-Italian border in panel 2. Each histogram bin represents a depth of 500 m was calculated, taking into account the 3D temperature and lithostatic pressure of the system as derived from the structural and thermal models (Spooner et al. 2019a(Spooner et al. , 2020a. Rheological parameters, also used to calculate structural model layer strengths, were assigned based on the compilation of laboratory measurements (Ranalli and Murphy 1987) for the dominant lithology interpreted for each layer from observed seismic velocities as well as the modelled density and thermal properties. Parameters used can be found in Table 1. Long term lithospheric strength of the Alps and both forelands were then calculated with a vertical resolution of 100 m, assuming grain-size insensitive secondary creep  Byerlee's law (Eq. 1), a temperature-independent function of confining pressure resulting from increasing density and depth (Byerlee 1978;Ranalli 1995) was used to calculate the brittle portion of the yield strength where Δf f is the static friction coefficient (set to a constant value of 3 to represent lithospheric stress as per Brace and Kohlstedt [1980]), b is the bulk rock density, g is the gravitational acceleration, z is the depth (below surface) and f p is the pore-fluid factor (the ratio of the pore-fluid pressure to the lithostatic pressure, set here to a constant value of 0.36, representing a fluid density of ∼1000 kg m −3 and a rock density of ∼2750 kg m −3 ).
Power-law creep (Eq. 2), representing non-Newtonian, temperature activated deformation of rocks at increasing depth (Karato and Wu 1993;Burov 2007), was used to calculate ductile strength ( Δ d ): where ̇ is the strain rate [set to a constant 10 −15 s −1 in a compressive regime, consistent with observed GPS measurements from the region (Houlié et al. 2018)], A is the powerlaw strain rate, n is the power-law exponent, Q is the activation enthalpy, R is the gas constant, T is the temperature. Tests were made to account for the onset of low temperature crystal plasticity in the lithospheric mantle layer [Peierls creep mechanism (Katayama and Karato 2008)]. However, these tests suggested that including Peierls creep does not affect the strength of the plate within our model for the range of depths relevant for our analysis on earthquake occurrence.
For the visualisation of maximum rock strength under a constant strain rate at every point in the model, yield strength envelopes (referred to as YSE hereon) showing the differential stress envelope (minimum between Δ b − Δ d ) versus depth were calculated (Brace and Kohlstedt 1980). The modelled strengths of both the crust and the entire lithosphere were then vertically integrated over their entire thicknesses to compare lateral changes in strength throughout the region. From the above stated relationships, it follows that in our model rates by which strain energy is dissipated by viscous creep can be expressed in terms of a non-linear with temperature effective solid viscosity as ( eff ): The calculated lithospheric strengths and viscosities were then compared to the seismic event catalogue to investigate the relationships between lithospheric strength and seismicity throughout the region.

Results
Across the Alps and their forelands, the pattern of variations in integrated strength of the entire lithosphere (shown in Fig. 7a) corresponds closely to the pattern of Moho depth (Fig. 3b). The weakest lithosphere (13 log 10 Pa m) occurs at the deepest portion of the Alpine crustal root, the largest Moho depth in the study area (55 km), west of the Giudicarie Fault. Similarly, around the Tauern Window, the eastern Alps are characterised by both a shallower Moho (45 km) and higher strengths (13.2 log 10 Pa m). In agreement with this correlation, both forelands exhibit significantly shallower Moho depths and higher lithospheric strengths than within the orogen. The lithosphere of the southern foreland was also found to be stronger than the northern foreland, with the highest strengths in the study area exhibited on the Adriatic plate in a region south of the Veneto-Friuli Plain (13.9 log10 Pa m). Similar results have been observed by previous works modelling the lithospheric strength of the eastern Alps along the Transalp profile (Willingshofer and Cloetingh 2003) and in the central Alps along the EGT profile (Okaya et al. 1996) and the entire orogen (Marotta and Splendore 2014).
Not all zones of high strength observed in the southern foreland correspond to shallow Moho depths. The strong Adriatic plate has some of the largest Moho depths encountered in either foreland at 40 km deep, south of the Veneto-Friuli Plain. Instead this region is observed to correspond to the deepest portion of the LAB (140 km, shown in Fig. 3d) and to significantly colder lithospheric temperatures in the Adriatic plate (shown in Fig. 4). The integrated crustal strength (shown in Fig. 7b) also positively correlates with temperature, being highest (13.175 log 10 Pa m) in the South East of the Adriatic plate where all isotherms are deepest (Fig. 4, around the annotation for the Adriatic Sea). North of the Periadriatic Lineament and West of the Giudicarie Fault, isotherm depths are consistently shallower than in portions of the orogen South and East of these locations. Consequently, the northern and western Alps feature lower crustal strengths (13 log 10 Pa m) than the Alpine crust in the south and east (13.125 log 10 Pa m) a finding also observed by Marotta and Splendore (2014). Where isotherms are raised below the Upper Rhine Graben (further referred to as URG), crustal strengths are also the lowest (12.925 log 10 Pa m).
The spatial distribution of seismic event epicentres in the southern foreland strongly correlates to the spatial trend of the computed integrated lithospheric strength (Fig. 7a) and not to the spatial trends of the computed crustal strengths (Fig. 7b), with most seismicity occurring in regions where the integrated lithospheric strength is weak, irrespective of the crustal strength. Within the northern foreland no correlation is observable between seismicity localisation and lithospheric strength, with epicentres instead corresponding closely to the weaker portions of the crust (Fig. 7b) around the URG in the west of the European foreland. Both integrated lithospheric and crustal strength maps portray lower strengths within the Alps proper, North and West of the Giudicarie Fault, corresponding to the location of the majority of Alpine seismicity.
Cross sections showing variations in strength differential (minimum between brittle strength and ductile strength) with depth and the relation to seismicity are plotted in Fig. 8. In line with previous work (Okaya et al. 1996;Willingshofer and Cloetingh 2003;Marotta and Splendore 2014), all cross sections show that most seismicity occurs within the strongest region of the upper crust (corresponding to ~ 1 GPa from our calculations), between 200 and 400 °C, with a largely aseismic and weaker lower crust mechanically decoupling the crust from the lithospheric mantle. Seismicity that occurs deeper than the upper crust of the model (~ 20 km depth) is however present in regions where the upper lithospheric mantle is cooler than 600 °C and strong (calculated as > 2 GPa), as shown in cross Sects. 1 and 3 (beneath the southern URG, the Po Basin and the Adriatic Sea). Regions seen in Sect. 1 (also beneath the southern URG) characterised by a stronger lower crust (corresponding to calculated values of ~ 1 GPa) also show seismicity to a greater depth. Section 4, which runs West to East through the centre of the orogen, does not portray a strong lower crust or upper lithospheric mantle, and shows Moho temperatures consistently higher than 600 °C, exhibiting no seismicity outside of the upper crust.
The pattern of variations in integrated effective viscosity of the lithospheric mantle (shown in Fig. 9) are the same as the integrated strength distribution of the whole lithosphere (Fig. 7a), with higher strengths analogous to higher viscosities, corresponding closely to the Moho depth (Fig. 3b). Seismicity is therefore observed in the orogen and southern foreland to spatially correlate with regions of lower viscosity, such as the Alpine root west of the Giudicarie Fault (trough of greatest Moho depth [20.6 log 10 Pa s]), that surround higher viscosity blocks like the Adriatic Plate (22.2 log 10 Pa s). In the Northern Foreland the opposite is observed, with seismicity occurring in the region of highest viscosity surrounding the URG (20.8 log 10 Pa s). Similar results, showing that viscosities below the orogenic root must be lower than beneath the forelands have also been modelled by Chéry et al. (2016), to fit the observed trends of glacial rebound.
Variations of effective viscosity with depth in relation to seismicity are shown in the cross sections on Fig. 10. As previously observed, most seismicity occurs in the upper crust, which largely corresponds here to effective viscosities of 23.5 log 10 Pa s or higher. Viscosities for the lithospheric mantle tend to be between 19 and 23 log 10 Pa s and for the lower crust between 21 and 23 log 10 Pa s. Across the majority of the region, both the lower crust and lithospheric mantle appear largely aseismic. Regions where seismicity does occur deeper than the upper crust also correspond to regions of the lower crust or lithospheric mantle where effective viscosities are also in excess of 23.5 log 10 Pa s, such as beneath the Po Basin and Adriatic Sea in Sect. 1 and beneath the southern URG in Sect. 3.

Mechanical strength
The strongest regions at depth in the study area correlate with the 600 °C isotherms within the upper lithospheric mantle, corresponding to the onset of crystal plasticity in mantle rocks (Boettcher et al. 2007;McKenzie et al. 2005). Additionally, the map of the integrated lithospheric strength Fig. 8 Cross Sects. 1, 2, 3 and 4 (Locations shown in Figs. 2, 3, 7 and 9) showing the variation of strength differential with depth throughout the lithosphere. The top of the upper crust, lower crust and Moho are shown as solid black lines. Isotherms up to 1000 °C are overlain as dashed blue lines and seismicity > M2 that lay within 20 km distance of the section has been marked as red dots. Abbreviations used to annotate the cross sections are shown in Fig. 1a caption, with the addition of: D: Dinarides ◂ (Fig. 7a) portrays a positive inverse correlation with respect to the Moho depth (Fig. 3b), allowing Moho depth throughout the region to be used as a first-order estimate for the relative strength of the whole lithosphere. Based on these findings we can conclude that the lithospheric mantle is both the largest contributor to the overall computed lithospheric strength variations and is also highly influenced by the temperature configuration across the entire orogen foreland system, thereby expanding upon previous findings derived from, but limited to, the central Alps (Okaya et al. 1996).
The crust accounts for only a small portion of the total lithospheric strength throughout the majority of the study area, as can be observed in Fig. 11. The foreland crust in particular only contributes to the overall lithospheric strength in a limited way (~ 30%). However, in locations where the crust is at its thickest, such as the crustal root of the orogen, west of the Giudicarie Fault (Fig. 3a), the pattern of crustal strength distribution can equal the whole plate strength distribution, similar to findings of Tesauro et al. (2009). The presence of a weak and thick crust North of the Periadriatic Lineament contributes to a significant weakening of the lithosphere underneath this domain, a feature that was also noticed by previous work (Marotta and Splendore 2014). Similar to results of our work, weaker lithospheric strengths are also noted to occur beneath thickened crust by Mareschal and Jaupart (2013). Indeed, the crustal thickness map shown in Fig. 3a corresponds very closely to the percentage of crustal strength, indicating that crustal thickness is important for the relative contribution of the crust to lithospheric strength. This relationship is to be expected as an increase of the crustal thickness means an increase of temperature at the base of the crust and thus in the mantle lithosphere. However, as shown in Fig. 7b, the thermal field is more important at determining the strength within the crust itself.
The temperature distribution throughout the region is primarily a function of lithospheric composition, with the relative contribution of model layers to the heat budget varying in response to their specific thermal properties and relative volume. Features such as a higher percentage of felsic (radiogenic) upper crust to mafic lower crust or the presence of thick insulating sediment deposits have been previously shown to raise the geothermal gradient within the study area (Spooner et al. 2020a), however the geothermal gradient is largely controlled by larger tectonic effects, that are reflected in the height of the thermal LAB, here derived from a global tomographic model (Schaeffer and Lebedev 2013). We therefore expect that these specific features would also exert a first-order control on the resulting mechanical configuration of the lithosphere and on the seismicity distribution.

Relation to seismicity
To further deepen the discussion about the implications derived from the thermorheological model on the seismicity distribution within the area we also note that within the Alpine domain, seismic events mostly occur North of the Periadriatic Lineament and West of the Giudicarie Fault, where the crust is both warmer and weaker. This relationship between lithospheric strength and seismicity distribution is part of a broader observed trend of West to East mechanical strengthening along the Alpine crust. This results in having significantly less seismicity East of the Giudicarie Fault, which we can explain by a deepening in the topography of the LAB moving eastward (Fig. 3d) and therefore a lower geothermal gradient.
The northern and southern forelands of the Alps display a markedly different pattern of seismicity in terms of their epicentre locations, potentially reflecting the different tectonic settings and driving mechanisms at play within each foreland. In the southern foreland seismicity primarily occurs at the boundaries of the European and Adriatic plates (e.g. Chiarabba et al. 2005). These locations mark a relatively sharp gradient in modelled lithospheric strength and effective viscosities from 13.9 log 10 Pa m and 22.2 log 10 Pa s within the plate to 13.2 log 10 Pa m and 20.8 log 10 Pa s along its edges, respectively. These mechanically stiff and rheologically strong lithospheric blocks move independent of one another (e.g. Nocquet and Calais, 2004). Therefore, it stands to reason that seismicity in such a setting would localise at the boundaries of these  (Figs. 7a and 9), as it has been also identified in previous works (Craig et al. 2011;Sloan et al. 2011;Tesauro et al. 2015).
The situation differs in the northern foreland, where seismicity occurs within an intraplate setting (e.g. Bonjer 1997) of the European plate and where the upper mantle is not seismogenic. It therefore seems logical to assume that under these tectono-thermal conditions the weaker regions of the crust would accommodate the majority of seismicity as depicted in Fig. 7b, where the lithospheric mantle remains relatively stable throughout the northern foreland.
These results are indicative of a structural correlation between the rheological configuration of the crust and the level of upper crustal seismicity. Weak crustal domains (northern Alpine foreland), or domains of rheological contrasts as in the southern foreland, represent zones of localised deformation, which in turn would promote rock failure at relatively low levels of accumulated stress. Depending on the local fault properties, such localised deformation would lead to a concentration of seismic activity in these areas, where faults would preferentially rupture with moderate magnitudes. It follows that, under such conditions, a rapid accommodation of permanent deformation by thermally activated viscous creep would efficiently limit the deepening of seismic energy release in these areas and therefore seismicity would be likely bounded to the shallower and weaker portions of the upper crust.
It is nowadays established that temperature is an important variable for determining the depth of the lithospheric seismogenic domain. This was noted in the seminal study of Sibson (1982) that demonstrated a correlation between intraplate seismicity focal depths and surface heat flow distribution, with shallow seismicity in areas of high surface heat flow and vice versa. The existence of an inverse correlation between heat flow and focal depth is easy to explain when considering the homogeneous configuration of oceanic plates. However, it is challenging to extrapolate these results to continental intraplate regions where the log-linear age-temperature relationship does not apply and the thermal state is a complex function of the history of a heterogeneous plate over geological times. In such cases, a conservative estimate for the lower bound of the seismogenic zone can be derived based on a description of the continental plate characterised by a non-linear creep law for its most abundant minerals. By relying on such a description, a maximum in the static strength would correspond to a particular value of the system's internal energy, which can then be described by the temperature at its maximum dissipation potential. Under this assumption, one would expect that to a first-order within the study area, the 600 °C isotherm, which represents the onset of crystal plasticity in ultramafic rocks (Boettcher et al. 2007;McKenzie et al. 2005), would mark the lower bound of seismicity, being associated with the highest strength in the upper lithospheric mantle and, possibly, in the lower crust.
Most observed seismicity occurs between the 200 °C and 400 °C isotherms, representing the strongest portion of the upper crust (up to ~ 1 GPa). Deeper in the crust, higher temperatures result in a gradual decrease in the plate strength and subsequent increase in its aseismic behaviour. These observations can be taken as indicative of the brittle-ductile transition within the crust, providing a conservative estimate to the lower bound of the seismogenic zone in that area. Willingshofer and Cloetingh (2003) performed an end-member sensitive analysis for the lithospheric strength along the Transalp section of the eastern Alps by considering either a strong or a weak crust. The main conclusion derived from their study was that a model with a strong crust provides a better fit with the seismicity. From their analysis they determine that the brittle-ductile transition occurs between 9 and 14 km.
By quantifying the maximum depth of seismicity across the region we are also able to comment on the likely depth of the brittle-ductile boundary. Figure 12a shows the density of seismic data available throughout the region, namely the number of seismic events within a 75 km radius around every grid point. Regions of very high event coverage such as the Po Basin and Apennines imply more robust estimates for the maximum depth of seismicity compared to regions of very low coverage such as the eastern portion of the European plate. By plotting the depth level above which 95% of seismic events occur (Fig. 12b), our results, based on a more up-to-date 3D structural and thermal model, estimate that the brittle-ductile boundary within the eastern Alps occurs on average at ~ 20 km depth, providing a better fit with the depths of recorded seismicity and adding validity to the model setup utilised here. Figure 12b also shows that the maximum depth of seismicity in the western Alps is slightly shallower, at ~ 15 km depth, whilst the maximum depth of seismicity in the northern foreland is slightly deeper at ~ 25 km. Maximum depths of seismicity within the southern foreland are much more variable, from ~ 25 km beneath the Veneto-Friuli Plain to ~ 60 km beneath the western Po Basin.
Additionally, whilst the majority of seismicity occurs within the upper crust of the model (~ 20 km depth), observations also show that seismicity in deeper layers occurs in both the northern and southern forelands (e.g. Bonjer 1995;Chiarabba et al. 2005). From the strength difference cross sections shown here (Fig. 8), it can be discerned that hypocentre depths of deeper seismicity vary with temperature, occurring where upper lithospheric mantle temperatures are ~ 600 °C or cooler. This thermal configuration leads to the presence of a relative weak lower crust, mechanically decoupled from and sandwiched between the strong upper crust and a strong upper mantle (> 2 GPa) visible in Fig. 8, thus providing a first-order explanation to the deepening of seismicity in the area (see also the discussion in the next paragraph in this regard). Seismicity within the lower crust is found to only occur in domains of higher strength (~ 1.4 GPa), though within the majority of the region the lower crust is estimated to be largely weak and therefore aseismic. This is also evident from the effective viscosity cross sections (Fig. 10) depicting how the lower crust mechanically decouples the upper crust from the lithospheric mantle, with effective viscosity values (22.5 log10 Pa s) indicative of a ductile regime and fitting well with the general lack of observed lower crustal seismicity.
The unimodal pattern of seismicity beneath the Alps being limited to the upper crust contrasts to the bimodal (crust and upper mantle) seismicity pattern found in other orogens worldwide, such as the Himalaya. Based on this observation, a recent modelling study by Dal Zilio et al. (2018) advanced the hypothesis of a structural correlation between plate convergence rates and seismicity distribution alone. The analysis of Dal Zilio and co-authors is based on a linear correlation between convergence rates and the resulting thermal configuration of the orogen, with faster rates resulting in a colder orogen and therefore in a more widespread seismogenic brittle domain. A major limitation of this reasoning is that it does not take into account the role of crustal inheritance. Whilst results presented here mostly agree with the unimodal nature of seismicity throughout the crust, they also suggest that the lithospheric makeup of the region such as crustal, and lithospheric thickness have a first-order effect on the location of seismicity in the region via their control on the lithospheric thermal field. In addition, relating the overall distribution of seismicity within the Alpine region to the background tectonic convergence rates cannot reconcile the diversity in the observed seismic style across the whole orogen. While the western Alps are characterised by seismicity showing mainly extensional faulting, the eastern and central Alps show mainly compressional seismicity. Therefore, it would be difficult to relate this difference in style to a common geodynamic process. Indeed, there is evidence from present-day uplift rates that the seismicity in the western Alpine domain could be instead related to still ongoing viscous relaxation from the waning of the last ice cap (e.g. Norton and Hampel 2010;Chéry et al. 2016;Mey et al. 2017;Sternai et al. 2019). In this regard a word of caution comes from the uncertainty in the mantle and (lower) crustal viscosities input in these studies. Past works have adopted values as low as ~ 20 log 10 Pa s for mantle viscosity (Norton and Hampel 2010;Chéry et al. 2016;Mey et al. 2017) whilst more recent work has made use of higher values ~ 22 log 10 Pa s for the Alpine lithospheric mantle (Sternai et al. 2019). The values derived in the present study, ranging from ~ 21 log 10 Pa s beneath the orogen to ~ 22 log 10 Pa s in the forelands compare favourably to those estimates. In addition, our model favours the presence of lower viscosities below the orogen domain proper in comparison to below the forelands, a result that agrees with the main conclusions derived from Chéry et al. (2016). This last result confirms how the pattern and style of seismicity within the complex Alpine area cannot be related to a single geodynamic parameter, such as convergence rates, but should be considered as a natural outcome of a rather complex crustal structure developed during the whole orogenic cycle in an ongoing plate tectonic setting.

Upper mantle seismicity-slab influence and other controlling mechanisms
The temperature at the maximum depth of seismicity (Fig. 12) is plotted in Fig. 13a, showing that most mechanically strong portions of the plate, whether within the crust or the upper lithospheric mantle, are effectively bounded by the depth of the 600 °C isotherm, previously discussed to represent the maximum temperature of seismicity. This trend is also visible in the various cross sections throughout the region (Fig. 8). We also note that there are regions where the maximum temperature of seismicity exceeds 600 °C, in domains where ductile conditions would be expected (North-East of the annotation for the central Alps, North of the annotation for Tauern Window, and running South-East from the Western Alps to the Adriatic Sea).
The majority of these domains correspond to the presence of actively subducting and/or previously subducted slabs (high velocity features at 100 km depth in Fig. 13b from El-Sharkawy et al. 2020), though some debate still persists in the community . The thermal field adopted in this study (Spooner et al. 2020a) neglects the cooling effect from active subduction dynamics, which remains largely unquantified in the study area. Therefore, the possibility of seismicity occurring in these regions at higher temperatures than expected is anticipated, as previously discussed in Spooner et al. (2020a). This effect is most pronounced at the location of the subducted Apennine slab where maximum temperatures of seismicity appear to approach 1000 °C. Seismicity above 600 °C is also mapped beneath the Alps, corresponding to the location of Alpine slabs, suggesting that even these older frozen-in subduction zones may still be having a thermal effect on the lithosphere, which remains unconsidered in our calculations. This could offer alternative explanations for observations of Fig. 10 Cross Sects. 1,2,3 and 4 (Locations shown in Figs. 2,3,7 and 9) showing the variation of effective viscosity with depth throughout the lithosphere. The top of the upper crust, lower crust and Moho are shown as solid black lines. Seismicity > M2 that lay within 20 km distance of the section has been marked as red dots. Abbreviations used to annotate the cross sections are shown in Fig. 1a caption, with the addition of: D: Dinarides ◂ lower crustal seismicity to high temperatures (> 600 °C) in northern Switzerland, having been previously correlated to active lower crustal fluid flow by Deichmann (1992) or to slab rollback by Singer et al. (2014).
Another possible explanation for the occurrence of the observed deep seismicity, can be due to the transfer of tectonic stresses into deeper domains through a lower crustal and upper mantle zone of weakness. For such a transient to grow into a full-scale instability it requires the stress drop to overcome the shear resistance of the surrounding shear zone, which could be easily attained if background stresses approach the local fault strength at depth. We thereby propose two additional mechanisms that can be responsible for such a process, that is: (1) ongoing lithosphere (re)strengthening from grain-size heterogeneity in concert with slab induced cooling promoting localised brittle instabilities which would subsequently propagate in the neighbouring ductile portions of the lower crust and mantle; and, (2) the occurrence of transient weakening mechanisms (e.g. thermal shear runaway and/or ongoing dehydration reactions along the developed shear band) in the ductile regime which would allow such deep earthquakes to initiate and propagate within a ductile regime. In the former case, seismicity would be caused by brittle failure within a cold (less than 600 °C) and possibly anhydrous mantle, where the presence of local heterogeneities within the fault would promote instabilities to nucleate and propagate into otherwise ductile mantle rocks. In the latter case, nucleation of earthquakes in the lower crust and/or lithospheric mantle would not be limited by brittle conditions, but would rather occur by ductile deformation. This would allow seismicity to occur at and above temperatures of 600 °C and would be consistent with the thermal configuration adopted in our study.

Workflow limitations
The thermal field exerts a first-order role in the long-term strength and deformation modes of the lithosphere. In this study, similar to other recent thermal models of the region (Limberger et al. 2018), we relied on the working assumption of steady-state thermal diffusion, thereby neglecting the effect that features such as sedimentary erosion and deposition within the basins could have on the thermal field. Unfortunately, whilst the omission of such features may influence the computed lithospheric strength, a lack of constraints from across the region through time, hinder the adoption of such a transient modelling analysis, which would be subjected to large extrapolation and uncertainties (please also refer to the discussion in Spooner et al. 2020a).
Lateral variations in the flow law parameters have also been not considered. Our choice stems from the still unresolved underlying issue on how to accurately constrain flow law parameters derived from laboratory experiments to the scale of the lithosphere, where uncertainties also exist on the effective boundary conditions at depths. Given the powerlaw dependence of differential stress on strain rate and the exponential dependence on temperature, values derived from monomineralic laboratory experiments cannot be simply upscaled to natural settings where the strain rates vary over ten orders of magnitude. Whilst some work has been published to attempt to address these issues (e.g. Handy et al. 1999Handy et al. , 2001Schrank et al. 2008), further work is still required to accurately constrain and generate higher resolution rheological models. Therefore, in an attempt to avoid introducing further uncertainties to the model, homogenous flow law parameters were used for each layer.
Lithospheric strain rates are known to be depth dependent and have been quantified in certain parts of the study area, such as in the URG where surface strain rates are shown to be lower and mantle strain rates higher (Heckenbach et al. 2021) than what we have assumed. However, accurate constraints on their values across the entire region are not presently available as no geodynamic model exists that can replicate the pattern of surface strain indicated by geodetic measurements. For this reason, in the present study we refrain from making use of variations in strain rates that would add bias on the final results, deciding instead to make use of an average distribution that still is considered representative of the bulk deformation velocities in the region.

Summary
In this work we have computed the long-term lithospheric yield strength for the Alpine regions and its forelands using available up-to-date structural, density and thermal input data. Variations in the strength of the upper lithospheric mantle exert the largest influence on the strength of the whole lithospheric column, with crustal strengths only contributing significantly to the whole plate integrated strength beneath the orogen proper, where the crust is thickest (55 km). The strengths, whether in the crust or mantle, are largely temperature dependant, with upper lithospheric mantle temperatures controlled by crustal and 1 3 lithospheric thickness, whilst crustal temperatures are controlled by thermal LAB depth, thickness of the radiogenic felsic upper crust and thickness and distribution of insulating sediments.
The results from the thermorheological modelling exercise has shed light on the relationship between background seismicity and resolved lithospheric strength variations in the region. We have been able to demonstrate how the occurrence of crustal seismicity in the study area is influenced by several factors acting at different scales with inherited geological crustal and upper mantle structures exerting a primary control in the seismicity distribution and style. The highest yield strengths within the crust (~ 1 GPa) and upper mantle (> 2 GPa) occur at temperatures interpreted as the onset of crystal plasticity (crust: 200-400 °C; mantle: ~ 600 °C) with almost all seismicity occurring in these regions. We also note the presence of a weak lower crust (< 1 GPa) that mechanically decouples the upper crust and lithospheric mantle across the entire region. Therefore, the lower crust appears largely aseismic, likely due to seismic energy being dissipated by ongoing creep in regions where effective viscosities are lower. Subducted slabs (Apennine and Adriatic) are also shown to significantly influence the maximum depth of seismicity possible above them, furthering the argument that lithospheric heterogeneity plays a key role in controlling the regional distribution of seismicity.
In the Alps, seismicity correlates spatially with a weaker crust and lithosphere, such as north west of the Periadriatic Lineament. Such a clear distinction cannot be derived uniquely for the northern and southern forelands, as each shows a different pattern of seismic distribution, likely reflecting their different tectonic settings. In the southern foreland seismicity preferentially occurs across boundaries between rigid lithospheric blocks, such as the strong Adriatic plate, whilst in the northern foreland seismicity localises beneath domains of crustal weakness as in the URG.
Author contributions MC and DA developed tools for calculations in the modelling workflow and plotting of results. MS-W advised on the entire workflow and the interpretation of results. CS carried out the modelling work and prepared the paper with contributions from all co-authors.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors would like to thank the Deutsche Forschungsgemeinschaft (DFG) for funding the Mountain Building Processes in Four Dimensions (4-D-MB) SPP that this work was produced as part of.
Code Availability GOLEM, the numerical calculator used throughout this work is publicly available and can be found at https:// doi. org/ 10. 5281/ zenodo. 999401.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.