Paleoenvironmental Evolution and Sea Level Change in Saronikos Gulf (Aegean Sea, Greece): Evidence from the Piraeus Coastal Plain and Elefsis Bay Sedimentary Records

Abstract

Thorough faunal (benthic foraminifera, ostracods, molluscs) and palynomorph analyses as well as magnetic susceptibility measurements performed on the Piraeus coastal plain sedimentary sequences have shed light on the paleoenvironmental evolution of the area since ca. 9000 cal BP. Benthic and palynomorph assemblages along with magnetic susceptibility suggest a typical lagoonal environment with significant freshwater inputs at the eastern part of the plain after 8700 cal BP. Between 7500 and 5400 cal BP, microfaunal assemblages, mollusc fauna and magnetic susceptibility suggest a shallow marine paleoenvironment, with Piraeus forming a tied island in the center of the bay. Since ca. 4800 cal BP a closed oligohaline lagoon is evidenced in the western part of the Piraeus plain further developed to a marsh after 2800 cal BP, while a coastal environment associated with the fluvio-deltaic system of Kifissos and Korydallos Rivers is continually developing to the west. Signs of cultivation and grazing activities in the area are evidenced since the Early Bronze Age, culminating during the Classical Period. A comparison with a well-dated marine record, recovered from the nearby shallow Elefsis Bay, provides a reasonable estimation of ~5 mm/yr for the absolute sea level rise rate in the inner Saronikos Gulf during the Mid-Holocene.

1. Introduction

In the Eastern Mediterranean coastal areas, the excessive population growth has been putting an unprecedented anthropogenic pressure on marine ecosystems, severely impacting their natural balance and resulting in extensive biogeographic shifts (e.g., [1,2]). A progress in the understanding of the paleoenvironmental changes in the Eastern Mediterranean ecosystems will be of great significance, since the region is characterized by great cultural, economical, political and demographic gradients, which are already affected by considerable environmental stress. In particular, in the NE Mediterranean, the semi-enclosed Aegean Sea is an ideal natural laboratory [3] for analyzing the complex interactions between climate, marine environment and humans in the present time and recent past, due to its unique physical and geographic configuration and the great development of human societies in its coastal zone.

In coastal marine and lagoonal environments, where environmental stress is inherently high, benthic foraminifera and ostracods are the most abundant shelled microorganisms, representing the most sensitive environmental indicators [4,5]. Changes in benthic faunal abundance and species composition offer good evidence for the variability of several environmental factors and can, therefore, be used as an efficient method for the determination of the ecosystems’ conditions [6,7,8,9]. In particular, changes in fossil benthic foraminifera, recorded by the study of coastal zone deposits, are a successful paleontological tool for recording past sea level changes and reconstructing detailed paleoclimatic and paleoenvironmental conditions (e.g., [8,9,10,11,12,13]). Likewise, the composition of the ostracod fauna captures in great detail any environmental change, especially in the transitional environments of coastal ecosystems, allowing the assessment of the potential impact of climate change [4,14]. The ostracod analysis constitutes a significant tool of the paleoenvironmental studies in coastal regions of Greece of both geoarcheological and geomorphological interest [8,12,15,16,17,18]. In addition, the analysis of mollusc assemblages is another good proxy for the provision of detailed records of paleoenvironmental and paleoclimatic changes [19], thus, being a valuable tool in the reconstruction of the Holocene paleoenvironments [20,21,22,23]. Further, coastal sedimentary deposits are rich in palynomorph content, namely pollen, dinoflagellates cysts and other “Non Pollen Palynomorphs” (NPPs) providing valuable additional insights on both the terrestrial and the aquatic realms (e.g., [24,25,26,27]).

Overall, microfossil studies, often complemented with palynological analyses, performed on several Holocene Aegean coastal plains and coastal lagoons [8,9,12,13,18,28,29], have provided a plethora of information concerning the environmental impacts and adaptation practices of human societies during the last few millennia (e.g., displacement of settlements, main agricultural and pastoral activities).

The application of environmental magnetism techniques performed on rock and minerals, is used to investigate the formation, transportation, deposition and post depositional alterations of magnetic minerals under the influence of a wide range of environmental processes. All iron-bearing minerals are sensitive to a range of environmental changes, which makes the magnetic measurements extremely valuable for detecting signals associated with environmental processes [30]. In particular, magnetic susceptibility measurements can be used as a rapid, simple and non-destructive proxy to indicate magnetic changes in the sediments corresponding to different depositional environments [31,32,33,34,35].

Saronikos Gulf is located at the western edge of the central Aegean Plateau, near the city of Athens. Its northern part, named inner Saronikos Gulf, receives the treated wastes of ~5 million people from the Wastewater Treatment Plant outfall that discharges south of Psittalia Island. The enclosed embayment of Elefsis Bay, connected to Saronikos Gulf by two narrow and shallow straits (an eastern and a western one), displays a physical setting with low freshwater inflows and limited water exchange, leading to strong seasonal density stratification, high nutrient accumulation and poor environmental status [36,37,38,39,40]. The sedimentary sequences deposited between the hill of Piraeus and the floodplain of Kifissos River (Attiki peninsula, Greece) constitute the Piraeus coastal plain in the inner Saronikos Gulf. The location of Piraeus coastal plain close to the metropolitan city of Athens makes it an excellent study area for tracing the human shaping of the landscape during Prehistory and Antiquity until the modern times. The history of Piraeus coastal landscape dates back to Early Holocene [28], with the main factors that featured its evolution being the low tidal range and the relative sea level rise configured by the glacio-hydro-isostatic changes, tectonic stability of the area (documented by the relative absence of earthquakes during the last few thousand years) and progradation of the deltaic fan of Kifissos River.

The paleoenvironmental evolution of the Piraeus coastal plain has already been presented in the study of [28], mostly based on the interpretation of benthic foraminiferal assemblages, and it has been implied that between 6800 and 5400 cal BP Piraeus was an island in the center of a shallow marine bay. However, the present study aims to provide additional interpretations concerning the paleoenvironmental evolution of the Piraeus coastal plain since ca. 9000 cal BP, using a multi-proxy analysis and synthesis of the benthic foraminiferal and ostracod assemblages, molluscan populations, palynological content and magnetic susceptibility in the local Holocene coastal deposits. Further, GIS tools (e.g., [41]) were used to define, analyze and visualize the spatial extent of the paleoenvironmental conditions in each stratigraphic interval and, finally, to reconstruct the environmental evolution of the area mainly after the Early Holocene. Finally, the comparison of the present study’s record with a well-dated marine sedimentary record recovered from the nearby shallow Elefsis Bay, allowed a reliable estimation of the sea level rise rate in the coastal zone of the metropolitan city of Athens during the Mid-Holocene. It should be emphasized that estimates of the sea level rise based on the analysis of local sedimentary archives corresponding to warm time intervals of the ‘recent’ past, analogous to the present time, are more reliable indicators for the projection of the future sea level rise in the local coastal areas of concern compared to the global predictions made by the modern climate models, because the resulting sea level rise will not be globally uniform, due to ocean dynamic processes and changes in gravity associated with water mass redistribution.

2. Study Area

The Piraeus coastal plain is part of the drainage basin of Saronikos Gulf, which is a semi-enclosed embayment located between the southwest coast of Attiki and the northeast coast of Peloponnese (Figure 1). The Piraeus peninsula, facing at its western part the Elefsis Bay, is part of the Athens Basin and is structured by Pliocene marine post-alpine formations with marls, sandstones, conglomerates and limestones. The alpine basement of the site is considered as one of the most aseismic areas in Greece with limited neotectonic activity [42,43]. The present geography of the area began to form in the Middle Pleistocene with the propagation of Kifissos River [44], while recent post-alpine sediments cover the foothills of the mountains and areas of lower altitude [43,44]. The lithostratigraphy of the Piraeus area [45] consists of the Upper Cretaceous limestones, the ‘Marls of Piraeus’ and the unconformably overlying Holocene deposits. The deltaic fans of Kifissos and Korydallos rivers were crucial for the formation of the Piraeus coastal plain [28]. The alluvial sedimentation, mainly from Kifissos River, led to the formation of Halipedon (salt field) in the area of Phaliron, a marshy area that covered a great part of the region, still existing in the first decades of the 20th century [46,47].

