From cylindrical to non‐cylindrical foreland basin: Pliocene–Pleistocene evolution of the Po Plain–Northern Adriatic basin (Italy)

The architecture of foreland basins and the resulting distribution of clastic sediments are related to the constant interplay between tectonics and sedimentation. Specifically, basin floor modifications strongly influence dimensions, continuity and connections of sand‐size and fine‐grained deposits. Given the increasing need to identify deep potential reservoir deposits, the large‐scale definition of clastic porous targets and their seals is a matter of interest for oil and gas industry. Here, we present the reconstruction of the Po Plain and Northern Adriatic Foreland Basin (with an extent of ca. 40,000 km2) and its Pliocene–Pleistocene evolution, as an example of a sedimentary clastic system controlled by strongly non‐cylindrical foreland geometry. The study is based on the basin‐scale mapping of six unconformity‐bounded sequences, performed by interpreting a dense network of seismic lines and correlating well‐log data. This provides a three‐dimensional model of the step‐by‐step evolution of the basin and a description of the sediment dispersal pattern. We found that the basin records the change from a continuous (cylindrical) to highly fragmented (non‐cylindrical) foredeep geometry during Late Pliocene. In the Northern Apennines case, the main factors driving the development of a non‐cylindrical geometry are mainly related to inherited inhomogeneity in the downgoing block linked to its Mesozoic extensional faulting, and the relative orientation of these lineaments with respect to the direction of orogen migration. During the late Pliocene–Pleistocene the two directions progressively became close to parallel, and the Northern Apennines system reacted changing from a cylindrical to a non‐cylindrical state.


| INTRODUCTION
The stratigraphic architecture of foreland basins provides the record of the interplay between tectonics and sediment supply in front of collisional belts. According to the broadest definitions (Allen, Homewood, & Williams, 1986), foreland-foredeep systems are described as sedimentary basins tectonically driven by the formation of a thrust belt and their longitudinal dimensions are almost pair to the length of the thrust belt (DeCelles & Giles, 1996). The 3D shape of the basins is in fact the result of the dynamic interaction among (a) thrust belt tectonic load, (b) distribution of subsidence in the foreland, and (c) sediment supply, influenced by both tectonic and climate, sourced by the hinterland and/or foreland. In many cases worldwide (e.g., in the Himalaya, Northern Alps, Taiwan, Betic Range, Southern Bolivia; see Bigi, Castellarin, Coli, Dal Piaz, & Sartori, 1990;Bigi, Castellarin, Coli, Dal Piaz & Vai, 1990;Burbank, Beck, & Mulder, 1996;Garcia-Castellanos, Fernandez, & Torne, 2002;Lin, Watts, & Hesselbo, 2003;Uba, Heubeck, & Hulka, 2006;Garzanti, 2019), foreland basins show a high lateral continuity and thickness homogeneity of sediment filling along-strike, over distances of hundreds to thousands of kms; this geometry can be defined as "cylindrical". In other cases, foreland basins show geometries far from that "cylindrical" state, as foredeep sediment thickness rapidly changes along-strike over tens of kilometres, recording the existence of strongly irregular depocentres. This type of geometry has been defined as "non-cylindrical", explained with the occurrence of along-strike basement heterogeneities in the downbending plate and/or along-strike tectonic variations in both lower and overriding plate (Catuneanu, 2004;DeCelles & Giles, 1996;Livani, Scrocca, Arecco, & Doglioni, 2018;Lucente & Speranza, 2001;Turrini, Toscani, Lacombe, & Roure, 2016). A non-cylindrical foredeep may reflect a break in the pattern of progressive outward propagation of the distal wedge, with jumps of deformation expressed as salient development across the foreland basin. Several studies approached the explanation of this "non-self-parallel" behaviour, though analogue and mathematical modelling, all highlighting the fundamental control of the stratigraphy on the deformation (Lickorish, Ford, Bürgisser, & Cobbold, 2002;Willett & Schlunegger, 2010). As most orogenic wedge-foredeep systems are exposed and variably eroded, tectonic phases and basin geometry are often poorly constrained from sedimentary record. This is a common issue during their reconstructions. By the contrary, subsurface investigations have the potential to overcome this limitation and develop a comprehensive understanding of the structures driving their geometric evolution. Specifically, a cylindrical versus non-cylindrical state of a foredeep basin strongly influences dimensions, continuity and connections of sand-size and fine-grained deposits. Given the increasing need to identify deep potential reservoir deposits, the large-scale characterization of clastic porous targets and their seals is vital for oil and gas industry. For all these reasons, detailed studies of well-constrained foredeep architectures and the understanding of if, when, and why they experienced the transition to a noncylindrical state, can improve both the general knowledge of orogen-basin dynamics and hydrocarbon exploration potential. The current architecture of the Northern Apennines foreland basin (Po Plain-Northern Adriatic foreland basin, hereinafter PPAF; Figure 1a) provides a classical example of non-cylindrical foreland basin. Due to the fairly continuous subsidence and high sediment supply related to the surrounding rising belts, the PPAF has been the site for the accumulation of a remarkable volume of sediments: up to 7-8 km-thick over thousands km 2 (Pieri & Groppi, 1981). This buried sedimentary framework records the time length of the complex evolution of the basin. The PPAF also hosts diverse hydrocarbon and water reservoirs (e.g., Cazzini et al., 2015 and references therein). For this reason, the basin has been the object of intensive geophysical surveys and drillings by oil companies, as well as scientific studies based on parts of these subsurface data (DISS Working Group, 2018;Fantoni & Franciosi, 2010;Ghielmi et al., 2010;Ghielmi, Minervini, Nini, Rogledi, & Rossi, 2013;ISPRA, 2015;Maesano, D'Ambrogi, Burrato, & Toscani, 2015;Maesano & D'Ambrogi, 2016;Toscani et al., 2014;Turrini, Lacombe, & Roure, 2014;Turrini et al., 2016 and references therein). Notwithstanding different projects spread over the study area and the extensive dataset acquired by industries, due to the fact that it remained reserved up to recent time, a homogeneous foredeep-scale reconstruction is still lacking. In this paper, we present the Pliocene-Pleistocene basin-scale (about 40,000 km 2 ) architecture of the PPAF, which can be considered a type example for non-cylindrical basin system by reason of the preservation of deformation and the extensive subsurface dataset available. Here, we describe the subsurface architecture through six maps of unconformity surfaces, obtained after the interpretation of seismic reflection profiles, calibrated and depth-converted with deep well stratigraphies and logs (Figure 1b). Thickness variation analysis (i.e., isochore maps) gives clues about the sediment distribution of the PPAF through time, describing the basin transition from a cylindrical to a non-cylindrical geometry during Pleistocene. This provides the opportunity to discuss controlling factors like the pre-orogenic settings, differential accommodation space availability and sediment distribution across the basin.
S a n t e r n o F m . EracleaFm.

