Late Holocene differential subsidence and relative sea level rise in the Tabasco Delta, Mexico

Coastal subsidence owing to compaction of Holocene strata and deeper-rooted components affects large delta plains such as the Tabasco delta in southern Mexico (Gulf coast). For this system, GNSS3-PPP ground-truthed LiDAR imagery of high-resolution dated beach-ridge series reveals considerable differential subsidence on either side of the present Usumacinta-Grijalva River mouth. Collected field-data allows for quantification of differential subsidence over several time windows and reconstruction of relative sea-level rise back to 5000 years ago. Observed differential subsidence of 1–1.5 m is regarded to be syn-sedimentary delta-subsurface compaction of buried strata in response to the accumulating overburden of the prograding beach-ridge complex.

1 Introduction

Coastal fronts and back-barrier plains of deltas tend to be areas with relatively high natural subsidence rates, owing to the stacked factors of Holocene auto-compaction, sediment and glaciohydro-isostatic loading and tectonic-basin subsidence (e.g. Törnqvist et al., 2008; Brain, 2016). In populated deltas, human activities can add to the natural subsidence rates (e.g. Jankowski et al., 2017; Koster et al., 2018; Minderhoud et al., 2018). In addition to these factors, ongoing climate change and associated sea-level rise, and reduced sediment input, substantiate threats of wetland loss and increasing flooding risk, call for monitoring and mitigating subsidence. Interpreting the human-induced components of modern-observed rates and patterns of subsidence requires reconstructions of natural background subsidence components and subsurface geological heterogeneity.

Figure 1 (a) Study area with location of GNSS points, and RSL index points. (b) Estimated subsidence along a beach-ridge formed around 3750 cal BP. Notice that the LiDAR surface elevation in the eastern part of the study area with mangrove swamp is ∼1 m overestimated.

In the Tabasco delta in southern Mexico, LiDAR imagery of the beach-ridge plain that stretches from the mouths of the Usumacinta and Grijalva rivers (Fig. 1a), reveals the system to have experienced considerable background subsidence. Central parts of the beach-barrier complex, stretching 3750–1800 cal BP (calibrated years before 1950), differ markedly in crest elevation west and east of the Grijalva branch, demonstrating differential subsidence between the two sides. On the eastern side, swales and all but the highest beach ridges are overgrown by peat (Fig. 1b). Knowing the age and topographic variation of each of the beach ridges (Nooren et al., 2017 and field data collected in 2019), allowed us to assess rates of west-east differential subsidence for older and younger parts of the complex and to derive data points for relative sea-level rise reconstruction: the double objectives of our research.

2 Methods

2.1 Ground-truthing LiDAR and borehole surface elevations

Revisiting key coring sites from 2011–2016 campaigns and while performing new corings in 2019, surface elevations were measured with a Trimble R8 GNSS3 receiver and groundtruthed with the 5×5 m LiDAR data from Mexican Authorities (INEGI) at 36 reference point locations (Fig. 1a). Surface elevation was measured by Precise Point Positioning (PPP) (Zumberge et al., 1997) using observations of at least 1 h. Post-processing used the Trimble CenterPoint RTX service (Leandro et al., 2011), yielding absolute position estimates in the ITRF2008 at cm-scale resolution. Ellipsoid heights were transformed to orthometric heights (NAVD88 datum) with the Geoidal Gravimétrico Mexicano GGM10 model (INEGI, 2010). The resulting elevations are reported as m + MSL. This is the same datum as used for the original LiDAR data product provided by INEGI.

2.2 Differential subsidence from beach ridge elevations

Differential subsidence estimates were derived in multiple ways. When we assume that past beach crests originally formed at the same heights above MSL as sub-recent ridges (for the eastern sector 0.5 to 1.5 m + MSL) we derive a minimum and a maximum subsidence magnitude for the ridges further in-land (Fig. 1b).

