“Continuous” Backstepping of Holocene Coastal Barrier Systems into Incised Valleys: Insights from the Ofanto and Carapelle-Cervaro Valleys

Two recently recognised incised valleys in the Manfredonia Gulf are described. The first (CCV) is correlated with the current Carapelle and Cervaro streams. The second (OSFV) is correlated mostly with the current Ofanto River. Six seismic facies and seven unconformity-bounded seismic units have been identified, which infilled CCV and OSFV. In CCV, during the sea-level ranges from −29 to −18 and from −18 to −4.7 m b.s.l., two barrier/spit-backbarrier systems formed in the most landward sector of the valley. The lower system was attributed to a time interval between 9.2 ka BP and ca. 8.3 ka BP, chronologically constrained by the ZS2 borehole. In OSFV, during the sea level ranges from −39 to −29, and from −29 and to −18 m b.s.l., two beach/spit-backbarrier systems, arranged in a “continuous” landward backstepping pattern, formed. The phase that contributed most to the beach/spit-backbarrier systems formation is that which is coeval with the formation of the sapropel S1 in the Mediterranean. The conservation of barrier/spit-backbarrier systems arranged in a “continuous” landward backstepping pattern, is due to a strong and continued sediment supply that occurred during the sapropel S1 formation, coupled with low-gradient settings and a regime of slow sea-level rise.


