Changing environments and human interaction during the Pleistocene–Early Holocene from the shallow coastal area of Dor, Israel

Abstract

The protected Tel-Dor coastal embayment in the eastern Mediterranean preserves an unusually complete stratigraphic record that reveals human–environmental interactions throughout the Holocene. Interpretation of new seismic profiles collected from shallow marine geophysical transects across the bay show five seismic units were correlated with stratigraphy and age dates obtained from coastal and shallow-marine sediment cores. This stratigraphic framework permits a detailed reconstruction of the coastal system over the last ca. 77 ka as well as an assessment of environmental factors that influenced some dimensions of past coastal societies. The base of the boreholes records lowstand aeolian deposits overlain by wetland sediments that were subsequently flooded by the mid-Holocene transgression. The earliest human settlements are submerged Pottery Neolithic (8.25–7 ka) structures and tools, found immediately above the wetland deposits landward of a submerged aeolianite ridge at the mouth of the bay. The wetland deposits and Pottery Neolithic settlement remains are buried by coastal sand that records a middle Holocene sea-level rise ca. 7.6–6.5 ka. Stratigraphic and geographic relationships suggest that these coastal communities were displaced by sea-level transgression. These findings demonstrate how robust integration of different data sets can be used to reconstruct the geomorphic evolution of coastal settings as well as provide an important addition to the nature of human–landscape interaction and cultural development.

INTRODUCTION

The Carmel coast of Israel is one of the main arteries that facilitated human movement from Africa to southwest Asia since the early Pleistocene (Olami and Ronen 1977; Zeder, 2008). Throughout prehistoric and historic periods, the importance of this corridor has been marked by some of the most famous Near Eastern sites (Rebollo et al., 2011; Meignen et al., 2017) that contain evidence of sociocultural change. Previously investigated sites located in the Carmel coastal zone, coastal plain, and adjacent Carmel ridge (Fig. 1) have established a detailed understanding of sociocultural evolution by identifying transitions from hunting and gathering nomadic societies (Geometric–Kebaran; Bar-Yosef, 2002) to semi-sedentary cultures (Natufian; Bar-Yosef, 1998) to settled farming communities (Pre-Pottery Neolithic–Neolithic; Kuijt and Goring-Morris, 2002; Goring-Morris and Belfer-Cohen, 2011) to the transition to complex urban societies (Bronze and Iron ages; Milevski, 2013) and later classical and recent historical periods. This makes the area a nexus for understanding cultural developments in the eastern Mediterranean (Naveh and Carmel, 2004; Nitschke et al., 2011; Galili et al., 2019).

Figure 1. Locality maps. (a) Direction of longshore transport along the eastern Mediterranean shown by black arrows, and marine and terrestrial paleoclimate proxy locations are shown by red stars. Dashed blue line represents the Nile; dashed red lines represent approximate international borders. Outlined rectangle shows location of inset (b), which illustrates surface lithology of Israel's Carmel coastal plain, stream location, and shelf bathymetry. Red and green triangles indicate location of Neolithic habitations (Galili et al., 2020). Numbered purple circles represent previous stratigraphic studies: (1) Kadosh et al. (2004); (2) Sivan et al. (2004a); (3) Cohen-Seffer et al. (2005); (4) Sivan et al. (2011); (5) Shtienberg et al. (2017); (6) Goodman-Tchernov et al. (2009); (7) Neev et al. (1978). Red dashed rectangle indicates locations of shallow marine geophysical surveys conducted by Shtienberg et al. (2016) and Goff et al. (2018).

Coastal environments protected by rocky headlands, embayment and submerged ridges along the Carmel littoral zone preserve key paleoenvironmental archives of conditions associated with human settlement. The subdued marine flow energies and protection from aeolian processes in these settings efficiently trap sediments and their associated faunal, floral, and isotopic signatures of late Pleistocene and Holocene paleoenvironments (Sneh, 1981; Sneh and Klein, 1984; Schattner et al., 2010; Shtienberg et al., 2016). These sedimentary records include evidence for climate change (Kadosh et al., 2004), sediment accumulation rates (Sivan et al., 2004a, 2011), hydrological events (Avnaim-Katav et al., 2017), and catastrophic tsunami events (Reinhardt et al., 1999, 2006; Tyuleneva et al., 2018; Shtienberg et al., 2020). Underwater excavations along the Carmel coast have identified 19 in-situ Stone Age settlements ranging from the Pre-Pottery Neolithic C (9–8.4 ka) to Pottery Neolithic (7.7–6.8 ka) periods (Fig. 1b; Galili et al., 2020). Investigations of the ancient constructions and remains found in these preserved underwater sites have revealed some of the earliest evidence of water-well construction (Galili and Nir, 1993), olive oil production (Galili and Sharvit, 1994–1995; Galili et al., 1997; 2018), table olive preparation (Galili et al., 2021), early fishing (Zohar et al., 2001), early developments in plant and animal domestication (Galili et al., 2002), as well as burial practices (Galili et al., 2005). However, the precise timing of paleoenvironmental change over the Holocene is still missing.

Here, we investigate the stratigraphy adjacent to Tel-Dor (Fig. 1), a key ancient urban center on the Carmel coast where an embayment (Figs. 2, 3) has protected the stratigraphic record sufficiently to provide a high-quality, integrated paleoenvironmental and archaeological record of the coastal plain. Habitation at Tel-Dor mound site began in the Middle Bronze Age II (3.9–3.5 ka), and continued to grow and develop throughout the Roman Period (2–1.7 ka), Byzantine (1.7–1.4 ka), and Crusader (0.9–0.75 ka) periods (e.g., Stern, 1994; Raban, 1995; Gilboa and Sharon, 2008; Nitschke et al., 2011). Recent excavations conducted in the south bay adjacent to the Tel identified an Iron Age II harbor (Fig. 3b; Arkin Shalev et al., 2019) and submerged Neolithic structures (Fig. 3b; Shtienberg et al., 2020), pointing towards Dor's complex and long pattern of human habitation, which is a focal point for our investigation.

Figure 2. Aeolianite ridge geometry and surface lithology from the coastal zone of Dor: (a) aerial view; (b) surface lithology map. The dashed red rectangle representing the approximate aeolianite ridge extent in the area outlines headlands, islands, and a wave-cut platform that are found in close proximity to the current shoreline. Purple stars mark the location of luminescence ages from Sivan and Porat (2004) and Mauz et al. (2013). Open rectangles mark approximate boundaries for close-up maps found in Figure 3.

Figure 3. (a, b) Inset maps (see Figure 2 for inset locations) representing surveying and sampling locations as well as geomorphological characteristics of the south bay of Dor and its pocket beach. (a) Aerial photograph with existing and new drilling locations as well as geophysical survey lines collected in the current study indicated (for fence diagrams, see Figure 6). (b) surface lithologies, archeological remains, and seismic lines (see Figures 7, 8, Supplemental Figures S5 and S6). (c) Modern coastal and seabed elevation surfaces from DEM and bathymetric data.

