Integrating Geomorphological Data, Geochronology and Archaeological Evidence for Coastal Landscape Reconstruction, the Case of Ammoudara Beach, Crete

Abstract

Understanding the processes that govern the transformation of the landscape through time is essential for exploring the evolution of a coastal area. Coastal landscapes are dynamic sites, with their evolution strongly linked with waves and sea level variations. Geomorphological features in the coastal area, such as beachrock formations and dune fields, can function as indicators of the coastal landscape evolution through time. However, our knowledge of the chronological framework of coastal deposits in the Aegean coasts is limited. Optically Stimulated Luminescence dating techniques are deemed to be very promising in direct dating of the coastal sediments, especially when they are linked with archaeological evidence. The dating of the sediments from different sediment core depths, determined by the method of luminosity, allowed us to calculate the rate of sediment deposition over time. More recent coastal evolution and stability were examined from 1945 to 2020 with the use of aerial photographs and satellite images. This paper presents the 6000 ka evolution of a coastal landscape based on geomorphological, archaeological, and radio-chronological data. Based on the results, the early stages of the Ammoudara beach dune field appears to have been formed ~9.0–9.6 ka BP, while the OSL ages from 6 m depth represented the timing of its stabilization (OSL ages ~5–6 ka). This indicates that the dune field appears to already have been formed long before the Bronze Age (5–10 ka BP). It became stabilized with only localized episodes of dune reactivation occurring. In contrast, while high coastal erosion rates were calculated for modern times.

1. Introduction

Understanding and modelling geomorphological and depositional processes, resulting in the transformation of the landscape from a diachronic perspective, are essential for exploring the cultural biography of a coastal area. Coastal landscapes are dynamic landscapes, with their evolution strongly linked with waves and sea level variations. There are several challenges in the study of the cultural landscape of coastal areas. In terms of environmental change, these challenges are principally related to the dynamic character of coastal areas. Wave and wind action result in processes such as coastal erosion and the formation of coastal dunes. Coastal erosion often leads to an irreparable loss of information used in the identification of maritime structures and activities. The Mediterranean, as a major example of a strictly coastal related inhabited area since the ancient times, is a hotspot of both opportunity and concern. Consequently, these changes in the coastline can generate new social relations as newly formed areas can be suitable for agriculture or already inhabited ones can be abandoned due to flood by seawater [1].

1.1. Tectonics of Crete

The island of Crete is located along the transition zone between the African and Eurasian plates in the fore arc of the Hellenic Subduction Zone. Differential tectonic behavior is observed in the continental area of Crete. Compression prevails in western Crete, while tectonic extension dominates in eastern Crete [2,3]. By the late Upper Miocene to early Pliocene, the conditions of deposition changed from open sea to coastal or subaerial environment [2,3]. The island uplifted gradually during the Quaternary, as evidenced by marine Pleistocene terraces and coastal deposits formed during the corresponding interglacial stages [4,5,6,7,8,9,10]. The fragmentation and uplift of the island continued during the Pleistocene [6]. During the Upper Holocene, subsidence affected the entire island of Crete and tectonic activity produced intense faulting [11,12,13,14]. Three subsequent phases of tectonic activity formed initially E–W-trending faults that affected Miocene sediments. Next, N–S-striking faults cut the E–W faults which had previously formed, while the younger phase included NE–SW- and NW–SE-dipping faults [15,16]. Arc-parallel (E–W) and arc-normal (N–S) extension faults are related to the Holocene tectonic activity and are possibly connected with large historical earthquakes and recent moderate seismicity (M = 6.0–6.5) [17,18,19,20]. An evolution scenario for the coast of central and eastern Crete [21] indicated that the eastern part of Crete has subsided gradually toward NE to a depth greater than 4.0 m in at least three subsiding phases over the last 4000 years [14].

Moreover, the eruption of Santorini in the late Bronze Age resulted in tsunami waves, based on the evidence of tsunami deposits along various coastal areas of Greece [22,23]. Numerical models have shown that tsunami wave amplitudes were significant, with inland inundation of up to 200 m along the coast of Crete [24]. More recently, evidence from the AD 365 earthquake, the largest known seismic event in the Eastern Mediterranean, has supported the idea that the respective tsunami not only reached the NW part of Crete but also severely affected the inhabited coastal areas in Heraklion [25], where vertical movements in the area reached −1 m [26]. The north coast of the island appears to be rather stable after the AD 365 earthquake event [27] (Figure 1).

1.2. Indicators of Sea Level Changes on Crete

Geomorphological features such as tidal notches, beachrock formations, and archeological evidence along the coastline of Crete are evidence of relative sea level changes over time (during the last 5000 years). In central and eastern Crete, the relative sea level change evidence from tidal notches and beachrock formations revealed five distinct sea level stands at −6.55 ± 0.55 m between 4200 ± 90 BP and 3930 ± 90 BP, −3.95 ± 0.35 m between 3900–3700 BP, −2.70 ± 0.15 m between 3600–3450 BP, −1.25 ± 0.05 m, and −0.55 ± 0.05 m. The lowest sea level stand can be identified by the oldest dated tidal notch of western Crete. Historical sources report a relative sea level rise of 0.70 m during the AD 1604 paroxysmal event. Over the last 400 years, the relative sea level rose by 0.55 m [14].

1.3. Beachrocks and Aeolian Deposits

Beachrock formations originate from the rapid cohesiveness of coastal loose sediments (sands, gravels) in the inland beach area. Their formation takes place either superficially and/or under a thin layer of sediments in the subtidal zone and includes sediments of various diameters and of classical and/or biogenic origin [28,29,30]. The dimensions of the beachrocks range from a few meters to a few kilometers in length, with a width from less than 1 meter to 300–400 m and thickness from 0.3 to 3 m. These formations present flat layers that slope toward the sea, either with the general slope of the beach [31] or with a greater slope than that of the rest of the beach [32,33,34]. Appearances have been reported on many coasts of the world [35,36,37] and in Greece [38,39]. Because they form in the tidal zone, beachrock formations are considered indicators of old coastlines and are often used as a fully accepted and valid method by the scientific community in research to detect sea level changes in earlier periods [36,40,41]. Beachrock formations in Crete are composed of sand and gravels (quartz, limestone lithoclasts, metamorphic, and volcanic) with generally less than 10% of bioclasts. The intertidal cement is cryptocrystalline to microcrystalline, mainly Mg-calcite and in some cases pure calcite [14,42,43]. The determination of the age and position of the respective past coastline and the corresponding sea level by the depth of occurrence of a beachrock formation and its age, using stratigraphic or radio-chronological methods, is now widely used [44,45,46,47,48]. On the other hand, aeolianite is any rock formed by the lithification of sediment deposited by aeolian processes. Coastal aeolianites are usually coastal formations consisting of carbonate sediment of shallow marine biogenic origin, formed into coastal dunes by the wind and subsequently lithified [49]. Quaternary aeolianites are typically coastal deposits composed of fine-to-medium-grained, well-sorted sand. The character of the sand grains significantly depends on the local environmental setting, although the dominant constituents of the aeolianites are quartz and feldspar grains and marine carbonate particles. Both the detrital mineral grains and the carbonate particles are often derived from local sources. Several studies have used dating methods to establish the time framework of aeolianites and beachrock formation and to assess their connection with sea level changes [50,51,52,53]. Aeolianites are mainly considered to be deposited around the world at sea level high stands [54,55], but several recent studies have indicated the cold and low sea level scenario as a genetic mechanism for the formation and accumulation aeolian deposits [56,57,58,59,60,61,62,63]. Differences in the timing of deposition and morphology of aeolian dunes during major sea level changes (high stands, low stands) among regions appear to be dependent on the interactions caused by minor-secondary sea level changes, local climatic conditions, sediment availability, and differences in the inland-near-shore topography [64]. However, our knowledge on the chronological framework of coastal deposits on the coasts of the Aegean is limited. Early studies have found sedimentary deposits on raised marine terraces along the coasts of Crete, which reveals the existence of characteristic Pleistocene marine fauna such as Cantharus viverratus, (Kien), Conus testudinarius (Mart), Polinices lacteus (Guild), Patella safiana, (Lmk), and Strobus bubonius (Lmk) [65,66,67,68,69].

