Factors controlling natural subsidence in the Po Plain

Understanding the causes and mechanisms of land subsidence is crucial, especially in densely populated coastal plains. In this work, we calculated subsidence rates (SR) in the Po coastal plain, averaged over the last 5.6 and 120 kyr, providing information about land movements on intermediate (103–105 years) time scales. The calculation of SR relied upon core-based correlation of two lagoon horizons over tens of km. Subsidence in the last 120 kyr appears to be controlled mainly by the location of buried tectonic structures, which in turn controlled sedimentation rates and location of highly compressible depositional facies. Numerical modelling shows that subsidence in the last 5.6 kyr is mainly due to compaction of the Late Pleistocene and Holocene deposits (uppermost 30 m).


Introduction
The effects of land subsidence could be devastating on heavily settled, low-gradient, coastal plains, in a scenario of 0.9 m of relative sea-level rise by 2100 (Rohling et al., 2013). Increasingly high costs are expected to protect coastal cities and touristic hotspots and to keep drained reclaimed lands (Syvitski et al., 2009;Erkens et al., 2016).
In the Po coastal plain, extending more than 100 km south of the Venice lagoon, decadal to centennial land movements in the order of several cm yr −1 (Teatini et al., 2011;Cenni et al., 2013), were strongly driven by human activities (mostly water and gas withdrawal, Teatini et al., 2006). To separate anthropogenic from natural subsidence, understanding the role of subsidence through time before human intervention is required. Long-term (10 6 years) natural subsidence has been calculated based on deep seismic analysis (Carminati and di Donato, 1999). On the contrary, spatial distribution of SR on intermediate (10 3 -10 5 years) time scales is poorly known.
In this work, we measured SR in the last 120 kyr along a 40 km-long transect in the Po coastal plain. We used as refer-ence marker beds two lagoon horizons identified in sediment cores, from the Last and Present Interglacials. To discuss the relations between subsidence and structural setting, we reconstructed stratal architecture down to 3 km through seismic analysis. We focused on the stratigraphic significance of marker beds and on their use for SR calculation. The contribution of sediment compaction to land subsidence was assessed through geotechnical analysis and decompaction modelling.

Geological setting
The Po Plain is a ∼ 40 000 km 2 -wide alluvial plain bounded by the Alps to the north and by the Apennines to the south. Frontal thrusts of both chains are sealed beneath the Po Plain (Amadori et al., 2019). The Apennine thrusts are north-verging and exhibit arcuate shapes. Thrusts SE of Ferrara ( Fig. 1) were active since the Early Pliocene. The Po Basin fill shows a shallowing-upward trend from deepmarine turbidites to coastal and continental units (Ghielmi et al., 2013). The rhythmical alternation of Middle-Late Pleis- tocene coastal and alluvial deposits reflects glacio-eustatic fluctuations at the Milankovitch scale (Amorosi et al., 2004). Two transgressive-regressive coastal wedges in the uppermost 130 m were assigned to the Last (MIS 5e) and to the Present (MIS 1) Interglacials based on pollen data, 14 C dates and electron-spin-resonance age determinations (Ferranti et al., 2006). Landwards, the maximum marine ingression is marked by a thin brackish-lagoon horizon, sandwiched between inner estuary (below) and upper delta plain freshwater deposits. Laterally extensive brackish water bodies settled in the Po coastal Plain around 7-5 cal kyr BP, after sea-level stabilization (Vacchi et al., 2016), in areas sheltered by fluvial sedimentary input .

Stratigraphy
Calculation of Late Pleistocene and Holocene SR relies on the elevation (meters above modern sea level -m a.m.s.l.) of two stratigraphic markers. The lower one is a lagoon horizon marking the top of the MIS 5e coastal wedge. The upper one is a thin lagoon horizon, dated to 5.6 ± 0.5 cal kyr BP, encountered at depths < 20 m. Chronological constraints for core correlation are provided by 76 14 C ages (see Supplement).
In order to assess the impact of buried faults and folds on lateral changes in SR, correlations were carried out perpen-dicular to these structures. A similarly oriented seismic profile was interpreted and time-depth converted through calibration with exploration-well logs (depth of 300-5000 m). Seismic interpretation relied on the identification and tracking of the main discontinuities (i.e. reverse faults) and unconformities.

Subsidence rates calculation
The calculation of SR for a selected time interval ( t) between the time of deposition (t d ) and the Present (t 0 = 0) was based on Eq. (1).
where Z td is the elevation of the stratigraphic marker at the time of deposition and Z t0 is the present elevation. An error of 0.15 m associated to Z t0 due to core stretchening/shortening was included in the calculation (Hijma et al., 2015). Sedimentological data  indicate that the 5.6 kyr BP lagoon horizon accumulated in a lower delta plain environment (∼ 0 m above coeval sea level). Estimation of sea level at 5.6 cal kyr BP was based on the prediction curve of Vacchi et al. (2016). As sediment deposition presumably took place in a water body < 1 m deep, an error of ±0.5 m was taken into account in SR calculation. The MIS 5e lagoon horizon likely accumulated in water depths < 2 m (associated error of ±1 m; Amorosi et al., 2004). An additional error of ±1 m was added to compensate possible inaccuracy in correlation due to low chronologic resolution. The lagoon horizon was generically assigned to the MIS 5e highstand (120-116 cal kyr BP, Rovere et al., 2016), 7 ± 2 m a.m.s.l.
Subsidence rates between 120-116 and 5.6-0 cal kyr BP were compared with SR averaged over the last 1.5 Myr. The 1.5 Myr unconformity (Amadori et al., 2019) is marked in exploration-well logs by the first occurrence of Hyalinea Baltica (Carminati and di Donato, 1999).