Our current study combines shallow-marine, high-resolution seismic reflection data, sedimentological data from the coastal zone, OSL dating chronology, and reference to existing archaeological records. The integration of these datasets reveals how regional processes of sea-level fluctuations, climate, and aeolian processes shaped the coastal area through time, as well as allowing assessment of the effect of these factors on the recently discovered early Holocene prehistoric coastal community from the south bay of Tel-Dor in northern Israel (Figs. 1–3).

REGIONAL SETTING

The coast of Israel in general and the Carmel coast in particular are ideal locations for studying long-term dynamics occurring between natural processes affecting coastal development and changing patterns of human adaptation. This is due to the microtidal nature of the Israeli coast (± 0.40 m, Emery and Neev, 1960; Davis and Hayes, 1984), its relatively tectonically stable character (e.g., Schattner et al., 2010; Sivan et al., 2010; 2016) with negligible isostatic variation based on the GIA model in Sivan et al. (2016) since the last interglacial (125 ka), and a constant, uniform supply of sediment from the Nile River since the Pliocene (Golik and Rosen, 1999), which eliminats the need to distinguish between sediment sources. Finally, humans have occupied the Carmel coastal plain and the (present-day) near-shore shelf since the late Pleistocene (Galili et al., 1993; Kuijt and Goring-Morris, 2002).

The width of the Israeli coastal plain varies from 15 km in the south to 3.5 km along the Carmel coast (Fig. 1b; Almagor and Hall, 1984). Up to eighteen calcareous sandstone aeolianite ridges are present along the coastal plain and continental shelf. These features are semi-parallel to the coastline and were formed during the late Pleistocene when sea level was lower and the shelf was exposed. However, a clear correlation between sea-level phases and aeolianite formation has yet to be established (Mauz et al., 2013). The number and size of these ridges diminishes northward, with only two ridges on land and two on the shelf found in the north-central to northern parts of the coast (Yaalon, 1967; Sivan et al., 1999). In these parts the shelf ridges are partly submerged and occasionally exposed in the form of small islands (Fig. 2). The troughs between the aeolianite ridges along the coast are filled with late Pleistocene–Holocene sequences up to 20 m thick. On land, these late Quaternary sediments are composed of several red-brown, sandy- to silty loam paleosol units that have been dated from ca. 60–12 ka. These paleosol units contain various subunits and hiatuses, which probably indicate long exposure to pedogenic processes, hence impeding a clear, lateral chronostratigraphical correlation to climate and sea-level variations (Sivan and Porat, 2004; Shtienberg et al., 2017).

Paludal/wetland deposits (constrained between ca. 21–9 ka; Kadosh et al., 2004; Sivan et al., 2011) and aeolian sand (< 6 ka; Cohen-Seffer et al., 2005; Roskin et al., 2015) rest on top of the late Pleistocene paleosol units. The Holocene wetland deposits consisting of organic-rich silty clay are found along the entire Israeli coast (Sivan et al., 2004a; Elyashiv et al., 2016; Shtienberg et al., 2017) and continental shelf (Avnaim-Katav et al., 2012; Shtienberg et al., 2016). Sedimentological, palynological, and faunal analyses show that these silty clay units were deposited in fresh- to brackish-water marshes (Kadosh et al., 2004; Cohen-Seffer et al., 2005), and were partly deposited during the rapid Holocene sea-level transgressive phase as well as during multiple wetter regional climate episodes.

Relative sea level (RSL) has risen from a global lowstand of ~135 m below mean sea-level (msl) during the last glacial maximum (LGM; 25 ka) to ~60 m below msl during the start of the Holocene, which closely correlates with eustatic sea-level records (Rohling et al., 2014).

The nearest RSL curve to the Israeli coast was reconstructed in the western Mediterranean along the relatively stable southern coasts of France by Lambeck and Bard (2000) from ca. 30 ka to present. Their curve indicates that sea level rose dramatically from the LGM, reaching ~35 m below present msl by the beginning of the Holocene (ca. 10 ka). Archaeological observations from the coast of Israel, which are the only available local indicators for this period, indicate that sea level continued to rise, and between ca. 9 ka to 7 ka, it transgressed from ~16 m to ~8 m below msl (Sivan et al., 2001; 2004b; Galili et al., 2019). From ca. 7 to 4 ka, RSL rose an additional ~7 m, and at ca. 4 ka, RSL was close to its present level, based on sedimentological indicators (Porat et al., 2008). From ca. 2 ka onwards, sea level was relatively stable, with sub-meter fluctuations (Dean et al., 2019).

Submerged Neolithic settlements inundated by post-glacial Holocene sea-level rise have been discovered along a 20 km stretch off the Carmel coast (Fig. 1b). The settlement patterns during this local transgressive phase and correlative landward shoreline migration suggest that sea-level rise resulted in the abandonment of coastal settlements, leading to their relocation landward (Fig. 1b; Galili et al., 1993, 2017, 2019). However, a detailed geochronological assessment that links sea-level transgression and the resulting environmental changes with cultural impact has yet to be conducted.

This paper provides a comprehensive overview and explanations of influences and processes operating at the land–sea interface following climate forcing and sea-level transgression while also considering how these changes influenced the distribution and timing of coastal settlements. The embayed setting of the coastal zone of Dor enables a detailed reconstruction of the geomorphic evolution history of the area through the last glacial–interglacial period, and allows us to assess natural influences on prehistoric coastal societies. We examine regional forcing factors (i.e., sea level, climate) that were active in the evolution of the coastal zone environment at Dor, the influence of these factors on the evolution of the study area leading to the formation of the coastal zone at Dor, and the effect of Holocene transgression on the prehistoric inhabitants at in the south bay of Dor.

METHODS

For this study, we combined surface mapping, detailed analyses of four new boreholes, OSL dating, and high-resolution, shallow-marine seismic reflection data with a large preexisting database of cores taken across the area for both academic and industrial purposes (See Fig. 1b for more detail). All elevations were referenced to local Israeli mean sea level to allow correlation across the datasets.

Geomorphological and topography mapping

The geomorphological characteristics and surface lithologies of the study area were mapped in detail through land and underwater surveys. A heavy-lift octocopter drone was used to collect photogrammetric data to create a digital surface model (DSM) of the terrestrial coast. Survey data were georeferenced using a South Galaxy G1 RTK-GPS with a vertical and horizontal error no greater than ± 7 cm, creating an orthorectified photomosaic.

New borehole drilling and laboratory analysis

New boreholes were cored during two expeditions (August 24, 2018; February 3, 2019) along the shoreline from the northern to southern end of the bay (Fig. 3a) using a Geo-probe 6620DT direct-push corer penetrating up to 8 m below the surface. Shallow marine drilling was conducted off an underwater platform at elevations of approximately −3 m relative to msl using a team of divers and a manually automated percussion system (Core AK994C1; Fig. 3a). Locations and surface elevations of the boreholes were measured using a Proflex 500 RTK-GPS with precision of ± 1 cm and ± 5 cm, respectively. The terrestrial and shallow marine cores were then split lengthwise for color analysis using a Munsell color chart, and lithological description.