The alternative is to compare beach ridge crest elevations between the western and eastern sectors. Such an analysis was performed for the 3750–1800 cal BP part of the complex (Fig. 1b), building on an earlier established relation of beach ridge elevation and progradation rate (Nooren et al., 2017). We refined this relation with a simplified river delta model (Nienhuis et al., 2016) to estimate the variability in progradation rates, and hence initial beach-crest elevation (Fig. 1b), at various distances to the former river mouth (SP y SP in Fig. 1a). The decreasing variability in beach-ridge elevation with distance from the former river mouth, apparent from the LiDAR-based DEM (Figs. 1a, 10; Nooren et al., 2017), is the result of a decreasing variability in progradation rate, and can be explained by the diffusive behaviour of wave-dominated delta growth (Nienhuis et al., 2015).

We constrain the river delta model to fit the 3750 and 1800 cal BP beach ridges delineated in Fig. 1a. The model was run with a WaveWatch III hindcast wave climate (Chawla et al., 2013), and a mean fluvial sediment flux of 3.1×10⁶ m³ yr⁻¹, varying with ±20 % over a 500 year cycle. Progradation rates were translated to beach-crest elevations assuming an aeolian accretion rate of 5 m² yr⁻¹, and an elevation of the top of the swash deposits of 0.5 m + RSL (Relative Sea Level at time of deposition). This method could not be applied to younger parts, because the Grijalva River began co-affecting promontory development (from 1800 cal BP), as did wave-erosive cannibalisation of the promontory of the former river mouth (after 900 cal BP), breaking the relation.

2.3 Establishing sea-level index points

We use lagoonal organic beds, organic debris layers incorporated in sandy beach deposits, and pedogenic decalcification depth in beach sands of the last 6000 years (Fig. 1b), to establish a series of sea-level index points of approximately known age and vertical indicative meaning (cf. Shennan et al., 2015). The ages of the organic material were obtained by directly Accelerator Mass Spectrometry (AMS) ¹⁴C dating on selected macrofossils. Their consistency was validated based on vertical and lateral position in the beach ridge complex. The ages of the decalcified beach ridges follow the Nooren et al. (2017) reconstruction (Fig. 1a).

Figure 2 Elevations measured in the field by Precise Point Positioning (PPP) compared with the LiDAR DEM elevation (mean value of a ø 20 m circle around the measure point).

The lagoonal organic beds from peaty environments overly muddy back-barrier deposits. Considering compaction, dates from the base of these beds are “intercalated” sea-level index points in the systematics of Brain (2015). Mangrove peat beds have been formed within 0 to 0.3 m + RSL (after Khan et al., 2017), considering a tidal range of 0.6 m for the area. Freshwater peats have been formed between −1 to 0 m + RSL, where in the dry season groundwater levels can drop below sea level (see Sect. 3.1).

Concentrations of organic debris within beach deposits were encountered at multiple instances at variable depth (down to −8 m + MSL; Nooren et al., 2017). Their preservation suggest deposition at levels up to RSL, buried shortly after deposition, and preserved below groundwater table (above this level: dry season oxidation). The shallowest dated organic macro-remains per core are used as sea level index point.

Decalcification depth is an unconservative proxy in sea-level reconstruction, but a promising method when applied on easily leachable CaCO₃-rich beach ridge sands. The decalcification depth marks the maximum depth of the repeated dry-season groundwater table, resulting from dissolution and leaching by percolating water. Especially the fine-grained beach sands east (updrift) of the Grijalva outlet are highly calcareous (5 %–10 % CaCO₃; Nooren et al., 2017). The decalcification depth was field-estimated while logging handcores (using 10 % HCl solution), and a detailed trajectory sampled at 5 cm increments for further analyses using a Scheibler calcimeter. We assume a maximum decalcification depth of −1 m + RSL, concordant with the lowest groundwater levels occurring during the dry season.

Figure 3 High water marks on trees growing in a low swale (GNSS point 11, Fig. 1). Notice the GNSS receiver on top of a 2 m long pole.

3 Results and discussion

3.1 Surface elevation accuracy, back-barrier water level seasonality

