Holocene Sea Level Recorded by Beach Rocks at Ionian Coasts of Apulia (Italy)

Abstract

Beach rocks are located along many coasts of the Mediterranean basin. The early diagenesis environment and the mean sea level along the shoreline make these landforms useful in the reconstruction of relative sea-level changes and, in particular, as SLIPs (sea-level index points). The beach rocks surveyed along the Ionian coast of Apulia were found to be well preserved at three specific depth ranges: 6–9 m, 3–4 m, and from the foreshore to about 1.20 m. Morpho-bathymetric and dive surveys were performed to assess both the geometries and the extension of the submerged beach rocks. Samples were collected at these different depths in the localities of Lido Torretta, Campomarino di Maruggio, San Pietro in Bevagna, and Porto Cesareo. Bivalve shells were identified and isolated from the beach rock samples collected at a depth of 7 m; AMS dating provided a calibrated age of about 7.8 ka BP. Their morphology and petrological features, along with the time constraints, enabled us to (i) reconstruct the local sea-level curve during the Holocene, (ii) corroborate acquired knowledge of the relative sea-level history, and (iii) identify possible local vertical land movement (VLM).

1. Introduction

Sea-level index points (SLIPs) are represented by geomorphological, biological, and archaeological markers of past sea-level positions. These markers are used to reconstruct the sea level from local to regional scales with a given accuracy [1,2,3]. One of these markers is provided by beach rocks (BRs), with an approximation conditioned by local physical and chemical parameters (i.e., beach slope and granulometry; tide; wave climate; water table elevation; pH, CO₂, CO₃⁻, and CaCO₃ concentrations; SST; air temperature; and pressure). BRs are sedimentary landforms deriving from the early cementation of the foreshore sands in a short time span, recording the position of the corresponding shoreline [4,5,6]. They are represented by cemented sandstone and are generally composed of many decametric levels, layered upon sequences of slabs [7,8,9]. They appear as collapsed “domino tiles”, following the beach slope. BRs are composed of medium-fine sands cemented by aragonite and high-magnesium calcite within the intertidal zone; the process of cementation is due to the permeation of fresh groundwater and seawater into the sediments following tidal phases [7,10]. The presence of BRs indicates the horizontal and vertical shoreline stability. Their current position above or below the mean sea level suggests past relative sea-level stands and related shoreline position. Furthermore, some critical questions on the study of BRs have been posed, mainly because their upper limit of cementation is not well defined [5,6,11,12]. For this reason, the use of a BR as a relative sea-level marker is affected by vertical errors, different from site to site, due mainly to the tidal range and wave climate (i.e., [12,13]). Most BRs occur on microtidal coasts with an average thickness of 2 m (e.g., [12]). Thus, the vertical errors typically fall between 2 m and 0 m. Every criticism in the identification of past sea levels using ancient BRs may be nulled by comparing them with the lithification environment of present BRs from the same area, thereby defining the modern analog [9].

BRs have been amply studied worldwide, and their geometry has allowed for making several considerations about past relative or eustatic sea-level changes (i.e., [7,13,14,15,16,17,18,19]). BR occurrence seems to be prevalent in the Mediterranean and Caribbean Seas, the tropical and subtropical Atlantic coasts, and the atolls of the Pacific and Indian Oceans (i.e., [4,14,20,21]). Their general features make them useful as robust SLIPs in those coastal areas of the Mediterranean basin (Figure 1) that are characterized by tidal range less than 1 m high as confirmed by geodetic data and tide gauge observations (i.e., [22,23] and references therein). This approximation can be reduced up to considering a data accuracy of +/−0.10 m owing to the performed survey and the comparison with the modern analog [9].

In particular, the Mediterranean climate and the microtidal regime allowed BR diagenesis in different times during the generalized Holocene transgression (see Supplementary Material S1 for references). BRs are widespread along the Ionian coast of Apulia, in southern Italy. They outcrop discontinuously from Taranto (NW) to Capo Santa Maria di Leuca (SE) at (A) Lido Torretta–Torre Sgarrata; (B) from Torre Ovo, Campomarino di Maruggio, and San Pietro in Bevagna; (C) from Porto Cesareo and Torre Squillace; (D) at Lido Conchiglie and Gallipoli; and (E) from Isola Pazze and Ugento for a total of about 110 km (Figure 2).