Elemental variation along the cores was measured at 1 cm resolution with an Avaatech X-ray fluorescence (XRF) system, with an excitation voltage of 10 kV and 35 kV, and a 2-cm-diameter beam. Raw element values (photon counts) were then normalized to silica values, a dominant element in Israeli coastal sediments, to enable relative difference assessment for each facies/unit sample (Löwemark et al., 2011). Additionally, magnetic susceptibility values were measured by a Bartington MS3/MS2 meter and point sensor, compatible with small-diameter cores, at 1-cm resolution. Particle-size distribution measurements were conducted with a Malvern Instruments Mastersizer 2000 laser particle-size analyzer on 39 samples from core D4 (Fig. 3a) collected at 10–15 cm increments.

Age control was obtained using optically stimulated luminescence (OSL) dating of sediments collected from unopened cores in darkroom conditions. OSL dating followed the single-aliquot regenerative-dose (SAR) technique of Murray and Wintle (2003) for single-grain and small-aliquot analysis of quartz sand. Optical measurements were performed on Risø TL/OSL Model DA-20 readers with blue–green light emitting diodes (LED) (470 nm, 36 mW/cm²) as the stimulation source for small-aliquot measurements and using a green laser (532 nm, 120 mW/cm²) for single-grain dating, both measured at 125°C through 7.5-mm UV filters (U-340). The signals were calculated by subtracting the average of the last 5 seconds (background signal) from the first 0.7s (4 channels) of the 40s signal decay curve for the small aliquot analyses and subtracting the average of the last 0.2s from the initial 0.05s of the 1s luminescence measurement for single grains. Dose response curves were fit within saturating exponential curves to calculate equivalent dose (D_E) values (Table 1; Supplemental Figures 1, 2).

Table 1. Optically stimulated luminescence (OSL) age information. OSL ages are presented as thousands of years (ka) before 2018. Age analysis using the single-aliquot regenerative-dose procedure of Murray and Wintle (2003) on 1-mm small aliquots of quartz sand unless otherwise noted. Number of aliquots used in age calculation and number of aliquots analyzed in parentheses. Datum for ages is AD 2018. Equivalent dose (D_E) calculated using the Central Age Model (CAM) unless otherwise noted.

* Age analysis using the single-aliquot regenerative-dose procedure of Murray and Wintle (2003) on single grains of quartz sand (Guérin et al., 2011). Number of grains used in age calculation and number of grains analyzed in parentheses. D_E calculated using the Minimum Age Model (MAM) of Galbraith and Roberts (2012).

** Ages presented in Shtienberg et al. (2020)

D_E distributions from samples show largely symmetrical distributions with over-dispersion (OD, scatter) below the 20% cutoff for scatter beyond instrumental error in all but two samples (Supplemental Figures 3, 4). D_E values were calculated using the central age model (CAM) for all samples except the uppermost sample (USU-2950), which has a positively skewed single-grain D_E population with high overdispersion (82%, Supplemental Figure 3). This sample was calculated using the 3-parameter minimum-age model (MAM-3) of Galbraith and Roberts (2012). Aliquots were rejected if they had evidence of feldspar contamination, recycling ratio <0.1 or >1.1, recuperation >10% of the natural signal, or natural D_E greater than the highest regenerative dose given. Environmental dose rate (Table 2) was determined using radio-elemental conversion factors (Guérin et al., 2011) and corrected for attenuation due to grain size (Brennan, 2003) and water attenuation (Aitken, 1998). Cosmic contribution to the dose rate was determined using sample depth, elevation, and longitude/latitude following Prescott and Hutton (1994). Errors on D_E values are reported at 2-sigma standard error and age estimates are reported at 1-sigma standard error (Table 1). Uncertainties include errors related to instrument calibration and dose rate and equivalent dose calculations, and were calculated in quadrature using methods in Aitken and Alldred (1972).

Table 2. Dose rate information for optically stimulated luminescence (OSL) samples selected from borehole D4, D6, and D11. Radioelemental concentrations (K, Rb, Th, U) determined using ICP-MS and ICP-AES techniques; dose rate is derived from concentrations by conversion factors from Guérin et al. (2011). Grain size for all samples is 125–212μm.

* Moisture content effect by drying of the core. 20% moisture content assumed over burial history and used in dose rate calculation.

** Beta DR uses 100% of chemistry from USU-2955a (flood layer), gamma DR uses 75% from USU-2955a and 25% from USU-2955b (underlying silt).

^† Radioelemental concentrations for samples were averaged in dose rate calculation: above sample (upper values) and below sample (lower values).

^‡ Ages presented in Shtienberg et al. (2020).

Geophysical survey

A shallow-marine geophysical survey was conducted to map bathymetry as well as sub-seafloor stratigraphy. The survey was conducted on board a catamaran with 0.5 m draft, equipped with an Innomar SES-2000 light plus Subbottom profiler operating in dual frequencies of 100 kHz and 12 kHz with a resolution of up to 5 cm. Location and navigation were established by two South Galaxy G1 RTK-GPS devices with a vertical and horizontal accuracy better than ± 7 cm. The measurements generated 35 km of seismic profiles extending from depths of −2 m to −13 m relative to msl, consisting of 45 shore perpendicular and 23 shore parallel lines spaced ~7 m apart (Fig. 3a).

The geophysical data set was uploaded to the Emerson/Paradigm software suite and all processing was performed in the ECHOS software. A vertical correction was applied to each trace based on RTK-GPS Z values in order to circumvent static vertical shifts created by swell and tide. Attenuation of short-path multiples (ghosts) and signal enhancement was performed with a minimum entropy deconvolution filter, and further signal enhancement was achieved with a lateral mixing filter. Interpretation-based amplitude scaling was used to enhance the appearance of the geological units.

RESULTS AND INTERPRETATION

Coastal-zone characteristics

The study area consists of a 280-m-long sandy pocket beach with a mean width of 230 m from the shoreline to the 4-m-high sand dune in the backshore. The pocket beach is bordered by an aeolianite headland to the north and an aeolianite-cored tombolo to the south, both of which project ~290 m seaward relative to the pocket-beach shoreline (Sneh, 1981; Sneh and Klein, 1984). On the northern headland, archeological constructions (wells, slipways) as well as the remains of an Iron Age and younger urban settlement are evident, while a rock-cut pool is apparent (Arkin Shalev et al., 2019) on the rocky headland of the tombolo (Fig. 3c).

The bathymetry of the study area ranges from −13 m to −2 m relative to msl, and the general slope gradient is 1:50, extending 760 m west of the shoreline. The entrance to the south bay is semi-constrained by the two edges of the northern headland and tombolo that are 140 m apart. Between the northern headland and tombolo, at water depths of approximately −3 m relative to msl, there is a N–S trending, submerged, elongated, calcareous sandstone platform that connects the headlands (Fig. 3b). A westward-trending, curvilinear depression 40 m wide and ~2.5 m deep truncates this platform at elevations of −7 m to −4 m relative to msl (Fig. 3c). Other than the submerged calcareous sandstone platform, the seabed is mostly covered by sand.

