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 (10³–10⁵ 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).

1 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⁻¹ (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⁶ 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³–10⁵ 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 reference 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.

2 Geological setting

The Po Plain is a ∼ 40 000 km²-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 deep-marine turbidites to coastal and continental units (Ghielmi et al., 2013). The rhythmical alternation of Middle-Late Pleistocene 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, ¹⁴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 (Amorosi et al., 2017).

Figure 1 Study area with location of cores, exploration wells and cross sections of Figs. 2 and 3. The projection of north-verging buried thrusts is depicted.

3 Methods

3.1 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 ¹⁴C ages (see Supplement).

In order to assess the impact of buried faults and folds on lateral changes in SR, correlations were carried out perpendicular 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.

3.2 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) was based on Eq. (1).

$$\begin{array}{}\text{(1)}& \mathrm{SR}\left(\mathrm{\Delta }t\right)={\displaystyle \frac{Z_{t\mathrm{d}}-Z_{t\mathrm{0}}}{\mathrm{\Delta }t}},\end{array}$$

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 (Amorosi et al., 2017) 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).

3.3 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).

Figure 2 Interpreted seismic profile (see Fig. 1 for location) perpendicular to the main buried tectonic structures.

4 Results and discussion

4.1 Deep stratal architecture

The subsurface of the study area, down to 3 km depth, consists of imbricated thrusts, locally associated with backthrusts (Fig. 2). Two major unconformities, dated to 5.5 and 1.5 Myr (Ghielmi et al., 2013), respectively, allow identification of three tectono-stratigraphic units (Fig. 2): (i) an intensely deformed and faulted pre-tectonic unit; (ii) a syn-tectonic wedge, made up of marine turbidites (Ghielmi et al., 2013); and (iii) a nearly undeformed, post-tectonic unit, sealing the main structures.

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.

4.2 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 (Bruno et al., 2017). A thin lagoon horizon, mostly sandwiched between freshwater deposits and dated to ∼5.6 cal kyr BP, is clearly recognizable along the whole cross-section.

Figure 3 (a) Subsidence rates averaged over the last 5.6, 120 and 1500 kyr; (b) core correlation of Late Pleistocene-Holocene deposits (see Fig. 1 for core location); (c) elevation of the 1.5 Myr BP unconformity.

4.3 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₁₂₀), and 1500 (SR₁₅₀₀) 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⁻¹, about twice the values measured for the last 120 kyr (0.4–0.8 mm yr⁻¹). 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 $\sim -\mathrm{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.

4.4 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.