In this study, all BRs were geomorphologically surveyed directly in-field. In particular,

i.

those in Campomarino di Maruggio were surveyed also by a multibeam echosounder (MBES);

ii.

the latter, as well as that of Porto Cesareo and Isola Pazze, were surveyed by scuba divers.

Moreover, (i) BR samples from Campo Marino and Porto Cesareo underwent petrographic analyses under an optical microscope, and (ii) samples of bivalve shells in BRs from San Pietro in Bevagna were collected and subjected to AMS analysis.

All the available data were analyzed in order to (i) assess the morphological features of BRs and their position compared with the present mean sea level, (ii) identify the past sea level and shoreline positions, and (iii) estimate the potential vertical land movements (VLMs) of the Ionian coastal area of the Apulian foreland.

2. Geomorphological Settings

BRs are widespread all along the mobile coastal system of the Ionian coast of Apulia. They are present in small pocket beaches as well as along mainland beaches [24,25,26]. BRs have been detected in correspondence to the present sea level and at different depths in the areas extending from Taranto to Capo Santa Maria di Leuca (Figure 2).

A detailed geomorphological field survey focused on the area extending between Campomarino di Maruggio and San Pietro in Bevagna as it is the longest stretch, about 10 km, along which BRs crop out continuously at different depths. This area is characterized by a mobile coastal system with transgressive dunes extending for about 15 km, having an elevation of up to 12 m, ascribed to the “Parco delle Dune di Campomarino”. Dune belts were chronologically attributed to the postglacial knickpoint of the sea-level transgression, which occurred about 7000 years BP ago [25,26,27,28,29]. In this area, beaches are characterized by medium-sized and well-sorted sands deriving from the Cretaceous and Pleistocene local basement and from modern biocoenosis [30,31].

The Salento Peninsula experienced a period of general subsidence during the Lower Pleistocene, interrupted by a rapid uplift during the Middle Pleistocene. The uplift ended at MIS 9.3 (about 330 ka) and was replaced by a phase of substantial stability, which lasted until the Late MIS 5.5 [32,33,34,35,36] in response to the recent doming of the region [37]. During the Middle Holocene, the Ionian coast showed significant beach accretion and dune belt growth [25,26]. Using archaeological evidence, recent papers [38,39,40] have stated that the latter, compared with the predicted sea level, highlights the vertical tectonic stability or the low-rate uplift of the Ionian and Adriatic sides of the Murge area. Some of them [38,39] seem to indicate the occurrence of local differential VLMs that affected some limited areas of the Ionian coast of the Salento Peninsula during the last four millennia.