Sedimentology of the south bay of Dor

The sedimentological sequence of the south bay of Dor (Figs. 1–3) consists of five lithological units (F1–F5) that unconformably overly the basal unit (BU) aeolianite. The lithological classification (Figs. 4, 5) of each unit and their spatial correlation (Fig. 6) is based on sedimentological characteristics, thickness, elevations, accompanying features, lithostratigraphical relations, and geochemical properties. Descriptions of the units from the oldest unit to the youngest modern unit is as follows:

Figure 4. Boreholes D8, D11, and AK944C1 with graphic log, lithological unit (see text for detailed lithological unit descriptions), lithology description, accompanying features, distinguishing Munsell color (dry), magnetic susceptibility, relative elemental variations, and optically stimulated luminescence (OSL) ages. See Fig. 3a for locations.

Figure 5. Borehole D4 (see Fig. 3a for location) with graphic log, lithological unit (see text for detailed lithological unit descriptions), Munsell color description (dry), relative elemental concentration, grain texture (percent grain size), magnetic susceptibility, lithostratigraphical interpretation, accompanying features (stars indicate OSL sampling), and OSL age results from this study and Shtienberg et al. (2020). Dashed lines indicate approximate facies boundaries.

Figure 6. Chronostratigraphic cross-sections from the coastal area of Dor based on OSL ages (white stars = this study) and sedimentological data. 14C results from lithologic samples results (green circles) and radiometric ages (red circles) taken from Kadosh et al. (2004). For map of profile locations and modern topography and bathymetry used in cross-sections, see Figure 3a, c. Because marine core AK994C1 did not extend deeper than −3.0 m msl, subsurface thickness and elevation is based on seismic profiles and seismic unit interpretation (see Supplemental Figure 5). BU = basal unit; F1, F2, F3, F4, F5a, F5b, F5c = lithologic facies; see text for full descriptions of unit/facies.

Basal Unit (BU): The upper ~0.5 m of the basal unit was reached in boreholes D8, D4, and D6 (Fig. 3, Supplemental Figure 3), and consists of yellowish white to pale gray (10YR6/2 to 2.5Y8/1), indurated to lightly cemented, calcareous-cemented to partly cemented sands. Ratios for Ti/Si, Zn/Si, Fe/Si, Rb/Si, and Al/Si are 0.15, 0.07–0.2, 0.15, 0.05, and 0.05, respectively, and magnetic susceptibility values ranged from 15–100 Si × 10⁻⁶ (Figs. 4, 5). This unit was dated in core D4 to 77.34 ± 7.34 ka (Table 1, Fig. 5), and resembles the sedimentology of headlands that border the south bay at Dor (Figs. 2, 3). It is interpreted as carbonate-cemented aeolian dune sands.

Unit F1: Unit F1 overlies the BU, and was identified in boreholes D8, D4, D12, and D6 (Figs. 3–5, Supplemental Figures 4, 5) at surface elevations ranging from −8.5 m to −3.6 m relative to msl, and is 0.5–3 m thick (Fig. 6). This unit consists of red to brown (10YR4/2 to 10YR5/8) silty loam with irregularly shaped, hard calcareous nodules and black manganese concretions, and lacks micro- and macrofaunal remains. The size and abundance of the nodules and concretions decreases upwards, and nodules are completely absent at the top of the unit. F1 has similar ratios for Ti/Si, Zn/Si, Fe/Si, Rb/Si when compared to BU, but has a strong signal for Al/Si. Magnetic susceptibility values for the unit vary from 100–300 Si×10⁻⁶ (Figs. 4, 5). Two OSL ages acquired from the bottom and middle parts of the unit date it to 59.76 ± 5.03 ka and 30.59 ± 3.20 ka respectively (Table 1, Figs. 5, 6; Supplemental Figures 4, 5). Based on these characteristics and lithostratigraphic correlation with previously collected cores (Fig. 6; Kadosh et al., 2004; Sivan et al., 2004a; Shtienberg et al., 2017), F1 represents silty, well-oxidized aeolian sediments of the paleo-coastal plain, with a Bt horizon interpreted as a paleosol that locally is known as Hamra.

Unit F2: Unit F2 was found in boreholes D11, D8, D6, D12, and D4, with surface elevations ranging from −6.8 m to −2.6 m relative to msl, and is 0.5 m to 1.5 m thick (Figs. 4–6, Supplemental Figures 4, 5). The unit is dark gray to dark brown (10YR4/1 to 10YR6/3) with a silty loam texture. F2 also contains brackish-water microfaunal remains as well as fresh- to brackish-water plant remains, such as charophytes (Supplemental Data Set 1). Its elemental composition includes higher Ti/Si, Rb/Si, Fe/Si ratios compared to the underlying unit and higher MS values (Figs. 4, 5). Two OSL ages acquired from the middle and upper parts of this unit date F2 to 14.94 ± 1.34 ka and 9.42 ± 0.85 ka, respectively (Table 1, Figs. 5, 6). Unit F2 is interpreted as a fresh- to brackish-water wetland deposit.

Unit F3 unconformably overlies unit F2 and was identified in boreholes D11, D8, D4, D12, and D6 (Figs. 4–6, Supplemental Figures 4, 5), with surface elevations ranging from −3.3 to −1 m relative to msl, thicknesses of 0.12 m to 1 m (Fig. 6), and a yellow to light-yellow color (10YR6/3 to 10YR5/6). The unit consists of poorly sorted sand composed of quartz grains, Glycymeris and Tucetona shells, limestone pebbles, gastropod shells, calcareous sandstone fragments and dark, silty clay rip-up clasts. Major element ratios for Al/Si, Ti/Si, Fe/Si, Rb/Si and MS values are lower compared to the underlying unit (F2), and an OSL age dates the unit to 10.19 ± 0.90 ka (Table 1). Unit F3 is interpreted as a high-energy coastal/shallow-marine sand. Previous work has interpreted this unit as a tsunami deposit (Shtienberg et al., 2020).

Unit F4: Unit F4 unconformably overlies unit F3 and was identified in boreholes D11, D8, D4, D6, D12, and AK994C1 (Figs. 4, 5, Supplemental Figures 4, 5), with surface elevation ranges of −6 to −0.5 m relative to msl and thicknesses of 0.4 m to 0.8 m (Fig. 6). This unit is dark gray to dark brown (10YR3/1) homogenous silty loam that contains brackish-water microfaunal remains, sea urchin spines, and gastropod shells in the upper 0.3 of the unit (Supplemental Data Set 1). Elemental ratios of Al/Si, Ti/Si, Zn/Si, Fe/Si, and Rb/Si, as well as MS values resemble those of F2 (Figs. 4, 5). Unit F4 was dated to 9.15 ± 0.78 ka by an age obtained from its base, and two ages of 7.78 ± 0.71 and 7.20 ± 0.70 from the upper part of this unit (Table 1, Figs. 4–6). F4 is interpreted as a dark, silty clay, fresh- to brackish-water wetland deposit similar to unit F2.