There is very good agreement between the 36 PPP-measured elevations and the elevations spotted from the LiDAR-based DEM (Fig. 2), except for heavily vegetated inland sites into back-barrier flood basins, where LiDAR surface elevation is generally overestimated by 0.1–1.5 m, due to the dense vegetation cover present (Fig. 1). In central parts of the beach-ridge complex, water-level marks seen on trees growing in the relatively low swales (Fig. 3) indicate water levels of 1 m + MSL during the wet season. Our GNSS2-PPP survey took place during a remarkably dry season (May 2019), with groundwater tables dropped to −0.78 m + MSL at the site of Fig. 3, and to −0.4 m + MSL just 200 m inland of the modern beach. We are not aware of any substantial human groundwater extraction at these sites and in the wider study area, and explain the deep groundwater level drop to the dry-season intensity (limited rainfall) and high evapotranspiration.

Figure 4 RSL index points for the Tabasco delta, and sites likely effected by differential subsidence. Subsidence is indicated by yellow arrow, representing the difference between modelled initial beach-crest elevation and measured elevation by Precise Point Positioning (PPP). A tentative RSLR curve for the study area is presented (this study), and compared with the RSLR curve for wider Mexico, reconstructed from compaction-corrected index-points (after Khan et al., 2017).

3.2 Relative sea-level rise (RSLR)

The sea-level index points, for the area west of the Grijalva river mouth, produces a decelerating RSLR signal, with a strong decrease in RSLR rate around 5000 years ago, and approaching a semi-linear background rate of 0.15 mm yr⁻¹ since 2500 years ago (Fig. 4), in broad agreement with Gulf-Caribbean regional RSLR reviews (Khan et al., 2017).

3.3 Absolute and differential subsidence

Western part – The barrier morphology and SL index points from the west of the study area indicate that this part of the beach ridge complex did not experienced substantial subregional subsidence, as the RSLR signal is shared with wider Mexico region. In the west, the beach-ridge deposits likely overly an incompressible substrate (Grijalva alluvial fan, extending in the subsurface offshore (Padilla and Sanchez, 2007). This said, fault-related differential subsidence up to 1.5 m, displacing beach ridges up to those from the 16th century AD (Fig. 1), did affect far western areas.

Eastern part – The barrier system in the east of the study area is affected by an additional subsidence term that could be quantified (Fig. 4). This term amounts for a maximum of 1.0–1.5 m lowering of sea-level index-points obtained from the beach ridge complex of the last 4000 years (Fig. 4), a number very similar to the west-east elevation differences of bundled beach ridge crests (1.0–1.6 m), with maximum values in the encircled area NE of Frontera (Fig. 1). Importantly, there is no increase of this value (length of yellow arrows in Fig. 4) with increasing age. Together with spatial patterns, this indicate that the cause of the subsidence is to be sought in the local delta substrate. We presume this sector to overly the Usumacinta low-stand palaeovalley with a muddy transgressive fill, over which foreshore and beach-ridge sands (some 8–12 m thick) began prograding in the last 6000 years (when RSLR inflected). This would cause syn-sedimentary compaction (accommodating extra beach-ridge thickness) of which the SL index points and beach crest displacements register the later phase.

We recommend to further test this attribution by establishing thickness of the beach ridge complex and its immediate substrate in next field campaigns. Cross-comparison with other Holocene delta-systems (Rhine system: De Groot and De Gans, 1996; Hijma et al., 2009; Mississippi system: Törnqvist et al., 2008) suggests that 4–6 m of excess thickness can be expected – of which our 1–1.5 m would be the residual part (delayed creep) and the difference the semi-instantaneous part. Such testing is essential for subsidence-impact evaluations of future activities (superficial drainage activities for agriculture, measures for lagoon-rim nature conservation, subsidence owing to hydrocarbon extractions at depth), and therefore highly needed for an adequate management of the Tabasco delta.

4 Conclusions

This study presents an integrated analysis of relative sea level rise and differential subsidence based on vertically quantified differences in coastal topography between western and eastern parts of the study area, a well-resolved beach-ridge chronology, and the collection of sea-level index points.

The results allow to separate a background RSLR mainly due to regional glacio-hydroisostasy, from local differential subsidence, interpreted to signal the presence of compressible strata immediately below the beach-ridge sands east of the Grijalva River mouth. This natural subsurface compaction has substantially contributed to the low elevation of the eastern part of the delta.