Data from the Taranto Tide Gauge Station, part of the official Italian Tide Gauge Network (Rete Mareografica Nazionale, https://www.mareografico.it/, accessed on 1 June 2021 Jume), indicate a maximum semidiurnal tidal range of about 0.4 mi [41].

3. Material, Methods, and Data

With the aim of reconstructing the morphology of the BRs and their relation to the present sea level, topographic, dive, and bathymetric surveys were performed.

The first was performed by means of a GNSS receiver with RTK methodologies. The GNSS position data were corrected in real time using an RTK correction service supplied by the Apulia Region GPS SIT RTK service. It consists of several GNSS station networks across Apulia (managed by the official regional administration; http://gps.sit.puglia.it/, accessed on 1 June 2021) that supply RTK correction in real time in an RTCM format through the NTRIP protocol via the internet. The last GNSS topographic position used in the survey achieved a centimetric RTK accuracy. Ellipsoid-measured elevations, which were reduced to the Official Italian High Precision Levelling Vertical Reference Surface (established by the Military Geographic Institute (Istituto Geografico Militare–IGM)) by applying an ellipsoid–geoid separation value, calculated for the study area, allowed us to achieve a final accuracy of +/−0.02 m

The underwater surveys were carried out by two divers each equipped with an air breathing apparatus (ARA), an underwater electric scooter, and a Scubapro Aladin 2G depth gauge; this water pressure apparatus has decimetric accuracy. The dives were made with an RHIB (rigid hull inflatable boat) as a support vessel equipped with a GNSS positioning system. The divers operated according to the semicircular relief technique, typical of the Italian Navy [42].

The bathymetric surveys were performed using the interferometric multibeam echosounder system (MBES) from the Italian Navy Hydrographic Institute (Istituto Idrografico della Marina—IIM) by means of a hydroboat of the Italian Navy ship ITS Galatea in August 2019 and an MBES R2Sonic equipment on board the vessel “DC7” owned by the Consorzio Nazionale Interuniversitario per le Scienze del Mare (CoNISMa) in September 2020, managed by DISTEGEO and ENSU (Figure 3).

As standard practice, our surveyor team adopted IHO (International Hydrographic Organization) procedures to execute all measurements performed for this study. In order to achieve the best determination of seabed bathymetry and morphology, an IHO S-44 Special Order standard was adopted and applied for the entire dataset, as described in the Manual of Hydrography issued by the Italian Navy Hydrographic Institute for IHO in 2016.

To obtain a comparable dataset, the survey vessel “DC7” was accurately measured in order to determine the relative position for each sensor installed on board (i.e., MRU, GNSS antennas, a multibeam echosounder, etc.) to minimize nonstochastic errors in seabed measurements. Offset sensors were carried out by means of total station and laser scanner measurement sessions under vessel dry-dock conditions. Once the center of gravity (COG) of the vessel had been determined, its value was inserted into the Data Acquisition System (professional hydrographic navigation and acquisition software installed on a PC with serial and Ethernet connections) for online correction of instantaneous measurements (GNSS position, sounding measurements, sound velocity correction, heave, pitch and roll values measured by MRU). During the data acquisition phase, a sound velocity probe was installed close to the ES transducer and connected to the Data Acquisition System to correct the water SV in real time. As good practice, at least once a day (or more, if seawater conditions abruptly changed), a sound velocity profile was performed by means of a CTD probe along the water column in the study area.

The GNSS position and attitude data were calculated by an INS (GNSS position + inertial platform + FO gyrocompass integrated) system connected to the Data Acquisition System and synchronized with the ES system through PPS sync.

A multibeam ES was integrated and connected to the Data Acquisition System, which performed water depth measurements at a 455 kHz frequency and up to a 40 Hz ping rate in order to obtain the maximum data density on the seafloor for the best determination of BR morphology (more than 100 soundings\sqm insonification). Such data density allows the production of the HR DTM of the seabed (10 cm cells).

The entire dataset was processed by Professional Hydrographic Data Processing Software. It allowed for filtering data outliers through manual, logical, and automatic filters and calculating the total propagated uncertainty (TPU, as defined in IHO S-44 publication) in order to remove non-SO data from the dataset as well as from the tide data application. Moreover, in order to remove low-frequency MRU data drifting, real-time heave data were replaced by delayed, more reliable heave data in the entire dataset. Finally, high-resolution DTM was produced to show detailed seabed features detected and to accurately evaluate the current depth of the BRs. A list of the instruments used can be found in Supplementary Material S2.

Such methodologies allowed the production of a high-resolution morpho-bathymetry with decimetric resolution.

3.1. Morphological Data

BRs situated in correspondence to the shoreline outcrop between about 0.25 m above the high tide level—perhaps as an effect of the mean run-up—and about 1.2 m below the mean sea level, showing a total thickness of about 1.65 m. They dip seaward with slopes ranging from 5° to 10°. They look like a mosaic of disjointed blocks/tiles with sharp edges (Figure 4). Fractures and generalized patterns derive from the recent collapse and fracturing of the BRs due to the differential adjustments of the underlying loose sands. BRs grade upwards into emerged beach/dune sediments that represent the foreshore/backshore area. Dive and bathymetric surveys jointly revealed that the BRs outcrop for about 9 km, while their vertical extensions are about 1.5 m. Moreover, they revealed the continuous occurrence of BRs at two depths, nearly parallel to the present shorelines (Figure 5 and Figure 6):

-

At a depth range of 3–5 (+/−0.10) m, there is a BR belt whose boulders are disjointed and well defined (Figure 6).

-

Another BR belt is situated at a depth range of 6–9 (+/−0.1) m; in this case, boulders are disjointed and well rounded, and their geometry tends to reproduce a “manmade artifact” (Figure 5D).

Submerged BRs are represented by the alignment of thousands of disjointed boulders, up to 7.0 × 5.0 × 1.5 m in size, placed, at least partially, as slablike debris [43] for a thickness of about 2 m. Considering the density values of 2.35/2.55 gr/cm³ measured, a maximum weight of about 90/140 tons can be calculated. Submerged boulders have different shapes.

BRs placed at a 6–9 m depth have a rectangular shape; they are characterized by well-rounded edges due to the continuous abrasion actions exerted by sands transported by waves and current (Figure 5). Their puzzle-like shape is reminiscent of a “manmade artifact” that has been lying there for a long time. Here, as well as in other places, they have been wrongly considered as harbor structures or submerged paths [44]. Their features differ from that of the boulders located on the foreshore or at a 3–5 m depth. Those found at the mean sea level have sharp edges, while those found at a depth of 3–5 m are sub-rectangular and less rounded. The very recent destruction of some boulders has not allowed for abrasion and consequent rounding of the boulder borders to occur (Figure 4C,D).

3.2. Thin Section Analyses

Analyses of thin sections using an optical microscope were performed. The samples studied from the Lido Torretta and Campomarino–Scorcialupi localities were emerged BRs, while those from San Pietro in Bevagna and Porto Cesareo were submerged BRs. Both emerged blocks and rounded submerged boulders are formed by medium-sized sand particles, well cemented and moderately sorted, classified as a packstone, composed of mollusk fragments, red coralline algae, echinoids, bryozoans, benthic foraminifers, intraclasts, and detrital grains. Isopachous carbonatic cement forms fringes of uniform crystals growing radially to grain surfaces. In the emerged boulders, pore spaces are filled with micrite rarely peloidal, in which small silty-sized skeletal fragments are dispersed (Figure 7).

Inside the Lido Torretta sample, a fragment of Sigillata Romana pottery, dating back to the II–VI CE century, has been preserved. These diagenetic features indicate early cementation of BRs in intertidal environments under the alternation of marine phreatic conditions. The silty matrix of the emerged boulders indicates an emergence and a second phase of diagenesis in a vadose environment [32]. The correlation with tide phases and diagenetic characteristics of the BRs indicates a fast cementation at the vadose–phreatic seabed interface, suggesting an approximation at a mean sea level identification of no more than 0.50 m [9,13].

3.3. Chronological Constraints

In the submerged BR samples, taken at a depth of 7+/−0.10 m, three bivalve shells were detected. With the aim of performing AMS age determinations, the shells were initially isolated and cleaned of clasts and cement using an electronic microhammer and then threated by means of an ultrasonic vibrator. Each sample was characterized for comparison with the mineralogical composition of present-day living ones before performing radiocarbon dating. The living shells and shell fragments were aragonitic. Indeed, in all the samples, an XRD qualitative analysis showed only the presence of aragonite within instrument detection limits of about 1%.

The shell samples were finely ground using an agate mortar. The powder obtained was homogenized, and an aliquot of about 1.5 g per sample was used for diffractometer analysis. Mineralogical analysis was performed with a Philips X’Pert Pro equipped with a spinner and a PANalytical X’Celerator solid-state detector. An X-ray tube with a Cu anode at 40 kV and 40 mA was used. The X-ray diffracted beam was filtered by a Ni foil to obtain only a Cu Kα1 radiation. The operating conditions were divergence slit, ¼°; secondary antidiffusion slit (antiscatter), ½°; sample rotation time (spinning), 2 s; secondary antidiffusion slit, 5.9 mm; step size, 0.023 °2θ; time for step size, 50 s; and angular range, 10 to 40 °2θ. The mineral phases were identified using a PANalytical High Score with the PAN ICSD database. Within instrument detection limits of about 1%, the qualitative analysis showed the presence of aragonite in all the samples, similar to that found in the present-day living ones.

Radiometric analyses provided an uncalibrated temporal range of 7.2 to 9.5 ka BP (

Table 1. Radiocarbon dating on the bivalve shells sampled in the beach rock.

Laboratory ID	Sample	Depth (m)	Type	Radiocarbon Age (BP)	δ13C (‰)	Calib Age BC (1σ)	Calib Age BC (2σ)
LTL21335	CM2	7 +/− 10	Cerastoderma spp.	7291 ± 40	−2.9 ± 0.3	5699–5559	5763–5481
LTL21336	CM3	7 +/− 10	Cardium spp.	7514 ± 40	−0.8 ± 0.3	5916–5758	5989–5695
LTL21334	CM1	7 +/− 10	Cardium spp.	9466± 65	0.4± 0.2	8294–8076	8418–7945

). Because the production of atmospheric radiocarbon has varied throughout geological time, the AMS ages were calibrated to provide ages in sidereal years with 1σ and 2σ ranges. The calibration performed with a CALIB 8.2 [45,46], using a Marine20 calibration curve, provided a temporal range of 5.5 to 8.4 ka BC.

Other chronological constraints derive from archaeological data. At Lido Torretta, near Torre Sgarrata, the uppermost tile of the BR, situated at roughly the level of the present intertidal zone, contains fragments of Roman pottery, identifiable as Sigillata Africana (II–VI century AD). The presence of these fragments hinders an absolute chronological attribution, but represents a “post quem” constraint in reference to the locality where the BR lithified. It is, however, plausible that the uppermost BR has an age dating back to at least the post-VI century along all the studied area.

4. Discussions

BRs can be detected along mobile coastal systems. They present a peculiar aspect owing to the early lithification of unconsolidated sediments. Beach rock lithification is a function of CO₃⁻² ion concentration in seawater, microbial activity, and degassing of CO₂ from seaward flowing groundwater [5,11,12]. The geographic distribution of BR deposits is mainly centered at midlatitudes, mostly between 20° and 40° N, and along microtidal coasts [4,11]. However, many authors underline that BRs can also form in supratidal environments since beach sediment lithification can be induced by cement precipitation from meteoric waters in relation to changes in water table elevation, temperature, and pressure (i.e., [18,47,48,49,50,51]).

The original morphology and texture of BRs vary from tiny patches—the collapsed “domino tiles”—to slablike debris of cemented sediments outcropping hundreds of meters. Their thickness varies from less than 0.5 m to more than 2.5 m, resulting in greater thickness in areas with more pronounced tidal fluctuations [5,6,11]. BRs dip mostly seaward with slopes of up to 15°, generally following the beach slope [8,9,47,48,49,50,51]. In other cases, BRs show markedly dissimilar slopes compared with those of the host beaches [16,20,51], and are characterized by sedimentary structures developing later due to the subsidence, breaking, and tilting of BR slabs [4].

According to Mauz et al. [12] and Strasser et al. [15], BR cement is indicative of the interface between seawater and meteoric water (marine and phreatic conditions), which ascribe both the marine vadose and marine phreatic environments. Under these conditions, the sediments are characterized by pore spaces through which the solutions can flow [8]. The solutions can have end members composed of calcite, high-magnesium calcite, and aragonite, which precipitate and fill the intergranular pore spaces. The presence of high-magnesium calcite and aragonite is common in foreshore zones, while the cement formed by low-magnesium calcite is common in the shoreface and backshore zones, facilitated by the presence of red algae and conditioned by the sediment composition [12,52,53]. The latter features suggest that the studied BRs are representative of past sea-level stands in the last 8 ka.

The validity of the proposed AMS age determination is supported by the XRPD analysis performed on the dated shell; it shows the quasi-total prevalence of aragonite and indicates that shells were not affected by recrystallization during cementation. The presence of an older sample compared with the others (c.ca 10 ka vs 7.8 ka) can be explained as the effects of the cementation of a reworked shell in a younger sediment. Moreover, the presence of Sigillata Africana pottery, dating back to the II–VI century AD, found in the uppermost BR, indicates that the three BRs formed during different phases in a span of time ranging from 8 ka to the present.

A more precise reconstruction of past sea-level stands during this time can be addressed by comparing the geomorphological evidence with the available modelled sea-level curves. We considered both the curve from Lambeck et al. [54] and the curve ICE-7G [55,56]; they present significant discrepancies between them, of even tens of meters, in the function of the considered glacial isostatic adjustments. Comparison between the two above-mentioned sea-level curves with our AMS ages showed good agreement with the ICE-7G curve, thereby corroborating the validity. Indeed, field data must be used to improve the modelled curve since the latter cannot discriminate local contributions, i.e., isostatic, tectonic, or sedimentary [57,58,59,60,61].

In order to draw conclusions, we must consider that (ì) the deepest BR is about 7.8 ka in age (samples CM2 and CM3) and that the youngest one is at the present sea level; (ìì) negligible vertical land movements on a long- [32,33,34] and short-term scale occurred in this area as indicated by permanent GPS stations [62,63,64]; and (ììì) archeological sea-level markers located in proximity to this area seem to suggest a limited and local downlift [38,39,40].

As a result, the present-day positions of the studied BRs should indicate three different relative sea-level stands. Following the rapid post-LGM rise, the sea-level rise decelerated sharply between 8.2 and 6.7 ka BP in correspondence to the final phase of the North American deglaciation, when the reduction of the meltwater input resulted in a significant deceleration on a global scale [54]. This is also recorded by BRs along several Mediterranean coasts, such as that of the Bonifacio Strait [65] and the coastal plains in Sardinia [66,67,68] and southern Tuscany [69]. BRs recognized along the Ionian coast of Calabria have been useful in reconstructing the local VLM [17,70,71]. In our case, the deepest BR, located at a depth range of 6–9 m, is chronologically correlated with a sea-level stand prior to the formation of the Holocene dune belt aged about 7.0 ka, occurring at the end of the rapid postglacial transgression. The BRs located at a depth range of 3–5 m could be chronologically related to a temporal range following the slowdown of the transgression, as shown in the ICE-7G curve (Figure 8).

Considering local conditions, these pieces of evidence together suggest that BRs are a useful tool in reconstructing the Holocene sea-level changes [1,65], but in the present case-study, they do not have the same effectiveness in the VLM assessment.

5. Conclusions

Due to their morphological extension and continuity, the BRs studied at Campomarino di Maruggio and San Pietro in Bevagna indicate that they could be considered representative of the Holocene sea-level history of the entire coastal area stretching between Taranto and Capo Santa Maria di Leuca. This means that they have regional significance and acquire a SLIP role. BRs provide data to assess the local sea-level changes during the Holocene as a result of the correlation between field, analytic, and modelling data.

In particular, bathymetric and dive surveys allowed the assessment of the continuity and the extension of the BR landform with the relative depth ranges, thereby confirming their regional significance as markers of relative sea-level changes. Furthermore, petrographic analysis and radiocarbon dating allowed the temporal constraints for the BR diagenesis to be obtained. The AMS data deriving from the deepest BR located at 6–9 m, as well as the archaeological contents found in the uppermost one, yielded two constraints, which led to the chronological attribution of the other BR.

The temporal phase, during which the relative sea level rose, conditioned the lithification of the three BRs. This corresponds to the reduction of the meltwater input from the ice sheets, which resulted in a significant deceleration of the global sea-level rise in the final phase of the North American deglaciation. The slow sea-level rise that followed the deceleration, about 7.8 BP, was characterized by a relative sea-level stand of 7.0 +/− 1 m b.m.p.s.l. that drove the lithification of the deeper BR. A second phase occurred when the sea-level stand was at 3–5 m b.m.p.s.l. about 6/5.5 ka BP. The last phase of the BR lithification continued up to the present, as evidenced by the study of the uppermost BR, as well as by the presence of Sigillata Africana pottery fragments found in it.

A comparison between our results and available isostatic models suggests that the ICE-7G model has proven to be the most appropriate to represent the Mid- to Late Holocene sea-level changes in this area of the Mediterranean basin. The slight discrepancies between the BR CM2 and CM3 ages and their positions on the modelled curves may be considered “physiological”.

Archaeological and geomorphological data obtained from this part of the Apulian foreland, relating to the last 8 ka, could suggest limited VLM in different directions. In any case, we have no evidence of their nature, and any attribution to one or another cause is, at the moment, only speculation. Actually, they could be considered evidence not only of local subsidence but also of a differential and local effect of the hydro-isostatic adjustment on the coastal part of the continental shelf as a consequence of 120 m in sea-level variation and the related increase in the water column load. Due to the absence of an accurate quantification of vertical displacements and of their fingerprints, at present, they cannot be used to correct the depth of beach rocks and, in turn, of the ICE-7G curve.