Unit F5: Unit F5 unconformably overlies unit F4, and generally has lower Al/Si, Ti/Si, Fe/Si, Rb/Si ratios (Figs. 4, 5) compared to the underlying unit. This unit has been broken down into three subfacies, described here from bottom to top:

1. 1. Facies F5a was identified in boreholes D11, D8, D12, D6, and D4, with surface elevations from −3.5 to −0.5 m relative to msl, and thicknesses of 1 m to 2.5 m (Figs. 4, 5). The sediment consists of sub-rounded, medium-grained, light brown to yellow (2.5YR6/4 to 10YR6/3) sand. Glycymeris shells are evident throughout the facies with abundance and size increasing upward. An OSL age acquired from the base of this subunit dates F5a to 7.10 ± 0.68 ka (Table 1, Figs. 5, 6) We regard this subunit as a coastal sand similar to modern beach deposits.

2. 2. Facies F5b occurs in boreholes D11, D8, D4, D6, and D12, with surface elevations of 0 to 0.3 m relative to msl, and a thickness of 0.4 m to 1 m (Figs. 4–6, Supplemental Figures 4, 5). The sediment is light gray to light yellow (2.5YR6/6 to 2.5YR3/0) coarse sand with abundant Glycymeris shells, aeolianite fragments, and chert and limestone pebbles. F5b also contains cultural material such as pottery sherds, pieces of iron, and charcoal fragments. (Figs. 4, 5). Trace-element ratios of Al/Si, Ti/Si, Fe/Si, Rb/Si as well as MS values are higher compared to the underlying facies (F5a), and OSL results from the base and upper portions of this facies produced ages of 2.80 ± 0.28 ka and 1.12 ± 0.20 ka, respectively (Table 1, Figs. 5, 6). Unit 5b is a coastal sand package associated with the Bronze Age–Crusader settlement at Tel-Dor.

3. 3. Facies F5c is found in boreholes D11, D8, D4, D6, and D12, with surface elevations of 1 m to 1.5 m above msl and a thickness of 0.8 m to 1.25 m (Fig. 5). The unit is light gray to light brown (10YR7/4 to 10YR6/2) medium sand with Glycymeris shells and iron and glass fragments in the upper 0.3 m of the unit (Figs. 4, 5). Similar characteristics were also identified in the upper few centimeters of the shallow marine core AK9444C1 (Fig. 3). Based on these sedimentological properties and accompanying features, facies F5c is interpreted as a coastal sand facies that is contemporaneous with the historic setting of the Carmel coast (Fig. 1).

Seismic stratigraphy

Geophysical analysis revealed four seismic facies within the bay and two seismic facies outside the south bay of Dor (Figs. 3b, 7). The acoustic basement (AB) is characterized as a high-amplitude unit with chaotic and highly diffracted reflections, and with limited seismic penetration. The surface of the unit consists of an irregular, uneven set of reflections (Figs. 7, 8), which extend over the entire study area from elevations of −1 m to −13 m relative to msl (Fig. 9). A north–south striking, elongate structural high is recognized at the entrance of the bay at elevations of −1 m to −4 m relative to msl, and ~0.2 km offshore of the current coastline. This morphological feature has a maximum width of 100 m, with crests that are up to 4 m above the surrounding surface topography. The main axis of this structure is perpendicularly dissected by a ~3-m-deep trough at water depths of −6 m to −3 m relative to msl.

Figure 7. (a) An E–W trending seismic section from Dor (see Figure 3a for location). Dashed rectangles indicate area of the enlargements presented in (b), (d), and (f); corresponding interpretations of these sections are presented in (c), (e), and (g). (h) A general description of the seismic units as described in the text, and accompanying appearance in the seismic data. TWT = two-way travel time. VE = vertical exaggeration.

Figure 8a. (i) A shore-perpendicular seismic section without interpretation, and (ii) the same seismic section with interpretation from seismic units to lithologies, which was achieved based on correlation to boreholes AK994C1, D4, D11, and D8, as well as surface lithologies. See Figure 3 for location of the seismic section and boreholes. Red line marks the seafloor. TWT = two-way travel time.

Figure 8b. (i) A shore-normal seismic section without interpretation, and (ii) the same seismic section with interpretation from seismic units to lithologies, which was achieved based on correlation to boreholes AK994C1, D4, D11, D8 as well as surface lithologies. See Figure 3 for location of the seismic section and boreholes. Red line marks the seafloor. TWT = two-way travel time.

Figure 8c. (i) A shore-perpendicular seismic section offshore from the south bay of Dor without interpretation, and (ii) with interpretation. See Figure 3 for location of the seismic section and boreholes. Red line marks the seafloor. TWT = two-way travel time.

Figure 9. Elevation maps of lithological facies and correlated seismic unit surfaces inside and adjacent to the south bay of Dor. (a) Elevation map of the submerged surface of the late Holocene sand unit (= Sb4, Sd2); (b) Elevation map of the submerged Early Holocene wetland unit surface (= Sb3, Sd1); location of Neolithic Pottery remains marked by a red rectangle. (c) Elevation map of the submerged Late Pleistocene wetland unit surface (= Sb2). (e) Elevation map of the submerged Late Pleistocene aeolianite unit surface (= AB). See Figure 8 for definition of seismic units.

Seismic unit Sb1 overlies the acoustic basement and was identified inside the south bay. It is categorized by a semitransparent unit that occasionally displays low-amplitude and hummocky reflections onlapping the acoustic basement surface on its lower surface and is bounded by a high-amplitude reflector at the unit's top (Figs. 7, 8a).

Seismic unit Sb2, which partly covers unit Sb1, is categorized by high-amplitude, hummocky to chaotic reflections that onlap the acoustic basement, and it has a toplap relationship with unit Sb3 (Fig. 7). The upper surface of unit Sb2 ranges from −5 m to −3 m relative to msl and dips westwards at ~2° (Fig. 9c). A sink-like depression, ~3 m lower than its surroundings, is evident southeast of the bay's entrance.

Seismic unit Sb3 is composed of parallel to subparallel, mid- to low-amplitude reflections that, for the most part, onlap the surface of Sb2 and are truncated by Sb4. The unit is up to 2.5 m thick (Fig. 10b), and onlaps the morphological high defined by the acoustic basement. The surface morphology consists of a sloping depression at −6 m to −3 m relative to msl, dipping westwards at ~1° (Fig. 9b). A rounded depression is identified inside the bay in a location and extent that resembles the underlying depression of Sb2.

