Tectonics and seismicity in the Northern Apennines driven by slab retreat and lithospheric delamination

subduction/collision system and to identify a set of best-fit model parameters allowing to capture the present-day spatial distribution of compressional and extensional normal stresses and of surface velocities. This is achieved by varying the rheology of the lower crust and the temperature of the slab, the asthenospheric wedge and of the lowermost crust of Adria. Once the location of present-day tectonic regimes is broadly reproduced, the resulting seismicity and surface velocities are analyzed in comparison with observations. Subsequently, we study how these short-term


Introduction
The geodynamic evolution and seismicity of many orogenic systems in the Mediterranean domain, such as the Apennines, Carpathians, Hellenides or the Betic-Rif systems, are driven by the rapid Cenozoic retreat of genetically associated slabs (e.g., Vergés and Fernàndez, 2012;van Hinsbergen et al., 2014;Faccenna et al., 2013;2014a,b;Horváth et al., 2015;Balázs et al., 2017;Jolivet and Burn, 2010;Jolivet et al., 2013). The evolution of such an orogen is typically characterized by the gradual migration towards the foreland of contraction in the external part and extension in the internal part, driven by slab retreat (Bertotti et al., 2006;Picotti and Pazzaglia, 2008;Leever et al., 2006). The extension is defined as back-arc, although often there is no stable magmatic arc observed (see discussion in Jolivet et al., 2013). Numerous reconstructions, geodynamic modelling studies of subduction dynamics and magmatism, as well as extensive analysis of present-day seismicity and surface to deep mantle observations are readily available in The link between lithospheric-scale structure and geodynamics and the resulting crustal tectonics has already been the subject of global observational studies. For instance, maximum earthquake size in subduction zones was suggested to correlate with convergence rate and lithospheric age (and thus temperature and slab pull) (e.g., Ruff and Kanamori, 1980). Furthermore, greater coupling between the slab and the rest of the subducting plate at the ocean floor was observed to correspond to less seismic moment release (Bilek et al., 2005). This is thought to reflect the effect of strength and structural integrity of the subducting plate on the localization of seismogenic brittle deformation. However, statistical correlations are often weak (e.g., Heuret et al., 2011). Furthermore, such studies are often limited by the short observational period compared to the recurrence times of earthquakes and by the difficulty of studying systems governed by multiple, interacting processes occurring at depths of tens to hundreds of kilometers and across various spatio-temporal scales. Modeling is thus needed to reveal the physical processes connecting deep and shallow deformation over different time scales (e.g., Dinther et al., 2013b;Dal Zilio et al., 2018).
To investigate the mechanisms controlling the regional tectonics and to address the limitations of existing research, this study considers long time periods and links small-scale tectonics and seismicity to large-scale dynamics and structure using a seismo-thermo-mechanical (STM) modeling approach (van Dinther et al., 2013b). This method can simulate long-term visco-elasto-plastic deformation of the lithosphere, dislocation creep and ma ntle flow together with short-term brittle-plastic failure and fault slip associated with accumulated elastic stress release. This ensures that faults can for the first time be loaded tectonically by velocities and stresses resulting from slab pull. Our model setup is thus based on the hypothesis that lithospheric dynamics driven by slab pull (i.e., with no imposed shortening or extension) can largely explain the active tectonics and seismicity of the modeled retreating subduction/collision system.
We performed a series of numerical experiments to both investigate the retreating subduction/collision system and to identify a set of best-fit model parameters allowing to capture the present-day spatial distribution of compressional and extensional normal stresses and of surface velocities. This is achieved by varying the rheology of the lower crust and the temperature of the slab, the asthenospheric wedge and of the lowermost crust of Adria. Once the location of present-day tectonic regimes is broadly reproduced, the resulting seismicity and surface velocities are analyzed in comparison with observations. Subsequently, we study how these short-term J o u r n a l P r e -p r o o f Journal Pre-proof features are affected by slab pull and lithospheric rheology.

