Articles | Volume 382
Pre-conference publication
22 Apr 2020
Pre-conference publication |  | 22 Apr 2020

Factors controlling natural subsidence in the Po Plain

Luigi Bruno, Bruno Campo, Bianca Costagli, Esther Stouthamer, Pietro Teatini, Claudia Zoccarato, and Alessandro Amorosi

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

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−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 (106 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 (103–105 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 km2-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, 14C 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 1Study 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 14C 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 (td) and the Present (t0=0) was based on Eq. (1).

(1) SR ( Δ t ) = Z t d - Z t 0 Δ t ,

where Ztd is the elevation of the stratigraphic marker at the time of deposition and Zt0 is the present elevation. An error of 0.15 m associated to Zt0 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 2Interpreted 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 (SR5.6), 120 (SR120), and 1500 (SR1500) 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 SR5.6. For example, SR5.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.

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.

Data availability

Exploration-well logs used in this work are available at (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 (, Geological, seismic and soil survey, Regione Emilia Romagna, 2019).


The supplement related to this article is available online at:

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 co-authors.

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


Amadori, C., Toscani, G., Maesano, F. E., D'Ambrogi, C., Ghielmi, M., and Fantoni, R.: From cylindrical to non-cylindrical foreland basin: Pliocene-Pleistocene evolution of the Po Plain-Northern Adriatic basin (Italy), Basin Res., 31, 991–1015,, 2019. 

Amorosi, A., Colalongo, M. L., Fiorini, F., Fusco, F., Pasini, G., Vaiani, S. C., and Sarti, G.: Palaeogeographic and palaeoclimatic evolution of the Po Plain from 150-ky core records, Global Planet. Change, 40, 55–78,, 2004. 

Amorosi, A., Bruno, L., Campo, B., Morelli, A., Rossi, V., Scarponi, D., Bohacs, K. M., and Drexler, T. M.: Global sea-level control on local parasequence architecture from the Holocene record of the Po Plain, Italy, Mar. Petrol. Geol., 87, 99–111,, 2017. 

Bruno, L., Bohacs, K. M., Campo, B., Drexler, T. M., Rossi, V., Sammartino, I., and Amorosi, A.: Early Holocene transgressive paleogeography in the Po coastal plain (Northern Italy), Sedimentology, 64, 1792–1816,, 2017. 

Carminati, E. and di Donato, G.: Separating natural and anthropogenic vertical movements in fast subsiding areas: The Po Plain (N.Italy) case, Geophys. Res. Lett., 26, 2291–2294,, 1999. 

Cenni, N., Viti, M., Baldi, P., Mantovani, E., Bacchetti, M., and Vannucchi, A.: Present vertical movements in Central and Northern Italy from GPS data: Possible role of natural and anthropogenic causes, J. Geodyn., 71, 74–85,, 2013. 

Erkens, G., Van der Meulen, M., and Middelkoop, H.: Double trouble: subsidence and CO2 respiration due to 1000 years of cultivation of the Dutch coastal peatlands, Hydrogeol. J., 24, 551–568,, 2016. 

Ferranti, L., Antonioli, F., Mauz, B., Amorosi, A., Dai Pra, G., Mastronuzzi, G., Monaco, C., Orru, P., Pappalardo, M., Radtke, U., Renda, P., Romano, P., Sanso, P., and Verrubbi, V.: Markers of the last interglacial sea-level high stand along the coast of Italy: Tectonic implications, Quaternary Int., 145–146, 30–54,, 2006. 

Gambolati, G., Giunta, G., and Teatini, P.: Numerical modeling of natural land subsidence over sedimentary basins undergoing large compaction, in: CENAS – Coastiline Evolution of the Upper Adriatic Sea due to Sea Level Rise and Natural and Anthropogenic Land Subsidence, edited by: Gambolati, G., Springer Netherlands, Heidelberg, the Netherlands, 77–102, 1998. 

Geological, seismic and soil survey, Regione Emilia Romagna: Webgis and interactive cartography, available at:, last access: 2 October 2019. 

Ghielmi, M., Minervini, M., Nini, C., Rogedi, S., and Rossi, M.: Late Miocene–Middle Pleistocene sequences in the Po Plain–Northern Adriatic Sea (Italy): The stratigraphic record of modification phases affecting a complex foreland basin, Mar. Petrol. Geol., 42, 50–81,, 2013. 

Hijma, M., Engelhart, S. E., Tornqvist, T. E., Horton, B. P., Hu, P., and Hill, D.: A protocol for a Geological Sea-level Database, in: Handbook of Sea Level Research, edited by: Shennan, I., Long, A., and Horton, B. P., Wiley, 536–553,, 2015. 

Ministry for Economic Development DGRME – Italian Geological Society – Assomineraria: ViDEPI Project Visibility of petroleum exploration data in Italy – Final wells logsl, available at:, last access: 2 October 2019. 

Rohling, E. J., Haigh, I. D., Foster, G. L., Roberts, A. P., and Grant, K. M.: A geological perspective on potential future sea-level rise, Sci. Rep.-UK, 3, 3461,, 2013. 

Rovere, A., Raymo, M. E., Vacchi, M., Lorscheid, T., Stocchi, P., Gómez-Pujol, L., Harris, D. L., Casella, E., O'Leary, M. J., and Hearty, P. J.: The analysis of Last Interglacial (MIS 5e) relative sea-level indicators: Reconstructing sea-level in a warmer world, Earth-Sci. Rev., 159, 404–427,, 2016. 

Syvitski, J. P. M., Kettner, A. J., Overeem, I., Hutton, E. W. H., Hannon, M. T., Brakenridge, G. R., Day, J., Vörösmarty, C., Saito, Y., Giosan, L., and Nicholls, R. J.: Sinking deltas due to human activities, Nat. Geosci., 2, 681–686,, 2009. 

Teatini, P., Gambolati, G., Ferronato, M., and Gonella, M.: Groundwater pumping and land subsidence in the Emilia-Romagna coastland, Italy: Modeling the past occurrence and the future trend, Water Resour. Res., 42,, 2006. 

Teatini, P., Tosi, L., and Strozzi, T.: Quantitative evidence that compaction of Holocene sediments drives the present land subsidence of the Po Delta, Italy, J. Geophys. Res., 116, B08407,, 2011. 

Vacchi, M., Marriner, N., Morhange, C., Spada, G., Fontana, A., and Rovere, A.: Multiproxy assessment of Holocene relative sealevel changes in the western Mediterranenan: Sea-level variability and improvements in the definition of the isostatic signal, Earth-Sci. Rev., 155, 172–197,, 2016. 

van Asselen, S., Stouthamer, E., and van Asch, T. W. J.: Effects of peat compaction on delta evolution: a review on processes, responses, measuring and modeling, Earth-Sci. Rev., 92, 35–51,, 2009.  

Zoccarato, C. and Teatini, P.: Numerical simulations of Holocene salt-marsh dynamics under the hypothesis of large soil deformations, Adv. Water Res., 110, 107–119,, 2017. 

Short summary
The effects of land subsidence could be devastating on heavily settled coastal plains. In a scenario of sea-level rise, high costs are expected to protect coastal cities and touristic hotspots and to keep drained reclaimed lands. In this work, we calculated subsidence rates (SR) in the Po coastal plain, over the last 5.6 and 120 thousand years, providing information about land movements before human intervention became the main driver of subsidence, through water and gas withdrawal.