1.4. Coastal Erosion

Coastlines in Greece are generally in a state of constant erosion due to anthropogenic interventions and climate change. Studies on the area [70,71] have estimated that about 25% of the Aegean coast consists mainly of low-lying coasts and deltas in a state of erosion. The northern coastline of Crete experiences extensive erosion, the highest severity in Greece, with erosion affecting 65.8% of its coastal zones. This large percentage of erosion is mainly due to the presence of extensive coastal zones with generally low beach slope, which are thus more vulnerable to erosion and the high frequency of large waves created by north and northeast winds in the Aegean [72].

1.5. Archaeological Setting

The coastal area of the Gazi-Ammoudara region is 8 km northwest of the palatial center of Knossos and not far from the Psiloritis mountain region. The region belonged to the administrative Minoan center of Tylissos, which is among the most poorly archaeologically investigated areas considering its proximity to some of the most important Minoan sites of Crete. This gap in archaeological research provides bias in regional and inter-regional archaeological interpretations and affects our understanding of the political, social, and economic networks in north-central Crete. Consequently, for heritage management, this gap is even more pronounced, since the lack of data and limited knowledge about existing archaeological sites and features implies difficulties in implementing effective strategies for their preservation. The archaeological importance of the region has been documented by early scholars [73,74]. Several archaeological remains in the area were destroyed by natural processes and modern human activity (agriculture, grazing, and construction). Evidence from rescue excavations, surface remains, and chance finds dating from the Neolithic to the Ottoman period have been found in many parts of the area, suggesting intensive patterns of habitation and other human activity. However, the only excavation in the coastal zone was carried out in the late 1970s. The excavation revealed a building, which was only partially excavated (Figure 2). It seems that it was a domestic unit, which most likely belonged to an extensive Neopalatial (MM III-LM IA) settlement that was part of the wider social contact network in the northern coastal zone of Minoan Crete. Pottery from the excavation is dated predominately to MM III-LM IA, though it includes Pre- and Postpalatial sherds [75,76].

1.6. Optically Stimulated Luminescence (OSL) Dating

Luminescence measurements are most frequently carried out on the quartz–rich fraction of aeolian dune deposits [77,78]. Environments rich in coastal and wind sediments [77] such as loess formations and dunes [79] have been shown to be suitable materials for Optically Stimulated Luminescence dating as they have been adequately exposed to sunlight before their deposition. However, in other types of sedimentary environments where there is doubt about their adequate exposure to sunlight, an investigation of the “zeroing” of the original luminance signal is required. Relatively recent applications of luminescence techniques in a variety of sedimentary environments [80], such as alluvial and colluvial sediments, glacial and peripheral deposits [81], riverbed sediments [82], sediments associated with cracks and seismicity [83,84], and terrace materials [85], have proven to be quite reliable, as the obtained ages were consistent with the chronological data found in existing reports. Optically Stimulated Luminescence (OSL) dating [86] has become a promising technique for dating fluvial, deltaic, and coastal sediments [27,87,88,89,90]. Thus, the Luminosity method is preferred as the most appropriate for investigating the location of the coastline [89,91,92,93,94,95]. The ages of the sediments from different depths, as they are determined by the method of luminosity, allow us to calculate the rate of sediment deposition over time and, therefore, the advance or erosion of an area such as the area of Ammoudara.

This paper presents the evolution of a coastal environment based on geomorphological, archaeological, and radio-chronological data that enabled the reconstruction of the coastal area of Ammoudara beach in northern Crete during the Bronze Age.

2. Study Area

The coastal area of Ammoudara (Figure 2) is located on the northern coast of Crete Island, approximately 2 km to the west of the city of Heraklion, Crete, Greece. The area is an alluvial plain and is characterized by low-lying regions, some parts of which have elevations either close to or even below the mean sea level. The entire coastal area has been shaped by significant geomorphological and tectonic processes as well as human activity [96]. The total coastline length is approximately 8 km, with a W to E trend, being slightly curved and having a subaerial width of up to 60 m. The beach is backed by a narrow dune field. where the dunes do not exceed 3 m in height. Most of the dune fields have been destroyed by anthropogenic activities. Three small rivers, the Gazanos (catchment area 192 km²), Xiropotamos (c.a. 35 km²), and Giofiros (c.a. 279 km²), discharge along the shoreline. At the west end of the beach, the Almiros stream drains the homonymous karstic spring (Figure 2). The Ammoudara coastal zone is microtidal with an astronomical tidal range of less than 10 cm [97], while the tide gauge installed in the harbor of Heraklion has recorded sea level fluctuations of more than 1 m over the period 1984–2020. These fluctuations incorporate the effects of astronomical tides and meteorological forcing. The Ammoudara beach zone is exposed primarily to NW, N, and NE wind-induced waves, with the NW direction being most frequent on an annual basis (23.6%), followed by the NW (12.4%) and the NE (3.7%). Wave heights are usually smaller than 2 m (88%), with 36% of them being less than 0.5 m and only 2% being over 4 m on an annual basis [98]. Likewise, 77.3% of the offshore wave periods are less than 5 s, while periods larger than 11 s have an annual frequency of only 0.74%. The maximum wave conditions are induced by NW and N winds during January and March, with a mean annual occurrence frequency of 0.01% and 0.03%, respectively. The maximum wave characteristics of the NW wind-induced waves are 2.3 m/4.6 s, while those of the prevailing N waves are 6 m/11 s [99].

A shallow submerged reef characterizes the shallow marine Ammoudara beach. The reef is ~4 km-long and has a mean width of around 35 m, while its rather planar crest lies, on average, at 0.6 m below the sea surface. Water depths exceed 3 m at its seaward side, while they are less than 2 m at its landward side. Previous studies in the area [99] have concluded that the reef represents a submerged beach rock formation, initially formed concurrently with the beach rocks along the western part of the present beach, based on the morphological and sedimentological characteristics of the reef. The reef consists of two distinct sedimentary layers: The upper layer, with a thickness of ~0.3, which consists of coarse-grained beach material having beackrock morphological characteristics (e.g., runnels, potholes, and grooves), and the lower layer, with a thickness greater than 0.5 m, which consists of fine-grained sandy material with cross bedding. The coastal zone of Ammoudara Beach can be subdivided into two sectors based on the presence or absence of the reef: The western sector, where the reef is absent, and the eastern sector, which is characterized by the continuous presence of the reef at less than 40 m from the shoreline.