| Seismic reflection profiles
The dataset we used in this study is part of the extensive subsurface dataset, composed of more than 5,000 reflection seismic profiles, collected over 60 years of hydrocarbon exploration activity in the PPAF (Cazzini et al., 2015). We selected and interpreted ~ 8,000 km of seismic reflection profiles from this dataset as the basis for the geological model presented here. In detail, our seismic grid is composed of 327 basin-length time-migrated 2D seismic reflection profiles, covering ca. 40,000 km 2 of the central-eastern Po Plain and the Northern Adriatic Sea (Figure 1b). Selection criteria for seismic grids are primarily based on quality, orientation of the seismic line and homogeneous covering of the basin, especially the offshore side where just few 2D seismic lines are currently available in digital format. Sedimentological interpretation of the seismic reflection profiles is mainly based on the description and interpretation of well logs coupled with the analysis of seismic geometries and facies (see also supplementary material). Sequence bounding surfaces of different order and importance were identified and correlated, at the basin scale, using well and seismic data and applying principles and procedures of seismic stratigraphy interpretation (Catuneanu et al., 2009;Vail et al., 1977 and references therein). The major unconformities are indicated in seismic profiles, along the basin margins and in correspondence of active tectonic structures, by seismic horizon terminations (i.e., onlap, downlap, erosional truncation, etc.). The unconformable nature of the surfaces was verified by detailed well-log correlations and by the biostratigraphic data of the available wells.