Figure 10. Isopach map of Dor's seismic litho-stratigraphy. These thickness maps were computed using the ArcGIS Minus Toolset that operated on the interpolated litho-surfaces that are presented in figure 9. The litho-facies are annotated from young (top left) to old (bottom right) (a) Thickness map of the middle to late Holocene sand unit (= Sd2, Sb4); (b) Thickness map of the early Holocene dark silty clay unit (Sd1, Sb3); (c) Accumulated thickness of the late Pleistocene–early Holocene dark silty clay unit and red-brown paleosol unit in the south bay of Dor (Sb1+Sb2+Sb3). See Figure 8 for definition of seismic units.

Seismic unit Sb4 displays semitransparent, wavy to subparallel mid- to high-amplitude west-dipping reflections (Fig. 7). This unit varies in thickness from 0.2 m to 2 m (Fig. 10a), and its surface is −3 m to −1 m relative to msl (Fig. 9a). Similar to underlying unit Sb3, a rounded depression was identified inside the bay ~2.5 m lower than its surroundings.

Seismic unit Sd1 is found west of the entrance to the south bay, and consists of high-amplitude, semicontinuous, sub-parallel west-trending reflections (Fig. 7). Unit Sd1 represents a prograded infill package that toplaps Sd2 and downlaps the surface of AB (Fig. 8c). The thickness of the unit varies from 1 m to 4 m. It is thicker in the southern parts of the unit, reaching elevations of −7 m relative to msl, 300 m from the present shoreline (Fig. 10b).

Seismic unit Sd2 was identified west of the entrance to the south bay and is characterized by continuous, medium-amplitude, medium-frequency, subhorizontal and subparallel reflectors onlapping the surface of AB (Fig. 8c). The unit was identified up to 250 meters from the entrance of the south bay at elevations of −5 m to −11 m relative to msl (Fig. 9a) and with thicknesses of 0.5–2 meters (Fig. 10a).

DISCUSSION

Stratigraphic architecture

The seismic units within the south bay boundaries (Sb1–Sb4) and west of its entrance (Sd1, Sd2) were correlated with the stratigraphic sequence of the Carmel coast based on geometric relations between the units, respective seismic facies, lithostratigraphical relations, and sedimentological correspondence with the on-land boreholes (Figs. 1b, 3, 5).

The morphology of the acoustic basement surface (AB) resembles the elevation differences, dipping angles, irregularity, and shore-parallel direction of the basal unit as observed on the adjacent coast (Figs. 3–6; Sneh, 1981; Sneh and Klein, 1984; Kadosh et al., 2004; Sivan et al., 2004a). The surface morphology and topography of AB (Fig. 7) also match those of the submerged aeolianite surface as interpreted from chirp surveys conducted on the shallow shelf between Hadera and Caesarea (~8 km south of Dor; Fig. 1b) at elevations of −5 m to −45 m msl (Shtienberg et al., 2016, 2017; Goff et al., 2018). Taken together, this leads to the identification of the acoustic basement as the top of the late Pleistocene calcareous sandstone surface that formed between 80–60 ka (Middle Paleolithic period).

The calcareous sandstone aeolianite surface is overlain by seismic unit Sb1 that was identified within the south bay. Boreholes drilled for this study (Fig. 6), as well as shallow seismic stratigraphy conducted ~8 km south of the study area (Fig. 1b), suggest that seismic unit Sb1is a late Pleistocene red-brown paleosol (Unit F1), which was deposited from ca. 60 ka to 15 ka (Upper Paleolithic to Natufian period).

The internal geometry and the strong bounding surface of seismic unit Sb2 that partly covers unit Sb1 suggest erosive contacts at the top and base of the seismic unit (Figs. 7, 8). Based on the coastal boreholes located a few tens of meters east of the seismic survey (Figs. 3a, 4, Supplemental Figure 5), as well as the erosive surface morphology (Fig. 5; Sneh, 1981; Kadosh et al., 2004), Sb2 is interpreted as a combination of a brackish wetland deposit (Unit F2) deposited between the late Pleistocene to early Holocene, and a coastal/marine sand (Unit F3) deposited between ca. 15 ka and 10 ka (Natufian to Pre-Pottery Neolithic B period; PPNB).

Seismic unit Sb3 overlaps unit Sb2 and consists of a strong bounding-surface reflector. Based on the lithology identified in core AK994C1 that penetrated this unit (Supplemental Figure 5), the core record of Dor (Figs. 3a, 6), and the unit's seismic geometry and erosional morphology, Sb3 is interpreted as the early Holocene brackish wetland unit (F4) that was deposited from ca. 9.15 ka to 7.20 ka (PPNB to end of the Pottery Neolithic period; PN).

Seismic unit Sb4 covers the coastal sequence in the south bay at Dor and is interpreted as the mid-Holocene coastal sand unit (correlative to F5 in the terrestrial borehole data) deposited between 7.1 ka–present (Chalcolithic–Modern). This interpretation is based on the pocket-beach core record (Figs. 3, 6, 8; Supplemental Figures 5, 6) as well as an underwater field survey conducted using a suction system that revealed the contact between the coastal sand and the underlying silty clay.

West of the bay, a prograded infill unit (Sd1; Figs. 7, 8) is identified above the submerged aeolianite (AB). Seismic unit Sd1 is interpreted as an early Holocene brackish wetland based on the similarities of this unit to three lens-shaped facies identified at similar depth beneath the seabed, with similar thicknesses (Shtienberg et al., 2016), 8 km south of the study area (Fig. 1b).

Finally, the topmost unit (Sd2) is identified as a mid- to late Holocene sand unit, based on an underwater survey conducted in the area (Figs. 1, 7–10) as well as similarities to the interpretation of Shtienberg et al. (2016).

Coastal paleoenvironmental reconstruction: a regional perspective

The combination of seismic stratigraphy, coastal coring, and dating reveals the Pleistocene–Holocene paleoenvironmental history of the Carmel coast. Between 77–8 ka, sea level was tens to several of meters lower than present (the last glacial maximum sea-level lowstand was 135 m below present msl) (Figs. 1b, 11m; Lambeck and Bard, 2000; Rohling et al., 2014). During maximum lowstand, the shoreline was ~13 km west of its current location; subsequent sea-level rise moved the shoreline to ~0.5 km west of present by ca. 8 ka, and it reached its current position ca. 4 ka (Sivan et al., 2004b; Galili et al., 2019). Hence, most of the stratigraphic sequence at Dor, other than the mid- to late Holocene sand unit (F5), was deposited in a terrestrial environment.