3. Materials and Methods

Ten rotational boreholes (P1–P10; 10 cm in diameter) were drilled in the Piraeus coastal plain (see Figure 1) and the Holocene stratigraphy of the recovered sediments (Units A–E) was recorded in detail, based on 42 calibrated accelerator mass spectrometry radiocarbon ages [28,48]. The most complete sequences (boreholes P2, P4 and P5; Figure 2) were selected for further investigation as their location is crucial for the understanding of the coastal plain paleoenvironmental evolution and the dominant processes that affected the landscape during Holocene. More than two hundred and eighty samples were studied from selected sedimentary layers of the P2, P4 and P5 borehole sequences, undertaking grain size analysis, benthic foraminiferal, ostracod, molluscan and palynomorph assemblage studies and magnetic susceptibility measurements.

Grain size analysis of samples from the boreholes P2, P4 and P5 was carried out using a Malvern Mastersizer 2000 laser scanner [48].

Magnetic susceptibility measurements were performed on two hundred and fifty samples from all selected boreholes. The samples were initially sieved in order to remove all impurities and, then, packed in cylindrical plastic boxes (of 8 cm³). The laboratory measurements of the volume-specific magnetic susceptibility relative to air (κ, SI units) were carried out using the Bartington MS2B sensor at a low frequency of 0.465 kHz. The samples were weighed prior to the measurements in order the acquired results to be expressed as mass-specific magnetic susceptibility (χ, 10⁻⁸ m³/kg). Each sample was measured at least 3 times and the average value was considered as the representative one for the sample.

Micropaleontological foraminiferal analyses of sixty-four samples from borehole P2, one hundred and eleven samples from borehole P4 and forty-two samples from borehole P5 were carried out. An amount of 10 g (dry weight) from each sample was treated with hydrogen peroxide (H₂O₂) to remove the organic matter and, subsequently, washed through a 125 μm sieve and dried at 70 °C. When feasible, a subset containing ~200 benthic foraminifera for each sample was obtained using an Otto microsplitter. The microfauna was identified under a Leica APO S8 stereoscope, following the generic classification of [49,50], while the foraminiferal density (FD; number of specimens/g) and the relative abundances in the benthic foraminiferal assemblage of each sample (%) were calculated. In addition, the PAST 2.12 software package [51] was used for the determination of the Shannon-Wiener diversity index (H′). Further, a set of proxy parameters for environmental interpretations were also calculated: (i) the ratio between large (L) and small (S) Ammonia tests, established as the A index [A = 100*L/(S + L)], used for the assessment of environmental conditions [9,52]; (ii) the P ratio (%) [P/(P + B −S)*100], representing the percentage of planktonic species (P) in the sum of planktonic and benthic (B) species in the foraminiferal assemblages discarding the stress-tolerant taxa (S), used for the identification of marine environments [53,54,55] and, finally, (iii) the BR ratio [(broken-reworked specimens/undamaged specimens + broken-reworked specimens)*100], based on the abundance of broken and abraded foraminiferal tests, used for the better evaluation of the high-energy paleoenvironments (e.g., coastal zone, shoreface), since, for example, a high number of broken foraminiferal specimens strongly reflect an upper shoreface environment (e.g., [56,57]).

A detailed quantitative and qualitative ostracod analysis was performed on fifty-nine samples from borehole P2, one hundred and eleven samples from borehole P4 and forty-one sample from borehole P5. Due to the observed great variability in the ostracod abundances among different samples, all ostracod specimens were collected from the >125 μm fraction. In cases where the specimen abundance was too high (thousands/g), aliquots were examined using a microsplitter and, subsequently, the total number of specimens was calculated for the purposes of the statistical analysis. Complete carapaces were counted as two valves, while all ostracods were identified in the genus and species level by using a stereomicroscope and a scanning electron microscope (Jeol JSM 6360). The ostracod specimens were classified according to Horne et al. [58] and their identification was based on several previous studies such as [17,59,60,61,62,63]. Finally, the relative abundances of each species were estimated for each sample and distribution diagrams per borehole were constructed for the most abundant taxa, while the PAST 2.12 software package [51] was used for the determination of the H′ index.

Molluscan analysis was performed on sixty samples from borehole P2, one hundred and eleven samples from borehole P4 and sixty-three samples from borehole P5. The residue of the 125 μm sieve was examined under a Leica APO S8 stereoscope. Molluscan assemblages were identified and counted in a semi-quantitative approach: (i) common (C) >30 specimens/10 g; (ii) rare (R) 10–30 specimens/10 g; and (iii) present (P) <10 specimens/10 g. In addition, other microscopic components (e.g., charophyta, spores, ceramic fragments, gypsum crystals) were recorded.

Palynological analysis was conducted using fifty-two samples from borehole P4, even though only in twenty-eight of them the pollen concentration was sufficient to be included in the present study. All examined samples were spiked with a known quantity of Lycopodium spores, chemically treated with hydrochloric (HCL 37%) and hydrofluoric (HF 40%) acids, acetolysed and, finally, sieved through a 10 μm sieve. The residues were mounted in silicon oil. Pollen identification was based on the key of Beug [64] and the Reille atlas [65], while for the identification of non-pollen palynomorphs (NPPs) the studies of van Geel et al. [66,67] were used. The pollen preservation was good, while the total pollen concentration ranged from 2700 to 32 grains per gram of dry sediment. The samples with concentrations lower than 150 grains/g were considered barren and excluded from this study. Finally, a percentage pollen diagram was constructed for borehole P4, based on a pollen sum of terrestrial pollen grains, excluding aquatic and hygrophilous pollen and spores. The fluvial input was determined based on the sum of the erosion indicating NPPs the mycorrhizal fungus Glomus and the algal remain Pseudoschizaea [68].

Eventually, in order to visualize the temporal and spatial evolution of the Piraeus coastal plain after ca. 9000 cal BP, the obtained paleoenvironmental results were grouped according to (i) specific stratigraphic intervals; (ii) borehole location; and (iii) main paleoenvironmental conditions. The broader coastal plain area, presenting slopes less than 2°, was mapped as an alluvial fan depositional system. The slope dataset was created using the ArcGIS Pro software, incorporating a 5-m-resolution digital elevation model derived from the area’s orthophoto map (National Cadastre and Mapping Agency, 2012), while the 1:50,000 geological map sheet “Athens-Piraeus” (Institute of Geology and Mineral Exploration, 1982) was used to digitize the extent of the alluvial deposits. The spatial extent of the paleoenvironments for each defined time interval was estimated using the ordinary kriging geostatistical interpolation method adopting the spherical semi-variogram model and the Jenks natural-breaks classification method. Finally, to better visualize the results, a 3D local scene for each time interval was created, combining the shaded relief of the area with the analysis results, the drainage network, the present-day coastline and the location of the boreholes.

4. Results

4.1. Borehole P2

4.1.1. Lithology and Magnetic Susceptibility

Borehole P2 is located near the present coastline, drilled to a maximum depth of ~20 m (Figure 1 and Figure 2). The most common sediment fractions of the sequence are sand and silt. Clay forms layers, often including shells and plant remains, in the depth interval of 17–10 m (see Unit A in Figure 2), while an evident sand layer is delimited between 7 and 6 m (see Unit B in Figure 2). Finally, medium to coarse sand layers (occasionally including pebbles) prevail up to 3 m core depth [48], associated with the lithostratigraphic Unit E of Goiran et al. [28].

