Probabilistic Assessment of Slip Rates and Their Variability Over Time of Offshore Buried Thrusts: A Case Study in the Northern Adriatic Sea

When sedimentation rates overtake tectonic rates, the detection of ongoing tectonic deformation signatures becomes particularly challenging. The Northern Apennines orogen is one such case where a thick Plio-Pleistocene foredeep sedimentary cover blankets the fold-and-thrust belt, straddling from onshore (Po Plain) to offshore (Adriatic Sea), leading to subtle or null topo-bathymetric expression of the buried structures. The seismic activity historically recorded in the region is moderate; nonetheless, seismic sequences nearing magnitude 6 punctuated the last century, and even some small tsunamis were reported in the coastal locations following the occurrence of offshore earthquakes. In this work, we tackled the problem of assessing the potential activity of buried thrusts by analyzing a rich dataset of 2D seismic reflection profiles and wells in a sector of the Northern Apennines chain located in the near-offshore of the Adriatic Sea. This analysis enabled us to reconstruct the 3D geometry of eleven buried thrusts. We then documented the last 4 Myr slip history of four of such thrusts intersected by two high-quality regional cross-sections that were depth converted and restored. Based on eight stratigraphic horizons with well-constrained age determinations (Zanclean to Middle Pleistocene), we determined the slip and slip rates necessary to recover the observed horizon deformation. The slip rates are presented through probability density functions that consider the uncertainties derived from the horizon ages and the restoration process. Our results show that the thrust activation proceeds from the inner to the outer position in the chain. The slip history reveals an exponential reduction over time, implying decelerating slip-rates spanning three orders of magnitudes (from a few millimeters to a few hundredths of millimeters per year) with a major slip-rate change around 1.5 Ma. In agreement with previous works, these findings confirm the slip rate deceleration as a widespread behavior of the Northern Apennines thrust faults.