| Well logs
We performed lithologic, sedimentological interpretation and the stratigraphic calibration of the seismic reflection profiles (in time domain) using more than 200 explorative wells. The well logs are stored into the ENI database and public ViDEPI Project database (available from http:// unmig.svilu ppoec onomi co.gov.it/videp i/) ( Figure 1b). The ENI well database contains updated primary information such as formation name, age, foraminifer and nannoplancton biozones, sedimentological core descriptions and well velocity data. In this work, the integrated analysis of all well data, including lithologic logs and associated petrophysical parameters (e.g., resistivity, spontaneous potential and sonic log), enabled a detailed correlation of the wells over the basin. The lithologic, stratigraphic and well-log markers of the wells constrained the horizon picking on the seismic profiles, by the use of a swift depth-time conversion. Location of the public wells and a direct link to log data are available in the supplementary material (file G).

| Seismic velocities
For this research, we have selected 65 wells to support the construction of a 3D velocity model used to convert the mapped surfaces from time to depth (Figures 4 and 5). The selected well logs are those which include velocity data (i.e., interval velocity and average velocity) from which velocity-depth curves and all the parameters needed as input values for time to depth conversion can be derived (see the following chapter).

| Workflow
The methodological approach is based on the interpretation of a dense grid of seismic lines calibrated and time-depth converted with a 3D instantaneous velocity model (Marsden, 1992). The designed workflow is articulated in time and depth domains.
Each domain is characterized by separate successive steps, sometimes connected ( Figure 4). The main steps of the workflow can be summarized as follows: (a) Seismic and well-log data selection from private and public database (ENI and ViDEPI Project, respectively), (b) Interpretation of the key horizons in 2D seismic lines and well-log analysis, (c) Construction of the 3D model in time domain, (d) Calculation of the 3D velocity model, (e) Time-depth conversion, (f) Construction of the first 3D depth-converted model, and (g) Refinement through comparison with available well stratigraphies and construction of the final 3D depth-converted model. This is a relatively standard procedure that does not require further description, also because it has been applied in similar contexts (Geomol Project; ISPRA, 2015). Details on procedure followed for passing from 2D seismic interpretation to 3D surface in the time domain, and for time-depth conversion of the 3D model in time domain, are described in the following section.

| 3D model in time domain and time to depth conversion
To create a three-dimensional time-surface from the 2D interpretation, we chose the Delaunay triangulation method because it honours all the input points, it optimizes the geometry of the triangles and, it fulfils the nearest neighbour relation. The 3D model in time includes also fault segments, but they were neither interpolated nor converted in depth domain. Dedicated structural analyses, like slip rates calculation, horizontal restoration, are not the aim of this research. After the 3D interpolation of all unconformities, the surfaces have been smoothed using a search radius of 3 km and resampled with 1 km cell size. The 3D architecture in time domain has been depth-converted using the Vel-IO 3D tool developed by Maesano and D'Ambrogi (2017) and tested in the central part of the Po Basin (ISPRA, 2015). We have analysed velocity data of the 65 wells (Figures 1b and 5a), extracting a linear velocity curve for each layer. This, provides a value of initial velocity (V 0 ) and a gradient (k) of increasing/decreasing velocity with depth (Figure 5c, d) for each well log. In the velocity model building, we imposed a barrier along the outer thrust front of the Northern Apennines. The barrier divides the study area into two sectors characterized by different tectonic settings: the foredeep-foreland area to the north-east, and the frontal part of the Apennines tectonic wedge with its thrust-top basins to the south (Figure 5a). The southern margin of the Po basin (i.e., Northern Apennines fronts), presents the most significant variability of velocity parameters V 0 and k due to irregular geometry, thickness variations and lithofacies distribution with respect to the more regular and homogeneous (in terms of reduced variability of the velocity parameters) northern and eastern depocentre zone.
The velocity model adopted in the time to depth conversion considers the PL1u as the bottom in the entire study area. In the onshore area, the velocity model is given by two layers considering the PS1u as the internal boundary whereas, in the offshore area, the model also includes the water column as additional third layer (Figure 5b). The bathymetry was obtained by converting the SRTM30 grid (Becker et al., 2009) falling within the study area, in time domain using the constant sound velocity in sea water of 1,500 m/s. We also used the "Velocity model optimizer" of Vel-IO 3D to improve the consistency of the velocity model obtained in the previous step. During the optimization, the velocity model includes the stratigraphic markers from wells without velocity information that locate far from the original velocity data. These wells work as "control points" for which the optimization process iteratively calculates the V 0 and k values. The model accepts these new velocity parameters only if they provide a calculated depth of the stratigraphic marker within a range of ±50 m from the measured depth in the well. After the optimization, the new velocity and gradient values, together with F I G U R E 5 (a) Map of the study area, well-log data with velocity time-depth table (red points) and imposed barriers for 3D time-to-depth conversion (bold black line). the existing velocity data, allow to work out the optimized grids of the velocity model (velocity maps of the key layers are available in the supplementary material, Fig, A). Finally, the Vel-IO 3D runs the time-depth conversion of the six mapped unconformities, automatically assigning each point to its correct layer in the 3D velocity model. We estimated the uncertainty in the time-depth conversion by analysing the variability of the initial velocity (V 0 ) across the entire area, and performing two additional depth conversions by which we added and subtracted, respectively, the standard deviation of the initial velocity. The resulting variability of the Z value is strongly dependent by the absolute depth of the data point itself. On average, the relative variability of depth is of 6% for PS1u and 9% for PL1u. The variability of depth is higher in the structural high (i.e., Ferrara arc), whereas the average absolute uncertainty is 100 m for PS1u and 163 m for PL1u, with the maxima located in the deeper depocentres where no constrain from well logs is available (see supplementary material, Figs. B, C). To reduce possible boundary effects, due to both 3D interpolation and time-depth conversion, the most external boundary sectors and the Northern Apennines piggyback depocentres, have been removed because seismic grid and/or well-log data were not dense enough to support outputs with the same quality of the foredeep sectors. 6 | RESULTS

| Seismic interpretation
This section describes the main tectono-stratigraphic characteristics of the PPAF, by means of some representative seismic reflection profiles over the study area (Figure 6a) (non-interpreted seismic images and facies analysis are shown in Figs. D, D2 into the supplementary material). Their most important features are briefly described following their geographic distribution. The north-eastern margin of the PPAF (i.e., Northern Adriatic Sea and Venetian region) represents the undeformed foreland, not involved into the propagation of Southern Alpine and Northern Apennines thrust and fold belts. In Figure 6b, the PL1s is made of a southwards prograding complex (Eraclea Formation) that overlies the Late Miocene succession or, locally, the Mesozoic carbonate substratum. An opposite, northwards, Late Pleistocene prograding complex (PS3s) is visible in the upper part of the section. Furthermore, along the S-dipping Venetian foreland ramp, all Plio-Pleistocene surfaces onlap on the PL1s (made of shales) or directly on PL1u. More in general, Pliocene turbidite foredeep deposits (PL1s, PL2s and PL3s) thin northwards, with a progressive NNE-ward migration of the pinch out on the foreland ramp (Figure 6b, c, g). The linedrawing in Figure 6c derives from composite deep seismic sections. From SSW to NNW, it crosses a thick piggyback basin and the entire Ferrara structural highs. In this transect, PL2s, PL3s and PS1s strongly thin towards NNE onto the anticlines crests of the Ferrara Arc. Here, stratigraphic relationships among the Pleistocene sequences outline an Early Gelasian deformation PS1 phase, responsible for the final isolation of the large NA piggyback from the foredeep axis. At the same time, in the western Po Plain, PS1 deformation enhanced the isolated thrust-related Piadena Anticline (P.A in Figure 1a), active since Intra-Zanclean PL2 phase (Figure 6i). In the foredeep, the PL1s is represented by the sand-rich Canopo Formation, identified by converging reflectors with high noisy signal contrasting in amplitude and seismic pattern. Diversely, the overlying PL3, PS1 and PS2 sequences are characterized by parallel reflectors and high lateral continuity (Figure 6d). These sediments, deposited over a limited area, form a complex of longitudinal highly efficient turbidite systems interfingered with coarse-grained turbidite systems fed by lateral entry points located in the northern margin (i.e., poorly efficient turbidite system, sensu Mutti, 1985;Mutti et al., 1999). In Figure 6e, seismic facies of PL1s and PL2s units show high-amplitude converging reflectors, reflecting their coarse-grained proximal lobe setting fed by some entry points along the Mantova Monocline (Figure 1a for location). By the contrary, Late Pliocene PL3 and Pleistocene PS1, PS2 and PS3 seismic sequences, are characterized by parallel, lateral continuous reflectors, due to turbidite sand and shale interbedded layers. An exception is made by the foresets of the prograding Po Plain Complex, visible in the upper part of the PS2s and bounded at the top by the PS3u. Figure 6f and g extend over the Parma recess, from the NA margin with reverse faults (Figure 6f), to the S-dipping foreland ramp, where pre-Pliocene extensional faults are visible (Figure 6g). NA thrusts were active before Pliocene time until Calabrian, when the top of the anticline was overlain by coastal, shelfal and slope deposits prograding from the south to the depocentre of the basin. Here (Figure 6g), PL1u, PL2u and PL3u mark an abrupt landward shift of the turbidite facies on the ramp; they can be defined as "drowning unconformities", sensu Rossi, Minervini, and Ghielmi (2018). The Pleistocene Prograding Complex (upper PS2s and PS3s) has opposite progradation direction with respect to the Figure 6f, due to a local southward clinoforms migration from the northern margin.
Moving along-strike to the south, in the central part of the foredeep (Figure 6h), we recognize the deformations related to the growth of an anticline with steep limbs, affecting by early-middle Pliocene (PL2 Intra-Zanclean phase) and active until recent times. Evidence of the still ongoing deformation is the PS3s (here represented by fully continental fluvial deposits) folded lateral onlaps on the crest, visible at 300-400 ms. The anticline is controlled by two reverse faults, evolved in different periods as evidenced by different unit thicknesses at the footwalls. First, the SE-dipping one migrated in the early Pleistocene (PS1s) whereas the NW-dipping back-thrust formed only in the early-middle Pleistocene (PS1s-PS2s) (Figure 6h). In both Figure 6h and i, syntectonic growth strata of the PS2s and slightly deformation of PS3s suggest some tectonic activity during middle-late Pleistocene. However, minor offsets turn out in a decreasing in morphological complexity over the entire study area since Calabrian time, with minimum thickness changes along Piadena Anticline. In the Emilia Arc (Figure 6i), the Plio-Pleistocene succession is separated by the structural high formed by the thrust front active since Late Miocene, creating the piggyback Cortemaggiore basin with prominent angular unconformities. Convergence of the PL2u, PL3u and PS1u with sediment condensation on the fold crest indicate that the Cortemaggiore thrust front acted like a submarine high during Pliocene and early Pleistocene. Moving towards NE, the profile cross over the isolated thrust-related Piadena Anticline. The last regional unconformity (PS3u) shows strongly erosive features along the margins of the basin. PS3u lies at the top of the PS2s (Figure 6d to i) and at the base of a forced regressive depositional systems, the Pleistocene Po Prograding Complex (PS3s), that strongly prograded eastward during the middle Pleistocene Alpine glaciations (Figure 6b, c).

| Isobath and isochore maps
Our main result consists of the basin-scale three-dimensional reconstruction of the architecture of the Po Plain-Northern Adriatic basin for the Plio-Pleistocene interval, and its evolution during the last ca. 5 Ma. In the following, the basin architecture and the tectonostratigraphic evolution are described step-by-step by isobath maps of traced unconformities (Figures 7 and 8) and by isochore maps of the present-day units' thickness (Figures 9 and 10). Isobath maps describe the depth-converted unconformities with respect to the present-day sea level, meaning the datum is 0 m above sea level. Isochore maps display lines of equal thickness in the sequences, where the thickness is measured vertically between two unconformities.
Isobath maps (Figures 7 and 8) are shown with a variable colour because maximum/minimum depth values significantly change and to appreciate the small physiography changes of the shallower Pleistocene PS2 and PS3 horizons, which is smoothed if compared with the Pliocene ones. Conversely, a common colour scale was used for isochore maps (Figures 9 and 10) in order to make them easily comparable.

| Base Pliocene PL1 phase (Figures 7a and 9a)
The base of the Pliocene (PL1u) corresponds to the base of Zanclean stage, dated 5.33 Ma, and it is the deepest and most deformed horizon included in our reconstruction. It was affected by all Plio-Pleistocene Northern Apennines tectonic phases. PL1u reaches the maximum depth of ca. 8 km in the depocentre located between Emilia and Ferrara Arcs (between the cities of Parma and Reggio Emilia, PR and RE in Figure 7a). Shallower depths (in the range of 0.5-2 km) correspond to the articulated fault-related anticlines near Ferrara, to the area of the Trento Plateau ( Figure 2) and to the northern Adriatic foreland ramp going towards NE. The thickness of PL1s, recorded in the foredeep areas close to the Apennines chain, in the recess between Emilia and Ferrara salient (the same area where PL1s shows its maximum depth) is ~ 700-900 m. Presently, the early Pliocene depocentres appear highly irregular due to multiple subsequent tectonic phases, but a narrow ca. 200 km elongated basin, NW-SEoriented, can be still recognized parallel and very close to the Northern Apennines front (Figure 9a). The higher thickness localized in the north-eastern Adriatic foreland (nearby the city of Venice, VE in Figure 9a) is the expression of the S-prograding delta complex of Eraclea Formation ( Figure  6b). In this sector, the delta complex evolved after an erosive phase responsible for the foreland exposure during late Messinian sea-level drop (Amadori, Garcia-Castellanos, et al., 2018). PL1u with erosive features (V-shaped canyons in Figure 6b) lies at the base of the sequence, evolved during the post-Messinian Salinity Crisis mega-flooding phase (Micallef et al., 2018) and supplied by the eastern Southern Alps.

| Intra-Zanclean PL2 phase (Figures 7b and 9b)
PL2s was deposited during an intra-Zanclean tectonic phase. This phase deformed the innermost NA Messinian front and caused the first structural highs that will progressively evolve to form Piadena Anticline to the west and the Ferrara Arc to the east. PL2u depocentres (Figure 7b) shifted few kilometres northward, maintaining similar dimensions and geometries with respect to the PL1s depocentres. The PL2u maximum depth is 6-6.5 km, about 1-1.5 km shallower with respect to the PL1u (present-day depth). The isochore map of the PL2s (Figure 9b) allows to recognize remnants of a narrow (less than 50 km wide) ca. 200-250 km long cylindrical foredeep, NW-SE-oriented parallel to the Northern Apennines front.

| Intra-Piacenzian PL3 phase (Figures 7c and 9c)
The intra-Piacenzian PL3 unconformity, dated 3.3 Ma, is represented in Figure 7c. During the Late Pliocene tectonic phase, the eastern Ferrara Arc continued to develop, progressively uplifting and migrating NE-ward. The growth of the Ferrara anticlines is constrained by the thickness variation of the sediments (Figure 7c). PL3s is deformed and thinner on the top of the structural highs with respect to the adjacent depocentres where the sequence is thicker and regular. The deepest depocentre, located between the cities of Parma and Reggio Emilia and limited by the Emilia and the Ferrara salients, reaches a depth of about 6 km, whereas the PL3s maximum thickness is of ~ 1 km (Figure 9c). This variation of the detrital distribution pattern indicates an incipient fragmentation of the foredeep and marks the transition to a noncylindrical stage.

| Early Gelasian PS1 phase (Figures 8a and 10a)
The PS1 unconformity is dated ~ 2.4 Ma and recorded the first tectonic event during Pleistocene time (Figure 8a). Isochore map (Figure 10a) highlights the non-cylindrical foredeep geometry, different from the early Pliocene stages.
In the western part of the basin, sediment thickness abruptly changes in correspondence of the Piadena Anticline (P.A. in Figure 1a) (Figure 10a). The well-structured Ferrara Arc anticlines have a hinge at depth ranging between 0 and 1.5 km and bound two deep independent depocentres up to 4.5 kmdeep, filled by PS1s 1.2-1.5 km-thick. The western depocentre is ~ 50 km wide and 70-80 km long, whereas the eastern one also extends to SE in the Adriatic area, outside of the study area.

| Early Calabrian PS2 phase (Figures 8b and 10b)
The Calabrian PS2u, dated ~ 1.5 Ma based on the first appearance of the benthic foraminifera H. baltica in Po Plain (Dondi & Papetti, 1968), is associated with the continuous non-cylindrical geometry evolution of the foredeep, coupled with the northwards shift of the Adriatic depocentre. Since Calabrian, the eastern and wider depocentre migrated NE-ward to the Adriatic side, far from the Apennines thrust fronts getting closer to the Venetian foreland. This migration was a consequence of a combination of the complete uplift of the Ferrara structural fronts (Figure 6c) with a fast shift towards NE of the lower Pleistocene turbidites pinch-out. This is evident by the extremely reduced sediment thickness in the area north of Ferrara city (FE in Figure 10b). By contrast, in the Po Plain area, the western depocentre appears stuck in the same position as it was during Pliocene (Figure 8a, b). In this area, the smoother morphology over the growing Piadena Anticline and the quasi-homogeneous sedimentary thickness above it, attest the decreasing tectonic activity during Calabrian time (see also Figure 6i). At basin scale, isochore map of the PS2 sequence ( Figure 10b) shows great difference in thickness between western and eastern depocentres, even if they both lie at maximum 2.5 km depth (Figure 8b). The western one is filled with more than 1.5 km of sediments on the contrary, the eastern one is only 800-900 m thick. This confirms that the Ferrara Arc acted as a barrier; the Adriatic depocentre could not be directly fed by the sediment flux from the NW until the accommodation space of the western depocentre was completely filled. An important environmental transition occurred during late Calabrian, from the marine to continental. As visible in seismics within the upper PS2 sequence, the fluvial Po Prograding Complex covers the Pleistocene turbidites (Figure 6d to i).

| Late Pleistocene PS3 phase (Figures 8c and 10c)
The PS3u is dated 0.87 Ma and reflects the onset of the major Pleistocene Alpine glaciations (Muttoni et al., 2003). Because of the decrease in tectonic activity since middle Pleistocene, this uppermost key horizon has regular smoother physiography (Figure 8c). The western sector is shallow, only 0-500 m deep, with large nearly flat areas. Here, the PS3 sequence records a complete continental sedimentation linked to the Po river. On the contrary in the Adriatic Sea, a foredeep basin with marine turbidite sedimentation was still active, at 1.5-1.8 km depth. PS3u isobath map shows its highest topographic gradient in correspondence with the shelf break whose trace follows the line among the cities of Bologna, Ferrara and Venice (Figure 8c). Regarding the thickness of the depocentres infill, the PS3s isochore map (Figure 10c) records an opposite trend with respect to the PS2s thickness map (Figure 10b). In fact, up to 1.2-1.3 km-thick of shelfal, deltaic and turbidite sediments occupy the eastern Adriatic side whereas, the western onshore sector accommodates no more than 600-700 m thickness of Po Prograding Complex sediments.

| Plio-Pleistocene evolution of the PPAF
The reconstruction of the 3D architecture of the Po Plain-Northern Adriatic foreland basin during last ~ 5 Ma revealed progressive changes of the basin floor shape. In detail, by means of thickness analysis it was possible to understand how the clastic systems reacted related to basin modifications; from cylindrical to a non-cylindrical foredeep passing through a "transition" phase. Also, we associate these results with the long-term dynamic stages a foredeep sensu Catuneanu (2004), from underfilled to overfilled depositional state. By integrating outcomes from this work with facies distribution maps (Ghielmi et al., 2010(Ghielmi et al., , 2013 and age of activity of Apennines thrust fronts (Livani et al., 2018;Maesano et al., 2015), we have drawn the PPAF evolution in six steps with synthetic profiles (Figure 11). Note that profiles from PL1 to PS2 do not describe the present-day subsurface situations, but sketch of paleo-morphologies inferred by the integrated outcomes listed above. The only profiles showing the actual subsurface situation are A-A' and C-C' in PS3 (top map in Figure 11). During early Pliocene PL1 and PL2 stages, only the innermost Northern Apennines thrust fronts were active and the thickness distribution of deep water turbidites defined a NW-SE elongated single foredeep. This is visible from isochore maps showing a narrow and along-strike continuous sediment thickness (Figure 9a, b and bottom Figure 11). As a result, the basin is defined as cylindrical, parallel to the advancing orogenic wedge. In consideration of the high base level of the Mediterranean Sea, restored after Messinian Salinity Crisis (Micallef et al., 2018), and the presence of deep-water deposits only next to the main thrusts, we infer an underfill phase of the PPAF during early Pliocene time (profile A-A' relative to PL1 and PL2 stages in Figure 11). At the intra-Piacenzian PL3 stage, the NA tectonic activity evolved with non-self-parallel arcuate thrust-folds, located in the central part of the foredeep basin (profiles A-A' relative to the PL3 map in Figure 11). These structural culminations lately developed in the Ferrara Arc. Turbidite sediments distribution of the PL3s is consequently affected by this basin floor changes. The longitudinal siliciclastic flux started to be slightly deviated and, locally, replaced by low thickness of submarine highly condensed fine-grained deposits that took place in correspondence to the morphological highs (PL3 map in Figure 11). Despite this local complexities, at basin-scale view, seismic images and well-log correlations of the clastic sequence demonstrate that western and eastern depocentres were still connected In consideration of these characteristics, this stage can be defined as a "transition" from cylindrical to a non-cylindrical foredeep. During PS1, due to the Gelasian tectonic phase, the structural culminations located in the foredeep grew up and evolved to create the Ferrara Arc and a closed inner piggyback (Figure 6c). PS1 isobath and isochore maps (Figures 8a and 10a, respectively) describe two independent depozones: the western one is deeper than the eastern one and accommodated a thicker volume of sediments. A large sector of lower sediment thickness separated the depocentres, breaking completely the connection between them. Submarine highly condensed sediments deposited over this wider area, corresponding to the Ferrara Arc's culminations (PS1 facies map and profile A-A' in Figure 11). At this stage, the sediments supplied from the NW (in competition with the subsidence) were insufficient to fill the western PPAF, pass over the Ferrara Arc and reach in the Adriatic area. The early Calabrian PS2 tectonic phase did not show substantial variation with respect to the PS1 phase in terms of horizontal shortening of the Ferrara salient. The weak deformations of the PS2u, shown in the Ferrara Arc inner sector ( Figure  6c) and the Piadena Anticline (Figure 6i), testify the reduction in tectonic activity of the Northern Apennines. Offsets of the PS2 unconformity can only be seen in the north-easternmost reverse faults of the external Ferrara Arc (Figure 6i). PS2s isochore map still highlights the remarkable difference of thickness between in the thicker western depocentre with respect to the eastern one, similarly to the PS1s (Figure 10b). The lack of sediment transfer along-strike to the south-western region confirms the disconnection of the Adriatic side and the non-cylindrical geometry of the PPAF. Despite the decrement in basin shape modifications, clinoforms (Figure 6d, e) and lithofacies distribution mark out a drastic change in the upper PS2s depositional system. The eastern margin PPAF completed a shallowing upward cycle due to the E-ward migration of the fluvial Po Plain complex; continental to shallow marine slope deposits replaced marine foredeep/foreland facies. In a broader view of the PPAF basin dynamics after Catuneanu (2004), PS1s and PS2s record the progressive switch from underfilled to filled phase of the foredeep during early-middle Pleistocene. The late Pleistocene last sequence, PS3s, deposited continental paleo-Po fluvial deposits since 0.87 Ma to recent in the eastern and central sector of the PPAF, covering the Ferrara Arc (PS3 map in Figure 11). When the sediment supply became higher than the accommodation space creation and the uplift of the tectonic barrier of the Ferrara Arc, the Po Prograding complex rapidly prograded eastwards and smoothed the basin morphology (Figure 8c). According with Catuneanu (2004), fluvial environment across the foreland testifies the overfilled phase and the end of the foreland cycle. Similarly to the present-day setting, marine turbidites deposited only in the Northern Adriatic Sea area, mainly fed by the Po Prograding Complex (profile C-C′ of PS3 stage, top Figure 11).

| Geodynamic implications from detrital pattern distributions
According to several authors (Marshak, 1988(Marshak, , 2005 and references therein), numerous interplaying factors lead to the growth of a non-cylindrical, asymmetrical mountain belt-foreland basin system. Among these, laterally variable rigidity of the downbending plate, stratigraphic variations along-strike in the foreland, structural indenters, differently oriented tectonic lineaments and foreland obstacles are the most commonly mentioned. Here, we discuss these possible causes applied on the study area. A cylindrical geometry of the Northern Apennines foredeep system is documented during Middle-Late Miocene. During Langhian-Tortonian age interval, deep-water siliciclastic turbidites of the Marnoso-Arenacea Formation filled an elongated foredeep basin NNW-SSE oriented (Figure 12a). The cylindrical and symmetric geometry of the basin is documented by turbidite key beds within the Marnoso-Arenacea Fm. They have been individually correlated for more than 200 km along the foredeep axis and their paleocurrent analysis exhibits along-strike opposite directions (Cibin et al., 2004;Di Giulio, 1999;Gandolfi, Paganelli, & Zuffa, 1983) ( Figure  12a, b). The coexistence of beds with opposite paleoflow directions into the same sedimentary body, is a clear proof of a general flatness and connection of the sea bottom. This foredeep basin shifted through time in response to the ENEward migration of the NA compressional front (Gandolfi et al., 1983;Ricci Lucchi, 1978, 1986Ricci Lucchi & Valmori, 1980;Zattin & Zuffa, 2004). Seismic profiles and stratigraphic analyses of the late Miocene succession (Rossi, Minervini, Ghielmi, & Rogledi, 2015), constrain the initial evolution of the Emilia Arc (Figure 1a) to the Messinian pre-evaporitic time. This means that the Emilia salient development broke this sedimentary continuum in the former Marnoso-Arenacea basin, leading to a non-cylindrical geometry. During the sequent Messinian post-evaporitic phase, Rossi et al. (2015) (Figure 10a). Mesozoic paleogeography of the study area shows NNW-SSE-oriented boundaries among carbonate platforms and adjacent basins, from W to the E they are: the Lombardian Basin, the Trento Plateau and the Belluno-Northern Adriatic basins (Carminati, 2009;Carminati et al., 2010;Fantoni & Franciosi, 2010;Livani et al., 2018;Masetti, Fantoni, Romano, Sartorio, & Trevisani, 2012;Turrini et al., 2016) ( Figure 2). Geodynamic and paleogeographic reconstructions describe the NNE-wards migration of the northernmost sector of Apennines belt progressively turning perpendicular to the Mesozoic inherited lineaments ( Figure  12c). The change in vergence was due to the anticlockwise rotation of the collisional system during the drifting of the Corsica-Sardinia block, with rotation pole located in the Ligurian Sea (Caricchi et al., 2014;Gattacceca et al., 2007). The emplacement of the Emilia and Ferrara arcs (during Messinian and Pleistocene, respectively) in the middle part of the PPAF, completely separated the foredeep into distinct depocentres (Figure 12c, d). It follows that the Apennine front rotation and foredeep migration encountered variable along-strike rheologic, structural and stratigraphic conditions, triggering the loss of cylindricity and symmetry of the basin floor.
To verify this hypothesis, a wider view of the Apennine orogen is needed. Towards the south, in the Central Italy the present-day Apennines thrust belt and foredeep system is ~NS-oriented (called Periadriatic basin, Artoni, 2013;Bigi et al., 2013) (Figure 12c). In fact, since Miocene time, this central sector experienced an eastward migration rather than counterclockwise rotation . Stratigraphic and tectonic evolution of the Periadriatic basin subsurface during Plio-Pleistocene time shows a general along-strike correlation of the sequences (Artoni, 2013). This is also in agreement with the detrital distribution pattern of Pliocene deposits inferred from outcrop analyses (Cellino Formation, grey box in Figure 12c). The Cellino Fm. is composed by well-developed and highly continuous thick-bedded turbidites (Casnedi, Moruzzi, & Mutti, 1978;Casnedi, 1983Casnedi, , 1991. Likely the middle Miocene deposits of Marnoso-Areneacea Fm., these features allow unravelling the cylindrical geometry of the central Apennine foredeep, where no severe rotation occurred.

| Main factors driving salient position
Within this general framework, the following question is: what is(are) the main factor(s) influencing the position of EAGE AMADORI et Al. salient and recesses in the Po Plain-Northern Adriatic thrust front? Mariotti and Doglioni (2000) suggest that structural salients usually develop in correspondence of the basins, whereas recesses are often located in correspondence of rigid, thick and relatively elevated sectors. In the Northern Apennines case, this correspondence between recess/salient position and the plateau-basin boundaries is only partial, as visible in Figure 2a. According to Livani et al. (2018), the mismatch is due to the widespread deposition of a highly plastic hemipelagic unit (Gallare Marls, Late  Gandolfi et al., 1983;Ricci Lucchi & Valmori, 1980). (c) Multiple depocentres and more complex detrital distribution pattern during non-cylindrical Early Pleistocene stage (box C, modified after Ghielmi et al., 2013). NA stepwise migration after . Grey area corresponds to the area where the Cellino Formation is exposed (Felletti et al., 2009) Eocene-Miocene) in the Meso-Cenozoic basins, whereas its deposition was limited or prevented on the top of structural high, like the Trento Plateau. This formation has been recognized as an effective décollement layer (Fantoni, Bersezio, & Forcella, 2004;Maesano et al., 2015;Massoli et al., 2006;Ravaglia et al., 2006), favouring the development of flat thrusts that partially override the basin-plateau boundary, causing the mentioned mismatch. Variability in the foreland slope is another driving factor for thrust spacing and nonself-parallel orogen evolution (Lickorish et al., 2002;Willett & Schlunegger, 2010). As dip increases, décollement depth, thrust wavelength and thrust-top basins dimensions increase (Doglioni, Merlini, & Cantarella, 1999;Mariotti & Doglioni, 2000) influencing also the seismicity distribution in the area (Vannoli, Burrato, & Valensise, 2015). The differential compaction rates, between the carbonate plateau and the pelagic basins (Figure 2b), may drive local subsidence response of the foreland and therefore, the southward dip of the regional monocline below the Apennines thrust front. Therefore, the low monocline dip (less than 10°) in correspondence to the plateau and its steeper attitude in correspondence of the basins can be considered as another controlling factor influencing the non-cylindricity of the PPAF. Other minor elements can emphasize the along-strike differences of the Northern Apennines thrust-foredeep system. For instance, inherited structures in the lower plate (Fantoni et al., 2004;Turrini et al., 2016) and volcanic-volcanoclastic deposits drilled in wells (Dalla et al., 1992;Di Giulio, Carrapa, Fantoni, Gorla, & Valdisturlo, 2001;Pieri & Groppi, 1981), located in correspondence to the Emilia Arcs, acted as obstacles during NA rotation.

| CONCLUSIONS
The Northern Apennines is a strongly arcuate orogen associated to a deep fragmented (non-cylindrical) pro-foreland basin. By virtue of its uncommon data coverage and preservation of tectonic phases in the sedimentary record, the Po Plain Northern Adriatic foreland can be considered a case history to investigate tectonic-sedimentation interaction in non-cylindrical orogen-basin systems. In fact, high-resolution studies on "when and why" the detrital distribution pattern (within potential hydrocarbon exploration targets) evolves during basin geometry changes, are particularly useful to both academia and industry. In detail, this work provides the three-dimensional reconstruction of the present-day PPAF subsurface architecture in depth domain. The model is composed of six basin-scale layers through Pliocene-Pleistocene time, bounding six sedimentary units. We have also described lithofacies and environmental variations of these sequences by seismic and well-log analyses, suggesting the PPAF experienced a complete foredeep cycle, from underfilled (deep marine) to overfilled (continental) stage, according to definitions by Catuneanu (2004). All these quantitative outcomes allowed to constrain the transition from cylindrical to non-cylindrical state of the PPAF during late Pliocene-Piacenzian time. To discuss and unravel the main factors controlling the non-cylindrical foredeep behaviour, we have compared our results with sedimentologic features from previous Apennines foredeep stages and geodynamic-paleogeographic reconstructions. We found that: 1. Change in Northern Apennines belt vergence with respect to pre-orogenic Mesozoic extensional faults in the lower plate. The foredeep system recorded a transition from a cylindrical to a non-cylindrical state, and the formation of a salient only when the thrust front became close to orthogonal to the inherited lineament. 2. Interrelated local geological features influenced the position of the salients, recesses and the related position of depocentres in the foredeep. They are: the shaley lithology on décollement level for thrust enucleation and foreland rigidity, expressed as monocline angle of dip under the orogenic wedge.
Moreover, the temporal resolution (~1 Ma) of the 3D stepwise evolution of the PPAF allows this work to be reconsidered for further quantitative estimations of the sediment flux through the last ~ 5 Ma. This would give new insights on climate versus tectonic implications for sediment generation as erosion rates of the surrounding mountain belts.

| DATA REPOSITORY
A: Key layers velocity maps. L1 and L2 correspond to PL1u and PS1u, respectively. B: Sensitivity analysis on the velocity model of PL1u. C: Sensitivity analysis on the velocity model of PS1u. D: Not interpreted seismic lines shown in Figure 6. D2: seismic units. E: A3 paper size of all isobath maps. F: A3 paper size of all isochore maps. G: wells dataset from the public ViDEPI database.