Figure 11. Comparison between a composite of sequences from the south bay of Dor (this study), general Carmel coast stratigraphy, and regional proxies. Columns (a–e) represent the following stratigraphic and age results: (a) generalized lithofacies of the south bay sequence; (b) sediment bed; (c) average thickness; (d) unit age range (from this study; Shtienberg et al., 2020); (e) sedimentological characteristics. Generalized Carmel coast sequence that includes: (f) generalized lithofacies of the Carmel coast; (g) bed/unit; (h) average thickness; (i) unit age range; (j) accompanying features. This compilation was computed through the integration of the thickness and ages: (1) Kadosh et al. (2004), (2) Sivan and Porat (2004), (3) Sivan et al. (2004a), (4) Cohen-Seffer et al. (2005), (5) Sivan et al. (2011), (6) Mauz et al. (2013), (7) Shtienberg et al. (2017) (see Figure 1 for locations of these studies); (k) correlative archaeological periods based on Galili et al. (2020) and Kuijt and Goring-Morris (2002); (l) paleoclimate reconstruction based on Core 9505 (Langgut et al., 2011), Soreq cave speleothem record as a proxy for precipitation and runoff (Bar-Matthews et al., 2003), and deep-sea record as a proxy for Levantine sea response to global ice accumulation (Almogi-Labin et al., 2009; Revel et al., 2010); (m) Global sea level stack (Spratt and Lisiecki, 2016) and relative sea-level curve (red polygon; Dean et al. 2019; Galili et al., 2019; Sivan et al., 2001) against approximated distance from present shoreline.

The south bay of Dor is protected by north–south trending islands, wave cut platforms, and headlands comprised of calcareous sandstone (Sneh, 1981; Sneh and Klein, 1984; Kadosh et al., 2004; Sivan et al., 2004a), and these sandstones are found above and below msl (Figs. 2, 3b, 9d). Based on luminescence age dates from the island headlands, 0.5 km south of the study area (Fig. 2; Sivan and Porat, 2004; Mauz et al., 2013), these sandstones represent aeolianites, and were initially deposited ca. 60–35 ka as parabolic sand dunes (Roskin et al., 2015, Shtienberg et al., 2017) on a ca. 77 ka aeolianite bedrock surface (Table 1, Figs. 5, 6, 9d). The dunes were then lithified, forming aeolianite mounds that were eroded over time into their current form. During the mid-Holocene, ca. 4 ka, when sea level rose to its current elevation (Dean et al., 2019), the lithified aeolianite headlands were further battered by wave energy (Fig. 2c), and were later mined for building stone (Figs. 2, 3c). As a result of these depositional and erosional processes, the aeolianite mounds consist of a rigid surface that is characterized by ridges, depressions, and E–W cutting troughs (Figs. 2, 3, 9d). The topography on top of the aeolianite ridges had a major impact on subsequent sedimentation (Fig. 9d) as well as the present coastal morphology of the area (Fig. 2).

From ca. 60–20 ka (Upper Paleolithic to Natufian period; Figs. 2, 11), the parent material of the paleosol unit (F1/Sb1) covered the aeolianite surface, with initial deposition occurring in topographic lows of 1–4 m relative to present msl (Figs. 6, 11). This relationship is in agreement with the findings of Sivan and Porat (2004), indicating that local settings (i.e., topography, hydrological conditions) influence differential processes of pedogenesis/lithogenesis (Shtienberg et al., 2017). Based on the chronology of this unit, it seems that F1/Sb1 mostly correlates to Israel's coastal plain Netanya paleosol (40–20 ka; Fig. 11) (Tsatskin et al., 2008), and represents extensive soil formation during the last glacial period along the coast of Israel (Fig. 11; Mauz et al., 2013; Shtienberg et al., 2017) as well as at other locations in the Mediterranean Basin (Muhs et al., 2010) and elsewhere (e.g., Canary Islands; von Suchodoletz et al., 2010). During this episode, climatic conditions were relatively stable (Almogi-Labin et al., 2009; Langgut et al., 2011) and global msl fluctuated between −40 to −135 m relative to present msl. The paleosol formation is associated with Heinrich stadials 6–5, which have been identified globally as periods of increased Saharan dust input (Mahowald et al., 2006; Torfstein et al., 2018) that affected the Mediterranean Basin and adjacent terrestrial areas. Wind-blown dust contributed ferromagnetic minerals to the porous parent material of the soil profile. These aeolian additions, along with contributions from weathering and leaching within the soil profile, account for the silt-clay lithology and elevated Al/Si, Fe/Si, and Ti/Si ratios and MS values (~200 Si × 10^-6) of unit F1 (Tsatskin et al., 2015; Shtienberg et al., 2017).

The upper parts of the red-brown paleosol within unit F1 gradually change into a ~1 m-thick, dark-gray loamy unit (F2/Sb2) that was deposited between ca. 15 ka and 9.9 ka (Natufian to PPNB; Figs. 6, 11). Elevated MS values (250–900 Si × 10⁻⁶) and high ratios of Fe/Si, Ti/Si, and Zr/Si, along with the faunal and floral remains, suggest that Unit F2/Sb2 is composed of brackish wetland sediment (Sivan et al., 2004a; Cohen-Seffer et al., 2005). The wetland may reflect climate conditions of the African Humid Period (ca. 15–7 ka), postulated to be the outcome of extreme insolation values at 65°N (Imbrie et al., 1984), and evident in a variety of local and regional proxies such as the Soreq cave speleothems, and southeastern Mediterranean deep-sea core record (Figs. 1a, 11l; Bar-Matthews et al., 2003; Almogi-Labin et al., 2009; Revel et al., 2010). Based on the global combined sea-level stack (Spratt and Lisiecki, 2016) and relative sea level of the Mediterranean (Fig. 11m; Rohling et al., 2014) at the time of deposition of Unit F2, the coastal wetland found along the northern Israeli coast would have been located several kilometers to the east of the paleo-coastline (Fig. 9).

Wetland deposition was interrupted by a coastal/marine sand deposit (Unit F3) that is bounded by the upper and lower wetland units (Units F2, F4). The sedimentological change of the wetland unit to marine/coastal sand is distinct, and is identified across the terrestrial and shallow marine parts of the bay. Unit F3 consists of poorly sorted sand containing clay clasts, marine gastropod and bivalve Glycymeris shells as well as limestone and calcareous sandstone pebbles (Figs. 4–6). Erosional features of the bounding reflector of Sb2 and evidence of an erosion surface at the base of Unit F3 in the core data indicate that this coastal/marine sand was deposited unconformably in the lower wetland deposit (Unit F2). The age of coastal/marine sand (F3) including its uncertainties (10.19 ± 0.90 ka), and ages from the underlying wetland surface from F2 (9.42 ± 0.85 ka) as well as the base of the overlying wetland unit F4 (9.15 ± 0.78 ka) brackets deposition to between 9.91 ka to 9.29 ka (Table 1). Relying on this age constraint (9.9–9.3 ka), sea level and projected paleoshoreline during the deposition of F3, as well as the marine faunal remains embedded in this unit, suggest that Unit F3 was deposited at a landscape elevation above the early Holocene sea level at the time. We suggest that this coastal sediment represents a tsunami deposit as discussed in Shtienberg et al. (2020). Similar deposits with similar ages have yet to be identified along the coast of Israel or in other coastal sites around the Levant basin.