INTRODUCTION
Buried and blind active faults pose one of the major challenges for understanding seismotectonics in many regions worldwide (Yeats and Lillie, 1991;Stein and Ekström, 1992;Berberian, 1995;Shaw and Shearer, 1999;Walker and Jackson, 2002). When such faults are located offshore, there is an additional challenge because they can be investigated only using geophysical subsurface data (Fisher et al., 2004;Sorlien et al., 2006;Hayes et al., 2010;Wolfe et al., 2019). Blind thrusts can be determined by inelastic deformation ahead of fault tip (Roering et al., 1997) that contributes to dissipate the strain in the overlying cover by backthrusting, coupling, and forethrusting (Dunne and Ferrill, 1988). Analog models also show that the upward propagation of thrust faults can be slowed down or even halted by the presence of pre-existing frictional discontinuity (e.g., decollements) above the fault propagating tip (Bonanno et al., 2017). However, they can also be determined by the relatively young inception of fault activity within the current regional stress state or by the interplay between fault activity and sedimentation rates in the cases in which the second one is significantly higher than the first one (Pueyo Anchuela et al., 2016). Blind thrusts are documented as to be the cause of several significant earthquakes and tsunami worldwide (Hauksson et al., 1995;Lettis et al., 1997;Borrero et al., 2001;Hayes et al., 2010;Burrato et al., 2012;McAuliffe et al., 2015) thus significantly contributing to seismic and tsunami hazard.
One of the key parameters for describing a fault system evolution is the rate of activity (slip-rate) and its variability through time. The slip-rate assessment is very sensitive to the scale and the methodological approach to its calculation (see Morell et al., 2020, for an overview). The different approaches to slip rate calculation and the variability of the results can be treated considering all the uncertainties involved and providing a probabilistic approach to the results (Zechar and Frankel, 2009;Styron, 2019). The determination of slip rates of blind faults is mostly based on the analysis of deformed syn-tectonic markers either on surface exposures (Ward and Valensise, 1996;Vannoli et al., 2004;Gunderson et al., 2013Gunderson et al., , 2018Livio et al., 2014;Tiberti et al., 2014) or in the subsurface (Dolan et al., 2003;Wang et al., 2014;Maesano et al., 2015Maesano et al., , 2020Bergen et al., 2017), and requires a good three-dimensional geometrical definition of the underlying fault (Shaw et al., 2002).
Inland, recent or ongoing tectonic activity and fault reactivation in individual seismic events can be studied by analyzing geologic exposures (Bresciani and Perotti, 2014;Livio et al., 2014;Amorosi et al., 2016Amorosi et al., , 2020Zuffetti and Bersezio, 2020), or by detecting river drainage anomalies (Burrato et al., 2003), or by inspecting GPS and remote-sensing data (Serpelloni et al., 2012;Pezzo et al., 2013Pezzo et al., , 2020Nespoli et al., 2016Nespoli et al., , 2017. Instead, the just mentioned approaches cannot be used in the marine environment, especially where the Pleistocene sedimentation rate exceeds the tectonic signal leading to a completely undisturbed and almost flat bathymetry. In such cases, the contribution of detailed geophysical investigations is always needed. Noteworthily, the tectonic activity in the Northern Adriatic Sea has recently been estimated from the analysis of Continuous GPS data combined with the analysis based on publicly available seismic reflection profiles and analytical models . Our study area (Figure 1) has a long record of both instrumental and historical seismicity; in particular, the offshore area near Rimini was affected by severe seismic sequences in 1875 (Mw 5.7) and 1916 (Mw 5.8) (Guidoboni et al., 2018;Rovida et al., 2020). The largest earthquake recorded in the instrumental era is the Mw 5.9, 1930 Senigallia earthquake, associated with a reverse fault running along the coast and controlling the deformation of Quaternary coastal terraces (Vannoli et al., 2004(Vannoli et al., , 2015. Small tsunami effects were also observed following the 1916 and 1930 earthquakes (Guidoboni et al., 2018(Guidoboni et al., , 2019. However, the tsunami potential associated with active faults in the region remains relatively small (Tiberti et al., 2008), and the tsunami hazard in the Northern Adriatic Sea is generally lower than the rest of the Mediterranean region (Basili et al., 2021). Nonetheless, considering the region's built environment, an advanced analysis of the tectonic activity is important to evaluate better and mitigate the risks derived from earthquake and tsunami hazards (Di Bucci et al., 2017;Antoncecchi et al., 2020).
This background information motivated us to investigate two key aspects of active faults that contribute to formulating better estimates of earthquake-related hazards. One aspect is the fault geometry that controls the ground-motion distribution with significant effects in the fault near field (Passone and Mai, 2017) and the tsunami wave height distribution in the fault near field and the far-field (Tonini et al., 2020). The other aspect is the slip rate that controls the maximum amount of tectonic moment converted into seismic moment release (Morell et al., 2020). To FIGURE 1 | (A) Simplified structural map of the Eastern Po Plain and Northern Adriatic Sea. Thrust fronts analyzed in this study are modified after Ghielmi et al. (2013). The color-coded bathymetry, also represented by dashed gray contour lines (every 20 m), is derived from the EMODnet Bathymetry Consortium (2020). Historical earthquakes (Mw > 4) are from Rovida et al. (2020), with labels indicating the year for the largest ones. The traces of the Northern Section and the Southern Section refer to two 2D seismic reflection profiles used in this study to compute the slip rates (see (C) Regional cross-section showing the structural arrangement of the main thrusts and detachment from Fantoni and Franciosi (2010). See the section trace in panel (A) for the location.
this end, we reconstructed the 3D geometries of the blind thrusts and a regional unconformity at the base of Pliocene deposits taking advantage of a wide dataset of 2D seismic reflection profiles (provided by ENI; Figure 2). Then, focusing on two high-quality regional seismic reflection profiles (Ghielmi et al., 2013), we performed a probabilistic assessment of the slip rates of four major thrusts by restoring the deformation of several key-horizons after removing the differential compaction effect. FIGURE 2 | Dataset of 2D seismic reflection profiles and wells (courtesy of ENI) used in this study to reconstruct the fault 3D geometries and the hor-PL1 regional unconformity ( Figure 5), the cross-sections (Figures 6-8), and to derive the slip rate values ( Table 3). The bathymetry, as in Figure 1, is derived from the EMODnet Bathymetry Consortium (2020). We adopted the trishear inversion method ; Figure 3) for recovering the fault-propagation folding and incorporated the uncertainty on both displacement and age in the slip-rate calculations using the formulation of Zechar and Frankel (2009). We calculated multiple slip rates and associated uncertainties on eight stratigraphic horizons, spanning the last 5.33 Ma, from the Pliocene to the present, supporting recent findings that suggest an ongoing activity of the Northern Apennines basal thrust (Barba et al., 2008;Carafa et al., 2015a;Pezzo et al., 2020) and detailing the previous work on this topic performed in adjacent areas but with coarser temporal resolution (Maesano et al., 2013(Maesano et al., , 2015. The reconstructed fault 3D geometries together with the computed slip-rates can be of help in better defining the tectonic evolution of coastal and nearoffshore active faults in the Northern Adriatic Sea and provide data that can be integrated into databases of seismogenic sources  and seismic and tsunami hazard studies.

THE PO PLAIN-NORTHERN ADRIATIC SEA FOREDEEP BASIN
The Northern Apennines is part of the system of fold-and-thrust belts (Northern Apennines, Southern Alps, and Dinarides) that wraps around the northern part of the Adria microplate, which is a promontory of the African plate (Malinverno and Ryan, 1986;Dewey et al., 1989), in the general framework of the Africa-Eurasia convergence. Starting from the late Mesozoic time, Adria was involved in multiple tectonic phases that led through times to the development of the Dinarides to the east (Handy et al., 2010), the Southern Alps to the north (Carminati and Siletto, 1997), and the Northern Apennines to the Southeast (Patacca et al., 1990;Carminati et al., 2012). Adria is presently overridden to the northeast by the outermost fronts of the Dinarides and the Southern Alps in the Veneto-Friulian Basin (Fantoni et al., 2002;Caputo et al., 2010;Toscani et al., 2016), to the west by the Southern Alps and Northern Apennines in the Po Plain (Fantoni and Franciosi, 2010;Toscani et al., 2014), and to the south by the Northern Apennines and Dinarides confronting to each other in the central Adriatic Sea (Argnani, 1998;Di Bucci and Mazzoli, 2002;Scrocca, 2006;Scrocca et al., 2007;Scisciani and Calamita, 2009;Kastelic et al., 2013;Finocchio et al., 2016).
The Po Plain and Northern Adriatic Foreland basin subsurface architecture resulted from a Mesozoic extensional tectonic phase developed in the western Tethys realm related to Upper Triassic-Lower Jurassic rifting, with evidence for Cretaceous-to-Paleogene structural inversion during Alpine pre-and syn-collisional tectonic phases (Bertotti et al., 1993;Masetti et al., 2012 and FIGURE 3 | Workflow adopted in this study to decompact and restore the regional cross-sections (Figures 7, 8) and derive the slip rate values (Table 3) using a combination of the fault-parallel flow algorithm for displaced horizons and the trishear algorithm for folded horizons. The trishear parameters have been defined through statistical analysis shown in Figure 9 and Supplementary Data 1. reference therein). Since the late Miocene, contraction affected the foreland and the surrounding orogenic belts due to Apennine north-to-northeast migration (Barchi et al., 1998;Basili and Barba, 2007;Franciosi, 2008, 2010;Ghielmi et al., 2013;Turrini et al., 2014). The buried thrusts and related fault-propagation folds have been investigated mainly using seismic reflection profiles and deep well logs taking advantage of hydrocarbon exploration activities in the region (Pieri and Groppi, 1981;Barchi et al., 1998;Regione Emilia-Romagna and ENI-AGIP, 1998;Regione Lombardia and ENI-AGIP, 2002;Fantoni et al., 2004;Fantoni and Franciosi, 2010;Ghielmi et al., 2010Ghielmi et al., , 2013Toscani et al., 2014;Turrini et al., 2014;Maesano et al., 2015;Amadori et al., 2019). The outermost fronts of the Northern Apennines buried in the southern Po Plain and the central Adriatic Sea are organized in structural arcs separated by recesses (Figure 1; Livani et al., 2018;Amadori et al., 2019). These arcs show an eastward increase in shortening and tectonic activity, consistent with their emplacement processes, while the Northern Apennines underwent a counterclockwise rotation within oblique convergence (Carminati et al., 2012;Maino et al., 2013). More specifically, the easternmost structural arcs were activated during the Pliocene, and their tectonic activity continued throughout the late Pliocene until the present, as demonstrated by the age of syn-tectonic deposits and the associated growth strata geometries (Boccaletti et al., 2011;Bonini et al., 2014, Maesano et al., 2015. The blind thrusts in our study area (Figure 1), which are part of the Romagna Arc and the eastern termination of the Ferrara Arc, are rooted on two regional decollement levels, at the bottom and on the top of the Mesozoic carbonates, respectively (Barchi et al., 1998;Fantoni and Franciosi, 2010), and deform Early Pliocene and Pleistocene deposits (Maesano et al., 2015). All thrust fronts and related anticlines show an arcuate shape with a general NW-SE trend and NE vergence (Bigi et al., 1992;Casero, 2004;Fantoni and Franciosi, 2010).
The Adriatic basin infill is organized in large-scale tectonoeustatic driven sequences defined as unconformity-bounded units whose boundaries correspond to abrupt major changes in the type and gross distribution of depositional systems, angular unconformities due to high-magnitude basin-wide tectonic phases or eustatic changes that occurred during the Plio-Pleistocene. From the bottom up (i.e., from the Zanclean base to the present-day seafloor), these sequences are referred to as PL1, PL2, and PL3 for the Pliocene, and PS1, PS2, PS3a, PS3b, and PS3c for the Pleistocene (Figure 4). To discriminate between seismic horizons (i.e., the sequence boundary unconformity and correlative conformity) and the stratigraphic unit or sequence they bound at the base, we use the prefixes "hor-" and "seq-" for horizon and sequence, respectively, before the age name abbreviation (Figure 4). Note that the unconformities always lie at the base of a sedimentary unit. The accurate determination of chronostratigraphic ages relies on micropaleontologic and magnetostratigraphic analyses ( Table 1; Muttoni et al., 2003;Ghielmi et al., 2010Ghielmi et al., , 2013, which were used to support and verify the correlations between well-logs and seismic horizons, as well as to aid the interpretation of the sedimentation environment. From a sedimentological perspective, the Plio-Pleistocene sedimentary infill of the Po Plain and Northern Adriatic Foreland represents an overall regressive cycle (sensu Catuneanu, 2004).
The Pliocene sequence is mostly represented by thick turbidites deposited in a deep marine environment (during the underfilled phase), upwardly and eastwardly evolving into shallow water (slope, shelfal, and coastal) and continental (fluvial) sequences during the Pleistocene (filled and overfilled phases).
After the abrupt transgression that ended the Messinian Salinity Crisis (Amadori et al., 2018;Micallef et al., 2018), Pliocene marine clays were deposited on the Adriatic foreland and foreland-ramp at the base of the turbiditic succession (Ghielmi et al., 2013;Rossi et al., 2015). These turbidites onlap the substantially undeformed SW-dipping Adriatic foreland ramp. In the study area, Plioceneto-Calabrian clastic deposits (seq-PL1 to seq-PS2) are present with thin sequences of marls intercalated with thin-to thickbedded fine-grained sandstones, interpreted as "highly-efficient" turbidite systems (according to Mutti et al., 1999;Ghielmi et al., 2010Ghielmi et al., , 2013Ghielmi et al., , 2019. The thickness of these sequences is highly variable due to the different depositional settings (from the Northern Adriatic foredeep depocenter to the shallower foreland ramp toward the northeast) and the coeval tectonic activity. Based on the works by Ghielmi et al. (2013) and Ghielmi et al. (2019), the maximum thickness for the early Pliocene sequence of the Canopo Formation (seq-PL1) is about 600 m, whereas for the first Pleistocene sequence of the Porto Corsini and Porto Garibaldi Formations (from seq-PL2 to seq-PS1) the maximum thickness is about 1,200 m each. These systems show remarkable lateral continuity in seismic reflection profiles. They are overlain by a maximum of 800 m of Middle-Late Pleistocene sand-rich turbidites of the Carola Formation (seq-PS2) and 1,500 m of the north-eastward progradational geometries of the Pleistocene Progradation Complex and Cyclothemes of the Ravenna Formation (seq-PS3).
Despite the tectonic uplift and exhumation of the Northern Apennines range to the southwest, the main area of sediment provenance for the Pliocene and Pleistocene sequences was located to the northwest, as shown by southeast-directed paleocurrents (Ghielmi et al., 2010(Ghielmi et al., , 2013. Therefore, the sediment supply's main direction during the Pliocene and Pleistocene was almost perpendicular to the fold-and-thrust vergence and propagation.

DATASETS AND METHODS
Here we briefly illustrate the subsurface dataset (Figure 2) used to reconstruct the fault planes and the Pliocene unconformity (hor-PL1) in the study area. Since the area is almost completely offshore, the ground-surface reference dataset includes the bathymetric data from the EMODnet Bathymetry Consortium (2020). We also present the methodological approach we followed to (1) depth convert and restore two regional seismic reflection profiles representative of the structural setting and (2) calculate the Probability Density Function (PDF) of the slip rate for each fault, for each deformed horizon within these sections.

Seismic Dataset
The subsurface datasets comprise public and confidential 2D seismic reflection profiles, all acquired in TWT units. Specifically, we collected (1) all the seismic reflection profiles available in FIGURE 4 | Simplified Plio-Pleistocene tectonostratigraphic scheme of the Northern Adriatic Sea (after Fantoni and Franciosi, 2010;Ghielmi et al., 2010Ghielmi et al., , 2013). The regional unconformities/stratigraphic horizons used in this work and related ages (see Table 3) are indicated. A brief description of the sedimentological characteristics and seismic facies is provided for each geologic formation (for details, see Ghielmi et al., 2013;Amadori et al., 2020). Vertical red arrows on the left-hand panel indicate the inception and duration of tectonic activity of the thrust faults detected on the Northern Section and Southern Section, as derived from decompaction and restoration (see Figure 1 for trace locations and the VIDEPI project, 1 (2) published sections and subsurface maps (Fantoni and Franciosi, 2010;Amadori et al., 2018Amadori et al., , 2019, and (3) confidential 2D seismic reflection profiles kindly provided by ENI in the framework of the SPOT project . This data collection resulted in a dense grid of seismic reflection profiles covering the entire study area (Figure 2).
Public data available in the VIDEPI project consist of 112 deep wells with associated stratigraphies and about 100 seismic reflection profiles, for a total length of almost 2,000 km. These well and seismic data were drilled and acquired at different times and with different techniques, resulting in an overall inhomogeneous dataset. The dataset quality and resolution permit to depict the general structural setting, but it is not suitable for detailed structural or stratigraphic reconstructions. 1 www.videpi.com Nevertheless, some of these data have been considered together with the confidential high-resolution dataset provided by ENI. This dataset includes more than 900 seismic reflection profiles for a total length of over 10,000 km. All data have the same datum plane and maintain a good quality up to 4,000-5,000 ms. Over an extended portion of the study area, the Pliocene succession and the youngest Messinian deposits are usually separated by a marked contrast of acoustic impedance, creating a strong and well recognizable high-amplitude reflector, also identified in neighboring onshore and offshore sectors of the Po Plain-Northern Adriatic basin (Ghielmi et al., 2013;Amadori et al., 2018Amadori et al., , 2020Mancinelli and Scisciani, 2020).
To detect possible evidence of recent/ongoing tectonic activity, one needs to analyze seismic reflection profiles with a good resolution in the first 4-500 ms. We thus selected two of the best seismic reflection profiles running SW-NE 1 | Ages and related uncertainties of the stratigraphic horizons that were used for restoring the cross-sections of Figure 6 and for calculating the slip rates ( through the study area to perform the structural restoration and slip rate analysis. These profiles originally belong to the ENI "Adria" 3D seismic survey and are interpreted and described in detail by Ghielmi et al. (2013).

Depth Conversion
The reconstruction of fault geometries, the structural restoration, and the following slip and slip-rate calculation necessitate depthconverted cross-sections. To this end, we took advantage of the well-log stratigraphies and TWT-velocity tables associated with the wells available in the study area, which survey the entire Plio-Pleistocene sequences (Figure 2). We derived the trendline and the related function for each time-depth table, selecting the one with the highest coefficient of determination (R 2 ), using an approach similar to that adopted for more complex 3D depth conversions (Maesano and D'Ambrogi, 2017). The obtained velocity function is linear [see also section "Fault surface and PL1 unconformity (hor-PL1)"], in agreement with velocity models recently carried out for neighboring sectors in the Adriatic Sea (Mancinelli and Scisciani, 2020), and fully comparable with other velocity models used for the same stratigraphic sequence in the Po Plain-Northern Adriatic Sea (Amadori et al., 2019).

Decompaction and Restoration
To carry out the complete restoration of the interpreted and depth-converted cross-sections, we used the MOVE Software (Petroleum Experts) (Figure 3). Different algorithms were adopted to appropriately restore geological structures, depending on the observed deformation style and related kinematics. For a correct restoration of the regional cross-sections, all deformation due to differential sediment compaction effects must be considered and excluded from the tectonic-related component before the slip calculation. Restoration can then effectively be carried out step by step after having decompacted the sedimentary body above the target stratigraphic horizon. In this view, decompaction is needed to remove the effect of sediment volume change over time-related with porosity reduction, and it represents a crucial operation to avoid an overestimation of fault slip and, consequently, the slip-rate. Sediment compaction can be driven both by natural (Teatini et al., 2011) and human activities (Stramondo et al., 2007). By removing the load associated with a layer of sediments, the residual deformation of the underlying levels can safely be ascribed to tectonic processes. The applied decompaction algorithm is based on principles described by Sclater and Christie (1980). The original sediment porosity and elastic properties were assumed homogeneous, and an equal mixture of shale and sand was also assumed (Scrocca et al., 2007;Maesano et al., 2013Maesano et al., , 2015. In the Northern Adriatic foreland basin sedimentary succession, the sediment thickness can change from 10 to 40%, with an average change of 25% (Maesano et al., 2015). Despite sediments can reduce their thickness up to 55%, the associate fault displacement loss varies between 5 and 15% (Taylor et al., 2008). After each decompaction step, we progressively recovered each key stratigraphic horizon residual deformation, from the youngest to the oldest, by calculating the fault slip amount needed to restore the selected stratigraphic horizons to the original horizontal depositional geometry. Different deformation styles associated with faults and folds require different restoration procedures. The study area is characterized by fault-propagation folding controlled by buried thrusts. The trishear algorithm (Erslev, 1991;Allmendinger, 1998) is suitable for studying complex fault propagation folds geometries and allows for reproducing some features commonly observed in natural folds, such as footwall synclines and thickness and dip variations in the forelimbs (Allmendinger et al., 2011). The trishear algorithm was thus adopted to restore the folding ahead of the fault tip. Since the trishear transfers the fault slip in a triangular shear zone emanating from the fault upper tip, the resulting fault and fold geometries depend on the fault tip location, the fault ramp dip, the amount of slip, the ratio between the propagation of the fault and the slip (p/s), and the angle of the trishear zone.
Trishear can be used for both forward and inverse modeling of fault displacement. In the case of inverse modeling, due to the number of variables in the trishear formulation, there is no univocal solution in the restoration Oakley and Fisher, 2015). We thus adopted the formulation of trishear proposed by Cardozo et al. (2011), which solves the inverse problem by minimizing an objective function, representing the least-squared linear regression of the restored horizons, by using an optimization method (simulated annealing). The optimization algorithm does not systematically test the parameter space (which is an inefficient approach from the computational side) but traverses the space searching for the best-fit model. Even if the optimization procedure proposed by Cardozo et al. (2011) is less affected by local minima compared to earlier versions (Cardozo and Aanonsen, 2009), we performed multiple runs for each model, particularly for the cases in which the amplitude of folding was small to have a more robust assessment of the stability of the inversion.
For each structure and each horizon, we performed between 10,000 and 50,000 inversions in each model run. In setting up the model, we decided to restrict the parameter space for the variables that are well constrained by the seismic data interpretation. We thus limited the boundary area for the tip position by the upper and lower coordinates of the fault ramp in the cross-section. According to seismic interpretation, the fault ramp dip was searched in a narrow range of variability (±5 • ) around the measured values. We left search limits for the other parameters (slip, trishear angle, and p/s) wider (slip between 0 and 2 km; trishear angle between 20 and 80 • ; p/s between 0 and 3). Since the parameter exploration space is a rectangular area between the upper and lower endpoints of the fault, we have several final solutions located at various distances from the fault trace. We assigned a weight to each solution obtained from the inversion to incorporate the epistemic uncertainty related to the restoration performed with the trishear method. The highest weight was assigned to the solutions with smaller residuals from the inversion and the fault tip position near the fault traces. The weight decreases linearly with increasing distance (inverse distance weighting method) and for the solutions with high residual. Note that the value of slip is independent of the criteria adopted for the weighting of the results. The slip of the best-fit solution was used to restore each horizon before the successive decompaction step.
The Fault Parallel Flow is an algorithm designed for kinematic modeling of the hanging wall blocks where deformation is accommodated by fault-parallel shear (Egan et al., 1997). This algorithm moves all the material within the hanging wall parallel to the fault trace along virtual flow paths and can be applied to a wide range of fault geometries, and it is better suited for compressional faults (Ziesch et al., 2014). Therefore, we restored the observed on-fault deformation with the Fault Parallel Flow method, where the displacement could be directly estimated from the cut-off of faulted horizons. However, the Fault Parallel Flow restoration does not completely recover the folding associated with off-fault deformation. This residual deformation was thus recovered using the Trishear, and the final slip amount associated with the fault is given by summing the results of both methods.

Probabilistic Slip Rate Assessment
Slip rate is a derived quantity that bears the epistemic uncertainty of both variables needed to calculate them: the amount of slip and the age of the reference horizon. To this end, we adopted the probabilistic approach for the slip rate calculation and reporting proposed by Zechar and Frankel (2009), which has been already applied in different tectonic contexts, continental and offshore faults, as well as at different spatial and temporal scales (Frankel et al., 2011;Gold and Cowgill, 2011;Chevalier et al., 2012;Maesano et al., 2020).
We first considered the age and uncertainty of eight Plio-Pleistocene horizons ( Table 1) and used a boxcar PDF with not-null probability in the time interval covered by the horizon age uncertainty.
If t 1 and t 2 are the minimum and maximum ages, respectively, for a given sequence, the probability for the slip to occur at any time during this interval, p(t), is expressed by Then, considering the slip obtained from the restoration of each horizon for each fault, we used a Gaussian PDF around the weighted mean displacement values.
If µ d is the mean displacement and σ d is the associated standard deviation, the Gaussian displacement uncertainty, p(d) is expressed by: Finally, we combined these PDFs to report the probability of the slip rate estimates. The maximum possible slip rate is given by the maximum displacement divided by the minimum age and vice versa for the minimum slip rate.

Stratigraphic Ages and Their Uncertainty
High-resolution geological information on the youngest sedimentary units is needed to assess details of recent tectonics. The seismo-stratigraphic, sedimentologic, and bio-chronostratigraphic characterization of the Plio-Pleistocene sequence reached in the Po Plain and Northern Adriatic Foreland area (Muttoni et al., 2003;Ghielmi et al., 2010Ghielmi et al., , 2013Amadori et al., 2019Amadori et al., , 2020 provides an effective tool to constrain the inception age of thrust faults and their activity age ranges, as well as to estimate displacement and slip rates for the Middle-Late Pleistocene time interval (Figure 4 and Table 1). Given the available temporal resolution, the probability of slip rate estimates is computed considering errors on the identified key-horizon age intervals. Here, we introduce the uncertainty associated with each of them based on bio-chronostratigraphic zones and subzones, paleomagnetism, or the tuning to Milankovitch astronomical cycles [e.g., astronomically tuned Neogene time scale (ATNTS2004) Lourens et al., 2004]. From the older to the younger, the first horizon in the Pliocene part of the sequence is hor-PL1, which corresponds to the beginning of the Zanclean Stage, calibrated at the base of the planktonic foraminiferal biozone MPl1, dated 5.33 Ma (Lirer et al., 2019). This reflector can be easily tracked in the seismic profiles because it was generated during an overall general reorganization of the Po Plain and Northern Adriatic Foreland basin at the end of the Messinian Salinity Crisis. During the Zanclean time, the study area was the foredeep deep-water depocenter of the Northern Apennines; therefore, we can exclude hiatuses within seq-PL1. The uncertainty associated with hor-PL1 age determination is ∼0.1 Myr, tied to the base of the Sphaeroidinellopsis s.l. biozone, and it cannot be older than the MPl1 biozone Lirer et al., 2019).
The third horizon is the Piacenzian hor-PL3, dated between 3.3 and 3.4 Ma and falling within the MPl4b biozone with a considered uncertainty of ±0.1 Myr (Ghielmi et al., 2010(Ghielmi et al., , 2013.
The first horizon in the Pleistocene part of the sequence, hor-PS1, is close to the Gelasian base, constrained in the lower part of MNN18 nannofossil biozones dated 2.4 Ma (Ghielmi et al., 2010(Ghielmi et al., , 2013. In the absence of specific studies on this age determination uncertainty but considering hor-PS1 to have similar stratigraphic detectability and importance as hor-PL3, we assigned an error of ±0.1 Myr even though the younger age would let suppose a smaller uncertainty. The second is the hor-PS2 unconformity, which is dated at 1.5-1.6 Ma because it falls within the lower MNN19c nannofossil biozones before the first occurrence of the Hyalinea balthica, dated at ∼1.5 Ma (Sprovieri et al., 1998). According to the age determination accuracy, we retain an error of ±0.1 Myr also in this case.
The Middle-Late Pleistocene is characterized by the regionwide Pleistocene Progradation Complex, within which we identified three main unconformities, namely hor-PS3a, hor-PS3b, and hor-PS3c. At the regional scale, the hor-PS3a unconformity is highly diachronic and falls into the MNN19f biozone, so that its age ranges in a large interval between 0.98 and 0.47 Ma (Lirer et al., 2019). However, in the Po Plain region, hor-PS3a is well tied at 0.87 Ma, which corresponds to the glacio-eustatic minimum peak MIS 22, through calibration using pollen and paleomagnetic inclinations (Muttoni et al., 2003). Based on the unconformity mapped by Amadori et al. (2019), we can safely attribute the age of 0.87 Ma to the hor-PS3a identified in our study area. The associated error can be estimated to be ±0.045 Myr based on the temporal range in which the MIS 22 peak can be defined according to the position of the two next midpoints of the sea-level curve.
The hor-PS3b unconformity is dated after the last occurrence of L. Gephyrocapsa sp. 3 at 0.597 Ma (Lourens et al., 1998). Based on the uncertainty usually attributed to the glacial-interglacial cycle absolute timing (less than ∼0.5 Myr; Lourens, 2004), and to keep consistency with the older horizon in the sequence, we assigned to hor-PS3b an error of ±0.045 Ma.
The uppermost unconformity, hor-PS3c, is the top lap of the Pleistocene Progradation Complex in the study area and does not correspond to any previously defined stratigraphic event. Hence, the hor-PS3c age determination had to be indirectly estimated by considering the sediment thickness deposited above PS3b and the average depositional rate, known to be slightly lower than 1.0 mm/year (Amadori et al., 2020). Given the average thickness of seq-PS3c and its average depositional rate, we estimate the age of hor-PS3c at 0.4 Ma. To validate this estimate, we compared the sedimentary setting and seismic facies association of the study area with the high-resolution chronostratigraphic correlation available for the neighboring Venetian-Northern Adriatic region (Massari et al., 2004;Amadori et al., 2020). In both cases, the first 500 ms of the seismic sequence corresponds to the topset portion of the Pleistocene Progradation Complex. According to the depth-converted position of hor-PS3c in the study area, the estimated age-depth relationship, and the correlation with dated horizons in the Po Delta area (Amadori et al., 2020), we assume an age of 0.4 and an error of ±0.05 Ma based on the age and scale of individual Milankovitch cycles.

Fault Surface and PL1 Unconformity (hor-PL1)
The interpretation of the seismic reflection profiles (Figure 2) allowed for a detailed 3D reconstruction of the hor-PL1 unconformity and 11 fault surfaces, referred to as NAF, for Northern Apennines Fault, followed by a progressive number (Figure 5). To produce this reconstruction in the depth domain, both 3D fault surfaces and the two interpreted seismic reflection profiles (Figure 6) have been depth converted following the same approach recently proposed for neighboring sectors of the Adriatic Sea by Mancinelli and Scisciani, 2020. Considering the Plio-Pleistocene sequence time-depth values, the function better representing the velocity increase with increasing depth is linear and is v = 0.3557d + 1945, where v is velocity and d is depth, with an R 2 = 0.881. Since all the faults mapped in the area are detached below hor-PL1, the velocity model validity domain includes the depth-converted fault geometries and the hor-PL1 surface entirely.
The hor-PL1 unconformity was reconstructed at a regional scale to provide a first-order geometrical constraint on the folding related to the underlying thrust activity and is particularly detailed in northeastern and foreland sectors. The hor-PL1 reconstruction was integrated with the available information from the literature (Bigi et al., 1992;Amadori et al., 2018) where the seismic dataset was not dense enough to provide an accurate interpolation (Figure 2). The depth contours depict a rather monotonous southwestward gently dipping monocline in the foreland area interrupted by a deformation pattern dominated by arc-shaped NW-SE oriented anticlines.
The fault geometries were initially obtained from interpreting the 2D seismic reflection profiles and interpolating the points sampled along the traces. Then, the final fault surfaces were obtained by an overall depth conversion. Coming from an industrial dataset, the seismic reflection profiles are not homogeneously distributed over the study area, but they are more concentrated in the correspondence of selected geological structures. For this reason, the fault reconstructions are not always laterally continuous and needed to be interrupted where the distribution of the seismic profiles did not permit their recognition. However, their continuity could be inferred by the anticline geometry detected in the hor-PL1, which aided the interpolation in zones of scarce data coverage or seismic profile poor quality.
All the faults strike NW-SE and dip SW coherently with the Northern Apennines thrust belt and are buried under Plio-Pleistocene deposits. For this reason, none of them has any clear seafloor evidence of activity that could be deduced from an inspection of the bathymetry. Instead, their overall structural arrangement could be characterized by analyzing the regional cross-sections (Figure 6), which show that they offset the Pliocene horizons and deform the Pleistocene horizons. They  Table 2.
also show that minor backthrusts contribute to accommodate the deformation produced by the identified major thrusts.
A summary of the main geometric characteristics of the mapped faults is presented in Table 2. Dealing with arc-shaped concave-up faults, the internal variability in strike, dip, and depth is represented by percentiles (10th, 50th, and 90th) of the encountered values. Overall, these faults range in size from 12 to over 400 km 2 and cut through the crust between 2 and 6 km depth. Note that the lower depth values do not necessarily correspond to the detachment level depth but only to the depth where a reliable reconstruction was possible.

Regional Cross-Sections and Chronology of Thrust Activity
The seismic reflection profiles that are shown in Figure 6 (location in Figure 5) have a high-enough resolution for observing growth strata on top of the main thrusts and recognizing the horizon geometries within the Pleistocene Progradation Complex (hor-PS3a, hor-PS3b, and hor-PS3c; Figure 4 and Table 1) in the first 500 ms. The Northern Section ( Figure 6A) is 27 km long, whereas the Southern Section ( Figure 6B) is an excerpt of a more than 70 km long seismic reflection profile and was cut where it intersects our study area (at 36 km in Figure 25 in Ghielmi et al., 2013). Overall, they intersect and provide insights for four of the 11 thrust ramps in the study area and their related secondary features. For some of them, the fault geometry can be followed down to their basal detachment (∼4,500 ms) and the bottom Pliocene unconformity. Detailed line-drawings of the tectonic structures and depositional sequence boundaries have been possible in the entire depth range spanned the profiles, allowing for estimating their inception age.
According to the structural map shown in Figure 5, the structures intersected by the two cross-sections are as follows. NAF1 is the innermost structure and is intersected by both cross-sections close to its structural culmination. NAF2 is intersected only by the Southern Section near its northern termination. NAF3 is composed of a series of NW-SE thrusts correlated from the onshore Ferrara Arc to the offshore Romagna Arc. The lack of seismic profiles between the two cross-sections prevents imaging the along-strike features; therefore, we keep the NAF3North and NAF3South geometric reconstructions separate. However, based on the structural trend deduced from other profiles, we can consider NAF3North and NAF3South as part of the same geological structure. NAF4, the outermost structure in this sector, is intersected only by the Northern Section, although it terminates just north of the Southern Section.
By analyzing these seismic reflection profiles, it was possible to infer the inception age of each structure by comparing the stratigraphic succession thickness in the hanging wall and footwall on each of them and verifying the presence of onlap geometries or growth strata in correspondence of fault-related anticlines. The presence of possible recent and ongoing tectonic deformation, albeit less evident than that affecting the Zanclean to Gelasian sequences, can be noted in the close-up views of the upper sectors of the seismic stratigraphy (Figures 6C-E).
NAF1 is made up of the main forethrust, with a deep splay visible only in the Northern Section and an associated backthrust. This structure clearly offsets hor-PL1 and hor-PL2. Above NAF1, seq-PL1 has the same thickness in both the hanging wall and footwall, whereas seq-PL2 is completely absent and seq-PL3 onlaps on the anticline forelimb. The main thrust activity can then be set at the beginning of seq-PL2. Later activity is testified by the warping of hor-PL3 and hor-PS1 on top of the anticline. The stratigraphic sequence above hor-PS1 shows only hints to elusive deformation, difficult to be traced. However, the stratigraphic markers are well detectable, and the weak signs of deformation appear more clearly through vertical exaggeration, as shown in the close-up view of Figure 6C. NAF2 includes a forethrust splayed into two branches and two separate backthrusts. The seq-PL2 thickness differences between the hanging wall and footwall suggest that NAF2 inception age can be set at the sequence base. This structure slightly deforms the bottom of the seq-PS1, whereas hor-PS2 is flat and undeformed, thereby marking the ceased activity.
NAF3North is constituted by a main forethrust and a few associated secondary backthrusts with very limited displacement. The main thrust shows a flat-ramp geometry and is rooted within the Plio-Pleistocene deposits. Above NAF3-North, seq-PL1 has the same thickness in the hanging wall and footwall. The main thrust inception age can be identified within the deposition of seq-PL3 that is strongly thinned on top of the related thrust anticline. Hor-PS2 is warped but not offset by the thrust, and signs of thrust activity are also detectable in the lower levels/reflectors of seq-PS2.
NAF3South is a rather complex structure in which the main forethrust with a flat-ramp geometry is associated with splay faults and backthrusts. This structure offsets and deforms both Pliocene and Pleistocene successions. Its tectonic activity starts at the base of seq-PL3, as evidenced by the sequence reduced thickness on the top of the related thrust anticline. The main fault slightly offset hor-PS2 at the base of the sequence, and growth strata near the main anticline demonstrate its activity carried on during the whole seq-PS2 deposition. Above the anticline, looking at the seismic reflection profiles, no clear evidence of tectonic activity is detectable for seq-PS3a and seq-PS3b due to the presence of clinoforms and prograding strata. However, the effects of tectonic activity can be detected for the overlying seq-PS3c (Figure 6D), suggesting that the tectonic process remained active during the entire sequence emplacement. NAF4 is a flat-ramp-shaped thrust with an associated backthrust. The ramp connects downward with a basal flat that merges with NAF3North. Seq-PL1 and seq-PL3 have the same thickness at the hanging wall and footwall of the main thrust, and reflectors are always parallel and coherent with the stratigraphic boundaries so that the hor-PS1 unconformity dates the thrust inception age. The thickness of seq-PS1 shows strong variations on the top of the anticline and significant differences relative to the undeformed layers in front of it. Hor-PS2, instead, is only slightly offset by the main thrust. On top of the anticline, seq-PS2 growth strata are present on both limbs, and the hor-PS3a and the reflectors immediately above it shows weak signs of deformation ( Figure 6E).

Restoration and Slip Rates
The two seismic reflection profiles shown in Figure 6 were depth converted to obtain two regional cross-sections suitable for restoration and concurrent calculation of the total slip and subsequent estimation of the slip rate on each fault. These two cross-sections intersect four structures: NAF1, NAF2, NAF3, NAF4. The depth-converted versions of the two cross-sections are shown in Figures 7, 8. The line drawing in the depthconverted sections is rather simplified relative to the original interpretation of the seismic profiles in Figure 6 to avoid numerical instabilities in the restoration process.
We performed the horizon restoration and fault-slip calculation from the younger to older horizons and from the outer to inner structures. The hor-PS3c is the shallowest horizon, placed around 300 m depth. Hor-PL1 reaches almost 6 km depth. We excluded from the restoration procedure the horizons located within the Pleistocene Progradation Complex since their clinoformal deposition prevents the correct discrimination of the tectonic component in determining their present geometry. The shallowest horizons were affected only by folding and were thus restored using the trishear method only. Deeper horizons were restored, combining the Fault Parallel Flow and the trishear methods. The step-by-step restoration progression is shown in Supplementary Figure 2.
We then selected the slip value that best represents the restored deformation by weighing each inversion component (p/s ratio, trishear angle, and ramp-dip angle). Figure 9 shows an example of the trishear inversion results, while the full model inversion results are provided in Supplementary Figures 3, 4 (numerical information in Supplementary Data 1). For the shallower horizons (hor-PS3a, hor-PS3b, and hor-PS3c), we made 50,000 trials to test better the stability of results associated with the very low slip values. This effect is represented by the high dispersion of the fault tips used for the trishear inversion. It is worth noting that for some of the shallow horizons (NAF1-South hor-PS3a, NAF1-South hor-PS3b, NAF3 hor-PS3a, NAF3 hor-PS3c, NAF4 hor-PS3a, NAF4 hor-PS3c), the resulting slip values are not null, and the position of the solution is in good agreement with the fault geometry even considering a large number of trials with a significant dispersion around the faults. We considered the weighted average and standard deviation of the slip values to incorporate these results in the slip rate calculations. The trishear restoration produced a larger uncertainty than the Fault Parallel Flow restoration because, in the latter, the uncertainty is limited to the accuracy of the offset measured at the fault cut-off. Overall, the slip values obtained from the trishear inversions have limited variability around the values with the highest weight, showing some dispersion only for the solutions with lower weights. In most cases, the solutions are located on the mapped faults or FIGURE 7 | (A) Depth-converted version of the Northern Section. Some minor faults (e.g., backthrusts) were not included in the restoration after having verified that their contribution to the total displacement was two orders of magnitude smaller. The main restoration steps are provided in Supplementary Figure 2. A summary of the slip-rate calculation results for the studied structures is shown by plots of the time-interval ages (B), the recovered slips (C), and the slip rates (D) for each fault. A boxcar distribution describes the interval-age uncertainty Probability Density Functions (PDF), whereas a Gaussian distribution describes the slip and slip-rate uncertainty PDFs. X-axes are in logarithmic scale for displacement and slip-rate values. Only the horizons considered for the slip rate calculations are represented (see Table 3 and Supplementary Data 2 for numerical information). within a 500 m distance from the faults. The slip values are coherent with the amount of deformation of the relevant horizon, being small (<100 m) for the younger and less deformed horizons and higher for older and more deformed horizons.
The slip rate calculations results, along with their uncertainty PDFs, for each fault and each restored horizon are graphically shown in Figures 7, 8 for the Northern Section and the Southern Section, respectively. Table 3 summarizes the horizon age, restoration method and restored slip, and slip rate results. In five cases, the displacement cannot be calculated because the relevant horizon is part of the Pleistocene Progradation Complex. For these five cases, we retained the slip associated with the underlying horizon. This missing information seems not to affect the general trend of the investigated structures.
Overall, the progressive restoration showed that the inception age is different for different faults in agreement with the inception age visual estimations based on the seismic stratigraphy.
Considering the structural relative positions, the inception age is progressively younger from the inner toward the outer structures. NAF1, NAF2, and NAF3 start in the Pliocene. The former two during seq-PL2 and the latter during seq-PL3. NAF2 remained active until the end of seq-PS1 or slightly afterward. NAF4 is the younger structure, starting its activity in the Pleistocene with the onset of seq-PS1 and has only a short period of activity overlap with NAF2 lasting about 1 Myr (2.5-1.5 Ma). The total cumulative slip is differently distributed on the different structures along the two regional cross-sections and varies from about 1,800 m for NAF4 to about 6,000 m for NAF3-South. This structure shows the highest structural elevation and amount of slip, probably due to the limited or almost null tectonic activity on the NAF2 structure after seq-PS1. The distribution of the cumulative slip over time is shown in Figure 10. Note that because of the adopted restoration methods (Figure 3), each discrete displacement increment is attained by each restored horizon at the bottom of each decompacted sequence or group of sequences (in the presence of clinoforms). Therefore, in Figure 10, the cumulative displacement distribution is attributed to the sequence termination age, regardless of the possibly finer internal displacement variations that could have occurred within the sequence or, in other words, between two successive horizons. This distribution can be approximated with an exponential function whose local slope and slope changes represent the slip rate and its decay over time. The fit parameters, along with some function features, are reported in Table 4.
The obtained slip rates for the various horizons span two orders of magnitude, from a few millimeters to a few hundredths of millimeters per year. The faster slip rates characterize the earlier stage of activity (Pliocene and Gelasian) and slower slip rates in the more recent times (Pleistocene), with the main rate change around 1.5 Ma within the deposition of seq-PS2. All faults retain some activity level until the present except for NAF2 that ceased its activity sometime after 1.5 Ma. Within the Pliocene (seq-PL1 to seq-PL3), NAF2 has rates of activity similar to those of the other structures. However, it must be noted that the Southern Section intersects NAF2 very close to its northern termination.
Further analyses are necessary to verify its activity to the south. Both cross-sections intersect NAF1 near the structural culmination; therefore, the deduced activity could be a good proxy to its maximum. NAF3 is a very long structure, and because of the lack of data between the two cross-sections, the two slip rate determinations NAF3North and NAF3South may represent the activity of two separate strands. NAF4, the youngest and outermost structure, has an onset slip rate slower than that of the older structures. However, its slip rate since hor-PS3c is almost twice as faster as the other faults in the same period.

DISCUSSION
We analyzed a large dataset of 2D seismic profiles to perform a structural reconstruction of an offshore sector of the Northern Apennines thrust front. The reconstructed faults generally agree with the tectonic lineaments and main faults identified in previous studies (Bigi et al., 1992;Fantoni and Franciosi, 2010;DISS Working Group, 2018;Pezzo et al., 2020). Nevertheless, the quality, density, and spatial distribution of the dataset used here allowed us to improve the 3D definition of the fault geometries. Relative to the DISS Working Group's (2018) fault compilation in the area, we characterized some new faults, such as NAF1, NAF3, NAF4, and NAF6, to be considered and provided data to improve some of the already mapped. Overall, our reconstructions provide a better definition, especially of the fault upper edges, whereas the deeper parts of the faults still require extrapolation unless data with deeper penetration become available. These improved fault geometries can help develop a new generation of seismotectonic models and earthquake hazards applications at the local scale.
As regards the seismotectonics of the area, the elusiveness of geomorphic signatures of tectonic activity in the bathymetry makes it difficult to assess the recent activity of the thrusts that has remained a matter of debate for years (e.g., Di Bucci and Mazzoli, 2002;Pezzo et al., 2020). The prevalence of sedimentary accumulation relative to the frontal anticline growth during the Pleistocene, as already observed within the Po Plain (Maesano and D'Ambrogi, 2016), causes the complete burial of the fault upper ends and the fault-related folds, thereby challenging the direct recognition of tectonic activity.
The present-day stress field acting in the Northern Adriatic Sea is also a matter of debate. Mazzoli et al. (2015) proposed maximum horizontal stress with orientation variable from NE-SW to NW-SE within a dominant transcurrent tectonic regime, based on minor seismic sequences in the area. Nonetheless, considering the distribution and clustering of stress data in the region, robust stress inversions indicate NE-SW to ENE-WSW directions of the maximum horizontal stress, which are almost orthogonal to the strike of the major contractional structures in the area (Carafa and Barba, 2013;Carafa et al., 2015b). The analysis of the seismological data of the largest past earthquake in the region (Senigallia, 1930, Mw 5.9) is consistent with a compressive stress field (Vannoli et al., 2015), as also confirmed by the focal mechanisms of the smaller seismic sequence (max Ml 5.0) recorded in the offshore of Porto San Giorgio (south of the study area) in 1987 (Riguzzi et al., 1989;Battimelli et al., 2019). This evidence suggests that the mechanism that produced the weak but not negligible deformation observed for the post-hor-PS3c (0.4 Ma) period is coherent with the current stress regime.
The structures with a 5-10 km wavelength of the outer Northern Apennines share a common basal thrust located on the top of the Meso-Cenozoic carbonatic succession (Barchi et al., 1998;Fantoni and Franciosi, 2010;Maesano et al., 2015). The flat thrust geometry is favored by the presence of a regional detachment level corresponding to the Gallare marls (upper Eocene to lower Miocene), a highly plastic formation locally FIGURE 10 | Cumulative displacement distribution versus age obtained from data listed in Table 3. Based on the adopted method, the displacement is attributed to each restored horizon at the bottom of each sequence, regardless of the displacement distribution during the sequence deposition. Error bars on the y-axis represent the error of interval displacement propagated from the present to older ages by adding in quadrature. Error bars on the x-axis represent the error on the age determination at the end of the interval. All datasets are fitted by an exponential function in the form y = a + be −x (solid line) shown along with their 95% confidence interval (dashed lines). See Table 4 for the fit parameters. characterized by overpressure zones (Livani et al., 2018). The presence of an Eocene-Miocene detachment level is a common characteristic of the entire Po Plain and Northern Adriatic Foreland (Fantoni et al., 2004;Massoli et al., 2006;Toscani et al., 2006;Maesano and D'Ambrogi, 2016;Turrini et al., 2016) and, together with a deeper decollement level located in the Triassic evaporites, is a key factor in controlling the structural style of the whole region. Basili and Barba (2007) estimated an average long-term (17 Myr) shortening rate of 2.9 mm/year of the structures above this basal thrust. Numerical modeling fitting GPS recordings and horizontal stress directions, both onshore and offshore, demonstrated the necessity of basal tractions to justify the observed ongoing deformation (Barba et al., 2008), which requires slip-rates between 0.6 and 2.5 mm/year on the basal thrust (Carafa et al., 2015a;Pezzo et al., 2020). The deformation along the common detachment should be transferred and distributed onto the outer frontal ramps. Our work focused on improving the geometric definition of these fault ramps to provide a numerical representation of how the deformation is distributed in both space and time.
As highlighted by the analysis of two regional cross-sections, we documented the asynchronous activation of four thrusts (NAF1, NAF2, NAF3, and NAF4), with a younger inception age for the outer faults as already observed in other sectors of the    Northern Apennines, the Po Plain, and the Adriatic foredeep (Barchi et al., 1998;Scrocca, 2006;Scisciani and Calamita, 2009;Fantoni and Franciosi, 2010). Our findings thus confirm a general trend that characterizes the Northern Apennines fold-and-thrust belt in its entirety. Importantly, the recent and ongoing tectonic activity can be effectively examined only after having decompacted and restored the depth-converted horizons. The removal of the differential compaction effect is necessary for the correct estimates of the The first derivatives at inception age, at 1.5 My, and at present provide a discretized evaluation of the fault slip rate as results from the statistical model. deformation of fault-related folds and growth strata. The impact of differential compaction in this type of structural setting could lead to an overestimation of the slip required to restore the fault geometry (Scrocca et al., 2007;Maesano et al., 2013Maesano et al., , 2015. After the decompaction and restoration procedure, all the analyzed structures showed evidence of activity throughout the Pleistocene. NAF2 is the only structure for which we did not find deformation evidence after the early Calabrian (hor-PS2). This result could either be due to the Southern Section position that crosses NAF2 near its northern lateral termination where its slip is minimum or to a migration of activity determining an accommodation of the deformation by the outermost structures (NAF3South). We underline the importance of two main requirements when dealing with inverse problems for fault slip assessment, implying a non-unique solution . The first requirement is the use of horizons that need to extend across both walls of the fault. The effect of using incomplete datasets has already been evaluated and tested by comparing results in natural and analog experiments, showing that a lower spread of the results is obtained when the restored horizon includes a significant portion located on both the hanging wall and footwall of the underlying fault . The second requirement is the use of well-constrained and accurate fault geometries. During the inverse restoration with the trishear method, we observed that some of the solutions with a good performance from a numerical viewpoint (i.e., low residuals) require a fault tip located far from the actual fault (Figure 9). Since seismic reflection profiles very well constrain the fault position and geometry, we adopted a solution that weighs both the distance from the fault and the residuals on the necessary slip. The solutions closer to the actual fault trace are preferred in the slip evaluation. One of our inversion results is that the spread of slip values is quite limited, showing some dispersion only for low-weighed solutions.
The treatment of slip-rate uncertainties in tectonic contexts involving slow-moving, blind or buried, faults is one of the major outcomes of this work. In our case study, we adopted age constraints provided by literature data in most cases (Table 1), and their uncertainties were related to the dating method. The age of the most recent horizon, hor-PS3c, instead was based on a detailed interpretation of the most recent reflectors, then correlated and compared with the interpretations and mean sedimentation rates calculated in neighboring sectors of the Adriatic Sea (Amadori et al., 2020). Although this horizon does not have a strong stratigraphic constraint from previous literature, it was important to confirm the recent tectonic activity. We thus recommend using more data on sedimentation rates because they can provide further constraints to the horizon definitions that could not be performed in this work. With the data at hand, we had to treat the age uncertainty as maximum and minimum values, in which all values within that range are equally likely (boxcar PDF). Combining the age boxcar PDF with the slip Gaussian PDF resulted in a probabilistic distribution of the slip rates, representing a significant improvement in treating this epistemic uncertainty relative to previous works in the area. However, the possibility of incorporating more horizons with peaked uncertainty distributions, like a Gaussian, also for the age term would lead to narrower uncertainty in the slip rate determination and a further potential improvement. When these improvements regard the most recent fault activity stages, they can have compelling implications when such data are used as input in earthquake hazard applications because the slip rate parameter controls the tectonic moment rate that is then converted into annual earthquake rates (see Morell et al., 2020, for an overview of these relationships).
In this work, we also estimated the slip rate variations over time. Slip rates calculated for different time intervals provide evidence for an exponential decrease over time for all the investigated structures. Over the Plio-Pleistocene, such variation spans two orders of magnitude (Table 3), from a few millimeters per year during the Pliocene to early Pleistocene, coincident with the early stages of fault activity depending on the inception age, to only a few hundredths of millimeters per year during the Middle-Late Pleistocene. For example, NAF3-South is the structure with the largest slip-rate variation through time, starting its activity at more than 3 mm/year during seq-PL3, then slowing down to 0.05 mm/year during seq-PS3. NAF4, which is much younger, starts its activity at a much slower rate of ∼1.3 mm/yr during seq-PS1 and then slows down to ∼0.2 mm/year while remaining the faster structure during seq-PS3 ( Table 3). The amount of restored deformation for seq-PS3c (post 0.4 Ma) is very small but not negligible in all cases. One must be aware, though, that these values are nearly in the same order of magnitude as the resolution limits of the adopted methodology, which is affected by the inherent limitations due to various factors, such as the manual picking of the seismic reflectors and the homogeneous decompaction parameters of the deposits.
For all the analyzed structures, we observed a marked slip rate reduction occurred during the Calabrian (seq-PS2), evidenced by the slope change of the fitted curves shown in Figure 10, in agreement with a general reduction of the Northern Apennines thrust fronts activity observed during the Pleistocene (Gunderson et al., 2013;Maesano et al., 2013Maesano et al., , 2015Bresciani and Perotti, 2014;Maesano and D'Ambrogi, 2016) coincident with the overfilling of the Po Plain -Northern Adriatic foredeep beginning in the Early Pleistocene (Gunderson et al., 2018). On average, slip rates obtained here are consistent with previous works and findings in adjacent areas (Maesano et al., 2013(Maesano et al., , 2015 and represent, in most cases, a finer temporal scan of thrust activity. However, although this trend seems to characterize most Northern Apennines thrusts, we must be aware of the exceptions such as the Mirandola Thrust, buried under the Po Plain, whose slip rates might have increased during the Middle-Late Pleistocene (Maesano et al., 2015) and be cautious in generalizing these results. Discriminating between decelerating and accelerating thrusts is key to understanding the recent Northern Apennines evolution and providing elements for mindful comparison with activity rates derived from modeling decadal geodetic (GPS and InSAR) observations and estimates of seismic efficiency (Carafa et al., 2017;Caporali et al., 2018). In this perspective, we note that the slip rates measured on the most recent horizon (hor-PS3c), and thus representing the last 0.4 Myr, for all the thrust ramps intersected by the two regional sections, collectively account for a fraction (25-50%) of the slip rates on the basal detachment (Carafa et al., 2015a;Pezzo et al., 2020). This apparent discrepancy can be explained considering that part of the basal detachment strain can be distributed on the innermost coastal or onshore active thrust (Maesano et al., 2015). The slip partition between disconnected faults and their common detachment may shed light on the mechanism that drives the slip-rate variability in the million-year timescale and contribute to better understanding the relative roles of the different factors controlling this process (Roy et al., 2016;Gunderson et al., 2018), including exogenous forcing due to the competition between erosion and deposition of the rock volume above the blind faults. In our study area, the main sediment supply comes from an axial system in an almost NW-SE elongated foredeep (Ghielmi et al., 2010(Ghielmi et al., , 2013Amadori et al., 2019) where the Pleistocene sedimentation rates are fast enough to conceal the tectonic signature. The kinematic relationships between a fold-and-thrust belt and syntectonic sedimentation are crucial in controlling the thrust geometries and activity (e.g., Storti and McClay, 1995;Duerto and McClay, 2009;Fillon et al., 2013). Syntectonic sedimentation could likely play an important role also in our case study by acting as an exogenous forcing in affecting the thrust kinematics. Indeed, the Pleistocene foredeep sediments mainly come from the northwest, i.e., from outside the thrust belt and foreland system. Looking at the regional cross-sections (Figure 6) and the slip-rate values ( Table 3 and Figure 10), one may notice that the clinoform deposition occurs while the tectonic signal dramatically decreases, suggesting a chronological correspondence between the high sedimentation rates and the low slip-rate values. While the increasing synkinematic sedimentation rate delays the development of frontal contractional structures and favors the formation and reactivation of more internal thrusts and backthrusts (Pla et al., 2019, and references therein), the presence of clinoforms typically indicates periods with reduced tectonic activity when new accommodation space in the basin is not created. Therefore, our data do not allow us to discriminate if the reduced tectonic activity after 1.5 Ma is the cause or the consequence of the sedimentation rate increase, but it is worth to be noted that these two processes have already been observed in large parts of the Northern Apennine outermost fronts (Maesano et al., 2013(Maesano et al., , 2015Amadori et al., 2020).

CONCLUSION
In this work, we emphasized the critical role played by subsurface data in characterizing the activity of buried faults, where sedimentation rates overtake tectonic rates, counteracting the geomorphic signature development at the ground surface or the seafloor.
Using a dense dataset of 2D seismic reflection profiles and wells, we performed a 3D structural reconstruction of eleven buried thrusts of the Northern Apennines chain in a nearoffshore sector of the northern Adriatic Sea. The investigated structures are consistent with previously published studies, but a well-constrained 3D reconstruction of the fault surfaces provided elements for improving the local structural setting description that can be useful for further seismotectonics and earthquake hazard analyses. The study area location is also a key factor in creating the opportunity to homogenize the results of thrust tectonics studies across the onshore-offshore transition of the Northern Apennines frontal part.
For four of the mapped thrusts, we documented the slip history over the last 4 Myr, revealing an exponential slip reduction implying a decelerating slip-rate pattern from a few millimeters to a few hundredths of millimeters per year with a major slip-rate change around 1.5 Ma. The slip histories revealed an asynchronous activation proceeding from the inner to the outer thrusts and that the youngest structure is also the fastest during the last 0.4 Myr interval. Comparing these findings with region-wide literature data confirms that, beyond the observed slip-rate variations within each fault, the strong slip-rate deceleration after 1.5 Ma is a widespread behavior of the Northern Apennines thrust faults. Yet, we are aware that also accelerating faults might exist in the thrust belt. Therefore, further studies should better address this occurrence to answer the remaining open questions on whether we can consider the deceleration due to a common regional forcing or just aleatory fluctuations of individual fault behavior.
The eight stratigraphic horizons with well-constrained age determinations (Zanclean to Middle Pleistocene) provided a slip-rate history reconstruction spanning a temporal scale that suits various geodynamic model needs. Slip rate values obtained through the analysis of older geologic markers could lead to apparent discrepancies when compared with deformation rates derived through geodetic data analyses. The slip-rate values in the most recent interval are the most appropriate to reconcile such possible discrepancies and would be useful to highlight how the geodetically derived basal thrust strain can be partitioned on different structures.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because the well-log and 2D seismic reflection profiles data analyzed in this study were obtained from ENI E&P, under a non-exclusive confidentiality agreement (https://www.eni.com/ en-IT/contacts.html). The VIDEPI dataset can be found online at www.videpi.com. Requests to access the datasets should be directed to https://www.eni.com/en-IT/contacts.html.

AUTHOR CONTRIBUTIONS
YP: restoration of regional cross sections, trishear analysis, slip rate calculations, and figure preparation. FM: conceptualization, performed seismic interpretation, 3D reconstructions, trishear analysis, and slip rate calculation. CA: stratigraphic horizon age definitions, reconstruction of the 3D model, stratigraphic setting, and preparation of related figures and tables. JF: seismic interpretation and 3D reconstructions. GT: conceptualization, performed seismic interpretation, 3D fault reconstruction, restoration of regional cross sections, and managed part of the funding. RB: conceptualization, slip rates calculation, supervised the work, revised the final version of the manuscript, and managed part of the funding. All authors actively contributed to writing various parts of the manuscript according to their individual specialty and agreed with the final content.

ACKNOWLEDGMENTS
We thank ENI for the permission to consult the reflection seismic profiles and wells datasets; Manlio Ghielmi for the fruitful discussions which contributed to constrain the seismic profile interpretation, the stratigraphic horizon age definitions, and the geological reconstruction, as well as for the feedback on an earlier version of the manuscript. We also thank Mauro Coltelli for coordinating the INGV team of the SPOT project and the SPOT consortium members for insightful discussions on fault activity and related earthquake hazards. Our views, opinions, and results expressed herein do not necessarily state or reflect MISE's and MIUR's or any agency thereof. Petroleum Experts are acknowledged for the donation of MOVE academic licenses to the University of Pavia. We are grateful to Nestor Cardozo for his suggestions and help on the use of his code. FM and RB acknowledge the resources made available by the SISMOLAB-3D at INGV.