Seismo-thermo-mechanical modeling
The STM modeling approach was developed, described, and first applied by (van Dinther et al. (2013b). This approach is based on a geodynamic numerical code that uses a fully staggered grid, conservative finite-difference scheme with marker-in-cell technique to solve for the conservation of mass, momentum and heat (Gerya and Yuen, 2007). To do so a Maxwell visco-elasto-plastic rheology is used. Drucker-Prager plasticity approximates both plastic and brittle yielding in a continuum mechanics framework and consists of plastic strain accumulation upon reaching a pressure-dependent yield strength yield  .
where G is the shear modulus, D Dt the objective co-rotational time derivative, the second invariant of the deviatoric stress tensor, and  a plastic multiplier which connects strain rates to stresses during yielding (Gerya and Yuen, 2007). Non-plastic effective viscosity  is computed as a function of stress according to a non-linear flow law with dependence on pressure P and temperature T : with R the universal gas constant, n the exponential coefficient, D A the pre-exponential factor, a E and a V the activation energy and volume (Ranalli, 1995 . To allow the spontaneous nucleation, propagation and arrest of seismic events, friction depends on local slip rate V according to a strongly slip rate dependent friction formulation (e.g., Cochard and Madariaga, 1994;Ampuero and Ben-Zion, 2008), in which a relation with V in the denominator is thought to represent friction at high seismic slip rates (e.g., Di Toro et al., 2011a)

Model setup and observational constraints
The reference model setup in Fig. 2 combines a suite of geological and geophysical observations.
We design the initial configuration on the basis of a geological and structural profile of Molli et al. (2010). We selected this profile because it synthesizes various geological and geophysical observations into a schematic yet detailed cross-section extending to mid-crustal depths. The faults in the profile are transposed into the model as thin zones with their own rock type. Their rheology is significantly weaker than the other rock types, and thus more prone to brittle-plastic yielding.
However, yielding and seismicity can and do also occur away from prescribed faults. At depth, the profile is supplemented by published information. In particular, we adopt the Moho geometry from a cut through the model of Spada et al. (2013), based on seismic reflection and refraction as well as receiver functions. In the setup, the lithospheric mantle and thus the Moho consist of two distinct and partly overlapping segments (Adriatic and Tyrrhenian), in accorda nce with the receiver function study of Bianchi et al. (2010). The Adriatic lower crust in the setup follows the Moho downwards and is shaped accordingly with the low-velocity zone tomographically imaged at moderate depths (35-55 km) beneath the middle and southwestern side of the range, above the Adriatic mantle (Di Stefano et al., 2009). The lithosphere-asthenosphere boundary is located 100-120 km beneath the southern Po Plain and northernmost Apennines and around 70 km beneath the eastern coast of the Ligurian sea, in accordance with the surface wave tomography of (Panza et al. (2003). We add a subvertical slab of lithospheric Adriatic mantle down to depths of 400 km, following the teleseismic tomography of Benoit et al. (2011). The uncertainty in the depth extent of the lithospheric mantle is implicitly addressed in the different temperature setups used ( Fig. 3 and Section 3.3), since it is temperature that distinguishes the lithospheric and asthenospheric mantle in terms of rheology and density (Tables 1 and 2). The reference set of material properties corresponding to each rock type is shown in Table   1. Three different lower crustal rheologies are tested, as their ductility is expected to control the style of deformation. Following Faccenda et al. (2009) and Ranalli (1995), these are mafic granulite, plagioclase, and wet quartzite.    From the corresponding generic rock types from Clauser and Huenges (1995).  conditions are also shown. The thin grey lines, shown as a spatial reference, are material contours.

Model specifications
The boundary conditions used are free-slip on all boundaries and zero lateral P gradient at the four corners (with =0 P at the top). Lateral boundaries are zero-heat-flux, while the top has fixed T (0 °C) and a fixed T gradient (2.413 °C km −1 ) at the bottom. A 25-km-thick weak "sticky air" layer ( =0 yield  , =1  kg m −3 , 17 = 10  Pa s) was placed at the top to approximate a free surface J o u r n a l P r e -p r o o f (Crameri et al., 2012). The models have a variable spatial resolution, from 250 m in the area of interest to 5 km at the edges and bottom. To minimize boundary effects, the model space extends vertically for 500 km and horizontally for 1000 km.

Modeling procedure
Following Dal Zilio et al. (2018), we run the models in two phases: (i) long-term, with timestep yr. The long-term phase allows the model to develop self-consistent stresses and velocities and attain isostatic equilibrium. Once the stress configuration and crustal velocities are broadly consistent with present-day observations, the timestep size is progressively reduced to 1 yr. We then switch to the short-term phase, during which the model is run for 20, 000 years and spontaneous events are simulated.

Evaluation of model seismicity
Due to computational challenges, we use a constant timestep of 1 yr during the short term modeling phase. Hence, seismic slip rates cannot be achieved. For the purpose of this paper, we focus on the type, distribution and size of events and do not attempt to realistically simulate the rupture process. To detect events, we record all occurrences of relative displacement ( 3.

Lithospheric and crustal dynamics
The present-day distribution of velocities and stresses in the study area of the northern Apennines can be broadly reproduced by our reference model through the combination of high temperatures in the mantle wedge and lowermost crust and the most ductile (wet quartzite) lower crust rheology ( Fig. 3 and 4). These parameters allow mantle material to protrude upward and northward at relatively high velocity ( 20 mm yr −1 ) into the lower crust. As a result, a hot, low-viscosity, increasingly weak channel forms above the Adriatic lithospheric mantle. This decouples the overlying upper crust, which includes the mountain range, from the slab and pushes upward and laterally with respect to the Tyrrhenian domain further to the southwest, causing uplift and extension. Conversely, the portion of the upper crust ahead of the low-viscosity channel, including the Po Plain, is compressed by the protruding material. It is also still mechanically coupled to the partly delaminated, sinking lithospheric mantle slab and therefore subsides. Two tectonic regimes thus form in the upper crust: extension and uplift throughout the range versus compression and subsidence in the foreland (Fig. 4). This is in agreement with present-day observations.
In contrast, the use of the strongest (mafic granulite) lower crust rheology, produces a rheological coupling of the slab to the upper crust, regardless of temperatures. This leads to tectonic regimes opposite to the observed in the study area, with thrusting and subsidence in the range and extension and uplift in the Po plain (Fig. 5).
On the other hand, the intermediate (plagioclase) rheology of the lower crust allows for wedge protrusion and lithospheric delamination only when the lower crust is brought to temperatures of over 450°C. However, normal stresses in the range are partially compressional, the deep area 20-35 km underneath the base of the mountain range is not entirely compressional, and the foreland beyond the external thrust front undergoes intense extension and uplift (Fig. 5).
With further modeling, we attempt to produce realistic tectonic regimes with less ductile lower crust, such as by imposing a zone of eclogitic material or a localized shear zone at the deep suture between the two lithospheric domains. These efforts did not lead to compression throughout the Po plain without significant changes in the defined initial crustal and lithospheric geometry.
Therefore, the simulations suggest that the simultaneous presence of the two tectonic regimes

Seismicity and surface deformation
We analyze the short-term behavior of the reference model in terms of surface displacement and seismicity.

Surface velocities
Horizontal model velocities show less than 2 mm yr −1 of extension distributed throughout the external half of the mountain range and 3 mm yr −1 of convergence across the whole of the Po plain. A similar trend is present in GPS measurements by Devoti et al. (2011) with respect to stable Eurasia (Fig. 6). In general, we found model velocities 1.5 2  mm yr −1 lower than observed, Vertical model velocities in the Apennines are compatible with some GPS observations, especially those by Bennett et al. (2012). In general, however, GPS data do not show consistent uplift in the mountain range like the models do. Furthermore, in the foreland the modeled subsidence peak is relatively low in amplitude (less than 5 mm yr −1 ) and extends across the whole plain.

Seismicity and stresses
Crustal stresses in the reference model are extensional underneath the Apennines and compressional throughout the Po Plain (Fig. 4). Average deviatoric shear stresses xy   reach their Modeled earthquakes also occur on thrusts around the base of the range, which produce fewer earthquakes than the external thrust front itself (Fig. 4). This is in agreement with both instrumentally observed and historical seismicity, which tend to cluster on the external thrust front complex. Extensional crustal seismicity in the reference model is focused in the middle of the mountain range and is absent in its most internal part (distance along the profile < 25 x km). This is in agreement with earthquake locations observed tens of kilometers to the southeast and with the proposed source of the w M 6.5 Lunigiana-Garfagnana earthquake of 1920 (Rovida et al., 2016). Nevertheless, several earthquakes have been observed instrumentally along our reference profile at < 25 x km (Chiarabba et al., 2005). The largest events in the model have a W M of 7.4. Reference model seismicity also includes events in the lithospheric mantle, at the top of the Adriatic slab. These occur on spontaneously formed normal faults located in the external hinge area of the bending slab (Fig. 4). Overall, the depth distribution of model earthquake hypocenters is bimodal, with a peak in the upper crust and one in the bend ing slab, and no events in the lower crust (Fig. 7). However, the cluster of instrumentally recorded earthquakes underneath the external thrust front only shows focal locations at middle-to -lower crustal depths (25-45 km), rather than in the lithospheric mantle.

Slab pull controlling the distribution of seismicity
In a novel development, we directly link slab pull and events in the same simulation. To explore the importance of slab pull in driving the dynamics of the northern Apennines, we perform two numerical experiments. In these experiments, we use alternative temperature setups (Section 2.2 and Fig. 3) with different isotherms in and around the slab, thus altering its buoyancy. This also allows us to explore the effect of different mantle temperatures on model dynamics, which is important given the lack of knowledge of the deep thermal conditions. Higher slab temperatures,

J o u r n a l P r e -p r o o f
Journal Pre-proof corresponding to a lower average density and less negative buoyancy and thus to less slab pull, decrease the intensity of crustal deformation (Fig. 8). Conversely, lower temperatures increase slab pull and the intensity of deformation. This occurs despite the absence of strong mechanical coupling between the upper crust and the slab (Fig. 3). The effect of slab temperatures on deformation is manifested in the total seismic moment released and in seismic rates (Fig. 8). In fact, the hotter slab reduces the moment released in both crustal regimes to around 20% of its value in the reference model (17% in the extensional crustal domain, 21% in the compressive one). The seismic rate in the crust is also reduced, though more unevenly in the two regimes. The seismic rate in the mountain range under extensional stresses is reduced by 68% and becomes 60% of that in the active thrust belt, where the reduction is by 50%. The different magnitude of the effect of slab temperatures on seismic moment release and on seismic rates is consistent with the distribution of model stresses and seismicity. In particular, in the reference model the compressional regime has a lower proportion of smaller earthquakes than the extensional one, as well as larger || xx   and II   . This is consistent with the inverse relationship between the Gutenberg-Richter b value and differential stress   observed in the laboratory (Amitrano, 2003), in earthquakes in continental areas (Scholz, 2015), and in STM modeling of orogenic belts (Dal Zilio et al., 2018). Therefore, the seismic rate in the extensional regime is more susceptible to weaker driving forces, as even a small reduction in crustal stresses stops the small extensional earthquakes from occurring. The influence of slab temperatures on crustal deformation is also evident in the surface velocities: the colder, more buoyant slab increases horizontal and vertical velocities, while the hotter slab has the opposite effect (Fig. S1). Based on our results, we can conclude that the lateral distribution of relative earthquake frequency and seismic moment release in the studied migrating orogen is critically controlled by the slab pull driving this system.

Effect of lithospheric mantle stiffness on seismicity
The impact of elastic properties of the lithospheric mantle on shallow tectonics is typically thought to be negligible. We evaluate this impact by varying its elastic stiffness. A higher shear modulus G in the lithospheric mantle leads to a stiffer, more rigid slab, which generally reduces slab bending and lithospheric delamination. As a consequence, less crustal deformation occurs.
Therefore, generally less seismic moment is released in the crust, especially in the compressional domain (Fig. 9). The vertical components of model surface velocities are also affected (Fig. S2). In particular, they decrease by up to 1 mm yr −1 in the external part of the mountain range (thus reducing uplift) and in the plain (which amounts to faster subsidence). Despite a general negative trend, the change in released seismic moment and vertical surface velocities with increasing slab G is variable. This could be due to a trade-off between two processes: on one hand, the lithospheric mantle undergoes less bending in response to the unchanged dynamic load. This forces the deformation due to hot mantle wedge and lower crust protrusion to localize to a greater extent in the crust, rather than in the bending and downwelling of lithospheric mantle. On the other hand, greater stiffness slows down slab bending and thus overall deformation through reducing crustal loading rates. Such influence of slab stiffness on seismicity demonstrates the importance of complex interactions between different subsurface regions. Furthermore, it shows that small changes in stresses and displacements due to preperties of material at sub-crustal depths (greater than 40 km) contribute to determine the lateral distribution and abundance of earthquakes.

Effect of rock strength on earthquake distribution
We also investigate the influence of the yield strength parameters on short-term tectonics in the models. Yield strength is expected to strongly affect local stresses and the specific distribution of

Surface deformation and GPS observations
The horizontal components Upward model velocities are in agreement with exhumation rates estimated near the southwestern end of our reference profile ( 0.7 mm yr −1 on average, 1.8 mm yr −1 maximum; Balestrieri et al. (2003)). However, present-day vertical velocities from GPS observations show a roughly even mix of uplift and subsidence in the mountain range (Fig. 6). This suggest that either localized deformation overprints the slow uplift signal, or that broad tectonic uplift has recently ceased.
The downward velocities produced in the models form a broad peak more than 100 km wide. Conversely, in GPS observations directly along the profile, subsidence in the Po plain is focused in the vicinity of the range and reaches rates of 8 mm yr −1 or more. However, results from the geodetic study of Serpelloni et al. (2013) are in better agreement with the subsidence rates produced in the model. In particular, their best-fit smoothed spline along a line 50 km to the southeast of our profile defines an 80-km-wide peak beginning with an amplitude of less than 4 mm yr −1 . This suggests that subsidence in the eastern part of the Po Plain might be compatible with a lithospheric-scale velocity field dominated by lithospheric delamination and slab flexure and retreat, as proposed by Carminati et al. (2003). The lack of a lithospheric delamination-retreat signal in the vertical velocities observed along our reference profile could result from slab retreat and flexure being significantly slower or currently absent altogether to the northwest of the lithospheric tear proposed by Piccinini et al. (2014). Alternatively, shallow crustal processes such as fault creep and sediment deposition might mask the signal of lithospheric delamination-retreat.

Geodynamics
Our results show that the characteristic tectonic configuration of the Northern Apennines can be explained by lithospheric and asthenospheric dynamics driven only by the buoyancy anomalies of the slab and mantle wedge through partial decoupling and retreat of the downgoing lithosphere, in agreement with inferences of previous studies (e.g., Ventura et al., 2007;Picotti and Pazzaglia, 2008). In particular, slab buoyancy and viscous flow in the hot lower crust and mantle wedge drive lithospheric delamination and slab retreat and reproduce the tectonic regimes of the orogen (Fig.   4). This geodynamic configuration fits well within the larger framework of the entire Western Mediterranean back-arc system being driven by upper mantle-scale convection cells associated with slab retreat, back-arc extension and dynamic topography (Faccenna et al., 2014a). In this context, our modeling results clarify how viscous flow and resulting dynamic topography relate to Furthermore, the rheology of Adriatic middle-lower lithosphere was inferred to control the tectonic style of the orogen by Chiarabba et al. (2014). They proposed a mechanism of propagation of the tip of the asthenospheric wedge that resembles the protrusion of asthenospheric and lower crustal material in our simulations, although they do not attempt to simulate such mechanism or its enabling conditions and tectonic effects.
Our approach imposes inferred present-day structure on the model setup and therefore cannot properly simulate the long-term geodynamic evolution that led to the current configuration where similar processes and features occur. This particularly holds where previous studies inferred the presence of highly ductile or molten lower crust associated with slab rollback, trench ret reat and back-arc extension, such as in the Aegean (e.g., Jolivet and Brun, 2010;Jolivet et al., 2013;Ersoy and Palmer, 2013;Menant et al., 2016;Kruckenberg et al., 2011).
The ductile lower crust and resulting asthenospheric wedge protrusion in the models could not reproduce the thrust-related seismicity observed beneath the external part of the mountain range, at horizontal distances along the profile around 60 km and at mid-crustal depths (Fig. 4).
This supports the idea that ductile deformation beneath the Tyrrhenian Moho is not dominant in the western sector of the Northern Apennines, while crustal underplating is likely taking place.
Such an interpretation of regional tectonics was proposed by Thomson et al. (2010) from thermochronological observations. It is also supported by Chiarabba et al. (2014), who in fact restrict their hypothesized mantle nose mechanism to the more southeastern part of the orogen.
Nevertheless, the models indicate that ductility of deep material at the Tyrrhenian-Adriatic suture and ongoing lithospheric mantle retreat are needed to obtain realistic stress distribution pattern.
This is the case, at least, given the assumptions of buoyancy-driven dynamics and predominantly axis-parallel deformation that underlie our model setup.

Seismicity
Earthquake magnitudes seem realistic, though the maximum values (up to W M 7.4) are significantly larger than those of any known historical event ( W M 6.5 maximum). In the reference model, 19 events larger than any observed earthquake ( w M in the 6.5-7.5 interval) occur in the 20,000 years simulated. If these simulations indeed reflect the true possible maximum event size, events larger than those known from the historical record of the last five centuries are long overdue and may occur in the near future. Alternatively, the models may overestimate the maximum possible event size in the area or its frequency of occurrence. For instance, a single event in the model might correspond in reality to multiple episodes of aseismic and seismic slip, or the J o u r n a l P r e -p r o o f rheology of crustal layers might be overly prone to yielding.
Extensional seismicity in the models occurs on major faults-including reactivated thrusts-underneath the crest of the range, while compressional crustal seismicity is focused on the external thrust front complex. This is in agreement with the observed hypocenter locations and focal mechanisms in an update version of the catalog of Pondrelli et al. (2006).
The moderately deep seismicity (deeper than 20 km) observed underneath both the base of the range and the external thrust front is not reproduced in the models (Figs. 4 and 7). This suggests that the models might not be fully capable to capture the rheological str ucture of the orogen. In particular, a higher-viscosity and more brittle lower crust in the Adriatic foreland may be present.
However, such brittle lower crustal material, if present, is likely to be localized to a relatively small area, since highly ductile material is needed at mid-crustal depths in the suture region between the two lithospheric domains (Sections 3.1 and 4.1.2).
Seismicity in the models is sensitive to multiple physical parameters, delineating a complex non-linear system. The different models highlight the importance of a sufficiently ductile lower crust and mantle wedge in allowing a realistic velocity and stress distribution (Figs 4 and 5).
They also show how significantly the temperature of the slab and its immediate surroundings affects crustal seismicity, particularly the seismic rate in the extensional regime and the seismic moment release in both the extensional and compressional regimes (Figs 8 and 9). Such influence of slab and asthenospheric wedge temperatures on model seismicity, together with the significant observed seismic activity in both regimes in the orogen, suggest that the deep structure of the Northern Apennines includes a significantly cold and negatively buoyant subducting slab and a distinctly hot mantle wedge. The lack of observed earthquakes in the mantle, produced in the models through slab bending (Figs 4 and 7), suggests that the lithospheric mantle in the models is either not strong or not ductile enough. As an alternative explanation, the lithospheric mantle might be undergoing aseismic creep or slow fault slip without generating earthquakes. However, it is also possible that the regional seismic velocity models used to locate earthquake hypocenters systematically underestimate deep velocities and thus the hypoce ntral depths. If so, at least some of the lower crustal earthquakes recorded there and not reproduced in our simulations could indeed be located immediately below the Adriatic Moho, thus eliminating some of the discrepancies between model results and observations. J o u r n a l P r e -p r o o f Journal Pre-proof

Implications of deeper processes affecting upper crustal seismicity
The significance of results and conclusions from our numerical study likely goes beyond the specific setting of the Northern Apennines. In fact, this study shows the importa nce of long-term, large-scale dynamics on determining short-term tectonic regimes and seismicity. It also shows the influence of deep temperatures, strength and stiffness on short-term tectonics. In particular, these parameters significantly affect the stress regime as well as the speed of tectonic loading in the various domains of a geodynamic system. They thus alter the spatial partitioning of earthquakes between different regions (Sections 3.3 to 3.5 and 4.2 and Figs. 8 and 9). The specific respone of the simulated system to changes in model parameters is complex and certainly depends on the details of the setup and rheological quantities. Nevertheless,there is no reason why a broadly similar setup representing another orogenic system dynamically driven b y its own buoyancy structure, especially if via lithospheric delamination-retreat, would lack any sensitivity of short-term crustal tectonics to deep rheology and the resulting dynamics. Certainly, our models indicate that the possibility exists for interplay between rheological complexity and geometrical structure resulting in the crustal seismicity being highly sensitive to various physical and rheological features of deep material. Furthermore, the models suggest that tectonic loading in systems with complex fault networks can be spatially complex and therefore cannot be modeled using a simple boundary condition with uniform far-field relative displacements. It thus needs to be simulated at least partly self-consistently and with the inclusion of realistic driving forces and tectonic loading to be able to correctly reproduce the distribution, frequency, size and type of seismicity in different areas. Overall, both long-and short-term dynamics, and therefore the lower crustal rheologies, deep thermal structure and material parameters that affect them, are needed for a complete physics-based seismic hazard assessment of a region with complex loading conditions resulting from tectonic forces. This importance of realistic tectonic loading and complex rheological structures on seismicity should be taken into account in physics-based seismic simulations (e.g., Dieterich and Richards-Dinger, 2010).
These results also imply that evidence of certain tectonic regimes, in terms of surface deformation and stress orientation, can be used to constrain the rheological structure and driving forces that cause them. Conversely, given a realistic rheological-structural model of the crust, present-day observations can be extended to also use seismicity to point towards plausible lithospheric and asthenospheric rheologies and long-term flow patterns and important parameters J o u r n a l P r e -p r o o f Journal Pre-proof such as mantle viscosity and material strength.

Limitations
This study shows that STM modeling can partly reproduce the observed seismically active tectonic regimes in a very complex system in which deformation is driven by buoyancy forces. This methodology can be applied to a variety of tectonic regimes in complex and diverse geodynamic settings. However, the current limitations of our numerical methodology need to be considered.
Two major limitations are the purely slip rate-dependent friction and the coarse temporal resolution. Both of these have been recently addressed in a newer version of the numerical code, which is yet to be tested on large-scale models (Herrendörfer et al., 2018). Short-term model characteristics such as the duration of events as well as their specific frequency and size may thus be affected, but the long-term stress regimes and resulting large-scale distribution of events in the extensional and compressional domains are likely to be robust. Another major limitation is the 2D model geometry. This allows a relatively simple model setup procedure and short run times, but implies the assumption of lateral homogeneity of the system, which is unrealistic for a real-world orogen.

Conclusions
We investigated the link between the lithospheric-scale geodynamic deformation and the relatively short-term tectonics and seismicity in the Northern Apennines orogen. We use seismo-thermo-mechanical numerical modeling and compare our results with earthquake catalogs and GPS observations. Results show that a large-scale dynamic regime driven by negative slab buoyancy can broadly reproduce the distinct, coupled extensional and compressional tectonic regimes of the Northern Apennines, in terms of both stress orientations and velocities. The latter, in particular, exhibit the coupling of uplift with extension and subsidence with compression that characterizes the region. Nevertheless, discrepancies rema in between modeled and observed vertical surface velocities. This misfit could indicate spatial and temporal changes in orogen behavior, deviating somewhat from the modeled geodynamic mechanism. However, it could also be due at least partly to shallow crustal processes.
The models can reproduce realistic tectonic regimes thanks to a highly ductile lower crust J o u r n a l P r e -p r o o f rheology (here, wet quartzite) and high temperatures. The resulting highly ductile lower crustal material, together with the hot asthenospheric wedge, protrudes upward and northwards as the negatively buoyant slab sinks and retreats. This flow of material produces uplift and extension in the overlying mountain range and decouples the upper crust from the lithospheric mantle, which delaminates. Ahead of the protruding material, where the upper crust is still coupled to the lithospheric mantle, subsidence and compression are produced. Our models thus indicate that thermally-and rheologically-controlled slab sinking and retreat and lower crustal deformatio n are needed to achieve the observed tectonic configuration of the orogen, as previously suggested by Chiarabba et al. (2014); Benoit et al. (2011).
The two tectonic regimes reproduced in the models generate seismicity, which is affected quite strongly in its cumulative released moment and spatial distribution by unexpected parameters. Slab pull controls critically both the distribution of seismic activity in the two tectonic regimes of the orogen. A colder and thus more negatively buoyant slab increases slab pull and thus the intensity of tectonic loading and the deformation and seismicity of the upper crust. Other material properties, such as slab stiffness and the strength of both bulk rock and faults, also influence the specific features of crustal seismicity. The shear modulus of the slab has a variable but broadly dampening effect on the release of seismic moment in the crust, which reflects the delamination-retreat mechanism that drives tectonics in the model. These results highlight the influence of the physical properties of deep material, at lithospheric mantle depths, on upper crustal seismicity. This influence is comparable with that of rock and fault strength, whose effect on earthquake ruptures is more direct and of widely acknowledged importance. This implies that research aiming to simulate the spatio-temporal distribution of earthquakes in complex fault systems should consider realistic forcing and rheologies to obtain appropriate loading of faults and surrounding rocks.
The success of STM modeling in reproducing some major features of the regional tectonics of the Northern Apennines and in describing the complex influence of key physical parameters on such tectonics paves the way for future applications, which could ultimately contribute to improved regional seismic hazard assessments.