A second wetland unit (F4), resembling unit F2 in sedimentological characteristics, caps unit F3, and based on its chronology, it was deposited from ca. 9.3–7 ka (Pre-Pottery Neolithic B to Pottery Neolithic period; Fig. 5). Faunal and floral remains identified in the unit as well as the chronostratigraphic correlation to previous studies from the Carmel coast (10.3–9 ka; Fig. 11) indicate that this wetland unit was deposited in a brackish-water environment (Sivan et al., 2004a; Cohen-Seffer et al., 2005). Relative (Rohling et al., 2014) and local (Sivan et al., 2001; Galili et al., 2019) sea-level curves suggest an abrupt sea level rise during the early Holocene, transgressing from 40 m to 16 m below present msl. As sea level rose, the shoreline migrated eastward from 4 km to 1.5 km west of Dor, elevating coastal aquifer levels that contributed water to the coastal wetland systems (Fig. 9; Sivan et al., 2011; Shtienberg et al., 2016). The flooding of the shelf, as well as prevailing wetter conditions during this part of the African Humid Period (Bar-Matthews et al., 2003, Kadosh et al., 2004; Grant et al., 2016), reduced coastal stream gradients and would have inundated topographic lows. OSL ages place deposition of the upper wetland deposit in core D11 at 7.2 ± 0.70 ka, and initial deposition of the overlying coastal sand at 7.10 ± 0.68 ka (Table 1, Fig. 6). These ages and their reported uncertainties constrain the termination of wetland deposition to sometime around 7.2 ka (Supplemental Figure 7). This age constraint suggests that this environmental transition from a wetland to a coastal site could be partly linked to the end of the African Humid Period (i.e., ca. 7 ka; Fig. 11l) (Bar-Matthews et al., 2003; Almogi-Labin et al., 2009; Revel et al., 2010; Langgut et al., 2011; Ehrmann et al., 2017).

The lowermost facies of the upper sand unit (F5a) consists of a 1–3 m-thick, medium-grained sand unit with subrounded grains (Figs. 4, 5) and abundant Glycymeris shells. The base of this unit was dated to 7.10 ± 0.68 ka in core D11 (Pottery Neolithic–Chalcolithic Period; Fig. 6) within the uncertainty bounds of the underlying wetland deposit (7.20 ± 0.70 ka), and is consistent with a similar age of 6.95 ± 0.68 ka acquired from the bottom of F5a in core D4 (Fig. 6). This overlap (Supplemental Figure 7) suggests that initial coastal sand deposition in unit F5a began sometime ca. 7.1 ka. During deposition of this facies, sea level was ~7–4 m lower than present, while the shoreline was located a few hundred meters westward of its current location (Fig. 11). Hence, at least initially, the coastal sand unit was probably deposited by storm wash-over and strong winds that moved sand into the adjacent wetland surface from the paleocoastline.

Overlying the coastal sand facies (F5a), facies 5b consists of poorly sorted, coarse sand and gravel, and Glycymeris shells. These sedimentological characteristics, two late Holocene ages (2.80 ± 0.28 and 1.14 ± 0.20 ka), and correlation with the relative sea-level curve of the Israeli coast, indicate that Unit F5a was deposited in a coastal environment when sea level was similar to present (Dean et al., 2019). Archeological remains (sherds, iron fragments), and high MS values as well as elemental ratios for Fe/Si, Zn/Si, (Figs. 4, 5) are the result of extensive human habitation from at least since Middle Bronze Age II and into the Roman period (Stern, 1994; Raban, 1995; Shahack-Gross et al., 2005; Sharon et al., 2005; Gilboa and Sharon, 2008).

The sand sequence is topped by ~1.5 m of facies F5c, which consists of medium- to coarse-grained sand distributed with shell material, and the facies also has low MS values. Given that the seasonal beach profile changes of the Carmel coast can fluctuate between 1–1.5 m (Shtienberg et al., 2014) and the presence of glass and other recent artifacts, this facies is interpreted as modern pocket-beach sands.

Pottery Neolithic settlement of the Carmel coast and environmental change

Underwater excavations have discovered a cluster of ten submerged Pottery Neolithic (PN) sites (8.25–7.8 ka; Galili et al., 2020) from Kfar Samir in the north to Habonim in the south (Fig. 1b). These sites were all found a few meters to a few hundred meters west of the present shoreline, in water depths ranging from 0 to 5 m below msl (Galili et al., 2017, 2020), and situated on top of the early Holocene wetland facies unit F4 (Fig. 11). Recent underwater excavations show that submerged settlements are also present 12 km south of Atlit, in the south bay of Dor (Fig. 3c). Excavations conducted in 2019 by the University of Haifa and University of California San-Diego uncovered PN sherds, a submerged rectangular building, and a circular feature 1.5 m below msl (Shtienberg et al., 2020) on the surface of the organic-rich, silty clay unit (F4/Sb3) that represents the top of the wetland facies. Relying on the youngest age constraint of the F4 surface dated to 7.20 ± 0.70 ka, we propose that the submerged rectangular building is associated with PN settlement.

A correlation between the PN period to the chronostratigraphy of the Carmel coast (Fig. 11) and Levantine climate patterns implies that the settlements of the Carmel coast were located on a dried wetland surface (Sivan et al., 2004a, 2011; Cohen-Seffer et al., 2005) close to the end of the African Humid Period (Bar-Matthews et al., 1999, 2003). Sea level rose rapidly during the PN period, transgressing from ~10 m to 5 m below present msl, at a mean annual rate of ~13–15mm/year (Stanford et al., 2011). During this period, the fast-transgressing shoreline would have been coupled with spring tides and extreme winter storms, resulting in wash-throughs that reached the entrance of the south bay of Dor. The time of abandonment of the coastal PN settlement is not known, but could be constrained based on the relationship of the relative sea-level change and corresponding PN site elevation. Taking these data into account, it is likely that the sea flooded the site and the south bay of Dor ca. 7.6–6.4 ka.

SUMMARY AND CONCLUSIONS

The present-day coastal landscape of Dor is the result of differential erosion of the north–south trending aeolianite ridge, which was deposited and lithified in two phases between 77–50 ka. The aeolianite bedrock is overlain by a terrestrial silty unit containing red-brown paleosols. Wetlands developed over the paleosols and lithified dunes during the African Humid Period. A high energy deposit, likely the product of a tsunami, interrupted wetland deposition. The upper surface of the wetland is associated with PN archaeological material, suggesting that coastal peoples settled on the surface of the dry marsh. These coastal communities were displaced by sea-level rise, which has contributed to the deposition of the upper coastal sand since ca. 7.1 ka. The preservation of late Pleistocene to Holocene sedimentary deposits and the associated archaeological record identified in our study was enabled due to the protective characteristics of the south bay of Dor and its coastal embayment created by the shore-perpendicular ridge.

The sedimentological and archaeological evidence identified in the stratigraphic record of the sheltered south bay of Dor demonstrate the importance of research-area selection. The morphological characteristics of the south bay of Dor has led to more complete preservation of the stratigraphical record, which helps to fill geochronological gaps in the evolution of the paleo-coastal plain over a glacial–interglacial period to the formation of the current coastline. This combined geological and archeological approach can help clarify the close relations between naturally changing environments and human settlement behavior in other coastal areas that exhibit a long historical settlement record and protective geomorphological conditions.