Introduction
The term coastal barrier refers to different geomorphological features at different times. Otvos [1] reviewed the nomenclature and classification of coastal barriers and considered two main types: (1) barrier islands, parallel to the shore and located between the inshore and open marine environments, and (2) barrier spits, which are linked to the mainland shore. Modern barrier islands are found worldwide, located mainly along wave-dominated coasts, and these islands were subjected to rising sea levels of the late Holocene [2]. Barrier-lagoon systems include the shoreface sector, the barrier, and the backbarrier sector, where wash-over fans and lagoons develop [1,3].
The formation of these systems is controlled by the interplay of wave energy, sediment supply, accommodation and sea-level oscillations [4][5][6]. In particular, barriers and backbarrier-lagoons can develop when the substrate gradient is less than 0.8 • [7].
The responses of coastal barrier and backbarrier systems to sea-level rise are crucial issues, in the light of the future sea level rise due to global warming. Studies of evolving barrier morphology indicate several responses to sea-level rise, all of which are strongly influenced by local factors. Most low-gradient coasts form transgressive barriers during sea-level rise. These barriers were composed almost entirely of tidal delta and wash-over deposits that retrograded into estuarine/lagoonal environments as sea-level rose [7]. The seabed exposed by the retreating barrier is a ravinement and of its depositional lobe, high-stand system tract (HST) Gargano sub-acqueous delta (gray tones) and bathymetry of the Manfredonia Gulf shelf with 10-m-spaced contours. Redrawn from literature data [29,31,32].
The Manfredonia Gulf was also affected by deposition of the so-called Gargano subaqueous delta [33]. This delta was formed on the eastern and southeastern sides of the Gargano promontory ( Figure 1). This subaqueous deposit represents the southernmost portion of the Holocene HST developing along the western side of the Adriatic. The HST rests above a regional down-lap surface (maximum flooding surface; mfs) that marks the time that the maximum landward shoreline shift occurred, approximately 5.5 cal. ka BP during the end of the late-Pleistocene to Holocene sea-level rise [30,[33][34][35].
Following the work of Maselli et al. [36] in the central Adriatic Sea, the TST is found to consist of a tripartite deposit, recording the step-wise nature of sea-level rise. The three distinct units, each deposited during a specific interval of the last sea-level rise, are the lower TST unit (lTST), middle TST (mTST) and upper TST (uTST). The lTST formed as sea-level rose between −120 and −95 m, the mTST between −95 and −60 m, and the uTST between −60 m and the position of the shoreline that was reached at the time of maximum marine ingression (5.5 cal. ka BP). The boundary between the lTST and mTST consists of the S1 surface, whereas the boundary between the mTST and the uTST consists of the S2 surface. S1 and S2 are regional erosional surfaces generated by reduced sediment inputs and enhanced marine reworking as a consequence of Meltwater Pulses 1A and 1B (MWP-1A and MWP-1B), respectively [37]. Within the mTST, another regional erosional surface (Si), can be detected, thus dividing the mTST into two different subunits (mTST-1 below and mTST-2 above the Si, respectively). The stratigraphic complexity of the mTST unit, recording the Bölling-Allerød and Younger Dryas events, can best be explained by sea-level oscillations, including a minor sea-level fall during the Younger Dryas event, which likely formed the Si surface. This conclusion is consistent with the formation of a Younger Dryas shoreline deposit ca. −75 m depth, which is within the mTST and of its depositional lobe, high-stand system tract (HST) Gargano sub-acqueous delta (gray tones) and bathymetry of the Manfredonia Gulf shelf with 10-m-spaced contours. Redrawn from literature data [29,31,32].

Seismic Survey
Seismic surveys were carried out in the Manfredonia Gulf (southern Adriatic) in 2016 and 2017 using the motorboat "ISSEL", which is owned by the Interuniversity Consortium for Sea Sciences (Con.I.S.Ma.). The following instruments and sensors were mobilised: (i) AppliedAcoustic Squid CSP-P 350J "sparker" source; and (ii) Geo-Resources Geo-Sense single channel seismic acquisition system with 8 hydrophone streamer and amplification/filtering stage.
The data was processed using the IXSEA DELPH software. Through the use of the system, all the necessary steps for filtering, gaining and interpreting the data have been performed.
Track-line positioning is based on D-GPS navigation, assuring a sub-metric position accuracy using the WGS84 datum.
The seismic survey comprised 15 sparker profiles oriented perpendicular to the coastline and 12 sparker profiles subparallel to the coastline for a total of ca. 580 km ( Figure 2). The minimum and maximum water depths were <8 m and ca. 97 m, respectively.
The data were subsequently converted into a digital format/image for interpretation. Each seismic profile was imported into TEI Delph Map with Seismic GIS, allowing for immediate georeferencing Water 2020, 12, 1799 4 of 32 and interpretation of acquired data. Finally, a buried erosional surface (ES1) was reconstructed using a triangular interpolation method in GIS.
The seismic survey comprised 15 sparker profiles oriented perpendicular to the coastline and 12 sparker profiles subparallel to the coastline for a total of ca. 580 km ( Figure 2). The minimum and maximum water depths were <8 m and ca. 97 m, respectively.
The data were subsequently converted into a digital format/image for interpretation. Each seismic profile was imported into TEI Delph Map with Seismic GIS, allowing for immediate georeferencing and interpretation of acquired data. Finally, a buried erosional surface (ES1) was reconstructed using a triangular interpolation method in GIS. In addition to the surveys described above, this work reports the seismic profile T1, which was acquired and processed within other surveys carried out between 2013 and 2014 using sub-bottom profiler Benthos CHIRP III in the Manfredonia Gulf, immediately NW of the area investigated in the present study ( Figure 2). For the technical details of that survey see [32,38].

Seismic Facies and Unconformity-Bounded Seismic Units
For this work, seismic facies and unconformity-bounded seismic units (UBSUs) belonging to the TST have been identified (Figure 3). In addition to the surveys described above, this work reports the seismic profile T1, which was acquired and processed within other surveys carried out between 2013 and 2014 using sub-bottom profiler Benthos CHIRP III in the Manfredonia Gulf, immediately NW of the area investigated in the present study ( Figure 2). For the technical details of that survey see [32,38].

Seismic Facies and Unconformity-Bounded Seismic Units
For this work, seismic facies and unconformity-bounded seismic units (UBSUs) belonging to the TST have been identified ( Figure 3).
By seismic facies we mean groups of reflections whose parameters (bedform geometry, lateral continuity, amplitude, frequency, and interval velocity) differ from those of adjacent groups of reflections [39]. We also used the spatial extension of seismic facies, along with the reciprocal position between facies, to infer depositional environments and formative processes. The seismic facies interpretation is also supported by previous study, e.g., some general seismic facies characteristics of the latest Quaternary transgressive deposits and incised valley infillings were established (i.e., [13,15,[40][41][42][43]). More specifically for the southern Adriatic, Maselli and Trincardi [30], Maselli et al. [36], and Maselli et al. [31] provide a complete description of the seismic facies that characterise the TST and HST. By seismic facies we mean groups of reflections whose parameters (bedform geometry, lateral continuity, amplitude, frequency, and interval velocity) differ from those of adjacent groups of reflections [39]. We also used the spatial extension of seismic facies, along with the reciprocal position between facies, to infer depositional environments and formative processes. The seismic facies interpretation is also supported by previous study, e.g., some general seismic facies characteristics of the latest Quaternary transgressive deposits and incised valley infillings were established (i.e., [13,15,[40][41][42][43]). More specifically for the southern Adriatic, Maselli and Trincardi [30], Maselli et al. [36], and Maselli et al. [31] provide a complete description of the seismic facies that characterise the TST and HST.
The seismic facies interpretations are also based on two vibro-cores drilled in 2015 in the Manfredonia Gulf: SP1_VC24 (intercepted by profile P2); and SP1_VC25 (close to profile PX2) on behalf of the Puglia Basin Authority (AdB) for a sand resource study of the Apulian continental shelf. The chronological constraint of the transgressive UBSUs is based on a borehole drilled in Zapponeta town (ZS2) in 2012 by the AdB during a seismic microzonation project; this borehole was directly analysed by the authors. Seismic facies, coupled with the identification of discontinuity surfaces, were used to separate the TST sediments into UBSUs. More precisely, UBSUs have been distinguished based on unconformities, interpretable as surfaces of erosion or discontinuity in sedimentation. The UBSUs in some cases are made up of a single seismic facies, in other cases of multiple seismic facies ( Figure 3) that transit into one another (i.e., the seismic units described in [17] and [44]).
In this work, the seismic facies have been recognised based only on seismic sparker profiles, because they are almost all of the profiles present in the study area. In the single sub-bottom profile T1, only the unconformity bounded seismic units have been recognised.

Age Range Evaluation and Chronological Constraints
For an approximate assessment of the age of the UBSUs, we have modified the method based on the current depth of the units, already applied in the southern Adriatic (see for example the age estimation made by Maselli and Trincardi [30] for the units 2C and 2D of area 3 of the MIV or the age estimation made by Pellegrini et al. [45] for the Palaeo Gargano compound delta).
One of the seismic facies (seismic facies 2) was found to be helpful in an approximated evaluation of the age range of the UBSUs. Seismic facies 2 has been interpreted in this study as a high-energy marine deposit and therefore deposited at (or close to) sea level; this interpretation is validated also because the seismic facies 2 is intercepted by vibro-cores SP1_VC24 and SP1_VC25, where it is present between the sea bottom and the bottom of the two vibro-cores. In addition, the interpretation of seismic facies 2 is in accordance with Trincardi et al. [28] (p. 86) which defined the The seismic facies interpretations are also based on two vibro-cores drilled in 2015 in the Manfredonia Gulf: SP1_VC24 (intercepted by profile P2); and SP1_VC25 (close to profile PX2) on behalf of the Puglia Basin Authority (AdB) for a sand resource study of the Apulian continental shelf. The chronological constraint of the transgressive UBSUs is based on a borehole drilled in Zapponeta town (ZS2) in 2012 by the AdB during a seismic microzonation project; this borehole was directly analysed by the authors. Seismic facies, coupled with the identification of discontinuity surfaces, were used to separate the TST sediments into UBSUs. More precisely, UBSUs have been distinguished based on unconformities, interpretable as surfaces of erosion or discontinuity in sedimentation. The UBSUs in some cases are made up of a single seismic facies, in other cases of multiple seismic facies ( Figure 3) that transit into one another (i.e., the seismic units described in [17] and [44]).
In this work, the seismic facies have been recognised based only on seismic sparker profiles, because they are almost all of the profiles present in the study area. In the single sub-bottom profile T1, only the unconformity bounded seismic units have been recognised.

Age Range Evaluation and Chronological Constraints
For an approximate assessment of the age of the UBSUs, we have modified the method based on the current depth of the units, already applied in the southern Adriatic (see for example the age estimation made by Maselli and Trincardi [30] for the units 2C and 2D of area 3 of the MIV or the age estimation made by Pellegrini et al. [45] for the Palaeo Gargano compound delta).
One of the seismic facies (seismic facies 2) was found to be helpful in an approximated evaluation of the age range of the UBSUs. Seismic facies 2 has been interpreted in this study as a high-energy marine deposit and therefore deposited at (or close to) sea level; this interpretation is validated also because the seismic facies 2 is intercepted by vibro-cores SP1_VC24 and SP1_VC25, where it is present between the sea bottom and the bottom of the two vibro-cores. In addition, the interpretation of seismic facies 2 is in accordance with Trincardi et al. [28] (p. 86) which defined the sediments cropping out at the sea bottom in the area were the two vibro-cores have been realised as "proximal transgressive sandy deposits of beach or paralic environment".
Having said that, because seismic facies 2 forms sedimentary bodies within different UBSUs arranged in a landward backstepping pattern [46] on the ES1, it follows that each of these bodies indicates an approximated sea-level range if we consider the base and the top of the body. Then, we considered the local sea-level curve (site 25 in [47]) to infer two ages corresponding to two sea levels that we associated with the base and top of each sedimentary body formed by seismic facies 2. These two sea levels have been corrected considering a possible subsidence calculated for the Holocene Water 2020, 12, 1799 6 of 32 deposits of our study area. In conclusion, we associate each sedimentary body formed by seismic facies 2 with two approximated palaeo sea-levels and therefore with an approximated age range ( Figure 4).
Having said that, because seismic facies 2 forms sedimentary bodies within different UBSUs arranged in a landward backstepping pattern [46] on the ES1, it follows that each of these bodies indicates an approximated sea-level range if we consider the base and the top of the body. Then, we considered the local sea-level curve (site 25 in [47]) to infer two ages corresponding to two sea levels that we associated with the base and top of each sedimentary body formed by seismic facies 2. These two sea levels have been corrected considering a possible subsidence calculated for the Holocene deposits of our study area. In conclusion, we associate each sedimentary body formed by seismic facies 2 with two approximated palaeo sea-levels and therefore with an approximated age range ( Figure 4). If the body formed by seismic facies 2 belongs to an UBSUcontaining multiple seismic facies, we extend to the whole unit the approximated age range obtained from the facies 2 body. This procedure allowed us to link the UBSUs at some sea-level ranges and therefore at corresponding age ranges.
The timing of the deposition of transgressive UBSUs was chronologically constrained by the ZS2 borehole, which intercepts the transgressive fill of the inland sector of an incised valley recognised in this study. Three samples of this perforation were dated using the 14C method at DirectAMS Laboratories, Bothell, Washington (laboratory ID: D-AMS). The radiocarbon dates were calibrated using the CALIB 7.1 program. The IntCal13, Marine/IntCal13, and Marine13 datasets [48] were used for terrestrial, transitional and marine samples, respectively ( Table 1). The percentage of marine carbon was calculated based on the δ13C value [49]. A ΔR correction (ΔR =121 ± 60) calculated for the southern Adriatic Sea off Barletta (http://intcal.qub.ac.uk/marine/) was used. Furthermore, the ZS2 borehole intercepted a pumice layer, which has been identified.

Data and Results
We present the following seismic and stratigraphic data regarding the most inner part of the Manfredonia Gulf continental shelf: (1) the erosional surface that formed during the last sea-level fall, sub-aerial exposure of the shelf and erosional reworking during subsequent sea level rise; (2) the main landforms described by the erosional surface; (3) borehole and vibro-cores; (4) seismic facies recognizable within uTST; and (5) UBSUs belonging to the uTST. If the body formed by seismic facies 2 belongs to an UBSUcontaining multiple seismic facies, we extend to the whole unit the approximated age range obtained from the facies 2 body. This procedure allowed us to link the UBSUs at some sea-level ranges and therefore at corresponding age ranges.

Erosional Surface ES1
The timing of the deposition of transgressive UBSUs was chronologically constrained by the ZS2 borehole, which intercepts the transgressive fill of the inland sector of an incised valley recognised in this study. Three samples of this perforation were dated using the 14C method at DirectAMS Laboratories, Bothell, Washington (laboratory ID: D-AMS). The radiocarbon dates were calibrated using the CALIB 7.1 program. The IntCal13, Marine/IntCal13, and Marine13 datasets [48] were used for terrestrial, transitional and marine samples, respectively ( Table 1). The percentage of marine carbon was calculated based on the δ13C value [49]. A ∆R correction (∆R =121 ± 60) calculated for the southern Adriatic Sea off Barletta (http://intcal.qub.ac.uk/marine/) was used. Furthermore, the ZS2 borehole intercepted a pumice layer, which has been identified.

Data and Results
We present the following seismic and stratigraphic data regarding the most inner part of the Manfredonia Gulf continental shelf: (1) the erosional surface that formed during the last sea-level fall, sub-aerial exposure of the shelf and erosional reworking during subsequent sea level rise; (2) the main landforms described by the erosional surface; (3) borehole and vibro-cores; (4) seismic facies recognizable within uTST; and (5) UBSUs belonging to the uTST.

Erosional Surface ES1
The seismic survey allowed us to identify a regional erosional surface (ES1). In some cases, the ES1 appears as a clearly visible reflector, with moderate to high amplitude, coinciding with a more or less evident angular unconformity between the underlying reflectors and the overlying reflectors; in other cases, the ES1 appears instead as an abrupt contact between very different seismic facies. However, in both cases, ES1 is identified as a surface that separates the substrate from younger sediments having a very different seismic architecture. The substrate is almost always characterised by seaward-inclined reflectors, referable to oblique-tangential clino-forms [39] which dip at low angles (usually fractions of a degree) and are truncated up-dip by the ES1. Sediments above ES1, on the other hand, are highly variable in that they are characterised by numerous seismic facies ( Figure 5).
Water 2020, 12, x FOR PEER REVIEW 7 of 32 The seismic survey allowed us to identify a regional erosional surface (ES1). In some cases, the ES1 appears as a clearly visible reflector, with moderate to high amplitude, coinciding with a more or less evident angular unconformity between the underlying reflectors and the overlying reflectors; in other cases, the ES1 appears instead as an abrupt contact between very different seismic facies. However, in both cases, ES1 is identified as a surface that separates the substrate from younger sediments having a very different seismic architecture. The substrate is almost always characterised by seaward-inclined reflectors, referable to oblique-tangential clino-forms [39] which dip at low angles (usually fractions of a degree) and are truncated up-dip by the ES1. Sediments above ES1, on the other hand, are highly variable in that they are characterised by numerous seismic facies ( Figure  5). , showing the erosional surface (ES1) on a regional scale. Note that ES1 is identified as a surface that separates the substrate from younger sediments having a different seismic architecture.

Erosive Landform: Incised Valley A and B
ES1 presents three main erosive landforms ( Figure 6). The north-western, not the subject of this work, is already known in the literature as Manfredonia incised valley (MIV [30,32]).
The central is an incised valley (hereinafter incised valley 1: IV1; Figure 6), that reaches a depth (in profiles T8 and P1) of ca. −37 m b.s.l. and ca. 18/22 m below the surrounding ES1. In the P2 profile, the maximum depth is ca. −38 m. b.s.l., which reaches −50 m. b.s.l. in profile P4.
This valley has a meandering course at the turn of the current coastline. The south-eastern landform is another incised valley (hereinafter incised valley 2: IV2; Figure  6), which appears, in the sector closest to the coast, as two different sub-incisions. These two , showing the erosional surface (ES1) on a regional scale. Note that ES1 is identified as a surface that separates the substrate from younger sediments having a different seismic architecture.

Erosive Landform: Incised Valley A and B
ES1 presents three main erosive landforms ( Figure 6). The north-western, not the subject of this work, is already known in the literature as Manfredonia incised valley (MIV [30,32]).
The central is an incised valley (hereinafter incised valley 1: IV1; Figure 6), that reaches a depth (in profiles T8 and P1) of ca. −37 m b.s.l. and ca. 18/22 m below the surrounding ES1. In the P2 profile, the maximum depth is ca. −38 m. b.s.l., which reaches −50 m. b.s.l. in profile P4.
This valley has a meandering course at the turn of the current coastline.
Water 2020, 12, 1799 8 of 32 The south-eastern landform is another incised valley (hereinafter incised valley 2: IV2; Figure 6), which appears, in the sector closest to the coast, as two different sub-incisions. These two sub-incisions merge in a single wide and shallow valley with a very rugged and irregular bottom. The IV2 reaches a maximum width of 11.9 km (profile P2) and narrows to approximately 2.5 km in correspondence with the Barletta strait ( Figure 6). The IV2 reaches depths between ca. Above the ES1 surface, a series of units were deposited during the last transgressive phase; they constitute the infilling of the IV1 and IV2.

Borehole and Vibrocores
The borehole and vibro-cores intersect only transgressive sediments lying above ES1.

Borehole and Vibrocores
The borehole and vibro-cores intersect only transgressive sediments lying above ES1.  Three samples from this borehole have been 14C dated (Table 1). In addition, the level of pumice present from 28.

Vibro-cores SP1_VC24 and SP1_VC25
Vibro-core SP1_VC24 was obtained 5 km offshore at a water depth of 20 m, and this core was intercepted by the seismic profile P2 ( Figure 2). Coring ceased at 5.13 m below the seabed, and from the bottom up the core consists of the following units ( Figure 7 from 0.79 m to the sea floor, bioclastic gravel of bivalves and gastropods in a clayey-sandy, dark grey matrix. Vibro-core SP1_VC25 was obtained 8 km offshore at a water depth of 24.6 m and is located very close to seismic profile PX2 (Figure 2). The coring ceased at 5.20 m below the seabed, and from the bottom upwards, the core consists of the following units ( Figure 7): from 0.78 m to the sea floor, fine and very fine ochre sand.

Seismic Facies
In this section, we describe the seismic facies recognised in the sediments belonging to the TST, which cover ES1 and infill the IV1 and IV2.
• Facies 1 ( Figure 8A,B). This facies is only present at the bottom of the incised valleys and is characterised by high amplitude down-lapping reflector or sets of concave and/or convex and/or wavy reflectors; reflectors also show frequent changes in dip direction; the lateral continuity is low-medium. In the seismic profiles normal to the narrow valley axes, the reflectors present a conformable channel infill pattern. • Facies 2. These seismic facies can be divided into two subtypes: 2a and 2b. Facies 2a ( Figure 8C): this facies, in general, has a chaotic, non-symmetrical internal configuration with reflectors of low or very low amplitude and very low lateral continuity. Commonly, these facies contains V and Λ-shaped reflectors. Alternatively, the reflectors are undulating, concave, convex or oblique-parallel, but always with low amplitude and very low lateral continuity. Facies 2b ( Figure 8D): the reflectors are sigmoidal and/or clinoform with moderate-high amplitude and medium-high lateral continuity, which alternate with low amplitude reflectors or with chaotic and non-symmetrical packages, and all are inclined landward. Facies 2 overlies ES1 on the shelf and previous infilling units in IV1 and IV2. From offshore to onshore, a transition between facies 2 and facies 3 or 4 is often observed. • Facies 3 ( Figure 8E,F). This facies is characterised by moderate amplitude and plane-parallel reflectors with medium lateral continuity; sometimes, reflectors are weakly undulating with low lateral continuity. This facies overlies ES1 or seismic facies 1. The relationship with seismic facies 2 is twofold: in some cases, facies 2 transitions landwards into facies 3; in other cases, facies 3 is surmounted by facies 2. • Facies 4 ( Figure 8G,H). This facies is characterised by moderate to high amplitude, medium-high lateral continuity plane-parallel reflectors alternating with wavy or gently inclined reflectors. The relationship of facies 4 with seismic facies 2 is twofold: in some cases, facies 2 transitions landwards into facies 4; in other cases, facies 5 is surmounted by facies 2. • Facies 5 ( Figure 8I). This facies is characterised by moderate amplitude and medium lateral continuity reflectors, with the presence of a topset, fore-set and bottomset; thus, it is a pro-gradational seismic facies. Downdip, this facies transitions into facies 3. • Facies 6 ( Figure 8). This facies is characterised by moderate to high amplitude and maximum observed lateral continuity reflectors that onlap onto inclined substrate and/or draping over uneven substrate. Reflectors are gently seaward-inclined ( Figure 8L) or sub-horizontal ( Figure 8M).  Facies 3 ( Figure 8E,F). This facies is characterised by moderate amplitude and plane-parallel reflectors with medium lateral continuity; sometimes, reflectors are weakly undulating with low lateral continuity. This facies overlies ES1 or seismic facies 1. The relationship with seismic facies 2 is twofold: in some cases, facies 2 transitions landwards into facies 3; in other cases, facies 3 is surmounted by facies 2.  Facies 4 ( Figure 8G,H). This facies is characterised by moderate to high amplitude, medium-high lateral continuity plane-parallel reflectors alternating with wavy or gently inclined reflectors. The relationship of facies 4 with seismic facies 2 is twofold: in some cases, facies 2 transitions landwards into facies 4; in other cases, facies 5 is surmounted by facies 2.  Facies 5 ( Figure 8I). This facies is characterised by moderate amplitude and medium lateral continuity reflectors, with the presence of a topset, fore-set and bottomset; thus, it is a pro-gradational seismic facies. Downdip, this facies transitions into facies 3.

Unconformity-Bounded Seismic Units
We identified seven UBSUs that lie above ES1 and/or infilled the IV1 and IV2.

Unit A u
Unit A u (the subscript letter u in the abbreviation of the units indicates that all units belong to upper part of TST) consists of seismic facies 1, occupies discontinuous portions of the IV1 and IV2 bottom and lies always above ES1 (Figures 9-12). Unit A u reaches a maximum thickness of ca. 9 m.

Unit B u
Unit B u extends between the offshore boundary of the study area and an ES1 depth of ca. −40 m b.s.l. and it consists exclusively of seismic facies 2a. This unit does not affect IV1 or IV2 and always lies above ES1 (Figures 9 and 12). Unit B u reaches maximum thickness of ca. 5 m.

Unit C u
Unit C u occupies only a small lower area of the IV2 at the intersection of profiles T5 and P3 (Figures 6 and 9), and it never exceeds a depth of −38 b.s.l. Overall, this unit has a lenticular shape, lies above unit A u and ES1, reaches a maximum thickness of ca. 5 m and consists exclusively of seismic facies 3.

Unit D u
Unit D u is present between depths of ca. −40 m b.s.l. and ca. −30 b.s.l., and lies on unit A u , unit B u , and ES1.
In the IV1 sector (profile T8, Figure 10), unit D u is present exclusively seawards and outside the IV1 and lies on the ES1. Here, unit D u is composed of seismic facies 2a, with some sets of landward-dipping reflectors, and the unit reaches a maximum thickness of ca. 7 m.
Passing to IV2 sector, unit D u changes proceeding from the NW to the SE. At the NW (profile T6), unit D u consists only of the seismic facies 2a, lies on ES1 and is present exclusively seawards and outside the IV2. Moving to the SE (profile T5; Figure 9), unit D u centres on the threshold that separates the IV2 from the continental shelf and consists of a central body formed by seismic facies 2a which transitions into seismic facies 3 inside the IV2 (landwards) and seismic facies 6 on the shelf (seawards). The central body is characterised by two sets of opposing inclined reflectors towards the NE and the SW (Figure 9). This internal organisation of the seismic facies within unit D u remains even further to the SE (profile T4; Figure 12), with the difference that the central body of seismic facies 2a appears narrower than that in profile T5 and transitions into seismic facies 4 inside the IV2.
Into the IV2, unit D u extends landwards much further in the NW area (where it appears with seismic facies 3) than in the SW area, where unit D u it is not present (profile PX2; Figure 11A-C).
The central body formed by seismic facies 2 stops at the Barletta strait. This is visible in profile PX3 which is elongated along the external threshold of the IV2 ( Figure 6) and then it runs across the central body of the unit D u .

Unit E u
Unit E u is present between depths of ca. −30 m and ca. −16 b.s.l., and lies on ES1, unit A u , unit B u , and unit D u .
In the IV1 sector, unit E u is formed by seismic 2a both on the shelf and in the seaward part of the valley (profile T8, Figure 10). In the landward part of the valley, unit E u consists of seismic facies 3 at the base, surmounted by seismic facies 2b ( Figure 10). In the IV1 sector, unit Eu is formed by seismic 2a both on the shelf and in the seaward part of the valley (profile T8, Figure 10). In the landward part of the valley, unit Eu consists of seismic facies 3 at the base, surmounted by seismic facies 2b ( Figure 10). Passing to IV2 sector, unit Eu appears to be the most complex system in terms of internal organisation, and it changes proceeding from the NW to the SE. To the NW (profile T6), unit Eu is composed of a very thin central body of seismic facies 2a, which seaward transitions into seismic facies 6 and landwards transitions first into facies 3 and then into seismic facies 5; this latter forms a body up to 10 m tick with an evident rollover point between the fore-set and the topset at ca. −19 m ( Figure 13). Moving to the SE (profile T5; Figure 9), the central body reaches a thickness of ca. 8 m and transitions landward into seismic facies 3 and seaward into seismic facies 6. Even further to the SE (profile T4; Figure 12), unit Eu is constituted by a body formed by seismic facies 2a, up to ca. 5 m thick, which transitions landwards into seismic facies 4. Passing to IV2 sector, unit E u appears to be the most complex system in terms of internal organisation, and it changes proceeding from the NW to the SE. To the NW (profile T6), unit E u is composed of a very thin central body of seismic facies 2a, which seaward transitions into seismic facies 6 and landwards transitions first into facies 3 and then into seismic facies 5; this latter forms a body up to 10 m tick with an evident rollover point between the fore-set and the topset at ca. −19 m ( Figure 13). Moving to the SE (profile T5; Figure 9), the central body reaches a thickness of ca. 8 m and transitions landward into seismic facies 3 and seaward into seismic facies 6. Even further to the SE (profile T4; Figure 12), unit E u is constituted by a body formed by seismic facies 2a, up to ca. 5 m thick, which transitions landwards into seismic facies 4. The central body of unit Eu, constituted by seismic facies 2a, is much wider in the NW sector of the IV2 compared to the SE sector. In fact, facies 2a reaches landwards the profile PX2 ( Figure  11A-C) only in the NW sector, while, in the SE sector, it gives way to facies 4 ( Figure 11A-C).
The deposition of unit Eu involves most of the IV2 infilling and its division into two residual sub-basins, which consist of two areas not completely infilled. The first is a smaller sub-basin located to the NW in front of the sedimentary body formed by seismic facies 5 visible in profile T6 ( Figure  13). The second is a larger sub-basin located to the SE in the sector behind the Barletta strait and is visible in profiles PX2 and T4 (Figures 11 and 12).  In the landward side of the profile, note the superimposition of unit E u and unit G u , which are two successive barrier/spit-backbarrier systems. In both systems, as the transgression progressed, the barrier/spit retreated landwards and covered the backbarrier sediments through over-wash and/or tidal inlet deposits, as evidenced by the reflectors in seismic facies 2b down-lapping landwards onto the plane-parallel reflector of seismic facies 3.
The central body of unit E u , constituted by seismic facies 2a, is much wider in the NW sector of the IV2 compared to the SE sector. In fact, facies 2a reaches landwards the profile PX2 ( Figure 11A-C) only in the NW sector, while, in the SE sector, it gives way to facies 4 ( Figure 11A-C).
The deposition of unit E u involves most of the IV2 infilling and its division into two residual sub-basins, which consist of two areas not completely infilled. The first is a smaller sub-basin located to the NW in front of the sedimentary body formed by seismic facies 5 visible in profile T6 ( Figure 13). The second is a larger sub-basin located to the SE in the sector behind the Barletta strait and is visible in profiles PX2 and T4 (Figures 11 and 12). 4.5.6. Unit F u Unit F u infills the NW sub-basin and lies always on unit E u . Unit F u formed after deposition of unit E u , in front of the sedimentary body formed by seismic facies 5 ( Figure 13). Unit F u consists of seismic facies 3 whose reflectors onlap both onto the down-lapping reflectors of the sedimentary body formed by seismic facies 5 ( Figure 13) on the landward side and on the chaotic reflectors (facies 2a) of unit E u on the seaward side of the sub-basin. The evidence that the unit Du is intercepted by this profile only in the NW sector of the OSFV indicates that this system is much more extended into the OSFV in this sector than that of the SE. Note also the position of the central body of unit Eu (seismic facies 2) above the backbarrier deposits of the previous unit Du (seismic facies 3) and the transition, within unit Eu, between the central body at the NW (seismic facies 2) and the seismic facies 4 in the SE sector.

Unit Gu
Unit Gu is present between depths of ca. −22 m b.s.l. and ca. −4.7 b.s.l., and lies on ES1, unit Fu, unit Du, and unit Eu. The unit Gu infills the last residual depressions of the IV1 and IV2 and cancels them completely.
In the IV1 sector, unit Gu is composed of seismic facies 2b, 3 and 6. Proceeding from the SW to the NE (profile T8; Figure 10), unit Gu shows a sedimentary body consisting of seismic facies 3 which is surmounted by seismic facies 2b and has a maximum thickness of ca. 10 m; seawards, seismic facies 2 transitions into seismic facies 6, which have a maximum thickness of ca. 7 m (profile T8; Figure 10).
In the IV2 sector, unit Gu is composed of seismic facies 2a and 6 and fills the residual accommodation within the sub-basin left in the SE sector in front of Barletta strait after deposition of unit Eu. In some cases, thin beds of seismic facies 2a, overlaid by the draping reflectors of seismic facies 6, occur (profiles T6 and T5; Figure 9). In other cases, only seismic facies 6, which drapes or onlaps onto underlying units and systems, occur (profile T4; Figure 12). A depo-centre of up to ca. 6 m thickness is observed at the turn of the intersection between profiles T4 and PX2 (Figures 11 and  12). The evidence that the unit D u is intercepted by this profile only in the NW sector of the OSFV indicates that this system is much more extended into the OSFV in this sector than that of the SE. Note also the position of the central body of unit E u (seismic facies 2) above the backbarrier deposits of the previous unit D u (seismic facies 3) and the transition, within unit E u , between the central body at the NW (seismic facies 2) and the seismic facies 4 in the SE sector.

Unit G u
Unit G u is present between depths of ca. −22 m b.s.l. and ca. −4.7 b.s.l., and lies on ES1, unit F u , unit D u , and unit E u . The unit G u infills the last residual depressions of the IV1 and IV2 and cancels them completely.
In the IV1 sector, unit G u is composed of seismic facies 2b, 3 and 6. Proceeding from the SW to the NE (profile T8; Figure 10), unit G u shows a sedimentary body consisting of seismic facies 3 which is surmounted by seismic facies 2b and has a maximum thickness of ca. 10 m; seawards, seismic facies 2 transitions into seismic facies 6, which have a maximum thickness of ca. 7 m (profile T8; Figure 10).
In the IV2 sector, unit G u is composed of seismic facies 2a and 6 and fills the residual accommodation within the sub-basin left in the SE sector in front of Barletta strait after deposition of unit E u . In some cases, thin beds of seismic facies 2a, overlaid by the draping reflectors of seismic facies 6, occur (profiles T6 and T5; Figure 9). In other cases, only seismic facies 6, which drapes or onlaps onto underlying units and systems, occur (profile T4; Figure 12). A depo-centre of up to ca. 6 m thickness is observed at the turn of the intersection between profiles T4 and PX2 (Figures 11 and 12).

Erosional Surface ES1
Based on data previously collected in the Manfredonia Gulf [30,32] and in the adjacent areas [36], the ES1 is interpreted as the surface of subaerial exposure formed during the last sea-level low stand attained during the LGM.
The ES1 describes the two incised valleys IV1 and IV2. Regarding the IV1, we believe that this valley is the seaward continuation of "incised valley 2" (Figure 6) from De Santis and Caldara [32]. Our data further support that it is unified valley of the Carapelle and Cervaro streams. Historical maps up to and including that of Marzolla [53] shows the two streams merged more or less at the current position of the Carapelle mouth. After 1851 the lower stretch of the Cervaro stream was diverted to the north, bringing its mouth into Salso Lake, and then the stream was diverted to its current position [54]. For this reason, we refer to this landform as the Cervaro-Carapelle valley (CCV).
The IV2 is an incised valley ( Figure 6) which appears, in the sector closest to the current coastline, as two different incisions. The main incision corresponds to the present Ofanto River mouth. The second incision is placed in correspondence with the probable mouth of an inactive stream; this mouth is masked today by the Margherita di Savoia saltworks. This stream is not indicated with its proper name in recent and historical maps; however, because of its proximity to

Erosional Surface ES1
Based on data previously collected in the Manfredonia Gulf [30,32] and in the adjacent areas [36], the ES1 is interpreted as the surface of subaerial exposure formed during the last sea-level low stand attained during the LGM.
The ES1 describes the two incised valleys IV1 and IV2. Regarding the IV1, we believe that this valley is the seaward continuation of "incised valley 2" (Figure 6) from De Santis and Caldara [32]. Our data further support that it is unified valley of the Carapelle and Cervaro streams. Historical maps up to and including that of Marzolla [53] shows the two streams merged more or less at the current position of the Carapelle mouth. After 1851 the lower stretch of the Cervaro stream was diverted to the north, bringing its mouth into Salso Lake, and then the stream was diverted to its current position [54]. For this reason, we refer to this landform as the Cervaro-Carapelle valley (CCV).
The IV2 is an incised valley ( Figure 6) which appears, in the sector closest to the current coastline, as two different incisions. The main incision corresponds to the present Ofanto River mouth. The second incision is placed in correspondence with the probable mouth of an inactive stream; this mouth is masked today by the Margherita di Savoia saltworks. This stream is not indicated with its proper name in recent and historical maps; however, because of its proximity to the town of San Ferdinando di Puglia, we named it San Ferdinando stream. As a consequence, we indicated the IV2 as Ofanto and San Ferdinando valley (OSFV), where the San Ferdinando stream is a tributary of the Ofanto river.

Borehole ZS2 and Vibrocores SP1_VC24 and SP1_VC52
Since borehole ZS2 intercepts the transgressive filling of the CCV (Figures 6), it allows the reconstruction of its sedimentary palaeoenvironments. From 30 to 29 m, the environment is continental and characterised by the presence of stagnant freshwater. From 29 to 27.5 m, a euryhaline and eurytherm lagoon (LEE, sensu [55]) is established, which was initially confined and then became more open. Lagoon sedimentation was interrupted by a pyroclastic episode. Starting from 27.5 m, there is a rapid marine transgression characterised by a strong sedimentary instability typical of the heterogeneous palaeo-community [56]. This palaeo-community developed at the passage between the infralittoral zone and the terrigenous mud biocoenosis (VTC, sensu [55]) of the circalittoral zone. After a transition phase (from 12.45 to 8.0 m), starting at 8.0 m, a regressive phase occurred characterised by fine well-sorted sand biocoenosis (SFBC, sensu [55]) with some beach-layers. Finally, the succession is closed by an anthropic deposit.
The borehole ZS2 also provides two important chronological constraints for the transgressive filling of CCV and allows for an assessment of the subsidence rate of the area during the Holocene.
The chronological constraints are: (i) sample ZS2 28.65, located at −27.45 m b.s.l., consisting of lagoonal environment clays rich in Cerastoderma glaucum, which provided an age of 9160 ± 137 cal year BP; and (ii) the level of coarse pumice, whose base is located at −27.30 m b.s.l., which is attributed to the Mercato pumice eruption that occurred at 8890 ± 90 cal year BP [51,52].
Above the pumice level, the lagoonal clays rich in Cerastoderma glaucum reoccur. When the Mercato eruption occurred, local sea level was ca.

Borehole ZS2 and Vibrocores SP1_VC24 and SP1_VC52
Since borehole ZS2 intercepts the transgressive filling of the CCV (Figure 6), it allows the reconstruction of its sedimentary palaeoenvironments. From 30 to 29 m, the environment is continental and characterised by the presence of stagnant freshwater. From 29 to 27.5 m, a euryhaline and eurytherm lagoon (LEE, sensu [55]) is established, which was initially confined and then became more open. Lagoon sedimentation was interrupted by a pyroclastic episode. Starting from 27.5 m, there is a rapid marine transgression characterised by a strong sedimentary instability typical of the heterogeneous palaeo-community [56]. This palaeo-community developed at the passage between the infralittoral zone and the terrigenous mud biocoenosis (VTC, sensu [55]) of the circalittoral zone. After a transition phase (from 12.45 to 8.0 m), starting at 8.0 m, a regressive phase occurred characterised by fine well-sorted sand biocoenosis (SFBC, sensu [55]) with some beach-layers. Finally, the succession is closed by an anthropic deposit.
The borehole ZS2 also provides two important chronological constraints for the transgressive filling of CCV and allows for an assessment of the subsidence rate of the area during the Holocene.
The chronological constraints are: (i) sample ZS2 28.65, located at −27.45 m b.s.l., consisting of lagoonal environment clays rich in Cerastoderma glaucum, which provided an age of 9160 ± 137 cal year BP; and (ii) the level of coarse pumice, whose base is located at −27.30 m b.s.l., which is attributed to the Mercato pumice eruption that occurred at 8890 ± 90 cal year BP [51,52].
Above the pumice level, the lagoonal clays rich in Cerastoderma glaucum reoccur. When the Mercato eruption occurred, local sea level was ca. −24.3 m b.s.l. [47]. If we consider the base of pumice level The subsidence value of one meter was used to correct the palaeo sea-levels obtained from the current position of the base and the top of the bodies formed by seismic facies 2.
The vibro-cores SP1_VC24 and SP1_VC25 show a predominantly sandy deposit rich in coarse bioclasts, with a minor content of silt and/or clay, which we interpret as a high energy beach deposited at (or close to) sea level (beach/shoreface).

Interpretation of Seismic Facies
Facies 1. The characteristics of this seismic facies, its position, and the similarity with other seismic facies of fluvial origin [40,42] suggest a fluvial lag which was most likely deposited under high-energy conditions. The deposition of this facies occurred when streams incised the substrate and developed during both the low-stand and transgression until the arrival of the sea.
Facies 2. The interpretation of the seismic facies 2 was based also on the vibro-cores SP1_VC24 and SP1_VC25 (Figure 2), which intercept only this facies from the sea bottom to borehole bottom ( Figure 14). Overall, the data allows us to interpret seismic facies 2a as a predominantly sandy deposit rich in coarse bioclasts, with a minor content of silt and/or clay. Based on this evidence, we interpret this facies as the expression of high-energy coastal sediments in a beach/shoreface environment. This interpretation is in accordance with Trincardi et al. [28], who defined the transgressive sediments cropping out at the sea bottom in the area where the two vibro-cores are located as "proximal transgressive sandy deposits of beach or paralic environment". Facies 2b is interpreted as a beach/shoreface deposit within a transgressive barrier or spit that migrated landward by wave erosion and over-washing. In fact, landward-dipping reflectors are dominant in transgressive barriers, and they can be interpreted as former flood-tidal deltas and wash-over fans, which generally overlie sandy-muddy deposits that formed in a backbarrier environment [61,62]. In conclusion, seismic facies 3 represents deposits that formed following arrival of the sea and wave action: this wave action reworked ES1 and created a coastal high-energy palimpsest deposit due to mixing among reworked sediments from the substrate with coeval sediments. Facies 3. The position of facies 3 with respect to facies 2 is a fundamental element for its interpretation: the transition between facies 2 and 3 is interpreted as a passage from a beach/barrier/spit environment to a low-energy environment; therefore, facies 3 is interpreted as a backbarrier deposit (perhaps a semi-enclosed bay/microtidal estuarine lagoon). The cases where facies 2 lies above facies 3 are interpreted as resulting from the migration of coastal environments landwards, when the beach passed over the backbarrier during progression of the marine transgression.
Facies 4. As in the case of facies 3, the transition between facies 2 and 4 is interpreted as a passage from a beach/barrier/spit environment to a lower energy environment. In addition, the position of facies 4 in the area of the IV2, near the opening towards the sea, allows us to interpret this seismic facies as a deposit of open backbarrier (perhaps bay/ microtidal estuarine lagoon), probably of silty nature. The proximity to the opening implies the presence of currents and channels, which create a moderate energy environment. Thus, the wavy undulations which characterise this facies can be interpreted in some cases as channel features, in other cases as deformation of the estuarine/lagoon bottom due to waves or currents.
Facies 5. This seismic facies is interpreted as a deltaic body. Thus, the body formed by seismic facies 5 within unit E u represents a small delta (D1; Figure 13) up to ca. 10 m thick. Facies 6. This seismic facies is interpreted as a lower shoreface/offshore, probably silty-clay deposit that lies around or below the wave action base. It is formed due to the rapid landward shift in the shoreline, which was in turn due to the low-gradient substrate.

Interpretation of UBSUs and Transgressive Architecture (Sea Level Ranges and Deposition of UBSUs)
We have linked the UBSUs at four sea-level ranges and therefore at four corresponding age ranges. The sea-levels which define the ranges are: −39 m, −29 m, −18 m, −4.7 m. The first two values are derived from upper and/or lower limit of bodies formed by seismic facies 2, corrected considering a subsidence of 1 m. The level −18 m is derived from the small coastal deltaic body (D1; seismic facies 5) visible within the unit Eu, by assuming its rollover point between the fore-set and the topset as a proxy of sea level (according to [30,45,63]), after a correction to remove the subsidence of 1m. The level −4.7 m is that coinciding with the maximum landward shift of the shoreline (ca. 5.5 ka BP; [47] and marks the end of the late-Pleistocene to Holocene transgression [30,[33][34][35] This phase coincides with a sea-level rise up to ca. −39 m b.s.l., i.e., until ca. 9.8 ka BP [47]. In this phase, the sea entered the OSFV occupying only the lowest reaches, corresponding to the Barletta strait and the northernmost corner (intersection between profiles P3 and T5; Figure 6). In this phase, unit Au, unit Bu, and unit Cu were deposited.
Unit Au contains only seismic facies 1, corresponding to a fluvial lag. In the OSFV, this unit extends both below and above −39 m b.s.l., while in the CCV, it extends only above the −39 m b.s.l. Thus, it was deposited in a subaerial environment on the bottom of the two valleys until the arrival of the sea.
Unit Bu consists exclusively of seismic facies 2a and is interpreted as a beach deposit formed on the shelf above the ES1 in conjunction with a sea-level rise of up to −39 m.
Unit Cu consists exclusively of seismic facies 3, and its deposition is correlated to the final phases of sea-level rise up to −39 m. We interpret this unit as an enclosed bay/lagoon deposit that formed when the sea began to enter the OSFV through the Barletta strait, without having yet submerged the IV2 threshold.

Interpretation of UBSUs and Transgressive Architecture (Sea Level Ranges and Deposition of UBSUs)
We have linked the UBSUs at four sea-level ranges and therefore at four corresponding age ranges. The sea-levels which define the ranges are: −39 m, −29 m, −18 m, −4.7 m. The first two values are derived from upper and/or lower limit of bodies formed by seismic facies 2, corrected considering a subsidence of 1 m. The level −18 m is derived from the small coastal deltaic body (D1; seismic facies 5) visible within the unit E u , by assuming its rollover point between the fore-set and the topset as a proxy of sea level (according to [30,45,63]), after a correction to remove the subsidence of 1m. The level −4.7 m is that coinciding with the maximum landward shift of the shoreline (ca. 5.5 ka BP; [47] and marks the end of the late-Pleistocene to Holocene transgression [30,[33][34][35] This phase coincides with a sea-level rise up to ca. −39 m b.s.l., i.e., until ca. 9.8 ka BP [47]. In this phase, the sea entered the OSFV occupying only the lowest reaches, corresponding to the Barletta strait and the northernmost corner (intersection between profiles P3 and T5; Figure 6). In this phase, unit A u , unit B u , and unit C u were deposited.
Unit A u contains only seismic facies 1, corresponding to a fluvial lag. In the OSFV, this unit extends both below and above −39 m b.s.l., while in the CCV, it extends only above the −39 m b.s.l. Thus, it was deposited in a subaerial environment on the bottom of the two valleys until the arrival of the sea.
Unit B u consists exclusively of seismic facies 2a and is interpreted as a beach deposit formed on the shelf above the ES1 in conjunction with a sea-level rise of up to −39 m. Unit C u consists exclusively of seismic facies 3, and its deposition is correlated to the final phases of sea-level rise up to −39 m. We interpret this unit as an enclosed bay/lagoon deposit that formed when the sea began to enter the OSFV through the Barletta strait, without having yet submerged the IV2 threshold. This range coincides with a sea-level rise from ca. −39 to ca. −29 m b.s.l., i.e., from ca. 9.8 ka BP to ca. 9.2 ka BP [47].
In CCV sector, this phase involves the deposition of part of unit A u and unit D u . Unit A u maintains the same characteristics as those described in the previous phase. Unit D u is constituted only of seismic facies 2a; thus, unit D u is interpreted as a simple beach deposit of high-energy environment deposited on the ES1.
In the OSFV sector, this phase involves the deposition of unit D u ; the OSFV threshold is not jet submerged in the NW sector of the valley (profile T6), while it is fully submerged in the central and south-eastern sector (profiles T5 and T4; Figures 9 and 12). Therefore, unit D u is confined to the outside of the valley in the first sector, while it centres on the threshold in the central and south-eastern sector. Considering the overall geometry of unit D u , as well as the seismic facies that constitute it, we believe that it consists of a central body composed of a beach/spit that, into the OSFV, transitions into backbarrier deposits (perhaps semi-enclosed bay/microtidal estuarine lagoon) and towards the shelf into the lower shoreface/offshore deposits. In particular, the presence of a spit is inferred from the two sets of opposing inclined reflectors that characterise the central body consisting of seismic facies 2a (Figure 9), from the progressive narrowing of the central body proceeding from the NW to the SE (compares profiles T5 and T4; Figures 9 and 12), and from its terminations at the correspondence of the Barletta strait (profile PX3). This range coincides with sea-level rising from ca. −29 m to ca. −18 m b.s.l., which occurred from ca. 9.2 ka BP to ca. 8.3 ka BP [47].

Sea-Level
In the CCV sector, deposition of unit E u only occurred. We interpret the infilling of the landward part by unit E u (profile T8; Figure 10) as a transgressive barrier/spit that migrated landward and covered the backbarrier deposits [64,65]. This system formed at the unified mouth of the Carapelle and Cervaro streams.
Passing to OSFV, deposition of units E u and F u occurred. In the NW sector (profile T6) rising sea levels submerged the threshold only in this phase, albeit slightly. The small delta D1 recognisable within unit E u is located seaward of the probable position of the San Ferdinando stream mouth before it was covered by the Margherita di Savoia saltworks. Thus, we named this feature the San Ferdinando delta (SFD).
The deposition of unit E u involves most of the OSFV infilling and its division into two residual sub-basins, the most northwest of which, placed in front of SFD, is rapidly filled by unit F u in the last part of this same phase, after the SFD formation ( Figure 13). Considering the overall geometry of unit E u , as well as the seismic facies that constitute it, we interpret unit E u as a complex coastal system consisting of lower shoreface/offshore deposits that pass landwards into to beach/spit deposits, then to backbarrier deposits (perhaps a semi-enclosed bay/ microtidal estuarine lagoon) in which a small fluvial delta flow. In the CCV sector, based on the seismic facies and their mutual positions, unit G u is interpreted as a transgressive spit-barrier sandy deposit that overlies backbarrier deposits (perhaps a semi-enclosed bay/microtidal estuarine lagoon), which passes seaward into the lower shoreface/offshore deposits.
In the OSFV sector, based on the seismic facies and their mutual positions, unit G u is interpreted as a sandy beach deposit of high energy at the bottom (seismic facies 2a) dating back only to the earliest moments of the phase. Then, the beach deposit was quickly replaced by lower shoreface/offshore deposits (seismic facies 6), which are probably silty-clayey, dating back to when the shoreline transgressed rapidly across the gentle slope of the shelf.

Focus on Beach/Spit-Backbarrier Systems in CCV and OSFV Infilling
In this section, we propose models for the formation and evolution dynamics of beach/spit-backbarrier systems in the infilling of CCV and OSFV.

CCV Infilling
We interpret the landward sector of CCV as the superimposition of two successive barrier/spit-backbarrier systems at the unified mouth of the Carapelle-Cervaro streams, with the lower system belonging to unit E u , and the upper systems belonging to unit G u ( Figure 10). In both systems, as the transgression progressed, the barrier/spit retreated landwards and covered the backbarrier sediments through over-wash and/or tidal inlet deposits, as evidenced by the reflectors in seismic facies 2b down-lapping landwards onto the plane-parallel reflector of seismic facies 3 ( Figure 10).
The ZS2 borehole (Figure 7) intercepts the infilling of CCV (Figures 6 and 15) and allows the correlation between sparker profile T8 and CHIRP profile T1, and a chronological validation of the age range attributed to unit E u . In particular, in profile T1 we recognise three seismic units correlated with units A u , E u , G u ( Figure 15B).
The lagoonal deposit with the enclosed Mercato pumice (dated at 9160 ± 137 and 8890 ± 90 cal year BP) in ZS2 borehole ( Figure 15B) is correlated with a series of very high-amplitude reflectors visible roughly at the same depth at the base of unit E u in the CCV sector shown by profile T1 (Figure 15A,B) and, in terms of both depth and sedimentary environment, with seismic facies 3 visible at the base of unit E u in the near profile T8 ( Figure 15C,D). This confirms the age range attributed to unit E u according to our method of age range evaluation (from ca. 9.2 ka BP to ca. 8.3 ka BP) which, in addition, is within the time of sapropel S1a formation in Mediterranean [66,67].
About the presence of two barrier/spit-backbarrier systems that characterise the inner part of the CCV (units E u and G u ), we hypothesise that it is due to the abrupt sea-level rise that occurred during the so-called 8.2 event [68][69][70]. This sea-level jump may have drowned the lower barrier/spit-backbarrier system (unit E u ); subsequently, the conditions were appropriate for construction of a new system over the previous one (unit G u ). This hypothesis is consistent with a scheme in which unit E u is attributed to the sapropel S1a event, and unit G u can be attributed to sapropel S1b, as these units are separated by a discontinuity dating back to the 8.2 event [67].

OSFV Infilling
With rising sea level, the inherited morphology of the OSFV favors the creation of beach/spit-backbarrier systems associated with the Ofanto river mouth. The first evidence of the creation of a real beach/spit-backbarrier system (unit D u ) occurs during sea-level range 2 ( Figure 16). Sea-level rise flooded the threshold separating the OSFV from the shelf. The threshold appears to play an important role in triggering the formation of the system. In fact, the first spit that forms (within unit D u ) is perched on the bedrock high of the threshold (Figures 9 and 12). The presence of the high provided material for the construction of the spit (as marine erosion affected the high) and promoted the deposition of material transported by the littoral drift. This feature appears to be like other barriers, such as barrier system of the Miquelon-Langlade island [71]. The landward expansion of the unit D u was enabled by the initial infilling of the OSFV (unit C u ).
Water 2020, 12, x FOR PEER REVIEW 23 of 32 relatively high energy due to the presence of channels, currents and/or waves in the sector closest to the opening to the sea.  The lagoonal level containing the Mercato pumice in borehole ZS2 has been correlated both to a series of very high-amplitude reflectors placed roughly at the same depth in the nearby profile T1 at the base of unit E u and to seismic facies 3 reflectors placed roughly at the same depth at the base of unit E u in profile T8.
The accretion of the unit D u spit left an opening in the Barletta strait (profile PX3). This landform, inherited from the time of the OSFV incision, is only modified by spit accretion but remains as an opening between the OSFV and the open sea. This opening influenced the different hydrodynamics that affect the backbarrier deposits. In the NW sector (profile T5, Figure 9), these deposits are composed of seismic facies 3, which indicates a low-energy environment. In the SE sector (profile T4, Figure 12), these deposits are composed of the wavy facies 4, which indicate relatively high energy due to the presence of channels, currents and/or waves in the sector closest to the opening to the sea.
During sea-level range 3, a general landward shift in the beach/spit-backbarrier system (unit E u ) occurs, together with an increase in sediment supply, resulting in a greater thickness and extension of the unit E u (Figure 16). Compared to the sea level range 2, the different hydrodynamics within the backbarrier is conserved between the NW and the SE sectors, with the same characteristics and same causes. In its final stages, this phase resulted in the near-complete infilling of the OSFV.
The small delta (SFD) that formed during sea level range 3 at the inner edge of the OSFV within unit E u is important to test our age model; in fact, the sea level of −18 m, reconstructed from the rollover point of SFD, corresponds to an age of ca. 8.3 ka BP [47], that is, within the formation period of sapropel S1a [66,72]. This data supports the attribution of the unit E u to the sapropel S1a, as in the CCV sector. During sea-level range 3, a general landward shift in the beach/spit-backbarrier system (unit Eu) occurs, together with an increase in sediment supply, resulting in a greater thickness and extension of the unit Eu ( Figure 16). Compared to the sea level range 2, the different hydrodynamics within the backbarrier is conserved between the NW and the SE sectors, with the same characteristics and same causes. In its final stages, this phase resulted in the near-complete infilling of the OSFV.
The small delta (SFD) that formed during sea level range 3 at the inner edge of the OSFV within unit Eu is important to test our age model; in fact, the sea level of −18 m, reconstructed from the rollover point of SFD, corresponds to an age of ca. 8.3 ka BP [47], that is, within the formation period of sapropel S1a [66,72]. This data supports the attribution of the unit Eu to the sapropel S1a, as in the CCV sector.

Palaeoenvironmental and Palaeoclimatic Considerations
The overall characteristics of the OSFV infilling suggest sediment movement proceeded from the NW to the SE. All the units decrease in size in this direction, and the OSFV infilled is faster in the NW sector than in the SE sector. In particular, the construction of a sandy spit points to a littoral drift from the N-NW to the E-SE; that is, in the opposite direction of the current littoral drift in the Manfredonia Gulf.
The direction of sediment transport reconstructed in this work suggests a scenario such as that described by De Santis and Caldara [73] for the period preceding the Holocene climatic transition (ca. 5.5-4.5 cal ka BP). For that period, the authors recognised a strong southward and

Palaeoenvironmental and Palaeoclimatic Considerations
The overall characteristics of the OSFV infilling suggest sediment movement proceeded from the NW to the SE. All the units decrease in size in this direction, and the OSFV infilled is faster in the NW sector than in the SE sector. In particular, the construction of a sandy spit points to a littoral drift from the N-NW to the E-SE; that is, in the opposite direction of the current littoral drift in the Manfredonia Gulf.
The direction of sediment transport reconstructed in this work suggests a scenario such as that described by De Santis and Caldara [73] for the period preceding the Holocene climatic transition (ca. 5.5-4.5 cal ka BP). For that period, the authors recognised a strong southward and south-eastward littoral drift in the Manfredonia Gulf, which was promoted by the prevalence and dominance of the N, E and NE winds, in turn promoted by more frequent cyclogenesis to the E, S and SE of the Italian Peninsula. Therefore, this littoral drift pattern in the Manfredonia Gulf seems to have been active from at least ca. 9.8 ka BP, when the first beach/spit-backbarrier system started to form (unit D u ).
We consider units D u , E u , and G u as beach/spit-backbarrier systems that were arranged in a landward backstepping pattern, which involved partial preservation of transgressive deposits or spit morphology. According to the literature, the preservation of relict barrier within a backstepping sequence is rare [12,21] and has been attributed to combinations of rapid sea-level rise, such as those associated with meltwater pulses [16], and accompanied by particularly favourable circumstances, such as early cementation (beach-rock and aeolianite formation) or gravelly sediments [12,14,17].
In our case, the barrier system overstepping occurred after MWP-1B and is therefore in a regime of sea-level rise slower than the previous period [13,74]. Furthermore, the overstepping affected loose and mostly sandy sediments, as demonstrated by the vibro-core collected in this study and by the previous studies [28,30,31]. Therefore, to justify the partial preservation of two barrier systems arranged in a landward backstepping pattern, we believe that a strong and continued supply of sediments attenuated or compensated for the erosive actions. In this context, subsequent units or systems may partially cover and protect previous systems ( Figure 17A-C).
Water 2020, 12, x FOR PEER REVIEW 25 of 32 south-eastward littoral drift in the Manfredonia Gulf, which was promoted by the prevalence and dominance of the N, E and NE winds, in turn promoted by more frequent cyclogenesis to the E, S and SE of the Italian Peninsula. Therefore, this littoral drift pattern in the Manfredonia Gulf seems to have been active from at least ca. 9.8 ka BP, when the first beach/spit-backbarrier system started to form (unit Du). We consider units Du, Eu, and Gu as beach/spit-backbarrier systems that were arranged in a landward backstepping pattern, which involved partial preservation of transgressive deposits or spit morphology. According to the literature, the preservation of relict barrier within a backstepping sequence is rare [12,21] and has been attributed to combinations of rapid sea-level rise, such as those associated with meltwater pulses [16], and accompanied by particularly favourable circumstances, such as early cementation (beach-rock and aeolianite formation) or gravelly sediments [12,14,17].
In our case, the barrier system overstepping occurred after MWP-1B and is therefore in a regime of sea-level rise slower than the previous period [13,74]. Furthermore, the overstepping affected loose and mostly sandy sediments, as demonstrated by the vibro-core collected in this study and by the previous studies [28,30,31]. Therefore, to justify the partial preservation of two barrier systems arranged in a landward backstepping pattern, we believe that a strong and continued supply of sediments attenuated or compensated for the erosive actions. In this context, subsequent units or systems may partially cover and protect previous systems ( Figure 17A-C). We correlate this strong and continued sediment supply with the sapropel S1 event. In fact, the chronological constraint provided by the borehole ZS2 and the sea level reconstructed thanks to SFD allows the consideration of unit Eu as coeval with formation of sapropel S1a. As this conclusion confirms and validates our age range evaluation method based on seismic facies 2, it appears that We correlate this strong and continued sediment supply with the sapropel S1 event. In fact, the chronological constraint provided by the borehole ZS2 and the sea level reconstructed thanks to SFD allows the consideration of unit E u as coeval with formation of sapropel S1a. As this conclusion confirms and validates our age range evaluation method based on seismic facies 2, it appears that unit D u and, probably, at least part of unit G u , were also formed simultaneously with the deposition of sapropel S1 in the Mediterranean. Thus, units D u , E u , and G u represent the shallow-water equivalent of the cm-thick sapropel layers that accumulated offshore in the deeper southern Adriatic basin.
During the minima of the precession of the equinoxes, the northward shift in the African monsoon led to increased precipitation and river runoff in northern African and eastern European regions, which favoured water column stratification and reduced deep-water ventilation of the eastern Mediterranean Sea and Adriatic. This allowed for the accumulation and preservation of organic matter in deep marine sediments and the formation of sapropels [77][78][79]. The youngest of such events, sapropel S1, occurred from ca. 10/9.0 to ca. 6.8 cal ka BP [66,67] and formed in two distinct pulses (S1a and S1b), which are separated by an interruption in sapropel formation, dating to ca. 8.3-7.8 cal ka BP in the Adriatic Sea [67] due to a short-lived episode of cooler and drier conditions [66].
Overall, units D u , and E u demonstrate a period of increased sediment influx to the sea, which has promoted the formation of complex coastal systems. In particular, the formation of the small delta within unit E u shows a period of increased precipitation during the S1a sapropel event in this area of southern Italy [80] and not only in the northern African and eastern European regions.
Therefore, the roles of the S1 event in the creation and partial preservation of our complex beach/spit-backbarrier systems (units D u , E u , and G u ) seem to be as follows: i) a direct role in the formation of the beach/spit-backbarrier systems of units D u , E u , and G u , due to the increased sediment supply; and ii) an indirect role in the partial conservation of previous beach/spit-backbarrier systems because the same increased sediment supply that created a subsequent system also allowed partial protection and conservation of a previous system from wave ravinement action.
The "in-place drowning" transgressive model for low-gradient shelves [75,76], although questioned later [11], predicts that a barrier island-lagoon system develops during conditions of slow relative sea-level rise (Figure 17A'); the barriers can be drowned and preserved due to a very rapid pulse of relative sea-level rise, such as the MWP-1A and MWP-1B, which are recognized in most of the sea-level curves. Later, a new barrier forms landwards during newly created conditions of slower relative sea-level rise ( Figure 17B'). This is the case, for example, in the barriers systems of the northern Adriatic [13], a low-gradient setting [81]; due to the episodic occurrence of very high rates of sea-level rise up to 60 mm/a [74], some barriers were overstepped and drowned. As a result, two barriers can be several kilometres apart [81]. A similar case is that described by Cooper et al. [15] in Tijucas Bay (southern Brazil), where a former sandy shoreface-barrier was overstepped during the 8.2 ka event and decoupled from the contemporary shoreline. Upon overstepping, the wave-influenced shoreline was shifted from the barrier to the landward margin of the former backbarrier (over 7 km cross-shore distance).
Our study instead suggests that, in situations of high and continued sedimentary inputs, low-gradient settings and sea-level rise lower than those recorded during MWP-1A and MWP-1B, a stacking pattern of the barrier systems that we define as "continuous" backstepping may be generated, in which the next barrier system partially covers the previous one ( Figure 17A-C). Continued sediment supply is a key factor in our model. In an extreme situation, a very strong and continuous sediment supply, along with a slow sea-level rise, can generate a prograding barrier system [82]. Instead, in our case, a fine balance between SLR and sediment supply causes a continuous landward backstepping.

Conclusions
Holocene sea-level rise has affected most of the LGM surface (ES1) in our study area between the last phases of the MWP−1B [13,74,83,84] and the subsequent slower rise until the maximum landward shift in the shoreline occurred. This implies that in our study area, the upper part of the TST was deposited (uTST [36]), as well as the subsequent HST, from 5.5 ka to present day.
In the innermost area of the Manfredonia Gulf, the post-Würm transgression resulted in the backfilling of two erosive landforms incised within the ES1 which we documented for the first time.
The first landform (CCV), located further to the NW, is correlated with the current Carapelle and Cervaro streams and consists of a meandering incised valley with stretches parallel to the current coastline in the most landward sector and a straighter route on the shelf, where it flows as a tributary into the so-called MIV.
The second erosive landform (OSFV), which is present farther to the SE, is correlated with the current San Ferdinando stream and Ofanto River. The OSFV consists of two separate incisions that join to form a very wide and shallow valley.
Six different seismic facies and seven UBSUs have been identified, which have covered the ES1 and/or infilled the CCV and OSFV. The UBSUs in some cases are made up of a single seismic facies, in other cases of multiple seismic facies that transit into one another.
We have identified some approximated age ranges for the different units, starting from the current position of one of the seismic facies identified, seismic facies 2. It has been interpreted as deposited at (or close to) contemporary sea level. Since seismic facies 2 forms several sedimentary bodies (belonging to different units) in a landward backstepping pattern on ES1, we associated two sea levels with the base and top of each sedimentary body. Then, we considered the local sea-level curve to infer two ages corresponding to these two sea levels. Thus, several sedimentary bodies formed by seismic facies 2 can be correlated to many sea-level ranges and therefore to age ranges.
Both in the CCV and OSFV, some barrier/spit-backbarrier systems, arranged in a landward backstepping pattern, were formed. One of these (unit E u ) was attributed to the sea level range between −29 to −18 m b.s.l. and, thus, to a time interval between 9.2 ka BP and ca. 8.3 ka BP. This age range attribution, obtained from the current position of sedimentary bodies formed by seismic facies 2 included in unit E u , was chronologically constrained both in CCV and OSFV.
In CCV the chronological constraint was obtained thanks to a dated lagoonal deposit with enclosed Mercato pumice present in the ZS2 borehole, which intercepts unit E u . In OSFV, the chronological constraint was obtained thanks to a small delta included in unit E u , in that the rollover point of the delta was assumed as a sea-level indicator, and thus as an age indicator.
The chronological constraint identified for unit E u validated the age range evaluation method based on the current position of sedimentary bodies formed by seismic facies 2. Thus, it results that most of the OSFV and CCV infilling is coeval with the formation of the sapropel S1a in the Mediterranean.
The role of the S1 event in the creation and partial preservation of our complex beach/spit-backbarrier systems seem to be both direct and indirect: direct in the formation of the systems, due to the increased sediment supply; indirect in the partial conservation of previous systems because the same increased sediment supply, coupled with the low gradient setting and slow relative sea-level rise, promoted what we call a "continuous" backstepping. This latter, in turn allowed the partial protection and conservation of the previous barrier system from wave ravinement action.
The progression of the infilling units into CCV and OSFV indicates a prevalent littoral drift from the N-NW to the E-SE; that is, in a direction contrary to the current littoral drift of the Manfredonia Gulf. This evidence allows us to extend the climatic phase of frequent cyclogenesis to the E, S and SE of the Italian Peninsula, identified by De Santis and Caldara [73] before the middle-late Holocene transition (5.5-4.5 ka cal BP), backwards in time to at least ca. 9.8 ka BP.