3. Materials and Methods

Geomorphological and sedimentological mapping of the area was conducted with the use of satellite images, drone imagery, and field surveys. For the geomorphological mapping, topographic diagrams from the Hellenic Army Geographical Service and Sentinel satellite images were used. The geomorphological surveys included elevation and slope measurements taken with a RTK GPS. Additionally, depth soundings were acquired using a single beam echo-sounder reaching 20 m depth within the water. Such data were synthesized to create current bathymetry and beach elevations. Furthermore, surficial sediment samples were collected along the beach and analyzed according to Folk’s procedure [100]. To determine the depositional environment of the coastal sediments, samples from four 10 m-deep research cores were utilized (Figure 2) at approximately 60 m from the present-day shoreline. From the resulting stratigraphic sections, the different sedimentary horizons were studied macroscopically, the cores were subsampled, and the organisms found within them were identified to determine the deposition environment. For OLS dating, lightproof plastic tubes were used to avoid light contamination.

3.1. Optically Stimulated Luminescence (OSL)

Dating and radioactivity measurements were implemented in 12 core samples using OSL and Equivalent Doses calculations (for locations and depth of the samples, see Figure 2; for stratigraphic characteristics, section). For the radioactivity estimations, the top layer (2–4 cm) of the core samples was used, and 4–8 cm of the inner part of the core was collected, dried at 50 °C, and passed through sieves sequentially decreasing mesh from 150–125 μm to 100–63 μm. Grain sizes from 100–150 μm were chemically treated for quartz purification. Parts of the samples used for radioactivity assessment were dried at 105 °C for 24 h and passed through a 500 μm mesh sieve. About 60–70 g of this sieved material was sealed in plastic containers and kept for radon equilibrium for 1 month. Afterward, the containers were measured for natural radioactivity. For chemical pretreatment, small portions of 100–150 μm grain size of the samples were placed in separate plastic tubes. All samples were treated successively with HCl (8% w/w) to remove carbonates, H₂O₂ (30% w/w organic residuals and HF (40% w/w), and feldspars. The outer parts of the quartz grains were etched off the samples. Several water washes were performed between chemical treatments.

3.2. OSL Measurements and Equivalent Dose (EDs) Calculations

Purified quartz grains (aliquots) were placed on small stainless disks for luminescence measurements. Typically, 18 disks in total were used for each sample. The SAR protocol was followed [78,101], and Paleodoses were measured after IR and blue-OSL signal recording. The IR protocol was used to check for feldspar impurities. Successively increasing doses were induced on disks to obtain “dose response curves,” and the same protocol was used for these doses as described for natural ones. To eliminate sensitivity differences of the aliquots, the same dose was induced, referred as the test dose, and the ratios of dose to test dose signals were used to obtain the “dose-response curves.” Recycling and recuperation tests, as well as dose recovery tests, were performed for quality assurance. All measurements were performed in the Archaeometry Center of the University of Ioannina. Approximately 100 g of dried material was ground to powder and sieved through a 500 μm mesh sieve at 105 °C until reaching constant weight. About 60–70 g of each sample was placed in a plastic container of standard geometry and hermetically closed with a silicone sealant. The samples could reach the equilibrium of radon daughters (for about 1 month) and were measured for natural radioactivity. The age of a sample is calculated in terms of the equivalent dose ratio or Dose Rate (Gy) [102]. The equivalent dose is the total accumulated radioactivity to which the precipitate is exposed during the time it is deposited and, in the laboratory, and is called the Paleodose (De) [102]. The unit of measurement of absorbed radioactivity in the International System of Units (S.I) is Gray (Gy) ≥ 1 Gy = 1 Joule/kg. The dose rate is the rate at which natural radioactivity is absorbed by the crystalline granules per year. The age was calculated by the ratio equivalent dose (De) per Dose Rate [102].

3.3. Coastline Evolution

Aerial photographs and satellite images were used to study the stability of the coastline and to measure its displacement over time. Based on these data, the coastlines were digitized with the utmost accuracy through ArcGIS 10.1 software, through which conclusions were drawn for the shoreline displacement over time. This was followed by processing with Digital Shoreline Analysis System DSAS v.5 [103] to quantify the changes. Three different types of remote sensing data—historical analogue panchromatic aerial photographs, orthophotos, and natural-color satellite images—were used for the identification of the beach areas alterations during the last 74 years. The acquisition of the aerial photographs took place during 1945, 1968, and 1998 from Hellenic Military Geographical Service (HMGS); the orthophotos were acquired during 2005 and 2010 from Ktimatologio S.A; and the satellite images during 2009, 2012, 2013, 2016, 2017, and 2019 (

Table 1. Date, type, scale, resolution, and source of the images used for coastal evolution assessment.

Date	Type	Scale	Pixel Size	Source
17/09/1945	Aerial photo	1:42,000	3 m	HMGS
1968	Aerial photo	1:15,000	1 m	HMGS
23/07/1998	Aerial photo	1:10,000	1 m	HMGS
2005	Orthophoto	1:15,000	0.6 m	Ktimatologio S.A.
14/06/2009	Satellite image	-	0.5 m	Google Earth (maxar technologies)
20/06/2010	Aerial photo	1:30,000	1.5 m	HMGS
2011	Orthophoto	1:15,000	0.6 m	Ktimatologio S.A.
06/11/2012	Satellite image	-	0.2 m	Google Earth (maxar technologies)
29/04/2013	Satellite image	-	0.5 m	Google Earth
27/06/2016	Satellite image	-	0.2 m	Google Earth (CNES/Airbus)
28/06/2017	Satellite image	-	0.3 m	Google Earth
24/03/2019	Satellite image	-	0.2 m	Google Earth (CNES/Airbus)

). The aerial photos were scanned with a high-resolution scanner (1200 dpi), and all image data wee georeferenced in ArcGIS 10.1, with the 2010 Ktimatologio S.A. aerial photo as reference dataset, using approximately 40 control points along the beach and 20 landmarks. For all spatial data, the WGS 84/UTM zone 35N Reference System was used. The beach area and shorelines of the georeferenced images were manually digitized with the maximum accuracy, determined by each image’s pixel size. The average accuracy of all images was about 0.5 m. Transects perpendicular to the shoreline were created through DSAS with 50 m intervals. The statistic parameter extracted by the DSAS and presented in this work is the Net Shoreline Movement (NSM), which is the total distance between the earliest and most recent shorelines for each transect, in meters. Table 1 shows the date, type, scale, and resolutions of each image used in the analysis.

3.4. Terrain Visualisation

The contemporary landscape was created by draping the 1968 georectified aerial photography onto the DEM using 1968 topographic diagrams from the Hellenic Military Geographical Service–HMGS with 1:5000 scale. The 1968 topographic diagrams and the 1968 aerial photograph were used, since they are the only combined data of the area before the first constructions in the coastal zone (Figure 3). After this, the correlation of sea level changes, beach erosion rates, sedimentological rates, and the geomorphological characteristics were considered for the possible reconstruction of the area, with the elevation ranging from −10 m to +10 m.

4. Results

4.1. Textural Characteristics

The subaerial part of Ammoudara Beach consists mostly of sand (S) and slightly gravely sand ((g)S). The inner nearshore zone (between the shoreline and the reef) is generally characterized as slightly gravely sand ((g)S), while seaward of the reef and in the nearshore zone of the western sector the seabed is predominantly covered by sands (S). At the vicinity of the mouth of the rivers Almiros, Gazanos, and Xiropotamos and toward the eastern end of Ammoudara beach (close to the channelized Giofiros River), the sediments become coarser, varying from gravely sand (gS) to sandy gravel (sG). Relatively coarser material (i.e., sandy gravel (sG) and gravel (G)) is also observed along the landward side of the reef. In the western sector of Ammoudara Beach, the sediments become progressively finer from W to E. The same trend exists in the central and eastern parts of the beach, with the exceptions of the small areas with coarser sediments (Figure 4)

4.2. Coastal Evolution

For the study of the stability of the coastline from 1945 until today, aerial photographs and satellite images were used. From the examination of the aerial photographs for the years 1945, 196,8 and 1998, a clear retreat of the coast can be observed due to the constant and permanent erosion induced by the waves. The retreat of the coastline was estimated to be 45 m on average. In the diagrams of all checkpoints, a displacement rate of ± 1.5 m was calculated, which was due to the errors of georeferencing and digitization of the coastline of aerial photographs and satellite images. Small shoreline displacements (< ± 5 m) were considered within the limits of the expected transient changes in the terrain of the coastal zone due to ephemerous changes in hydrodynamic conditions and displacements due to transitions of the beach from cross section to storm cross section and, therefore, were not rated as significant. The results of the comparisons are presented in Figure 5 and

Table 2. Changes in the position of the broader beach coastline in meters, on selected transects (for location see Figure 5), during the period 1945–2020.

Transects	1945–1998	1998–2002	2002–2009	2009–2015	2015–2020	1945–2020
20	0.0	−1.9	3.4	−3.3	−4	−3.4
37	−17.8	−11.0	3.6	14.7	−12.4	−23.1
50	−33.0	−15.8	15.2	−7.3	−0.6	−41.4
75	−40.9	−8.7	0.6	−0.3	−4.7	−53.9
90	−49.2	−10.2	3.3	−2.9	3.3	−55.6
110	−53.7	−7.2	2.6	−3.5	11.2	−50.5
125	−42.4	−7.4	−2.5	−3.9	9.8	−31.45

.

From the 1945 aerial photograph (Figure 5), it is seen that the dune field in 1945 was more extensive, in respect to current state, with the wider coastal area having more accommodation space and availability of sediments. This indicates that during past periods, the area also had an extensive wide terrestrial coastal zone and an abundance of fine-grained sediment. The 1945 aerial photograph shows the boundary of the stabilized field of dunes (by sparse vegetation). Furthermore, the boundary of the fully developed field of dunes with denser vegetation is clearly visible. The analysis of the 1945 aerial photograph shows that, at that time, the area of the transect 37 and Drill G3 was not affected by the waves. Instead, the affected area was located at the marsh of Almiros karstic spring delta, which is in line with the reports of early residents of the area.

4.3. Stratigraphy

To determine the stratigraphy of the area, four research cores were drilled. The geographical characteristics of these drills are presented in

Table 3. Topographical characteristics of the research drills.

No	Lat	Lon	Altitude	Distance from the Coastline	Depth
G1	35°20′15.62″ N	25°6′7.79″ E	2.45 m	67.5 m	10 m
G2	35°20′14.73″ N	25°4′53.74″ E	2.33 m	55 m	10 m
G2′	35°20′14.40″ N	25°04′53.57″ E	2.37 m	51 m	8 m
G3	35°20′22.75″ N	25°3′33.88″ E	2,62 m	76.6 m	10 m

while their locations are shown in Figure 5.

The stratigraphy of the first drill (G1) indicates that the soil layer extends from the surface down to 3 m depth. It consisted of red sand soils and modern anthropogenic materials (e.g., building material, plastics). The next layer (3–6 m) was about 3 m-thick and consisted of sand with a small presence of clay (Mz: In the next 1.5 m (6–7.5 m), there were sandy deposits (Mz: 164.84 μm) in the upper part until 7.5 m. At a depth of 7.5 m, a layer of sandy gravel (Mz: 270.20 μm) appeared which extended, at least, down to 10 m. At a depth of 3.5 m, 5.8–6 m, and 8–8.3 m, blind samples with codes G1CAD1, G1CAD2, and G1CAD3 were collected in a plastic tube for optical luminosity analysis.

In the second drill (G2), the soil horizon was found at the 0–2 m drilling depth, which consisted of red sands and anthropogenic materials. Then, at the depths of 2–4 m, the sediments were sandy (Mz: 156.62 μm). At a depth of 3.8 m and for 20 cm (3.8 –4 m), a layer of slightly coarse-grained material was found in the presence porous igneous volcanic rock that was identified as Pumice. After that, and at least until the depth of 10 m, the sediments were gravelly sandy deposits (Mz: 266.51 μm). Five blind samples were taken at depths of 4.1–4.4 m, 4.9–5.2 m, 6.2–6.4 m, 7.2–7.4 m, and 9.5–10 m, with sample codes G2CBD1, G2CBD2, G2CBD3, G2CBD4, and G2CBD5. In the same area in 2016, an initial research drill (G2′) was performed for the collection of soil samples. The drilling data are in line with the latest data. The soil horizon, which consisted of red sands and carried anthropogenic materials, was found at the drilling depth of 0–1.8 m. In the next about 20 cm (1.8–2.0 m), the material was sand with a small presence of clay (Mz: 109.73 μm). In the next 1.5 m (2.0–3.5 m), there wee sandy deposites (Mz: 132.87 μm). At 3.5 m depth and for 20 cm (3.5–3.7 m), a layer with slightly coarse material and the presence of pumice was found. From 3.7 m depth to 4.5 m, the material was gravelly sand (Mz: 261.42 μm), while at 4.5 m depth, the material was aquifer, and a small horizon with clay sand at appeared the depth of 4.5–4.7 m. From this depth until the end of the drilling (4.7–8 m), the material was slightly gravelly sand (Mz: 259.72 μm). At a depth of 6.5 m (4.2 m below sea level), a remnant of pottery was found, which, after a conversation with the archaeologist who has the license to study the area, we discovered a possibility that it is of Minoan origin, as it has similarities with those found on the adjacent plot and at the same altitude [104].

In the third drilling (G3), the stratigraphy was slightly different due to the significant presence of organic materials. The soil horizon, which consisted of grey sand soils and anthropogenic materials, was visible at a drilling depth 0–0.4 m. In the next 1.6 m (approximately 0.4–2.0 m), the material was sandy (Mz: 124.57 μm), with a small presence of roots. In the next 4.2 m (2.0–6.m), there were grey sandy (Mz: 119.36 μm) and gravelly sandy (Mz: 221.31 μm) deposits, similar to coastal land dunes, in the vicinity of a river system. From 6 m to the end of the drilling depth of 10 m, there were grey sands (Mz: 216.12 μm) with the presence of cobbles (0.5–1 cm). At depths of 3.3 m, 6–6.2 m, and 9.8–10 m, blind samples with codes G3CARD1, G3CARD2, and G3CARD3 were sampled in plastic tubes for optical luminosity analyses. G3CARD1 and G3CARD2 consisted of slightly gravelly sand, while G1CAD3 was grey sandy gravel. The stratigraphy of all drills is shown in Figure 6.

In all three locations, each core showed that the deposits in the study area had the typical characteristics of subaerial coastal land and aeolian deposits in the upper sections, and no evidence of marine deposits was found. However, the coarser material in the lower section indicated that other coastal processes may affected the area. The stratigraphy of Ammoudara beach indicates that the dune system in the area is very old, since sand deposits were found in all cores with ages reaching 11.5 ± 3.15 ka. Sedimentological finds between 1.8 m and 8 m depth consisted of sandy and gravely sandy deposits, which indicate typical coastal deposits. Few different layers were recognized, e.g., in G2 and G2′, a layer between 3.5 and 3.7 m was found with slightly coarse material and the presence of pumice.

4.3.1. OLS Results

The equivalent doses (De) for the samples were estimated by interpolation of the natural dose in the dose response curve produced following the above protocols. The average values, with the corresponding standard deviation, are given in

Table 4. Equivalent doses estimated with OSL dose response curve lines.

		G1				G2					G3
Sample Code		G1CA D1	G1CAD2	G1CAD3	G1CA D4	G2CBD1	G2CBD2	G2CBD3	G2CBD4	G2CB D5	G3CARD1	G3CARD2	G3CARD3
Equivalent Dose (Gy)	Mean ± 1sd (N = 18)	3.4 ± 1.8	9.2 ± 3.9	9.9 ± 2.7	10.4 ± 3.3	3.7 ± 1.6	6.5 ± 1.7	5.6 ± 1.1	7.0 ± 2.3	6.9 ± 1.8	5.2 ± 1.7	8.0 ± 3.9	6.3 ± 1.0

.

Quality assurance parameters of the OSL measurements of the samples were also calculated. The recycle ratio was well below 10%, and the recovery (signal of 0 Gy-induced dose compared to the natural dose) was less than 7.0% in all east borehole samples, 5.5% in central drill samples, and 3% in western drill samples. Finally, the dose recovery ratio (the estimated dose from a new dose response curve to a known given dose, i.e., ~ 15 Gy) was also less than 10%. All these values ensure the quality of the estimated equivalent doses for the samples. Recycling and dose recovery values are presented in

Table 5. Quality assurance parameters of the OSL measurements.

	Sample Code	Recycling Ratio	Recuperation (%)	Dose Recovery Ratio
G1	G1CAD1	1.06 ± 0.09	3.89 ± 3.29	0.93 ± 0.17
	G1CAD2	1.04 ± 0.10	2.18 ± 1.58	0.96 ± 0.13
	G1CAD3	1.02 ± 0.25	6.97 ± 6.02	0.91 ± 0.21
	G1CAD4	0.98 ± 0.13	2.15 ± 2.44	0.98 ± 0.19
G2	G2CBD1	1.08 ± 0.09	5.24 ± 5.88	0.92 ± 0.15
	G2CBD2	1.11 ± 0.11	1.54 ± 3.35	0.94 ± 0.24
	G2CBD3	1.08 ± 0.10	3.08 ± 2.72	0.97 ± 0.17
	G2CBD4	1.01 ± 0.14	0.67 ± 2.97	0.91 ± 0.13
	G2CBD5	1.01 ± 0.14	2.88 ± 5.41	0.95 ± 0.19
G3	G3CARD1	1.06 ± 0.05	2.76 ± 2.27	0.95 ± 0.14
	G3CARD2	1.02 ± 0.10	2.12 ± 2.00	0.93 ± 0.11
	G3CARD3	1.04 ± 0.05	2.92 ± 1.68	0.98 ± 0.19

.

4.3.2. Radioactivity Results

Measurements of the activity of isotopes of ²³⁸U, ²³⁵U, the ²³²Th series, and ⁴⁰K and ⁸⁷RB in Bq/kg and the calculated dose rate in Gy/ka are presented in

Table 6. Radioactivity measurements and calculated dose rates.

		Isotope	238U	235U	232TH	40K	87RB	Total
G1	G1CAD1	Radioactivity (Bq/kg)	17.5 ± 1	1.52 ± 0.59	11.9 ± 0.9	185.9 ± 7.9
		Dose rate (Gy/ka)	0.337 ± 0.022	0.012 ± 0.005	0.175 ± 0.013	0.494 ± 0.021	0.019	1.037 ± 0.033
	G1CAD2	Radioactivity (Bq/kg)	18.4 ± 1.1	0.87 ± 0.63	15.7 ± 1	232.1 ± 8.5
		Dose rate (Gy/ka)	0.341 ± 0.022	0.007 ± 0.005	0.25 ± 0.016	0.668 ± 0.025	0.019	1.285 ± 0.037
	G1CAD3	Radioactivity (Bq/kg)	14 ± 0.8	0.91 ± 0.44	9.3 ± 0.7	141.8 ± 5.9
		Dose rate (Gy/ka)	0.274 ± 0.018	0.008 ± 0.004	0.148 ± 0.011	0.41 ± 0.017	0.019	0.860 ± 0.027
	G1CAD4	Radioactivity (Bq/kg)	14.6 ± 0.9	0.66 ± 0.33	9.4 ± 0.7	133.8 ± 6.9
		Dose rate (Gy/ka)	0.283 ± 0.019	0.006 ± 0.003	0.156 ± 0.012	0.402 ± 0.021	0.019	0.866 ± 0.031
G2	G2CBD1	Radioactivity (Bq/kg)	12.8 ± 0.7	0.97 ± 0.36	6.8 ± 0.5	104.9 ± 4.5
		Dose rate (Gy/ka)	0.244 ± 0.015	0.008 ± 0.003	0.106 ± 0.008	0.297 ± 0.013	0.019	0.675 ± 0.021
	G2CBD2	Radioactivity (Bq/kg)	13.8 ± 0.8	1.17 ± 0.51	7.8 ± 0.6	117.9 ± 6.6
		Dose rate (Gy/ka)	0.272 ± 0.018	0.01 ± 0.004	0.127 ± 0.01	0.344 ± 0.019	0.019	0.771 ± 0.028
	G2CBD3	Radioactivity (Bq/kg)	9.9 ± 0.7	0.37 ± 0.29	8.9 ± 0.6	168.1 ± 11.2
		Dose rate (Gy/ka)	0.176 ± 0.013	0.003 ± 0.002	0.14 ± 0.01	0.477 ± 0.032	0.019	0.815 ± 0.036
	G2CBD4	Radioactivity (Bq/kg)	11.3 ± 0.7	0.36 ± 0.3	8 ± 0.6	146.9 ± 6.6
		Dose rate (Gy/ka)	0.212 ±0.015	0.003 ± 0.002	0.125 ± 0.01	0.413 ± 0.019	0.019	0.772 ± 0.026
	G2CBD5	Radioactivity (Bq/kg)	11 ± 0.7	0.41 ± 0.29	8.1 ± 0.6	148.6 ± 16.1
		Dose rate (Gy/ka)	0.189 ± 0.014	0.003 ± 0.002	0.117 ± 0.009	0.389 ± 0.042	0.019	0.717 ± 0.045
G3	G3CARD1	Radioactivity (Bq/kg)	16.6 ± 0.9	1.07 ± 0.56	8 ± 0.6	136.8 ± 6.7
		Dose rate (Gy/ka)	0.243 ± 0.016	0.008 ± 0.004	0.114 ± 0.009	0.355 ± 0.017	0.019	0.739 ± 0.026
	G3CARD2	Radioactivity (Bq/kg)	16.7 ± 0.9	0.84 ± 0.42	11.3 ± 0.8	190.7 ± 7.2
		Dose rate (Gy/ka)	0.301 ± 0.019	0.007 ± 0.004	0.181 ± 0.012	0.55 ± 0.021	0.019	1.058 ± 0.031
	G3CARD3	Radioactivity (Bq/kg)	13.2 ± 0.8	0.68 ± 0.46	9.3 ± 0.7	155.7 ± 6.4
		Dose rate (Gy/ka)	0.255 ± 0.017	0.006 ± 0.004	0.151 ± 0.011	0.457 ± 0.019	0.019	0.887 ± 0.028

.

Dose rates were corrected for the moisture content of the raw samples. Using the equivalent doses given in Table 6 (averages ± 1sd), the calculated ages and the range (min–max) are given in

Table 7. Equivalent doses, dose rates, and calculated ages in y before present.

	Sample Code	Equivalent Dose (GY)	Dose Rates (GY/KA)	Ages (Y)	Min–Max (Y)	Moisture (%)
G1	G1CAD1	3.4 ± 1.8	1.037 ± 0.033	3300 ± 1750	1500–7400	21.280
	G1CAD2	9.2 ± 3.9	1.285 ± 0.037	7150 ± 3050	3600–17,550	13.120
	G1CAD3	9.9 ± 2.7	0.860 ± 0.027	11,500 ± 3150	7800–21,000	12.523
	G1CAD4	10.4 ± 3.3	0.866 ± 0.031	12,050 ± 3800	7900–25,300	8.893
G2	G2CBD1	3.7 ± 1.6	0.675 ± 0.021	5500 ± 2400	2400–12,900	14.591
	G2CBD2	5.6 ± 1.1	0.815 ± 0.036	6850 ± 1300	4450–9500	14.519
	G2CBD3	6.5 ± 1.7	0.771 ± 0.028	8450 ± 2150	5300–13,050	11.732
	G2CBD4	7.0 ± 2.3	0.772 ± 0.026	9100 ± 2950	6150–16,450	15.477
	G2CBD5	6.9 ± 1.8	0.717 ± 0.045	9650 ± 2500	5250–13,850	22.991
G3	G3CARD1	5.2 ± 1.7	0.739 ± 0.026	6550 ± 1100	4850–8800	23.923
	G3CARD2	8.0 ± 3.9	1.058 ± 0.031	6900 ± 1350	3800–9400	12.839
	G3CARD3	6.3 ± 1.0	0.887 ± 0.028	7150 ± 1150	5600–9500	11.221

.

Distributions of the calculated ages for the samples are given in Figure 7, while individual ED values are given in

Table 8. Individual OSL ages (y).

	OSL	G1				G2					G3
Age (Y)N		CAD1	CAD2	CAD3	CAD4	CBD1	CBD2	CBD3	CBD4	CBD5	CARD1	CARD2	CARD3
1		1511	3604	7781	7892	2377	5313	4447	6170	5246	4854	3808	5586
2		1551	4862	8010	8666	3505	5671	5070	6702	5632	5416	5251	5638
3		1796	4886	8506	9271	3842	6728	5816	7009	7083	5457	5681	6080
4		2072	5014	8775	9345	4000	6873	5939	7054	7547	5675	6153	6153
5		2167	5671	9312	9980	4114	6878	6258	7064	8322	5757	6176	6205
6		2192	5738	9913	10,644	4340	7102	6355	7123	8659	5776	6309	6217
7		2195	5830	10,219	11,085	4495	7107	6362	7618	8762	5982	6361	6333
8		2211	5888	10,344	11,101	4527	7461	6433	7878	8960	6098	6406	6630
9		2373	6316	10,969	11,250	4562	7901	6443	7948	9083	6237	6506	7023
10		2709	6647	11,319	11,336	4637	8065	6446	8282	9265	6694	6855	7273
11		3069	6943	11,440	11,367	4933	8493	6816	8969	10,100	6710	7002	7456
12		3397	7233	11,463	11,979	4988	8612	7071	9009	10,264	6765	7123	7487
13		3689	7260	12,293	11,998	5749	8831	7082	9113	11,037	6891	8080	7775
14		4144	7854	12,718	12,655	6670	10,016	7286	9284	11,380	7194	8247	7825
15		5191	8376	13,543	13,949	7180	10,232	8410	11,092	12,131	7789	8319	7965
16		5202	8832	14,363	14,390	8048	11,706	8478	11,095	12,743	7792	8392	8266
17		6410	9786	14,877	15,033	8203	11,899	8878	16,066	13,644	8281	8448	9137
18		7425	17,558	20,977	25,282	12,904	13,062	9492	16,427	13,835	8799	9401	9516

.

Considering the values of the total of the disks used for the OSL measurements, the averages and the associated 1sd are shown in the 3rd, 7th, 11th, and 15th rows of Table 8. The methodology used to calculate the mean age and relative standard deviation considers the total number of tablets (i.e., pellet samples, placed in suitable tablets for passage through the machine), which is N = 18. Then, the methodology applies appropriate exclusion criteria and rejects some noncontributing tablets, namely N = 6, and from the rest N = 12, and gives values that are rounded to multiples of 50 years. Considering the values of the total number of tablets used for the OSL measurements, the averages and the relative sd were obtained.

Excluding the mean calculations, the top three (G1CAD1, G2CBD1, G3CARD1) and bottom three ages (G1CAD3 G1CAD4 G2CBD5), as well as the calculated averages, were still in the 1sd range. Samples G1CAD1, G1CAD2, G1CAD3, and G1CAD4 presented relative errors of 33.9%, 15.3%, 12.3%, and 10.5%, respectively, which were lower than 53.0%, 42.7%, 27.7%, and 31.5% of the initial values (

Table 9. Calculated age averages with the associated standard deviation considering the total number of disks (N = 18) and excluding some of the values (N = 12). Values are rounded to multiples of 50 y.

	Sample Code		18	12
G1	G1CAD1	Average ± 1sd (y)	3300 ± 1750	2950 ± 1000
		Relative error (%)	53.0	33.9
	G1CAD2	Average ± 1sd (y)	7150 ± 3050	6550 ± 1000
		Relative error (%)	42.7	15.3
	G1CAD3	Average ± 1sd (y)	11,500 ± 3150	11,050 ± 1400
		Relative error (%)	27.4	12.7
	G1CAD4	Average ± 1sd (y)	12,050 ± 3800	11,400 ± 1200
		Relative error (%)	31.5	10.5
G2	G2CBD1	Average ± 1sd (y)	5500 ± 2400	5000 ± 1000
		Relative error (%)	43.6	20.0
	G2CBD2	Average ± 1sd (y)	6850 ± 1300	6750 ± 650
		Relative error (%)	19.0	9.6
	G2CBD3	Average ± 1sd (y)	8450 ± 2150	8150 ± 1150
		Relative error (%)	25.4	14.1
	G2CBD4	Average ± 1sd (y)	9100 ± 2950	8350 ± 1200
		Relative error (%)	32.4	14.4
	G2CBD5	Average ± 1sd (y)	9650 ± 2500	9650 ± 1350
		Relative error (%)	25.9	14.0
G3	G3CARD1	Average ± 1sd (y)	6550 ± 1100	6450 ± 650
		Relative error (%)	16.8	10.1
	G3CARD2	Average ± 1sd (y)	6900 ± 1350	6950 ± 800
		Relative error (%)	19.6	11.5
	G3CARD3	Average ± 1sd (y)	7150 ± 1150	7050 ± 700
		Relative error (%)	16.1	9.9

). For samples G2CBD1, G2CBD2, G2CBD3, G2CBD4, and G2CBD5, the calculated averages were still in the 1sd range and showed relative errors of 20.0%, 14.1%, 9.6%, 14.4%, and 14.0%, respectively, which were lower than 43.6%, 25.4%, 19.0%, 32.4%, and 25.9% of the initial values (Table 9). For samples G3CARD1, G3CARD2, and G3CARD3, the calculated averages were still in the 1sd range and showed relative errors of 10.1%, 11.5%, and 9.9%, which were lower than the 16.8%, 19.6%, and 16.1% of the initial values (Table 9). Finally, typical OSL curves of the samples are given in Figure 8.

5. Discussion

In the Eastern part, core samples consistently produced slightly younger ages, mainly around 6ka. In the central part, a clearer stratigraphy is present. Ages started from around 5 ka ago at 4 m depth and 6.5ka ago at 6.5 m depth, where archaeological remains dating back to the Bronze age were found. Note that the Minoan civilization was a Bronze Age Aegean civilization on the island of Crete, and other Aegean Islands were flourishing from c. 3000 BC to c. 1450 BC. At the depth of 7.2 m, the age of the sediments appeared to be around 9ka, while at the depth of 9.5 m, the age was around 10 ka. The older age cluster of sediment ages fell between 10.5 and 12.0 ka and at a depth of 7.5–8.5 m in the western area of the beach. This period encompassed the period of the initial growth of the dune field while it was affected by the Almiros karstic spring delta. Based on those findings, the sedimentological rate for the formation of the dune fields in each area was 0.84 m/1000 y in the east part, 1.04 m / 1000 y in the central part, and 1.4 m / 1000 y in the western part, which indicates that a faster sedimentological process that may have been related to the availability of river sediments (Figure 9). To examine the evolution of the coastal landscape of the area, this work aimed to determine the rate and extent of beach evolution and how it might affect the discoverability of archaeological sites. If the sediment transfer and coastal erosion were rapid enough, they may not have been easily discernible using OSL dating as the progression of ages would fall within the ages’ standard errors, especially at specific locations. However, it appears that any dune evolution in the area, represented by the wider range of ages for the samples from 4 m depth, as well as the younger age cluster from shallower depth, did not provide any indications of age inconsistencies. At all drill sites, OSL ages indicate a correlation between depth and ages. Notably, the 4.1 m depth sample in Drill G2 is consistent with the depth of the Minoan site excavation, and the pumice layer indicates a volcanic event [105].

Based on the above findings, a potential DEM of the area was generated by inputting the tectonic movements and the sea level rise rates for the last 6000 years and combining these with archaeological data from the same area. The proposed DEM can be used as the future basis for a more detailed reconstruction to better understand the landscape evolution of the area [106]. However, there are several methodological limitations that need to be addressed. The topographic diagram of 1960, although of relatively large scale (1:5000), requires careful interpretation of boundaries and features. The map is also a record of property and agricultural land. This raises the important question of whether the mapping should be constrained by the actual features shown on a particular map or whether other sources of documentary evidence and physical factors (soil type, topography) should be used to interpret the past and create a representative rather than an actual landscape. Based on the topographic and bathymetric field measurements, the current elevations and bathymetry were created. The Ammoudara beach zone of 6ka BP appears wider by 80 m in the reconstructed DEM, and the dune field is more developed (Figure 10). The area behind the dune system had a higher altitude, while the area where the Minoan building was found at a distance from the dune system. During the last decades, there has been a noticeable loss of beach sediments, as indicated by the degradation of the dunes and the overall decline of the coastline.

Erosion is also associated with the increased distance between the reef and the coastline in the central and eastern parts of the coastal zone and the presence of coastal rocks in its western part. The loss of large amounts of sediment is intricately linked to the circulation of water on land and the relative transport of sediment, as well as human intervention and the relative rise in sea level. The former includes reductions in land sediment flows because of infrastructure development, coastal urbanization, and arbitrary disturbance of the dunes, while the latter may be related to climate change and sea level rise of >3 mm/year since 1980 [107], as well as a possible natural reduction of the coastal zone through the condensation of river deposits.

This gradual loss of marine sediments in recent decades has caused changes in the topography of the sea (depth, reef development) to which the hydrodynamic conditions have been dynamically adapted, continuing the removal of marine sediments. This loss of sediment is generally associated with sediment transport to the east, while the loss of sediment in the offshore zone is favored by reef openings in front of river estuaries. These findings agree with the conceptual model of Alexandrakis et al. [100]. In this model, it is reported that as the sea level rose, newer dunes eroded, and stable wind reserves came to the surface. Then, as the sea level reached its current height, coastal erosion removed the sediment from the beach, and finally, the beachrock separated from the beach, sank, and formed the current reef. In this model, there were multiple lines of dunes, with the older being quite stable and hosting an alluvial area behind it. As the sea level rose, younger dunes were eroded, and stable aeolian deposits were brought to the surface. Concurrently, a beachrock formation was developing along the coastline. In the areas of Giofiros and Almiros, the low topographic relief and muddy sediments act as indicators of a potential lagoon system behind the dune field.

The proposed DEM and shoreline position of the Ammoudara beach 6 ka ago is in agreement with findings from archeological excavations in the central Northern Crete. The sea level rise of 1.20 m between the Minoan Protopalatial and Neopalatial phase could be associated with a large seismic event that most likely occurred around 1700 BC [105]. The change of 1.45 m between the sea level stand of the Neopalatial phase (−2.70 ± 0.15 m) and the sea level stand with an upper age dated to the fourth century BC (−1.25 ± 0.25 m) would be either contemporary with the second destruction of the Minoan centers (~1450 BC) or in one or more paroxysmal events between 1450 BC and 400 BC. After the Roman conquest of Crete and during the subsequent centuries, 35 earthquakes from moderate-to-large magnitude (MO 6.5 to 8.2) have been reported in this region [108,109,110,111,112,113]. Particularly, the earthquake of 365 (MO 8.2) split the island and caused crustal uplift in western Crete. The tsunami triggered by this event struck coasts of most parts of the Mediterranean basin [25,26,114]. Sea level indicators depict the submerged ancient coasts of Crete. The presence of submerged tidal notches and beachrocks observed in the nearby coasts provide evidence for the existence of old shores [115]. The proximity of installations such as breakwaters, jetties, ship sheds, and fish tanks provide sufficient sea level indicators. Measurements of the emerged part with respect to the average sea level, the estimated time of construction and abandonment, and the evaluation of both the height and functional depth to the mean sea level play important roles in understanding mean sea level changes in correlation to archaeological remains [116]. In Sicily, the sluice gate of a sunken fish tank helped archaeologists deduce that the functional height of the Roman quarry’s floor (nowadays sunken) would be estimated at around 0.60 m above sea level [117]. The reported rock-cut ship sheds and fish tanks from Crete are in the vicinity of coastal quarries. The fish tanks at Chersonissos, 35 km away from the study area, could be correlated with the beachrock quarry. The floors of three fish tanks at Chersonissos, arranged in a row, now lie between 1.95m and 2.20m below modern sea level [69,118]. To the north at Agia Pelagia, 10km north from the study area, hollow channels and some holes in the shape of tubes have been spotted, which sometimes protruded up to 0,15 m from the surface of the sandstone. Those were interpreted as remnants of previous dunes that were covered with plants [119]. The same applies to the coastal quarry of Malia 40 km from Ammoudara [120,121]. Some of these holes could have been opened by the Minoans quarrymen for the insertion of poles and wooden beams [121]. Moreover, 10 km east of Amnisos, archeological findings of coastal infrastructures indicate a −2m sea level at 3500 BP [122] (for locations, see Figure 11).

Similar approaches of proposing conceptual geomorphological evolution models were proposed for the coastal region of the Sylt Island, a barrier island in the southern North Sea of North China that presented the sediment material sources during glacial and interglacial periods [123]. Moreover, in a case study in the Itapicuru coastal plain, north coast of Bahia State, paleoenvironmental information and OSL dating were used to define the possible geomorphological evolution scenario during the Upper-Middle Holocene of marine terraces and floodplain, where it indicated a slow and gradual regressive phase that was responsible for the highest marine terrace’s formation [124]. Additionally, Naples coast since the Middle Holocene was reconstructed using geomorphological analysis based on historical sources, ancient pictures and maps, and high-resolution digital terrain model (DTM) combined with iconographic research. The resconstruction led researchers to propose three new paleoenvironmental scenarios of the study area during the Holocene, which were derived from the evaluation of the relative sea level changes and vertical ground movements of volcano-tectonic origin [125].

When applying OSL dating in archaeology, one challenge is how to properly interpret the ages and their contexts. OSL could be used to date deposits that are related to otherwise defined cultural layers [126,127,128] or deposits that are directly associated with cultural material, where it essentially dates human disturbance of that sediment. Consequently, to assess the reliability of an OSL age, a detailed understanding of the archaeological context is especially critical to be able to evaluate whether enough of that disturbance had occurred to fully reset the luminescence signal. The archaeological site in the area near drills G2 and G2′ is a well-dated site which contains numerous archaeological features and ceramic artifacts [86]. Moreover, in the research drill G2′, a small ceramic was found at the 6 m depth (~4 m from current sea level) [104]. The OSL samples, which came from drill G2 above and below the dated layer, were dated to 1.4 ± 0.1 ka at both 0.3 and 0.6 m depth. The OSL results date the disturbed sand deposit. It is not necessarily surprising that OSL ages are slightly older. The assumption that the samples were disturbed by human activity, combined with the potential lack of complete sunlight exposure of the sand, make it likely that some residual signal may have been retained in the grains, causing some overestimation of the OSL age [129]. Also, it is possible that there was no more extensive human presence there and that the ~2 ka OSL ages may indicate some natural sand remobilization that occurred at that time. Whereas localized dune formation cannot be discounted, it is not possible to fully interpret these dates without further archaeological testing, as well as further OSL dating at the same depth in other parts of the Ammoudara dune field. This illustrates how important a broader archaeological and sedimentary context is for the interpretation of OSL data. Often, sparse artifact scatters are incorporated into sand deposits without clear indication of whether they are in their original context or perhaps within a buried lag surface, making the interpretation of OSL dates and their relationship with the original occupation time extremely difficult. Even if that is the case, careful examination of OSL data can provide insight into the conceptual development of a coastal landscape [95,130,131,132].

6. Summary and Conclusions

The puzzling scarcity of archaeological sites on the Ammoudara beach prompted an investigation into the chronology of the development of this dune field. The collected soil samples were dated using Optically Stimulated Luminescence to refine the chronology of the area and to assess the variability of ages with depth. The use of drilling in the dune field enabled the collection of samples from a 10 m-long stratigraphic column and the capturing of the early stages of the dune formation (OSL ages ~9.0–9.6 ka), while the OSL ages from 6 m depth represented the timing of its stabilization (OSL ages ~5–6 ka). The dune field appears to have formed long before the Bronze Age and became stabilized with only localized episodes of dune reactivation occurring thereafter. Based on the sediment texture, these localized episodes occurred mainly near river mouths. The shore-parallel reef in the Ammoudara coastal zone is a submerged beachrock that was formed during the Upper Holocene [99,133]. Its current position, about 0.5 m below sea level, indicates a relative sea level rise of more than 0.5 m, whereas the latter has caused a coastline retreat of approximately 50 m. During the last decade, the erosion rates have decreased and the beach zone of Ammoudara seems to be more stable, reaching a morphodynamics equilibrium with the nearshore hydrodynamic conditions. Luminescence dating of the aeolian deposits in the coastal area of Ammoudara yielded ages ranging from the Bronze Age to before the Bronze Age, suggesting that the coastal dunes in the area accumulated during the past 10,000 years. The study also revealed that, at a depth of 6 m below the ground surface, sediments dated back as old as the Bronze Age. Whatever the extent of the more recent dune reactivation in the study area, it was probably localized and not substantial enough to cause significant sand accumulation that would place archaeological deposits outside of the reach of standard proportionate methods. Since the geologic setting does not appear to preclude the archaeologists’ access to potential resources, the sparse distribution pattern of human occupation on the Ammoudara Beach might be better explained by exploring cultural habitation choices of past populations. Overall, the landscape reconstruction suggests that, geomorphologically, the area was rather different in the Bronze Age. Based on the archaeological and geomorphological findings, alluvial plain sediments seem to have been appropriate for agriculture since they were higher than sea level, while the location of shoreline was indicated approximately 200 m away from its current position. This can be also supported by the agricultural activities that are visible in the 1945 aerial photograph. Future work will need to address the depth to which OSL ages are adversely impacted by turbation due to wave and river processes.