Magnetic susceptibility ranges from 1.6 to 37.4 × 10⁻⁸ m³/kg^. in the lower part of the sequence (17–13 m), displaying enhanced values. Within 13–6 m depth, the magnetic signal recorded is up to 10 × 10⁻⁸ m³/kg, while the lowest value throughout the borehole length appears at the depth of 6 m. Finally, the higher values are observed in the upper part of the borehole in the 4–3 m depth interval (Figure 2).

4.1.2. Benthic Foraminifera

Throughout the P2 sequence, FDs vary between 0 and 448 specimens/g. Very low values are exhibited in the lower portion of sandy and clayey layers and in layers of coarse sand and pebbles within the upper part of the sequence (up to 20 specimens/g), while the rest of the studied samples are characterized by an average FD of 100 specimens/g (Figure 3a). In total, 57 foraminiferal species were identified, including 30 hyaline and 27 porcellaneous specimens. The diversity index H′ varies between 0.45 and 3.06, with the higher values recorded within the 10–7 m depth interval, coupled with maxima (>4%) of the P ratio distribution. Finally, the A ratio profile displays prominent peaks at 13.5 m and 6 m depth, while the maximum ΒR ratio (100%) is marked within the upper part of the P2 sequence (4.5–4 m) (Figure 3a).

From the bottom of the borehole up to ~10 m, the benthic foraminiferal assemblages are marked by the presence of Ammonia tepida (Figure 3b). In particular, within 17–13.5 m depth, the foraminiferal tests are practically absent, with only sporadic specimens belonging to the tolerant to restricted conditions species Haynesina germanica. In between 13.5–10 m core depth the foraminiferal abundance increases (Figure 3a) and the assemblages are primarily composed of A. tepida, which gradually dominates (up to 88%), H. germanica (max 50%), A. beccarii (max 75%) and Bolivina spp. (max 50%), accompanied by Cribrelphidium gunteri (max 9.2%), Aubignyna perlucida (up to 8.6%) and miliolids (Figure 3b). Within the overlying sediment (10–7 m), the assemblages demonstrate their highest diversity, with their composition considerably differentiated from that of the underlying interval. They mainly consist of Rosalina bradyi (max 13.3%), Bolivina spp. (max 15.4%), Elphidium complanatum (max 14.8%) as well as miliolids (up to 40.7%) that are represented by Quinqueloculina seminula, Q. parvula, Q. berthelotiana and Q. bicarinata. Other species such as Asterigerinata mamilla, C. gunteri and Buccella frigida generally represent less than 10% of the foraminiferal fauna. However, A. tepida, H. germanica and A. beccarii, continue to constitute part of the occurring assemblages. Eventually, in the upper part of the borehole (7–3 m) (Figure 3b), a decrease in the foraminiferal abundance is recorded, with the populations mostly represented by A. beccarii (21% on average) and Q. seminula (6% on average).

4.1.3. Ostracods

A total of 15 species were identified, while twenty-one samples, mostly retrieved from the 6.4–3.35 m depth interval, were barren (Figure 3b). The ostracod assemblages present a low number of taxa, usually 3–6 species per sample, that is normal for shallow marginal marine environments [4], which, subsequently, affects the diversity index H′ showing very low values (0.2–1.6) (Figure 3a). The ostracod populations from the bottom of the sequence up to 10.15 m are characterized by the dominance of Cyprideis torosa (up to 100%), being along with Loxoconcha elliptica (up to 50%), Cyprinotus salinus (0.8–50%) and Xestoleberis communis (1.8–21.4%). In between 17–13.5 m depth, the ostracods are featured by few specimens, whereas the upper part demonstrates the highest ostracod species density in terms of valves per sample (Figure 3a). However, upwards a decline of the ostracod density is recorded, marked by the maximum abundances of C. salinus (50%) and Cytherois fischeri (8.3%). Within the overlying sediment layer (10.05–8.7 m), C. torosa still remains the dominant species, forming in most cases monospecific assemblages (Figure 3b). Between 8.6 and 7.1 m, the ostracod fauna mainly comprises Pontocythere turbida (28.6–50%), L. rubritincta (14.3–50%), L. elliptica and L. affinis. Finally, the upper part of the sequence (7–3 m), consisting of medium sand and coarse sand, is characterized by the negligible occurrence of ostracods. Only one sample at 4.7 m contained a poor ostracod assemblage, including C. torosa, C. salinus and L. elliptica.

4.1.4. Mollusc Assemblages

Τhe lower part of the sequence, below 13.6 m, is barren of molluscan specimens (Figure 3c). The major part of the mollusc fauna was recorded in the 13.6–7.1 m borehole depth, including several shells principally belonging to Cerastoderma glaucum (small-sized and juvenile forms), Abra sp. and Hydrobiidae, followed by Loripes lacteus and Bittium reticulatum. In contrast, the topmost portion of this interval (8.3–7.1 m) displays limited mollusc fauna marked by the presence of a few thin-shelled Cerastoderma spp. fragments. Finally, the upper 7 m of the sequence, in particular the 6.8–3.24 m depth interval, is characterized by scattered highly reworked marine mollusc shell fragments of Cerithiidae and Cerastoderma spp. as well as by indeterminable fragments of possible land snail shells.

4.2. Borehole P4

4.2.1. Lithology and Magnetic Susceptibility

Borehole P4 is located north of borehole P2, drilled into the mainland of Piraeus peninsula and being further away from the present coastline. It corresponds to a more complete sedimentary archive compared to the P2 sequence (Figure 1 and Figure 2), since it includes the majority of the lithostratigraphic units identified in the Piraeus coastal plain by Goiran et al. [28].

In Unit A, conglomerates, occurring at the 17–16.75 m depth level, are mainly overlain by mainly mixtures of silt and clay up to 10 m (Figure 2). In Unit B, silty sand with shells prevails between 9.5 and 6 m (Figure 2), while silt with sand and clay prevails from 5.42 to 5.12 m in Unit C. Finally, within Unit D, pure clay forms a layer extending up to 3 m, overlain by a deposit of rubbles [48].

The magnetic susceptibility profile of borehole P4 shows a great variability throughout the sequence, ranging from almost a zero value to 100 × 10⁻⁸ m³/kg (Figure 2). From the bottom of the borehole up to 12 m, magnetic susceptibility displays an enhanced mean value of 41 × 10⁻⁸ m³/kg, while moderate values characterize the 12–9.5 m interval. From the previous depth level to 6 m, a reduction of the magnetic susceptibility is observed reaching the lowest values throughout the P4 sequence (mean of ~4 × 10⁻⁸ m³/kg). Finally, up to 4.15 m, magnetic susceptibility remains at low levels and, then, it increases again up to the top, demonstrating the highest values in the sequence (mean of 68 × 10⁻⁸ m³/kg).

4.2.2. Benthic Foraminifera

Even though relatively low FDs characterize the P4 sequence (max 341 specimens/g), higher values with an average of 65 specimens/g occur in the sandy interval between 9 and 6 m (Figure 4a). In total, 63 foraminiferal species were identified, including 41 hyaline, 20 porcellaneous and 2 agglutinated specimens. The diversity index H′ varies between 0.31 and 2.85, with the higher values mainly observed in the 9.5–6 m interval (Figure 4a). Within the previous depth range, the P (up to 4%) and A ratios (50% on average) display their maximum values, while the BR ratio shows relatively low values, not exceeding 20% (Figure 4a).

In the lower part of the P4 sequence (17–12 m), the foraminiferal fauna is dominated by A. tepida, displaying abundances higher than 50%, and H. germanica, occasionally reaching abundances of 60% (Figure 4b). In the overlying interval (12–9.5 m), A. tepida is still the prevailing species (70% on average), however, H. germanica decreases (max 30%) in respect to the previous interval, whereas C. gunteri (max 30%) and A. perlucida (up to 20%) are displaying increased values. Within the sandy sediment layer between 9.5 and 6 m, the foraminiferal content changes with A. beccarii and miliolids comprising up to 70% of the assemblage. Among miliolids, Q. seminula is the most abundant species (max 50%), followed by Q. berthelotiana (max 15%), Q. bicarinata (up to 15%) and Q. parvula (max 8%). Moreover, other typical marine infralittoral species like E. complanatum, R. bradyi, A. mamilla and B. frigida also occur, representing less than 25% of the foraminiferal fauna. In the upper part of the sediment occurring between 6 and 4 m, the strong dominance of A. tepida (70% on average) along with the presence of Q. seminula and H. germanica document an oligospecific assemblage, whereas from 4 m to the borehole top, the foraminiferal species are completely absent (Figure 4b).