Sediment compaction
In order to assess the contribution of sediment compaction to natural subsidence, a finite-element 1-D decompaction model (Gambolati et al., 1998;Zoccarato and Teatini, 2017) was applied to core B4 (see Supplement). The mechanical characterization of the main lithofacies associations was based on oedometer tests, bulk-density and loss on ignition tests (van Asselen et al., 2009).
The 1.5 Myr unconformity depicts a major fold, with a wavelength of ∼ 30 km resulting from faults propagation and imbrication. Lower-rank folds (wavelength of ∼ 3-5 km) are associated with single-fault propagation. The thickness of the post-tectonic unit ranges between 300 m at the culmination of the main anticline and 2000 in the backlimb and forelimb depocentres. Post-1.5 Myr turbidites accumulated in synclinal areas (Ghielmi et al., 2013), whereas truncation of Pliocene strata is observed above the anticline. Deformation of Pleistocene strata was observed close to the anticline culmination, suggesting local post-1.5 Myr tectonic activity.

Shallow stratigraphy
The stratigraphy of the uppermost 140 m is depicted in the correlation panel of Fig. 3 (see Supplement for facies description). The MIS 5e sediment wedge, at 60-100 m depth, consists of coastal sands, replaced by lagoon deposits at inland locations (core 204S4). A thin (0.5-3 m) lagoon horizon drapes the coastal sands and is overlain by a shallowing upward succession of coastal plain and alluvial deposits, 35-65 m thick, deposited between ∼ 110 and 10 kyr BP. Delta plain estuarine and swamp deposits, pinch out toward the anticline culmination. Fluvial-channel deposits are abundant around MIS 4 and MIS 2, whereas closely spaced paleosols mark MIS 3 floodplain deposits. A prominent paleosol marks the Late Pleistocene-Holocene boundary and is overlain by estuarine deposits . A thin lagoon hori-zon, mostly sandwiched between freshwater deposits and dated to ∼ 5.6 cal kyr BP, is clearly recognizable along the whole cross-section.

Factors controlling subsidence rates
Subsidence rates in the Po coastal plain appear to be strongly influenced by the buried tectonic structures. Particularly, SR averaged over the last 5.6 (SR 5.6 ), 120 (SR 120 ), and 1500 (SR 1500 ) kyr decrease towards the anticline culmination (Fig. 3). In syncline areas, subsidence may have been enhanced by high sedimentation rates and by the preferential accumulation of highly compressible prodelta, delta plain and swamp muds (Fig. 3).
Subsidence rates in the last 5.6 kyr are in the range of 0.8-2.0 mm yr −1 , about twice the values measured for the last 120 kyr (0.4-0.8 mm yr −1 ). Superimposed to the structurally controlled trend, Holocene stratigraphic architecture influenced lateral changes in SR 5.6 . For example, SR 5.6 are higher in EM1, where the 5.6 kyr BP horizon overlies soft estuarine muds, and lower in EM8, where the same marker bed lies onto the less compressible, 7 m-thick coastal sand body. The 5.6 kyr lagoon horizon, encountered at ∼ −9 m a.m.s.l. in core B4, was placed at −3.5 m after decompaction, ∼ 1 m below sea-level predictions (Vacchi et al., 2016). Therefore, subsidence at this location is almost entirely accommodated by the compaction of the uppermost 30 m-thick sediment package. Particularly, the Late Pleistocene succession did not act as an "incompressible substratum" but compacted of ∼ 23 % (see Supplement), likely due to soft swamp peaty clays.

Conclusion
Subsidence rates for the last 5.6 and 120 kyr were assessed in the Po coastal plain within a chronologically constrained stratigraphic framework. Two laterally extensive lagoon horizons were used as reference marker beds.
Subsidence rates in the last 120 kyr appear to be mainly controlled by the location of buried Apennines thrust-related folds, which in turn influenced sedimentation rates and distribution of highly compressible depositional facies. Compaction of Late Pleistocene and Holocene sediments strongly impacted SR in the last 5.6 kyr. Data availability. Exploration-well logs used in this work are available at https://www.videpi.com/videpi/pozzi/pozzi.asp (Ministry for Economic Development DGRME -Italian Geological Society -Assomineraria, 2019). Oedometer tests were downloaded from the database of the Geological Seismic and Soil Survey of Regione Emilia Romagna (https: //ambiente.regione.emilia-romagna.it/en/geologia/cartografia/ webgis-banchedati/webgis-e-banche-dati?set_language=en, Geological, seismic and soil survey, Regione Emilia Romagna, 2019).
Author contributions. This work relies on a stratigraphic study carried out by LB, BCa and BCo and coordinated by AA. Laboratory analyses were carried out by BCo under the supervision of ES. Decompaction modelling was carried out by CZ and PT. LB prepared manuscript and figures, with the contribution of all coauthors.
Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "TISOLS: the Tenth International Symposium On Land Subsidence -living with subsidence". It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17-21 May 2021.
Financial support. This research has been supported by the RFO grant to Alessandro Amorosi. This paper is also a contribution to the International Geoscience Programme (IGCP) Project 663 "Impact, Mechanism, Monitoring of Land Subsidence in Coastal Cities".