4.2.3. Ostracods

A total of 21 species were identified, whereas 23 samples were barren of specimens. Generally, the ostracod assemblages present low number of species, usually composed of 2–8 specimens per sample, with the diversity index H′ displaying low values, ranging from 0.1 to 1.6 (Figure 4a). The micropalaeontological analysis reveals that the changes in the composition of the ostracod populations follows the lithological variation along the P4 sequence. C. torosa is the dominant species over a major part of the sequence, i.e., 17–9.5 m (Figure 4b). In the majority of the samples retrieved from the previous interval, the frequencies of the prevalent C. torosa fluctuate between 60% and 99%. The rest of the fauna occurring, in general, in the 17–12 m depth range, consists of Ilyocypris spp. (I. bradyi and I. gibba; max 64.4%), Limnocythere inopinata (max 27.8%), C. salinus (0.6–33.3%) and Sarscypridopsis aculeata, excluding, however, the 15.5–15.6 m and 14.1–14.4 m intervals, where the previous species are replaced by X. communis, L. elliptica, C. fischeri, C. salinus and Leptocythere spp. Upwards, between 12 and 9.5 m, C. torosa is mainly accompanied by X. communis (up to 86%), L. elliptica (up to 37%), C. fischeri (up to 7.3%), C. salinus (up to %) and rare Leptocythere lagunae. The middle part of the borehole P4 (9.5–6 m) is a sandy deposit, characterized by a great reduction in the number of specimens per sample, while many samples were barren of specimens. At the lower half of this interval, ostracod assemblages are composed of C. torosa (max 100%), L. elliptica (max 67%), X. communis (max 100%) and rare Cytheromorpha fuscata. Upcore, sporadic samples contained differentiated assemblages including Semicytherura incogruens, L. affinis, L. rubritincta, Xestoleberis spp. and Aurila woodwardii. In contrast, ostracods are found in extremely high numbers within the overlying sediment (6–4 m), where the assemblages are dominated by C. torosa (more than 90%) and complemented by Ilyocypris spp. (I. bradyi and I. gibba; up to 25%), L. inopinata (max ~6%), C. salinus (1.7–33.3%) and rare S. aculeata and Neglecandona neglecta. Finally, the assemblages show a striking change in the uppermost part of the P4 sequence (from 4 m to top), characterized by Ilyocypris spp. (I. bradyi and I. gibba; up to 60%), N. neglecta (5.5–37.5%), C. salinus (up to 45%) and rare L. inopinata.

4.2.4. Mollusc Assemblages

C. glaucum, Abra sp., Hydrobia sp. and a few Bittium reticulatum and Parvicardium sp. constitute the macrofaunal assemblage in the lower part of the sequence (17–12 m). C. glaucum, Abra sp. and Hydrobia sp. along with B. reticulatum and a few Rissoa sp. and Loripes lacteus constitute the characteristic population of the 12–10 m depth range (Figure 4c), while from 10 to 6.7 m, B. reticulatum, L. lacteus, small Cerithium sp., fragments of small C. glaucum, Hydrobiids, Rissoa sp., calcareous tubes of Spirorbis sp. worm and sponge spicules as well as very few Donax sp., Parvicardium sp., Gibbula sp., Tricolea sp. and Retusa sp. are present. Upwards (6.7–4 m) the sandy layers contain plenty of shells of Hydrobiids, juveniles and fragments of C. glaucum, Abra sp. and a few Cerithium sp. and Bittium sp. fragments. Finally, the clay layer between 4 and 3 m shows a lack of molluscan specimens (Figure 4c).

4.2.5. Palynological Analysis

Within the lower part of P4 (17.00–9.50 m), Arboreal Pollen (AP), comprise the ~30% of the pollen assemblage (Figure 4d). Pinus is the most common tree, followed by deciduous Quercus and other deciduous tree taxa. Abies shows a profound short-lived peak that coincides with an increase in Artemisia. Poaceae (20%) and Amaranthaceae (15%) are the main herb taxa recorded, while Amaranthaceae fluctuate exhibiting increased abundances in the base and around 12m depth. Increased Cerealia-type is recorded between 12.00 and 9.50 m. The occurrence of several dinoflagellate cyst species and foraminifera test linings, as well as, Pseudoschizaea and Glomus has been documented (Figure 4d). The interval between 9.50–6.50 m is devoid of palynological content. From 6.5 to 4.00 m maxima abundances of the halophytic Amaranthaceae (55%) are recorded, while other herb taxa such as Asteroideae, Cichorieae and Poaceae abundances appear abridged. All deciduous tree taxa are also reduced, while Pinus and Olea are the main tree species during this interval. An increasing trend in Cerealia-type as well as in the coprophilous fungal remains is observed. Dinoflagellate cysts are scarce and foraminifera linings are absent. The upper part of the core (4.00 m-top) the maxima values in Olea, Cerealia-type and increased values in Cichorieae are recorded. Coprophilous fungi, and Zygnemataceae record their maxima.

4.3. Borehole P5

4.3.1. Lithology and Magnetic Susceptibility

Likewise with borehole P4, the drilling site of borehole P5 is located on the mainland of Piraeus peninsula, being northwest of borehole P2 (Figure 1). In Unit A, the sediment layers between 15 and 10.5 m comprise clay and silt (Figure 2). In Unit B, silty sand and silt prevail from 10.5 to 5.8 m, while in Unit C, up to 4.5 m, a silt-dominated deposit occurs, including shells and shell fragments (Figure 2). Finally, in Unit D, sandy clay layers prevail, overlain by rubbles.

The magnetic susceptibility measurements performed on the P5 sequence reveal a moderate variation throughout the borehole length. The most enhanced values are recorded from the bottom to 9.7 m, particularly below the 12 m depth level, and the observed peak of 120 × 10⁻⁸ m³/kg at 13.5 m (Figure 2) is rather related to the presence of ceramics sherds. The interval between 9.7 and 4.5 m is characterized by very low to null magnetic susceptibility, a sudden upcore increase is noticed within the overlying sediment (Unit D), with the maximum intensity recorded at 4 m (see Figure 2).

4.3.2. Benthic Foraminifera

The highest FDs (up to 1843 specimens/g), among the three investigated boreholes, occur in the P5 sequence (Figure 5a). A total of 73 benthic foraminiferal species were identified, including 40 hyaline and 32 porcellaneous specimens, and 1 agglutinated species. The diversity index H′ varies between 0.26 and 3.12, with the higher values observed, in particular, between 11.5 and 5.9 m. Finally, the P ratio shows maximum values (4–6%) at 8.4 and 7.4 m, the A ratio presents relatively low values (generally less than 40%) and the BR ratio profile demonstrates a prominent peak of 50% at 14.2 m (Figure 5a).

In the lower part of the P5 sequence, between 15 and 11.5 m (Figure 5b), the foraminiferal fauna is dominated by A. tepida (50–100%), followed by C. gunteri (max 24%), H. germanica (max 15%), A. perlucida (max 11%) and Q. seminula (up to 10%). Other epiphytic species, such as E. complanatum and R. bradyi, display discontinuous distributions and, generally, represent less than 10% of the foraminiferal fauna (Figure 5b). Above the previous interval, between 11.5 and 5.9 m, A. beccarii and miliolids (represented mainly by Q. seminula and Q. bicarinata) show increased relative abundances (up to 50%). Within the same interval, several taxa like E. complanatum (up to 15%), B. frigida (max 10%), R. bradyi (max 9%) and Bolivina spp. (max 9%) display their maximum percentages. Upwards, between 5.9 and 1.5 m, A. tepida increases again (being more than 80%), followed by other common species such as C. gunteri (max 10%) and A. perlucida (max 10%).

4.3.3. Ostracods

A total of 31 species were identified, whereas three samples that contained a very low number of valves (3–5; see Figure 5a), belonging to juveniles, were considered as barren of ostracod specimens. In general, the ostracod assemblages displayed a low number of species, usually 3–6 species per sample. The ostracod diversity, based on the H′ index, appears higher than that determined for the P2 and P4 sequences, ranging from 0.112 to 2.091, presenting, however, its lower values between 5.8 and 4.5 m (Figure 5a).

The micropalaeontological analysis revealed that C. torosa is the most abundant species over the largest part of the P5 sequence, but the composition of the total ostracod assemblages varies in terms of the accompanying species (Figure 5b). In particular, the ostracod fauna from the bottom to ~11.44 m (where the highest density of specimens in terms of valves per sample was observed) is characterized by high frequencies of C. torosa (43–96%). It is followed by L. elliptica (occurring in all examined samples, up to 38%) and C. salinus (up to 27%), as well as, by X. communis (0.1–33.3%), mainly being present in the uppermost part of the interval. In the 11.34–5.85 m depth range, the ostracod populations become more diverse, with the greatest diversity and highest number of species observed in the 9.96–9.66 m and 7.8–5.85 m intervals. C. torosa is the most abundant species with relative abundances of 22–94%, occurring along with C. salinus (up to 33.3%), L. elliptica (up to 26%), X. communis (up to 67%), P. turbida (1–21.2%), L. rubritincta (0.3–16.9%), L. affinis (0.3–7.1%) and Leptocythere spp. (0.4–7.7%). Within the previous interval, at the level of 8.52–8.43 m, the ostracod assemblage temporarily changes, presenting an oligospecific character, while at the top the presence of S. aculeata, Ilyocypris spp. and L. inopinata is significant (Figure 5b). The 5.8–4.5 m interval presents almost monospecific ostracod assemblages, having as dominant species C. torosa (more than 95%), accompanied by C. salinus, S. aculeata and Ilyocypris spp. occurring at low frequencies. Finally, the upper part of the P5 sequence is characterized by oligospecific and low abundance communities, excluding the depth intervals of 3.95–4 m and 1.74–1.81 m, where there is an increase in ostracod density. These assemblages are composed of C. torosa accompanied by C. salinus, L. elliptica and X. communis.

4.3.4. Mollusc Assemblages

C. glaucum, Abra sp., Hydrobia sp. characterize the lowest part of the sequence (15–13.08 m), while the presence of Nucula sp., B. reticulatum, Parvicardium sp., Rissoa sp., Tricolea sp., Retusa sp., Spirorbis sp., a few Donax sp., Gibulla sp. and fragments of Tapes decusatus is worth mentioning (Figure 5c). A comparable assemblage continues up to 13.08 m, with a small increase of Nucula sp., B. reticulatum, Rissoa sp., Tricolea sp., Retusa sp. and Spirorbis sp. observed at 11.3–10.9 m. C. glaucum, Abra sp. and Hydrobia sp. were identified in the fine-medium sandy layers of the 9.84–5.85 m interval, accompanied by B. reticulatum, Bornia sebetia, Rissoa sp. and very few Tricolea sp., Turbonilla sp. and Eulimella sp. The mud of the 5.8–5.5 m interval is featured by C. glaucum, Abra sp. and Hydrobia sp., but also contains Planorbis sp. and charophytes. Finally, very small fragments of little thin-shelled gastropods, possibly land snails, are included within the sandy clay of the 4.55–3.95 m interval, while upwards, in the 3–1.7 m depth range, C. glaucum, Abra sp., Hydrobia sp. and a few Mytilaster sp., B. reticulatum and Retusa sp. are worth mentioning (Figure 5c).

5. Discussion

5.1. Paleoenvironmental Interpretations

The boreholes P2, P4 and P5, drilled into the Piraeus coastal plain of Athens Basin in the fluvio-deltaic area of Kifissos and Korydallos rivers (Figure 6), unravelled the paleoenvironmental conditions prevailing during the deposition of the lithostratigraphic units (Units A–E) of the relevant sedimentary sequences. The multi-proxy analysis provided the opportunity to reconstruct the paleogeographic evolution of the Piraeus coastal plain and the associated inner Saronikos Gulf during Holocene. It should be mentioned that dating and lithostratigraphic correlations referred to in the present study are presented after Goiran et al. [28].

5.1.1. The First Stage of the Piraeus Lagoon (Unit A)

The deposition of Unit A occurred between 8700 and 7500 cal BP [28]. Even though this deposit has been identified within the lower part of all analyzed sedimentary sequences, its paleoenvironmental characteristics do not demonstrate a spatial consistency. Hence, Unit A of [28] can be further divided into two subunits with different paleoenvironmental characteristics ([70]; present study).

The older subunit is identified in the borehole P4, including benthic foraminiferal assemblages dominated by the euryhaline species A. tepida, which reflects wide-ranging salinity and temperature conditions met in nearshore marine environments, lagoon systems and deltaic zones [52,71]. This together with the simultaneous enhanced presence of the tolerant to restricted conditions H. germanica, ([72,73], within the clayey silt layers formed between 17 and 12 m, assigns the older subunit of Unit A to mesohaline to oligohaline biofacies. This is actually the common biofacies recorded in the closed lagoons of the Aegean region [9,22,36]. The previous interpretation is further supported by the characteristics of the abundant ostracod species, constituting the assemblages of this stratigraphic interval. These are dominated by the euryhaline C. torosa, which occurs together with I. bradyi and I. gibba. The latter two species are widespread in all freshwater environments [74], mostly preferring permanent shallow water bodies with clayey, muddy or sandy substrates [61], with I. gibba also recorded in shallow oligohaline and brackish coastal environments [13,75,76]. Moreover, the accompanying freshwater to oligohaline species such as L. inopinata inhabit mainly shallow, fresh to oligohaline water bodies but found also in oligohaline inland coastal waters [61] whereas S. aculeata, is common in coastal rockpools and tolerates salinities up to 17 [61,77]. Finally, in addition to the previous evidence, the associated molluscan assemblages, consisting mainly of C. glaucum, Abra spp. and a few Hydrobiidae, further reflect the establishment of a typical meso-oligohaline lagoonal environment [22,78,79,80] during the deposition of the older subunit of Unit A.

The dominance of the coastal euryhaline dinoflagellate Lingulodinium machaeropholum, and the scarcity of other species supports the deposition in a lagoon environment with terrigenous influx [81,82]. Between the bottom of the P4 sequence and 12 m (Figure 4d and Figure 7), the observed enhanced presence of Pseudoschizaea and Glomus palynomorphs, both regarded as indices of soil erosion [27,68] and increased fluvial runoff, and the simultaneous relative increase in magnetic susceptibility (~40 × 10⁻⁸ m³/kg on average), further evidence the existence of closed lagoon conditions [30,33], implying increased freshwater inputs. Nevertheless, the occurrence of numerous gypsum crystals at the 10.65–10.6 m, 10.45–10.4 m and 10.25–10.1 m depth intervals may be linked with short-term desiccation episodes in the lagoon.

Therefore, based on the above interpretations, the lower part of Unit A can be defined as the result of a closed lagoonal environment, named Paleoenvironmental Unit Aa (see Figure 7). This depositional system was restricted in the eastern part of the Piraeus coastal plain (Figure 8a). As the paleo-relief comprising the Pliocene substrate had formed a ~15 m topographic high with reference to the deepest part of the plain at −17 m (recorded at the bottom of borehole P2, most probably corresponding to the riverbed of Korydallos River; see Figure 8c and [28]), it hinders a further expansion of the lagoon to the west. This is in agreement with the predicted from the glacio-hydro-isostatic model of Lambeck [84] sea level stand at ca. 8000 yr BP (−17 m) for the broader area of the Piraeus coastal plain. All evidence imply that the established lagoon and the associated marine impact on the eastern part of the plain might have been achieved either due to restricted connection with the open sea via channels in the fluvio-deltaic area of Kifissos and Korydallos rivers backwater transition zone (see Figure 6 and Figure 8, [28,69]) or due to salinization processes of the unconfined aquifer [85].

Unit A displays relatively different features in the P2 and P5 sequences and at its upper part in the P4 sequence (Figure 7). The benthic foraminiferal analysis reveals once more the dominance of A. tepida, together with the concomitant presence of H. germanica, C. gunteri and A. perlucida, which are typical species of estuarine and shallow marine environments [9,10]. This assemblage marks the mesohaline biofacies, which has already been identified in the open lagoons of the Aegean region [9,22,36]. This is in accordance with the domination of the ostracod C. torosa in, which is a highly euryhaline species that is primarily associated with shallow (<30 m) environments and salinities ranging from 0.4 to 15 [86,87], even thriving under salinity fluctuations ranging from 2 to 17 [88]. Furthermore, the accompanying species L. elliptica, C. salinus and the mud-dweller C. fisheri, which represent common fauna of the brackish waters [75,76,86,89,90,91], provide additional evidence for the establishment of mesohaline conditions [8,12,92]. In agreement with the previous interpretations, the molluscan fauna of this deposit includes similar species with those occurring in the paleoenvironmental Unit Aa, but their considerably increased abundance is evidencing a better communication with the sea. Eventually, the presence of diverse dinoflagellate cysts and foraminifera linings (Figure 4d) indicates the marine influence on the depositional system [24].

The documented mesohaline open lagoon conditions for the sediment of the upper part of Unit A, named Paleoenvironmental Unit Ab (see Figure 7), are also supported by the evident presence of clay, organic material and the limited sand content [45,48]. Further, Unit Ab exhibits lower magnetic susceptibility in relation to Unit Aa, signaling the marine influence on the established open lagoon conditions. Interestingly, the low frequencies of shallow marine foraminiferal species, such as A. beccarii, Bolivina spp. and E. complanatum [93], combined with the scarce presence of the mostly marine, epiphytal littoral ostracod X. communis [89,94] and the limited occurrence of the mainly epiphytal L. affinis that inhabits shallow marine environments [62,94,95], provide evidence of increased marine impact, while the presence of gypsum crystals (e.g., between 11.7 and 10.4 m in the P2 sequence) may reflect partial desiccation conditions. However, the topographic high of the paleo-relief (see above) was still obstructing the lagoon to extend towards the western part of the coastal plain. Hence, likewise with Unit Aa, Unit Ab was restricted in the eastern part of the plain (Figure 8b).

Despite the generally low pollen concentrations of the deposits, which resulted in a fragmented record, the vicinity of Piraeus coastal plain to the City of Athens offers valuable insights on the plant landscape and human activities in the area during Prehistory and in the Antiquity. Pollen assemblages of P4 denote the predominance of open vegetation around the site. In a regional scale a mosaic plant landscape with open pine Mediterranean evergreen woodlands and mixed deciduous oak forests further inland is evidenced. The pollen record displays a short-lived peak of Abies and increase in Artemisia abundances in between 14.50 and 14.00 m, that may be related to a climatic deterioration event during the Early Holocene. Following this interval, the thermophilous deciduous woodlands, mainly composed of deciduous Quercus and Carpinus/Ostrya further developed in agreement with other vegetation records from south Greece [96,97,98]. Unit A corresponds to the Early Neolithic period, yet there are no clear indications about human impact on vegetation. However, the sporadic microscopic ceramic fragments found out at the depth intervals of 16.6–16.63 m in P2, 14.9–14.94 m in P4 and 15–13 m in P5 sequences, as well as the upward increase of Cerealia-type pollen near the upper part of Unit A might signal the presence and the activities of the first farming communities in the area.

5.1.2. The Shallow Marine Bay and the Island of Piraeus (Unit B)

Unit B (ca. 6800–5400 cal BP; see [28]) extends between 10 and 7 m, 9.5 and 6 m, and 11.34 and 5.8 m in the P2, P4 and P5 sequences, respectively. The benthic foraminiferal assemblages within this deposit present increased FDs and a high A index (Figure 3a, Figure 4a and Figure 5a), implying the establishment of full marine conditions [9]. This is further enhanced by the composition of the foraminiferal assemblages, comprising primarily marine species characterized by the significant presence of A. beccarii as well as the lower frequencies of other typical marine infralittoral species such as E. complanatum, R. bradyi, A. mamilla, B. frigida and Bolivina spp. In addition extra evidence for the establishment of marine conditions during the deposition of Unit B is provided by the remarkable presence of miliolid species, such as Q. seminula, Q. berthelotiana, Q. bicarinata and Q. parvula (being typical representatives of the infralittoral and upper circalittoral zones [99], which is comparable to the significant presence of planktonic foraminifera within the deposit (P ratio obtaining its maximum value). The minor presence, however, of the euryhaline A. tepida and H. germanica also indicates a constant, albeit low, freshwater input. In parallel, the ostracod C. torosa remains highly dominant in the relevant assemblages of Unit B, forming in most cases monospecific concentrations, thus, indicating the persistence of shallow brackish waters especially for the lower part of the deposit. However, the radical change in the abundances of the ostracod populations, i.e., C. torosa, L. elliptica, X. communis and C. fischeri in the lower portion of Unit B versus Xestoleberis spp., S. incogruens, L. affinis, L. rubritincta and A. woodwardii (marine association; e.g., [60,95]) in the upper portion of the deposit (Figure 3b, Figure 4b and Figure 5b), suggest the gradual establishment of a shallow coastal marine environment. The previous interpretation is in agreement with the prevailing lithology (sand alternations with clay and silt) in the upper part of Unit B in the P4 sequence, which marks the termination of the lagoonal environment. Further, the mollusc fauna of the deposit, mainly represented by Cerithiidae, Rissoa spp. and Tricolia spp., strengthens the interpretation for the occurrence of shallow-water marine conditions, even though the common presence of C. glaucum implies lagoonal conditions [78,79,80,100]. Eventually, the occurrence of B. reticulatum shells and some calcareous worm tubes of Spirorbis spp. within Unit B infers potential sea grass vegetation [101], while the deposit’s moderate to low magnetic susceptibility (see Figure 2) signals marine sedimentation.

Therefore, based on the evidence presented above, the development of a shallow marine bay in the Piraeus coastal plain during 6800–5400 cal BP is a reasonable paleoenvironmental interpretation, with Piraeus forming a tied island in the center of the bay (after [28]), probably connected to the mainland via a sandy isthmus (Figure 9a,b).

5.1.3. The Second Stage of the Piraeus Lagoon (Unit C)

The microfaunal content of Unit C within the 6–4 m and 5.8–4.5 m depth intervals of the P4 and P5 sequences, respectively, (ca. 4800–3500 cal BP; see [28]) reveals the relative increase of A. tepida and H. germanica, along with Q. seminula, while the molluscan fauna consists of C. glaucum, Abra spp. and numerous Hydrobiidae, suggesting the reestablishment of the closed lagoon conditions (Figure 7). The previous interpretation is further supported by the relatively increased magnetic susceptibility within Unit C (Figure 2). Even though, a few charophyte gyrogonites indicate freshwater inputs, the sporadic presence of gypsum crystals in Unit C can be associated with short-term desiccation episodes in the lagoon. Comparably, the ostracod abundance is highly enhanced in the deposit, with high numbers of C. torosa and oligohaline to freshwater species, reflecting the transition to a restricted oligohaline lagoonal environment. This oligohaline closed lagoon, featured in the pollen diagrams by the expansion of the Amaranthaceae halophytes (Figure 4d), was used for grazing, as indicated by the occurrence of the coprophilous fungal remains. Furthermore, the human presence is clearly detected during the deposition of Unit C by the increase of cultivars like Cerealia-type and Olea featuring Early Bronze Age cultivation practices.

According to Goiran et al. [28], the closed lagoon was established during 4800–3500 cal BP and was separated from the sea by beach barriers. The paleoenvironmental reconstruction suggested by the present study delimits this restricted depositional system in the western part of the Piraeus plain, while at the same time a coastal environment associated with the fluvio-deltaic system of Kifissos and Korydallos Rivers, was being formed to the east (Figure 9c,d). The transition to closed lagoon conditions, already documented in several coastal plains of the Aegean region [9,12,13,15,18], is most probably related not only to the fluvial inputs and the subsequent deltaic progradation, but also to an initial deceleration of the sea level rise since 6000 cal BP, followed by a slowdown recorded at ~4000 cal BP [13,102]. Obviously, the warm and humid climate conditions associated with the ongoing Mid-Holocene African monsoon forcing during 5400–4300 yr BP [103,104], affected the relatively higher sea level rise recorded before ca. 4000 cal BP.

5.1.4. The Piraeus Marsh (Unit D)

Throughout the upper 4 m of the P4 sequence and between 4.5 and 2.5 m in the P5 sequence (Unit D, being younger than 2800 cal BP; see [28]), the molluscs are practically absent, while the foraminiferal and ostracod microfauna indicates an oligohaline to brackish marshy environment (Figure 7). Oligohaline to freshwater ostracod species gradually predominate in the deposit, mainly with Ilyocypris spp., N. neglecta (a freshwater species occurring in temporary and permanent freshwater bodies; see [61]) and C. salinus. Therefore, the formation of an oligohaline to freshwater marsh can be suggested. This is also supported by the maximum values of the erosion indicating palynomorphs, revealing enhanced fluvial inputs (Figure 7). Furthermore, the occurrence of a plethora of numerous Zygnemataceae (Figure 4d), is indicative of the freshwater depositional environment of Unit D. Finally, the magnetic susceptibility of the clay layer occurring in the P4 sequence demonstrates its higher values throughout the sequence, reaching a maximum of 100 × 10⁻⁸ m³/kg (Figure 2), thus, strengthening the interpretation for the high impact of freshwater inflows [33].

The Piraeus marsh, formed after 2800 cal BP, was restricted in the western part of the Piraeus plain (Figure 9e,f). A coastal paleoenvironment was being developed to the east (i.e., the Piraeus Coast, Unit E), mostly due to the Kifissos and Korydallos water and sediment discharges, and to the gradual expansion of the deltaic floodplain, in an area of prominent tectonic stability.

During this interval, corresponding to Archaic/Classical Periods, the maxima of both Olea and Cerealia-type, as well as of coprophilous fungi attest the intensification of human activities in the area during the Antiquity. However, the Piraeus olive/cereal cultivation pattern does not present the shift from cereal cultivation during the Bronze Age, to olive cultivation in the Antiquity already recorded in the eastern Attica (e.g., Vravron site), [27]. Olea and Cerealia-type are emblematic cultivars in the Mediterranean culture (e.g., [105,106], therefore common in most pollen records from south Greece (e.g., [27,107,108,109]); while preference on either crop plants is associated with local conditions [110], as well as with societal choices and trade practices [111].

5.1.5. The Piraeus Coast (Unit E; ca. 3400–2500 cal BP)

In the P2 sequence, between 7 and 3 m, the developed medium to coarse sand layers display an enhanced BR ratio, thus, indicating an impact from hydrodynamic processes [57]. At the same time foraminiferal assemblages dramatically decrease and are principally composed of the marine species A. beccarii and Q. seminula. In addition, ostracods are practically absent, since only one sediment sample at 4.7 m contained a poor assemblage with brackish water species, while the rest of samples mostly included scattered, reworked marine mollusc fragments. Hence, all evidence points to a high-energy depositional environment, i.e., a supralittoral zone (Figure 7), formed at the eastern edge of the Piraeus lagoon/marsh after 3400 cal BP in the vicinity of the Kifissos River floodplain (Figure 9c–f).

5.2. Absolute Sea Level Rise Scenario for the Inner Saronikos Gulf during Mid-Holocene

A high-resolution marine record derived from a sediment core (named S2P) recovered from the deepest part of Elefsis Bay (see the coring location in Figure 1) provides good evidence for the paleoenvironmental conditions identified in the area. In particular, the establishment of full marine conditions in the Elefsis Bay, with a sufficient water column depth, i.e., at least 20 m, to sustain the pelagic coccolithophore community [112] documented by the dominance of Emiliania huxleyi in the S2P record, can be dated at 7500 cal BP according to the age model presented in the study of Kouli et al. [83] (Figure 7). Elefsis Bay is connected to the inner Saronikos Gulf via two narrow straits, a western and an eastern one, both having sills occurring at water depths of ~8 and ~12–13 m, respectively. The establishment of marine conditions in the bay at 7500 cal BP suggests that during this time, when the sea level position was at −13 m [84,113], seawater masses started spilling over the eastern sill, and subsequently, reached the western part of the Piraeus coastal plain. Hence, the dating at the base of Unit B of the Piraeus plain sequences should be reevaluated at ~7500 cal BP (Figure 7), instead of 6800 cal BP, thus, improving the age constraint (i.e., 6800–5400 cal BP) suggested by [28] for the shallow marine depositional environment of Unit B.

As shown from the paleoenvironmental evolution of the Piraeus coastal plain, an 8-m-deep lagoon had been developed in the eastern part of the plain at least till 7500 cal BP (see Figure 8). Based on the previously provided interpretations (see above), during 7500–5400 cal BP seawater inundated the Piraeus plain through its western part since the eastern part was an uninterrupted coastal area in the vicinity of the Korydallos and Kifissos River deltaic plain (see Figure 9c–f). Considering the paleo-relief of the western part of the Piraeus plain (see Figure 8c and [28]), the sea level rose 10 m during a period of 2100 years, thus, demonstrating an absolute sea level rise rate of ~5 mm/yr for the Mid-Holocene. The estimated rate appears excessive compared to the average current sea level rise rate of ~1 mm/yr in the Aegean Sea as derived from Late Holocene data (see [114]). However, as the Mid-Holocene climate was warmer than today by 2–3 °C in average [103,104,115] offers a reasonable explanation for this deviation. Indeed, the proposed by the IPCC RCP 4.5 scenario, for an average increase of the global temperature by 2–3 °C, predicts a sea level rise rate of 5.2 mm/yr [116].

6. Conclusions

The multi-proxy investigation of the Piraeus coastal plain paleoenvironmental evolution during the last ca. 9000 cal years BP showed a very good correlation of the data derived from the benthic foraminiferal analysis with the data obtained from ostracod, molluscan and palynological analyses as well as magnetic susceptibility measurements. Hence, the use of such a multi-proxy approach for the paleogeographic reconstructions of the dynamic coastal environments shows great potential to be established as a very promising paleoenvironmental tool. Regarding the study area, which is part of the metropolitan city of Athens, the following evolutionary stages can be defined:

• During 8700–7500 cal BP the first stage of the Piraeus Lagoon was established, reflected in the Unit A deposit of the local sedimentary archive. However, based on the data of this study, two sub-stages in the evolution of the lagoonal environment are identified:

  • (i). a closed lagoon system, resulting in the formation of the Unit Aa Paleoenvironment (the older part of Unit A) and restricted to the eastern part of the coastal plain; and

  • (ii). a mesohaline open lagoon system, resulting in the deposition of the Unit Ab Paleoenvironment and still remaining confined to the eastern part of the plain, due to a topographic high formed by the Pliocene substrate. Since the sea level position at ca. 8000 yr BP was at −17 m, the marine influence should have been established either via channels incised in the fluvio-deltaic area of Kifissos and Korydallos rivers or due to the salinization of the unconfined aquifer. The increased occurrence of Cerealia-type pollen in Unit Ab may be indicative of the first activities of the farming Neolithic communities in the area.

• During 7500–5400 cal BP a shallow marine bay was developed (reflected in the Unit B deposit of the local sedimentary archive), being the result of the Holocene sea level rise. At that time, Piraeus was a tied island connected to the mainland through a sandy isthmus.

• During 4800–3500 cal BP the second stage in the Piraeus Lagoon evolution occurred (reflected in the Unit C deposit of the local sedimentary archive), featured by a restricted environment with oligohaline conditions. The lagoon was delimited in the western part of the Piraeus plain, partly affected by the warm and humid climate conditions associated with the ongoing Mid-Holocene African monsoon forcing during 5400–4300 yr BP. To the east, the fluvio-deltaic deposits of Kifissos and Korydallos rivers were configuring the Piraeus Coast (reflected in the Unit E deposit of the local sedimentary archive; ca. 3400–2500 cal BP). During this stage, the borderlands of Piraeus Lagoon were used for cultivation and grazing activities, as evidenced by the palynological data.

• After 2800 cal BP the Piraeus Marsh was formed (reflected in the Unit D deposit of the local sedimentary archive), being confined to the western part of the coastal plain, while the Piraeus Coast continued to expand to the east. The increased occurrence of cultivars, such as Cerealia-type and Olea, in Unit D provides good evidence for human activities during the Archaic and Classical Periods.

Finally, the establishment of a shallow marine paleoenvironment in the Piraeus coastal plain at 7500 cal BP is in agreement with the pure marine conditions documented in a high-resolution sedimentary record obtained from the nearby Elefsis Bay. Based on the previous finding, an absolute sea level rise rate of ~5 mm/yr can be estimated for the broader study area during the 7500–5400 cal BP interval of the warm and humid Mid-Holocene. Provided that the current climate trend will approach- in the future the climate conditions prevailing in the Mid-Holocene, an analogous rate for the sea level rise in the inner Saronikos Gulf (being adjacent to the metropolitan city of Athens) should be expected.

Author Contributions

Conceptualization, M.V.T., T.T., I.P.P., J.-P.G., E.F. and K.P.P.; methodology, M.V.T., T.T., K.K., O.K., M.D.D., E.V.A., G.S., V.A.; investigation, M.V.T., T.T., K.K., O.K., J.-P.G., E.F. and K.P.P.; data curation, O.K., M.D.D., D.V., A.P., S.C., E.S.; writing—original draft preparation, M.V.T., T.T., K.K., O.K. and M.D.D.; writing—review and editing, M.V.T., T.T., K.K., O.K., M.D.D., E.V.A., G.S., I.P.P., K.P.P.; visualization, T.T., O.K., M.D.D., E.S.; project administration, M.V.T., K.P.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Centre National de la Recherche Scientifique (CNRS), CONCLIMA: Constraining large-scale climate feedbacks in the Earth system using paleo data, HORIZON2020, grant number 56.90.13459 and the Greek National Project CLIMPACT: Flagship Initiative for Climate Change and its Impact by the Hellenic Network of Agencies for Climate Impact Mitigation and Adaptation.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Generate image description

Figure 1. Map of the study area: (a) map showing the Saronikos Gulf; (b) Map of Greece where Saronikos Gulf is indicated; (c) Piraeus coastal plain showing the location of boreholes P2, P4, P5; (d) Elefsis bay showing the S2P coring location (source Hellenic Navy Hydrographic Service).

Generate image description

Figure 2. Lithostratigraphic columns of the studied boreholes and the results of the performed sedimentological and magnetic susceptibility analyses. Lithostratigraphic units (A–E, separated by dotted lines) and the chronostratigraphic framework are based on [28].

Generate image description

Figure 3. (a) Distribution of foraminiferal and ostracod indices along borehole P2 in respect to lithostratigraphic units; (b) Relative abundances of the major foraminiferal and ostracod species in borehole P2; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P2. Lithostratigraphic units (A,B,E) are separated by dotted lines.

Generate image description

Figure 3. (a) Distribution of foraminiferal and ostracod indices along borehole P2 in respect to lithostratigraphic units; (b) Relative abundances of the major foraminiferal and ostracod species in borehole P2; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P2. Lithostratigraphic units (A,B,E) are separated by dotted lines.

Generate image description

Figure 4. (a) Distribution of foraminiferal and ostracod indices along borehole P4 in respect to lithostratigraphic units; (b) Relative abundances of the major benthic foraminiferal and ostracod species in borehole P4; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P4; (d) Relative abundances (exaggeration ×10) of the most important pollen and other palynomorph assemblages in borehole P4. Lithostratigraphic units (A–D) are separated by dotted lines.

Generate image description

Figure 4. (a) Distribution of foraminiferal and ostracod indices along borehole P4 in respect to lithostratigraphic units; (b) Relative abundances of the major benthic foraminiferal and ostracod species in borehole P4; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P4; (d) Relative abundances (exaggeration ×10) of the most important pollen and other palynomorph assemblages in borehole P4. Lithostratigraphic units (A–D) are separated by dotted lines.

Generate image description

Figure 4. (a) Distribution of foraminiferal and ostracod indices along borehole P4 in respect to lithostratigraphic units; (b) Relative abundances of the major benthic foraminiferal and ostracod species in borehole P4; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P4; (d) Relative abundances (exaggeration ×10) of the most important pollen and other palynomorph assemblages in borehole P4. Lithostratigraphic units (A–D) are separated by dotted lines.

Generate image description

Figure 5. (a) Distribution of foraminiferal and ostracod indices along borehole P4 in respect to lithostratigraphic units; (b) Relative abundances of the major benthic foraminiferal and ostracod species in borehole P5; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P5. Lithostratigraphic units (A–D) are separated by dotted lines.

Generate image description

Figure 5. (a) Distribution of foraminiferal and ostracod indices along borehole P4 in respect to lithostratigraphic units; (b) Relative abundances of the major benthic foraminiferal and ostracod species in borehole P5; (c) Distribution and semiquantitative analysis of the main mollusc taxa and other microscopic findings in borehole P5. Lithostratigraphic units (A–D) are separated by dotted lines.

Generate image description

Figure 6. (a) The present landscape of Piraeus coastal plain and the location of the drilled boreholes; (b) A map of the late 19th century [69], projected on the present landscape, showing the city of Piraeus and the Korydallos River flowing within the Piraeus coastal plain; Korydallos was an active element of the drainage system until the urbanization of the area.

Generate image description

Figure 7. Distribution of foraminiferal and ostracod assemblages along the studied P2, P4, P5 boreholes and intercomparison of the recorded paleoenvironmental changes. Correlation with the paleoenvironmental interpretation of the marine core S2P from Elefsis Bay, is indicated [83]. Paleoenvironmental units Aa–E are separated by dotted lines.

Generate image description

Figure 7. Distribution of foraminiferal and ostracod assemblages along the studied P2, P4, P5 boreholes and intercomparison of the recorded paleoenvironmental changes. Correlation with the paleoenvironmental interpretation of the marine core S2P from Elefsis Bay, is indicated [83]. Paleoenvironmental units Aa–E are separated by dotted lines.

Generate image description

Figure 8. Paleogeographic reconstruction of the First Stage of the Piraeus Lagoon in respect to the present coastline: (a) The closed lagoonal environment of Unit Aa; (b) the mesohaline open lagoonal environment of Unit Ab; (c) Simplified sketch-transect showing the development of the First Stage of Piraeus Lagoon (8700–7500 cal BP), restricted at the eastern part of the plain (depths according Goiran et al. [28]).

Generate image description

Figure 9. Paleogeographic reconstruction in respect to the present coastline and simplified sketch-transects (depths according Goiran et al. [28]), showing: (a,b) The Shallow Marine Bay of Piraeus (Unit B; 7500–5400 cal BP); (c,d) The Second Stage of the Piraeus Lagoon (Unit C; ca. 4800–3500 cal BP); (e,f) The Piraeus Marsh (Unit D; younger than 2800 cal BP).