Long and Short-Term Coastal Changes Assessment Using Earth Observation Data and GIS Analysis: The Case of Sperchios River Delta

Abstract

The present study provides information about the evolution of the Sperchios River deltaic area over the last 6500 years. Coastal changes, due to natural phenomena and anthropogenic activities, were analyzed utilizing a variety of geospatial data such as historic records, topographic maps, aerial photos, and satellite images, covering a period from 4500 BC to 2020. A qualitative approach for the period, from 4500 BC to 1852, and a quantitative analysis, from 1852 to the present day, were employed. Considering their scale and overall quality, the data were processed and georeferenced in detail based on the very high-resolution orthophoto datasets of the area. Then, the multitemporal shorelines were delineated in a geographical information system platform. Two different methods were utilized for the estimation of the shoreline changes and trends, namely the coastal change area method and the cross-section analysis, by implementing the digital shoreline analysis system with two statistical approaches, the end point rate and the linear regression rate. Significant river flow and coastline changes were observed with the overall increase in the delta area throughout the study period reaching 135 km² (mean annual growth of 0.02 km²/yr) and the higher accretion rates to be detected during the periods 1805–1852, 1908–1945 and 1960–1986, especially at the central and north part of the gulf. During the last three decades, the coastline has remained relatively stable with a decreasing tendency, which, along with the expected sea-level rise due to climate change, can infer significant threats for the coastal zone in the near future.

1. Introduction

Coastal zones are some of the most complicated, highly productive, and extremely diverse Earth ecosystems and are considered as areas of major socio-economic importance, worldwide, occupying 18% of the land surface, while approximately 60% of the human population (and 1/3 of the population of Greece) lives no farther than 2 km from the coastline [1,2,3]. They constitute an extremely dynamic environment that is subjected to changes due to natural processes, such as sea-level rise, sediment supply, wave energy, tidal inundation, tectonic setting, and floods, and human actions, such as industrial, agriculture, residential and leisure or touristic activities [4,5,6]. Those processes play an important role in shoreline change and coastal landscape evolution [7]. The shoreline change rate is the most essential characteristic of the coastal area because it reveals future harmful impacts driven by long-term climate changes and short-term episodic events of natural hazards [8].

Between the various morphological types of the coastal zones, river deltas have been identified and described by historians (e.g., Herodotus) since ancient times, for their triangular shape, which is similar to the inverted Greek letter Δ [9]. Deltas represent a highly active land-water interface that is also exposed to the strong pressure of several natural and human-caused processes interacting in a complex way [10,11,12]. In general, deltaic systems generally respond promptly to these controlling factors, which determine, according to their balance, whether the coastline is expanding or declining in the delta [13].

The eustatic, isostatic, and tectonic movements of the relative sea level, as well as the sediment volume delivered through river outflow or due to severe flood events, are the most important parameters guiding the morphodynamics of coastal systems [14]. The two most essential processes affecting coastal morphology are erosion and accretion dynamics. Erosion and accretion can occur because of human-induced or natural causes [15]. Riverine sediment amount is related to the regional geologic structure, weathering processes, river flow characteristics, land use/cover, and catchment shape and topography (i.e., slopes, slope aspect, curvature, etc.) [16].

Today’s coastal geomorphology has developed, roughly, over a time span of a few thousand years (less than 7000), when the sea-surface rise was reduced considerably, reaching its current level in some areas and ranging within a few meters of it in others. This lower rising rate allowed coastal climatic, sedimentary and tectonic processes to control its accretion, shaping various coastal morphologies. Analogous coastal geomorphologies have also appeared several times throughout the Earth’s geologic history, especially during the Quaternary [11]. In recent decades, human activities (e.g., construction of dams, distributary excavations, spillways deforestation, changes of rivers’ position or status, etc.) have become the most important controlling factor of the coastal zone evolution by inducing considerable pressures to the coastal areas, modifying their natural evolution, and causing considerable alterations in the morphology of deltaic environments [12,13,17,18]. Direct or indirect human modifications of distributary networks substantially impact river deltas evolution [19,20]. As a result, a worldwide decrease in riverine sediment supply to the coastal areas has been noted over recent decades [16,21].

Throughout history, researchers were driven by the need to assess short- and long-term shoreline change rates, which are used to define coastal development and protection zones [22]. Long-term shoreline changes exploration has tended to estimate shoreline change rates by utilizing old aerial photographs, old maps, historical documentation, and data deriving from pioneer geographers [13,23]. Aerial photographs available are dating back to the 1920s, and existing georeferenced topographic maps that extend the shoreline changes record to the mid-to-late 1800s [24,25,26]. On the other hand, the exploration of short-term changes is managed in general by field-based surveys, aerial photographs, and satellite images that offer detailed and accurate data for coastal change detection [22,27,28].

Earth observation (EO) data, passive and active, are a unique tool that can deliver efficient information for coastal change mapping and monitoring, by offering continuous, multitemporal, and low-cost airborne and satellite data in a variety of spatial, spectral, and temporal analysis [29]. Starting from the early 1970s, the National Aeronautics and Space Administration (NASA) and the Department of Interior United States Geological Survey (USGS) mission of Landsat satellites, offers the longest continuous space-based record of Earth’s surface [30]. Moreover, the European Space Agency (ESA) Copernicus Program, with the optical satellite Sentinel-2 (S2), made available a new land monitoring and classification mission [31]. These passive and freely available systems of the Landsat and Sentinel programs offer multitemporal, high-resolution multispectral data that, along with very high-resolution aerial imagery, constitute a powerful toolbox for the mapping of coastal changes [32,33]. Numerous satellite image processing algorithms and techniques have been used to overcome the difficulties in detecting shoreline position and change analysis, such as classification (supervised, unsupervised, etc.), single or multiple bands, spectral indices (normalized difference water index, tasseled cap wetness, etc.), image enhancement, principal component analysis, etc. [34,35,36]. Similarly, active microwave radar systems or synthetic aperture radar (SAR), such as ERS-2, Sentinel-1, etc., can be useful for coastal shallow water and delta mapping since their wavelengths are much longer than those of the passive sensors that can penetrate water [37,38].

Furthermore, the synergy of EO imagery with geographic information systems (GIS) can provide a significant framework for data processing and an efficient set of tools for evaluating and obtaining spatial information, which is essential in coastline changes studies [39,40].

An equally important tool in coastline mapping is the use of multitemporal satellite data, as existing maps are often neither accurate nor updated, hence unable to follow rapid coastline changes. Numerous studies, e.g., Blanco-Chao et al. [14], Aquino da Silva et al. [16], Jabaloy-Sánchez et al. [13], Petropoulos et al. [41] and Ford [22], have been conducted worldwide utilizing old maps, topographical data and EO information with different methods, techniques, and datasets, in order to monitor and map long- and short-term coastal changes.

The area under investigation was the eastern coastal part of the Sperchios River basin, and the highly dynamic delta formed at the river mouth. During the last 4500 years, the coastal area has undergone many and significant changes driven by the multiple alterations in the Sperchios River route and delta shape, due to the synergistic influence of tectonic and sediment-depositional procedures [42,43,44].

In the present study, an attempt was made to map and understand long-term and short-term shoreline changes, as well as erosion/deposition patterns due to natural phenomena such as earthquakes, flash floods, etc., further including also sea-level rise, and anthropogenic activities along the coastal area [45]. For that reason, and due to the different characteristics of the available datasets, the coastal changes that took place in Sperchios River coastal deltaic plain were divided into two distinct periods, one that spans from 4500 BC until the beginning of the 18th century and another one from 1852 to 2020. In the first period, a qualitative assessment of coastal changes with relative accuracy was conducted, while in the second period, a quantitative and qualitative evaluation of changes was made utilizing georeferenced data having accurate coordinate systems. The identification and the comparative analysis for the determination of the shoreline changes and their change rates over time were supported by a GIS platform. Two methods were utilized for the assessment of the coastal and shoreline changes and trends: (i) the coastal change area method (CCA) of erosion or accretion zones and change rates were applied for both the first and second periods, (ii) the cross-section analysis (CSA) utilizing the digital shoreline analysis system (toolbox (DSAS; a GIS Software tool for estimation of shoreline change areas and change rates), was implemented only for the second quantitative period due to the difficulty to create representative cross-sections for the complicated shoreline changes. Regarding the CSA method, two statistical techniques, i.e., the end point rate (EPR) and the linear regression rate (LRR) were implemented and their results were employed. Finally, the causes of the observed changes were analyzed and investigated.

2. Study Area

The Sperchios River basin is located between 38°44′ to 39°05′ N and 21°50′ to 22°45′ E, in the central Greece district, covering an area of approximately 1830 km² (Figure 1a). The basin has an elongated shape, orienting from west to east, delimited by the mountains of Tymfristos and Oxia on the west, Othrys on the north, Oiti, and Vardousia and Kallidromo on the south. The basin is a graben-like asymmetrical depression as a part of a tectonic trough, controlled by major NW–SE and E–W trending faults parallel to the Atalanti normal fault zone [46,47]. The tectonic activity of this graben separates the catchment into a northern and a southern part, according to the theory of tectonic dipoles [48], with the latter being lifted and the former sinking [49,50]. The northern part is characterized by mild relief and lower mean elevation, and the southern displays strong relief, steep slopes, and diverse topographic features. The river, approximately 85 km long, originates from Mt. Tymfristos (2315 m), following a west–east course [28,43]. It crosses the basin towards the flat plain area located at its eastern part, which is characterized by very gentle slopes (0–2%) while creating an exceedingly long system of meanders where finally forms an active delta, discharging into Maliakos Gulf. It has a mean annual water discharge of approximately 62 m³/s, varying between 110 m³/s (in January) and 22 m³/s (in August) [51], while the total annual sediment load of the Sperchios River has been estimated to be more than 1.5 × 10⁶ tons/year [52].

The Maliakos Gulf is a semi-enclosed embayment with water depths less than 30 m [37,43,52] (Figure 1b). The average length of the Gulf, in an east–west direction, is about 11 km, whilst its north–south width is almost 9 km [28,43,44]. It experiences low to moderate winds and low wave activity due to restricted fetches; with the presence of a weak coastal current flowing from east to west, mainly along the northern part of the gulf which is associated with distinctly weak coastal erosional processes. The tidal range varies between 0.95 m (during spring tides) and 0.15 m (during neap tides) [44].

3. Materials and Methods

3.1. Datasets

The geospatial data exploited in the present study incorporate products of different accuracy comprising old maps, comprising mostly drawings and gravures without geo-registration, and modern topographic maps of high accuracy. They also include several different Earth observation products, such as aerial photographs and satellite (optical and radar) images (

Table 1. Analytical description, temporal distribution, and origin of the datasets used in the present study, and their separation to qualitative and quantitative mode according to their accuracy characteristics.

	Type/Scale	Reference Date	Source
Qualitative Data	Paleogeographic evolution Map/no scale	4500 BC	Mariolakos et al. [53]
	Battle of Thermopylae Map/no scale	480 BC	Herodotos—History Book [54]
	Map of Strabo/no scale	1st Century	Strabo [55]
	Map of Olivier/no scale	17th Century	Based on the map of B. Du Bocage (17th century)
	Charta of Rigas Feraios (Map of Greece)/no scale	17th Century/480 BC	Web portal [56]
	Map of Martin Leake/no scale	1805	Zamani and Maroukian [44]
	Map of Laurengebrio/no scale	1825	Davanellos and Stavropoulos [57]
Quantitative Data	French Geographic Map (Imprime chez Kappelin Quai Voltaire 17 Paris)/1:200,000	1852	Laboratory of Mineralogy and Geology (Agricultural University of Athens) Archive
	Topographic Map/1:75,000	1908	Hellenic Military Geographical Service
	Aerial photographs/1:42,000	1945/1960	Hellenic Military Geographical Service
	Topographic Map/1:50,000	1971	Hellenic Military Geographical Service
	USGS EROS Archive—Declassified Data	26 July 1965, 12 July 1975	Earth Explorer/USGS
	Landsat 5 archive/1:75,000	26 July 1986	Earth Explorer/USGS
	Topographic Map/1:50,000	1990	Hellenic Military Geographical Service
	Orthophotos/1:10,000	1997	Hellenic Military Geographical Service
	Landsat 7 archive/1:40,000	28 July 1999	Earth Explorer/USGS
	ERS-2 radar images/1:75,000	18 June 1998, 14 January 1999, 21 October 1999	Laboratory of Mineralogy and Geology (AUA) Archive
	Orthophotos/1:5000	2007	Hellenic Mapping and Cadastral Organization
	Orthophotos/1:5000	2014	Hellenic Mapping and Cadastral Organization
	Landsat 8/1:40,000	2 August 2016	Earth Explorer/USGS
	Sentinel-2/1:30,000	2 August 2016	OpenAccessHub/ESA
	Sentinel-2/1:30,000	12 July 2020	OpenAccessHub/ESA

).

The cartographic data were collected from the Hellenic Military Geographical Service and several other sources, such as web portals, literature and databases of various organizations, universities, etc. (Table 1, Figure 2a).

Concerning the optical satellite data, three images of Landsat 5, 7 and 8 of 1984, 1999 and 2016, Level 1 geometrically corrected to the UTM 34N WGS84 ellipsoid with standard terrain correction applied were obtained free of charge (Path/Row: 184/033 and 183/033) through the United States Geological Survey Earth explorer portal (https://earthexplorer.usgs.gov/, accessed on 27 January 2021) (Table 1). From the same portal two digital declassified military intelligence photographs were acquired, one from the Corona/Argon/Nanyard satellite systems collection (acquisition date: 1965, DS1022-2104DA076, coordinates: 38.92/22.305, stereo medium) and the second from the KH-7/KH-9 satellite systems (acquisition date: 1975, DZB1210-500130L016001, coordinates: 38.892/23.182, stereo medium) (Figure 2b). Additionally, two Sentinel-2 cloud-free and atmospherically corrected (LEVEL 2A) images to provide bottom-of-atmosphere reflectance values, in cartographic geometry (UTM/WGS84 projection), were acquired free of charge via the ESA portal (https://scihub.copernicus.eu/, accessed on 27 January 2021).

For the detailed delineation of the basin’s downstream area, and the shoreline, multi-date very high-resolution and accuracy orthophotos and aerial photographs datasets were acquired for the time points of 1945, 1960, 1986, 2007 and 2014 from the National Cadastre and Mapping Agency SA. Furthermore, three SAR.PRI radar images of ERS-2 satellite (ascending orbit, C band, λ = 5, 6 cm and VV polarization system) were obtained with acquisition dates 18/6/1998 (orbit 16530, frame 0765, 16 bits, Lat.: N 38.587, Long.: E 22.896), 14/1/1999 (orbit 19536, frame 0765, 16 bits, Lat.: N 38.587, Long.: E 22.896), and 21/10/1999 (orbit 23544, frame 0765, 16 bits, Lat.: N 38.587, Long.: E 22.896).

The satellite image processing was accomplished utilizing ENVI (v.5.5) and SNAP (v.7.0) software (L3Harris Geospatial Solutions, Pearl East Circle Boulder Co., Boulder, CO, USA and European Space Agency, respectively), while the GIS-based analysis was made utilizing ArcGIS (version 10.7, Environmental Systems Research Institute-ESRI, Redlands, CA, USA).

3.2. Data Processing

The complete study workflow is summarized in Figure 3, demonstrating the datasets used, the processing steps and the methodology followed in the present study, which is analyzed in the following subsections.

3.2.1. Aerial Photos—Orthophoto Images

The aerial photographs were scanned and radiometrically enhanced by using their histogram breakpoint editor and were geometrically corrected, using image-to-image registration. The reference image for the geometric correction was a high-accuracy collection of 1945, 1997 and 2007 orthorectified aerial photographs (orthophotos) at 1 m/pixel resolution. The very high resolution of these orthophotos assisted in the accuracy enhancement of the rest aerial photographs of lower resolution. At first, the photographs of 1986 were registered to the 1997 orthorectified photos using the nearest neighbor resampling. Then, the aerial photos of 1960 were ortho-corrected, using the images of the former period utilizing a detailed Digital Elevation Model (20 m resolution), that derived from the topographic maps of the area (scale 1:50.000, contours interval 20 m), for the elevation information. Finally, the histograms of the final orthophotos were matched, and a mosaic image was created for each period.

3.2.2. Optical Earth Observation Data

Optical data acquired by Landsat (5, 7 and 8) satellites have relative sensor observation attributes, due to their similar technical and orbital characteristics. Nevertheless, some basic pre-processing steps were necessary for each image, to enhance the quality and accuracy of the finally used products. For instance, as atmospheric conditions can differ both spatially and temporally, standardized atmospheric models—such as the one provided by the Quick Atmospheric Model of ENVI (QUAC; an empirical atmospheric correction code that produces reflectance spectra (results) of extreme accuracy, and enables the retrieval of reasonably accurate reflectance spectra even when the sensor does not have a proper radiometric or wavelength calibration, and especially for continental areas)—were essential for the atmospheric correction of the images [58]. The preprocessing of the Landsat Thematic Mapper, Enhanced Thematic Mapper + and Operational Land Imager (OLI) data concerned the spatial enhancement (at 15 m using the panchromatic band of Landsat 7 and 8). Sentinel-2 data preprocessing comprised image resampling to 10 m since S2 spectral bands operate on the different spatial resolution of 10 m (4 bands, B2, B3, B4, and B8), 20 m (6 bands, B5, B6, B7, B8A, B11, and B12), and 60 m (3 bands, B1, B9, and B10). Furthermore, the collocation of 10 spectral bands [by excluding bands 1 (coastal aerosol), 9 (water vapor), and 10 (cirrus)] was made [59,60].

Towards compensating for the divergence in the calibration characteristics of the different satellite sensors, all data were radiometrically calibrated and converted to radiance values. The radiometric correction steps comprised the subtraction of the atmospheric contribution, reduction in illumination, view angles and terrain effects, and sensor calibration [61]. The dark object subtraction method was applied, adopting the principal that if there are areas in an image with very low actual reflectance values, any apparent reflectance should be due to atmospheric scattering, and this information can be used to calibrate the rest of the image [62,63]. The radiometric correction was implemented by using the appropriate parameters (offset/gain, sun elevation and satellite viewing angles, etc.) which are included with the satellites’ metadata documentation.

Several, transformation models, such as principal component analysis [28,43,64], and spectral indices, such as the modified normalized difference water index (MNDWI) provided by Xu [65] (an improved form of the simple normalized difference water index (NDWI) provided by McFeeters [66]) and the normalized difference turbidity index (NDTI) provided by Lacaux et al. [67] were applied to all-optical data, to better delineate the mainland from the shallow waters and the sea. The calculation of MNDWI (Figure 4a) and the NDTI and the associated spectral bands are presented in Equations (1) and (2).

$$\mathrm{MNDWI} =\left({\mathrm{B}}_{Green}-{\mathrm{B}}_{SWIR}\right)/\left({\mathrm{B}}_{Green}+{ \mathrm{B}}_{SWIR}\right)$$

(1)

$$\mathrm{NDTI} =\left({\mathrm{B}}_{Red}-{ \mathrm{B}}_{Green}\right)/\left({\mathrm{B}}_{Red}+{ \mathrm{B}}_{Green}\right)$$

(2)

where B_Green, B_Red and B_SWIR are the reflectance of the Green, Red and short-wave infrared (SWIR) bands of Landsat and Sentinel satellites.

3.2.3. Radar Data

The SAR.PRI radar images were processed, trying initially to reduce the speckle noise (salt and pepper texture) that is caused by the random effect of many small individual reflectors within a given pixel [37]. Thus, the enhanced Lee adaptive filter was applied iteratively, using two continual passes of the filter with different parameters each time (sigma multiplier and window size) [68]. Then, the radar images were georeferenced to a common map projection, using the Transverse Mercator coordinate system. Finally, a temporal differentiate image was created using the image of October 1999 as Red, the image of January 1999 as Green, and the image of June 1998 as Blue. Changes in the backscattering coefficients between the three images provoke a multicolor image that reveals, significantly, the shallow parts of the Gulf and, consequently, the delineation of the shoreline (Figure 4b).

3.3. Methodological Approach for Shoreline Change Detection

The present study attempts to examine the shoreline changes of the estuary (delta) area of the Sperchios river and the internal part of the Maliakos gulf. The analysis is based on the comparison of the current position (2020) against the shoreline positions of ten different historical periods (4500 BC, 480 BC, 1805, 1852, 1908, 1945, 1960, 1984, 1997, 2007).

Considering the topographic maps analysis and the ancient settlements of the area, Lamia (old name Zitouni), Anthili, Komma, Agia Marina, Thermopylae, and Heraklia were selected as reference points of the coastal region (Figure 1b). Τhe Old riverbed or channel of Sperchios, the old estuary between the bay “Agia Triada” and the bay “Kakolaspi”, the New riverbed or channel of Sperchios at the central part of the delta, in the position of “Bouka”, north of the Bay “Spilidi” and the Main riverbed or Spillway, the artificial estuary of the river at the northern part of the delta in the position of “Tsoumari” were used as placenames (Figure 1b).

The datasets were divided into two time periods (Table 1, Figure 5). In the first period, a qualitative assessment of coastal changes with relative accuracy was conducted from 4500 BC to 1852, by considering the old drawings and gravures, and by reviewing historical documentation and literature. These data provided useful but less accurate and not verifiable and information and thus the investigation of sea-level rise, the geological and geomorphological evolution, and related impacts on coastal processes were quite rough and of lower accuracy. In the second period, a quantitative and qualitative evaluation of changes was made by considering geolocated data of coordinated topographic maps, aerial photographs (orthophotos), and satellite images for eight different dates in 1852, 1908, 1945, 1960, 1984, 1997, 2007, and 2020. In this period a very detailed designation of the shoreline was achieved, making possible the accurate quantitative assessment of the coastal changes. The qualitative and quantitative approach of the changes, spanning from 1852 until 2020, is based on very high-resolution aerial photographs (orthophotos), topographic maps of high accuracy with a precise coordinate system, and satellite images of high radiometric, spectral, and spatial resolution. Landsat satellite imagery also was used to complement the aerial photographs in the case they did not provide complete cover of the coastal area. Finally, radar images were used to better delineate the shoreline, especially in shallow water areas.

The geometric correction of the unregistered data was made utilizing mainly the very high precision orthophotos of 1945, 1997, and 2007 and the recent high precision topographic maps. The image to image (or map to image) co-registration was made by using the closer in time data and moving from the most recent to the oldest datasets. Then, all the corrected and mosaiced (where it was necessary) raster data were handled in the geographical information system and different coverages were created following the visual interpretation of the final products. The final step regarded the digitization of the shorelines of the selected periods and the calculation of statistical attributes, such as their length, areas, etc.

For the estimation of the coastal and shoreline changes and trends, the CCA method of erosion or accretion zones and change rates was applied for both the first and second periods. The CSA analysis, utilizing the DSAS toolbox, was implemented only for the second quantitative period, due to the uneven time intervals of the first period and the big distances of the shorelines, which led to the creation of complex cross-sections and uncertain results.

The methodology used to estimate the erosion-accretion change areas along the coastline for the two time periods included the following: (a) Acquisition of shoreline position from the eleven different dates. These information layers derived from several sources such as the digitalization of cartography, the photogrammetric restitution of aerial photographs and orthophoto maps photointerpretation, and the appropriate processing and analysis of the satellite images. (b) Overlapping of shorelines in the GIS environment, corresponding to different years, to obtain erosion and accretion areas. This allowed us to distinguish and calculate the erosion and accretion areas. (c) Analysis of the evolutionary trends, defining the behavior of the system [40] (Figure 5).

The calculation of the shoreline changes rate was based on the measurement of the areas that have augmented (accretion) or reduced (erosion) over a period, according to the shoreline length. The rate of the shoreline changes was calculated based on Equation (3) developed by Doukakis [69]:

$$R=\frac{\Sigma {\rm E}a-\Sigma {\rm E}e}{L} \ast T$$

(3)

where ΣΕa is the total area of accretion, ΣΕe the total area of erosion, L is the length of an average, between two time points (approximately can be considered as the shortest of its two shoreline lengths in the respective periods) and T is the period between the two time points.

The shoreline change analysis with the cross-section method was performed by using the DSAS toolbox (v.5), derived by the United States Geological Survey (USGS) as an extension (add-in) of ArcGIS [70,71]. All the shorelines were added in a single file in a personal geodatabase. A baseline, from which the changes are to be calculated, is required in another shapefile in the same Database. The actual baseline taken is 1852. DSAS casts several transects perpendicular to the baseline and records the intersection position between the transect and each shoreline.

Then, DSAS automatically generated several statistical methods. In the present study, shoreline changes were estimated using two statistical approaches, such as the LRR and the EPR. The EPR was computed by dividing the distance of shoreline movement by the time elapsed between two consecutive time points at each transect. The LRR was used to express the long-term rates of shoreline change [13,23,72,73]. A linear regression rate-of-change statistic can be determined by fitting a least-squares regression line to all shoreline points for a particular transect. The regression line is placed so that the sum of the squared residuals (determined by squaring the offset distance of each data point from the regression line and adding the squared residuals together) is minimized. The linear regression rate is the slope of the line [73].

3.4. Fieldwork

During July and November of 2016, several field visits took place across the coastal area, and measurements of the sea-bottom depth of Maliakos Gulf were made using a boat. Sampling was performed at 14 points, which are presented in Figure 6. More specifically, the samples 3, 4, 5, 6, 7, and 8 were collected close to the new riverbed area, the samples 1, 2, and 14 in the same straight line but in greater depths, while samples 9, 10, 11, 12, and 13 were close to the old riverbed and more specifically in the area between the old and the new riverbed.

4. Results and Discussion

4.1. Overview of Long and Short-Term Coastal Changes

Long-term coastal changes, from a geological point of view, can be related to changes in the mean sea level. These changes can occur either due to tectonic shifts, which are associated with the expansion of the seabed or the movement of the crust and to climate changes, which are connected to glacier periods. Tectonic movements of the earth’s crust have caused in the geological past huge changes in the position of the coastlines. Such tectonic micro-movements still occur to the present day. These movements are also related to the seismic activity of the area where they are observed. About two million years ago, during the Quaternary, there were several glacial and intermediate glacial periods, when the ice expanded and receded. In the first case, the sea level was lowered, while in the second it was raised by the melting of the glaciers. During the last interglacial period (120,000 years ago), when the average temperature was one to two degrees Celsius higher than today, the sea surface was about 6 m higher. Conversely, during the last glacial period, 80,000 years ago, when the Earth was 5 degrees Celsius colder than today, much of the ocean water was icy and the sea level was 100 m below its current level. The initial rise in sea level and possibly the initial melting of continental ice appears to have accelerated over the last 18,000 years, and especially over the last 6000 years, at a rate of 0.2 mm per year. In light of the above, it becomes clear that sea level is an indicator of global climate change occurring on the planet because it responds immediately to these changes. Therefore, it is obvious that the formation of the coasts and the coastal zone, as well as the adaptation of the environment to the changes of the sea level and all the geological processes that followed, lasted for several thousand years and, in some cases, they continue to the present day.

Short-term changes are related to both tectonic movements and slower geological and geomorphological processes, that lead to significant changes in coastal morphology. The serious seismic events and the accretion-erosion processes due mainly to the deposition of the river sediments and several severe flood events, the wave action, and coastal currents, along with human interventions affect the evolution of the coastal zone. Phenomena causing periodical changes such as the tide, the constant blowing of strong winds or waves, and barometric changes have catastrophic consequences especially in shallow coastal areas but have a minor effect in the Sperchios coastal area. The increase in the volume of seawater due to the global temperature growth and the melting of glaciers can result in a fast and significant change in the sea level and the morphology of the coastal areas globally.

4.2. Historical Evolution of the Area

The inhabiting of Sperchios Valley dates from the Early Neolithic period (about 6000 BC) due to its strategic location. The name Sperchios derives from the Greek word “σπέρχω”, which means «to drive forward, or to move rapidly». On the other hand, the Maliakos Gulf was named after the Greek tribe of Malians that resided in the coastal area since the tenth century BC [74,75].

Several historians, such as Herodotus, geographers and travelers such as Strabo and Pausanias and other contemporary scientists have been dealing, from ancient times till nowadays, with the processes and the evolution of the Sperchios river and the coastal delta area [28,42,44,54,55,76,77,78,79]. According to their descriptions and scientific analysis, followed by information retrieved from paleoenvironmental studies conducted in the coastal plain area of the basin, it is concluded that the Maliakos gulf was submerged by seawater, about 8000 years ago and the deltaic sediments were deposited during the Holocene [74,77]. Local tectonic movements and river sedimentation processes play the most important role in the coastal area development, mostly due to periodical flood flows and less to erosion caused by sea currents or sea-level changes [42,49,80,81]. The basin has undergone several serious earthquakes and flood events throughout the last couple of millenniums [43,48,50,82]. Flood events occur frequently due to the geomorphological characteristics of the basin (very mountainous and with steep slopes, especially at its southern part) and the coastal part (intense meandering and very gentle slopes) [43,48]. The most significant ones were those of 1889, 1939, 1954, 1984, 1987, 1994, 1997, 2001, 2003, 2012, and 2015, in response primarily to sudden rainstorms or snowmelt, covering usually areas from roughly 20,000 to 100,000 ha [44,59].

Moreover, the Sperchios coastal area is associated with several important historical events such as the historic Battle of Thermopylae in 480 BC between the Greeks and the Persians, the battle between the Byzantine and Bulgarian armies in 997, and the Battle of Alamana between the Greeks and Ottomans during the Greek War of Independence in 1821. The most famous of them, according to the Historians, took place there because of the narrowness of the passage between the sea (Maliakos gulf) and the steep slope of Mt. Kallidromo, allowing King Leonidas and his three hundred Spartans warriors to stand against the massive forces of the Persian Empire. The description of the exact location and the conditions of the battle terrain were documented by Herodotus nearly half a century (or more) after the battle [83]. Back then, the historical narrow passage was then approximately 90 m wide, while today is approximately several kilometers larger due to the coastal expansion [40]. Both catastrophic and progressive geological processes occurred, leading to massive changes in the morphology of the Sperchios river and its estuary, and consequently to the “passage” of Thermopylae and the head of the Malian plain. The shoreline has prograde eastward into the Gulf of Maliakos at least 15 km within the last 4500 years [43,74].

According to the Davis theory for Earth basins’ evolution, each river is initially in its infancy and forms its bed through backward erosion [84]. The youth stage takes place mainly in mountainous areas and is characterized by erosion. It then goes through stages of maturity when the relief forms progressively lose their sharp character and appear rounded. A flood plain reflects the stage of maturity of a river. The loss of river transport energy/capacity due to climate and tectonic changes are key factors in this trend. Eventually, the river system flows through wide valleys, creating meanders. This final form of the basin and river system is called the ageing stage [76,85,86]. Today the Sperchios river is in a stage of maturity, while some tributary rivers at the basin’s southern part, such as the Gorgopotamos and the Assopos, are still in the youth stage. In such conditions, the accretion or erosion of the coastline is principally due to the neotectonic activity, the sedimentation processes, climate change, and the direct and indirect human actions at the coastal area [48,87].

4.3. Coastal Area Status, Shoreline Position, and Drainage Network Development from 4500 BC to 1805

The qualitative analysis of the first period concerns the span from 4500 BC till 1805 where the data were obtained from historical documents references and old maps or gravures. Based on them and utilizing significant descriptive information deriving by scientists, historians, and geographers the changes of the coastal area were approximately accessed.

The sea level elevation that took place during the melting of glaciers, has been attested by many researchers [81,88,89,90]. They were estimated that the valley of Sperchios was inhabited from the 6th millennium BC, during the Neolithic era, and they claim that the Neolithic settlement of Lianokladi—Amouri (very close to the city of Lamia)—must have been very close to the sea.

As well, in the descriptions of Herodotus in 480 BC [42,54], it is mentioned that the river Sperchios flowed to the north side of the Maliakos Gulf, while about 3.2 km (20 stadia; the stadion was an ancient Greek unit of length that, according to Lev Vasilevich Firsov, who made an empirical determination of its length by comparing 81 distances given by the Eratosthenes and Strabo, was about 157.7 m [91,92]). South of the river Dyra (its name today is Gorgopotamos) was the river Melas (now called Xerias). About 1 km (5 stadia) ahead was the village of Trachina, located at the hills around the Asopos River, on the southeastern edges of the valley [93] where the Persian army encamped in 480 BC, since it was the valley’s widest area. These three rivers (Dyras, Melas, and Asopos), at that time, flowed directly into the Maliakos Gulf in contrast to nowadays, where they are tributaries of the Sperchios river. After the Asopos River, in a distance of 4 km, the passage of Thermopylae was situated. Additionally, in the general description given by Herodotus, for Thermopylae, it was mentioned that to the west there was a mountain, inaccessible, high and steep (Oiti-Kallidromo), while to the east the road was almost in contact with the sea which had very shallow waters [54,76]. Furthermore, Strabo (who lived from 67 BC until 23 AD) [55], mentioned the same in his descriptions, saying that the Sperchios river flowed between Lamia and Thermopylae and was about 6 km (30 stadia) away from Lamia crossing a plain that outflows into the Maliakos Gulf, stating also that Thermopylae was seventy stadia away from the estuary of Sperchios. The mapping of the area for this period (480 BC) was attempted by Strabo [55], as well as the map (Charta) of Rigas Feraios [56], which was created around 1796, based on older descriptions and depictions of the area.

For the period that followed and until the beginning of the 18th century, references to the area were made around the 10th century with the battle between the Byzantines of Emperor Basil II and the Bulgarians of Tsar Samuel (997 AD) [94]. According to the historical descriptions, maps, calculations [95,96], and the report of Vortselas [97], it was estimated that from the period of the Persian wars until the battle of 997 AD, the shoreline of the Maliakos Gulf shifted more than 6 km to the east. At this time, the Sperchios River still flew into Maliakos Gulf at the north part of the basin, while a lake of considerable size (Eropoli) existed along a tributary of Mavroneri River.

Pococke [98], in 1745, mentioned that in 1740 Sperchios River was situated at the north part of the basin and that Zitouni (Lamia), was 4 miles northwest of the MaliakosGulf. He also observed an earthquake that occurred in 1758 [74,99], estimated as magnitude ML = 6.5 on the Richter scale that could be responsible for the alteration of the river course from the northern to the southern part of the basin, near Thermopylae [43]. Moreover, Gel in 1805 [57] writes that the river Sperchios had made a great turn to the south and flowed into the southern part of the bay, and that there are two artificial canals (one wider) from the river to the sea at the northern part that drains the coastal area. Additionally, Francis Puckeville (1806) [57] referring to older records and says that although in the time of Strabo the distance of Lamia from Sperchios was estimated at 5–6 km (30 stadia), due to the change in the flow of Sperchios, the distance had become 16 km. Due to this change in the watercourse flow, the south side has started to be covered with sediments, the first protrusions of the delta front began to form and the river appeared to have meanders, which radically change the shape and properties of the riverbed and the coastal area, provoking frequent flood events [59,100].

From these descriptions, it is obvious that at the mid-end of the 18th century, the Sperchios River flowed at the southern part of the basin close to Thermopylae, and consequently the rivers Asopos, Melas and Dyras were then tributaries of Sperchios. The analysis of these coastal changes was cross-evaluated with several other studies of the past, such as those of Tziavos [76], Zamani and Maroukian [44], Sigalos and Alexouli-Livaditi [101], Kraft et al. [42] which provided similar results in shorter period studies. According to the above, the possible positions of the shoreline and the Sperchios riverbed from 4500 BC to 1805, are presented in Figure 7a–d.

4.4. Coastal Area Status, Shoreline Position, and Drainage Network Development from 1852 to 2020

Starting in 1852, and taking as a basis the first known topographic map with an established coordinate system (French Geographical Map; Imprime chez Kappelin Quai Voltaire 17 Paris), the shoreline along with and the Sperchios river and its tributaries situations were digitized (Figure 8a). Once again, this map shows the new course of Sperchios, which follows a major turn, starting after the village of Komma and forming several meanders to the south, discharging in the southern part of the Maliakos Gulf. In this area of the delta, at least three different protrusions (not of the same size) of the delta can be observed in the south part [102]. These delta estuaries, which were quite close to each other, indicated that probably the river, due to a variety of severe flood phenomena or tectonic movements, was constantly changed its outflow point several times. The river’s position, according to the topographic map of the Hellenic Military Geographical Service has shifted again, discharging at the central part of the gulf. In the period from 1852 to 1908, some very significant natural phenomena took place, such as a massive flood event that occurred in 1889 causing the failure of the river’s embankments, as well as the strong earthquake (ML = 7.0 magnitude on the Richter scale) in the area of Atalanti (about 45 km southeast of the study area) in 1894 [99,103] that is mentioned by several researchers as an important historical event (Figure 8b) [104,105,106,107]. One of these two occasions (or both for that matter) played a decisive role in the evolution of the delta area at this time [102,108]. In the same map, the various attempts to construct an artificial drainage network in the northern part of the gulf to reduce the flood phenomena generated mainly by the multiple meandering of the river and the very gentle (almost flat) slopes of the coastal part, were easy to notice.

For the periods from 1908 to 1945 and 1960, the mosaics of the orthophotos and aerial photographs of 1945 and 1960, respectively, were utilized to examine the coastal area progression. In 1944, a ditch (named German) close to the city of Lamia was constructed to drain, the Xerias stream (a big tributary river of Sperchios River) basin driving its flow to the north part of the gulf. The great sedimentation process and the extended accretion of the delta mouth in the new position of the river in the central part are apparent in the datasets of 1945 (Figure 8c). This strong embankment is based both on a large amount of sediment transportation from the river as well as on the fact that this area of the gulf was quite shallow, especially near the shoreline and therefore the filling of this part with sediments was easier and shorter in time. Additionally, more intense use of the coastal zone by human activities, such as cultivations, creation of a modern and denser road and drainage network, etc., can be observed. Between 1945 and 1960, the Sperchios river main flow was diverted at the northern part of the gulf (Figure 8d) A major technical drainage-flood relief project related to the construction of a large Spillway canal, having a width of almost 20 m, was accomplished in 1957–1958. The creation of the spillway was started a few years before having its gate barrier next to the village Komma. The canal direction followed the old riverbed of Sperchios, at the north part of the valley, and was outflowed at the north part of the gulf, receiving large amounts of water and aiming to reduce the flood phenomena that occurred in the coastal area. As a consequence, more sediments were once again transported by the river once again to the northern part of the gulf and less to the central part which caused the accretion of the north part (0.31 km²) and the erosion of the central part (0.13 km²). At the point of the new estuary, two river mouths were remained active, flowing into the Gulf, instead of the three mouths that existed in 1945. Furthermore, in 1960, due to the extended construction of new irrigation and drainage canals, the coastal zone has shifted drastically to highly productive agricultural land with the expansion and dominance of irrigated arable crops, a fact that led also to the increase in the rural population at the coastal zone.

During the period from 1960 to 1986, mainly satellite images and photographs deriving from satellite and aerial platforms were mainly utilized for coastal change monitoring. In 1986, the protrusion of the spillway mouth started to become particularly distinct in the northern part of the gulf where a large section of land was created (Figure 8e). The wave-originated erosion at the new riverbed mouth, at the central part of the gulf, was pronounced particularly intense, also due to the reduction in the sedimentation, since a large amount of Sperchios river water supply, along with its sediments are now led at the northern part of the gulf. Hence, from the two previous river mouths, the northern one was completely eroded, and the river now only had the outflow branch of the southern part. At the same time in the two old estuaries, on the south side of the gulf, there are signs of a slight decrease due to the erosion processes, but the general picture gives a clear indication that in the coming years there will probably be a greater retreat. Likewise, there is a greater intensification of the cultivations (summer crops), accompanied by an extended improvement of farm conditions (more well-maintained irrigation systems, drainage canals, roads) as well as the increase in constructions, such as buildings, warehouses, agricultural units, etc.

For the coastal condition monitoring of 1997, the highly accurate dataset of orthophoto maps was utilized. Slight changes in the shoreline were observed, from 1986 to 1997, concerning the widening of the spillway, low accretion of the river’s mouth area and a minor erosion mainly on the sides of the new riverbed estuary (Figure 8f). For the description of the coastal zone concerning the periods from 1997 to 2007 and 2020, orthophoto maps and Sentinel-2 images were employed, respectively. In the period from 1997 to 2007, several constructions of significant size and importance occurred, which affected, directly or indirectly, the coastal area (Figure 8g). More specifically, the new National highway, in the part from Agios Konstantinos to Stylis, was constructed [43] and at the same time, the new modern railway line was formed having an almost parallel direction to the highway, crossing the coastal part and the delta area. Furthermore, trying to support these works, multiple drainage culverts were created or improved. The most significant was the formation of the new artificial distributor next to the Komma village (Figure 8h). The design of this distributor was intended to divide the water flow primarily to the spillway and less to the new riverbed (ecological supply) and to radically decrease the river’s water velocity to reduce the flood phenomena that can cause problems both to the New highway and railway constructions and also to the adjacent villages of the coastal area (Anthili, Moschochori, Komma, etc.) [109]. Consequently, this new distributor drastically interrupts the deposition of sediments in the delta area of the river, since it retains most amounts of them.

This condition, along with the fact that the outflows of river Sperchios (new riverbed and spillway) have reached deeper parts of the Maliakos Gulf, the reduction in atmospheric precipitation due to climate change as was presented by Psomiadis [43], as well as the increase in irrigation and fresh-water needs (due to the intensification of the cultivations, the increase in summer water-consuming crops and the increase in population, tourism and industrial activities) [43,110], has led to the considerable decrease in water and sediments amounts towards the delta area [111,112,113,114]. As a result, during the period from 1997 to 2007 only a few negligible changes were noticed in the shoreline, concerning a very small accretion of the shoreline in the spillway and the new riverbed outflow area. On the other hand, in the span from 2007 to 2020, a noteworthy stabilization of the shoreline was observed, showing at the same time a tendency for the shoreline to recede in the coming years.

4.5. Coastal Changes Rates and Qualitative Analysis Results from 4500 BC to 1805

As it was previously mentioned regarding the changes and the change rates assessment for the first period (4500 BC to 1805) only the CCA method was applied due to the uncertainty of the available historical data. The areas of accretion and erosion that took place among each time-pair of the shoreline position, which corresponds to two different and consecutive time points, were assessed and the calculation results of shoreline changes are presented in

Table 2. Shoreline area changes and annual change rates for the second quantitative period from 4500 BC to 2020.

Periods	4500–480 BC	480 BC–1805 AD	1805–1852	1852–1908	1908–1945	1945–1960	1960–1986	1986–1997	1997–2007	2007–2020
	1st Period			2nd Period
Years—difference	4020	2285	48	57	38	16	27	11	10	13
Accretion (+) (km2)	43.93	69.72	12.79	6.19	4.69	2.13	3.99	0.21	0.08	0.02
Erosion (−) (km2)	0	0	2.16	3.95	0.83	0.66	0.79	0.14	0.07	0.03
Area—Difference = E (km2)	43.93	69.72	10.64	2.24	3.86	1.47	3.20	0.07	0.01	−0.01
Coastline length = L (km)	34.5	24.78	24.78	35.44	42.63	48.69	42.43	41.32	41.32	41.32
Coastline Change = E/L (m)	1273.3	2813.6	429.38	63.21	90.55	30.19	75.41	1.69	0.24	−0.24
Annual Rate of Change	0.31	1.23	8.95	1.11	2.38	1.89	2.79	0.14	0.02	−0.02

.

In the period from 4500 to 480 BC and 1805, the accretion covers an area of 43.93 km² and 69.72 km², respectively (Table 2; Figure 9a,b), but the change rate appears quite low being a little higher during the second shoreline pair, from 480 BC to 1805 (1.23). These enormous accretion areas are explained by the fact that the basin evolution was probably at the first stages of its development, and the erosional and sedimentation processes, were very high, without any disturbances from human activities (Figure 10a–c).

4.6. Coastal Changes Rates and Quantitative-Qualitative Analysis Results from 1852 to 2020

As it was previously mentioned in the methodology section, the approach of the shoreline changes, covering the period from 1852 to 2020, was based on maps having coordinate systems, aerial photographs and satellite images, which enabled the quantitative, except the qualitative, assessment of the shoreline variations, as well as the calculation of the rates of changes. Therefore, both CCA and CSA methods were applied.

The results of the CCA method (Table 2), indicated a continuous alteration of the coastal area and the shoreline in different parts of the delta. These changes were related mainly to natural phenomena (several flood events and the strong earthquake of 1894) and less to human interventions (spillway construction) and the corresponding shift of the river outflow during these years, demonstrating an overall accretion of the delta with the change rate values ranging from 1.1 to 2.79. Therefore, the higher accretion of the delta appeared during the periods of 1852–1908 and 1960–1986 (3.86 and 3.20 km², respectively) where the river shifted from the south to the middle part of the gulf (1894), and the new spillway started to bring greater amounts of water and sediments to the northern part of the gulf (1958; Figure 10d–j).

As far as it concerns the outcomes of the CSA approach, the dynamics of the shorelines, over the second study period, were calculated using the LRR and EPR methods. The shoreline of 1852 was utilized as a baseline, setting random intervals for the transects not to have intersections, due to the very complicated shape of the shorelines. DSAS generated transects along the shoreline of the study area, oriented perpendicularly to the baseline recording in parallel the statistical methods of LRR and EPR. The results presented in Figure 11 demonstrate the overall shoreline change rates calculated from their analysis. The positive values indicate the accretion of the shoreline, while the negative values are related to coastal erosion.

The LRR method provided 45 transects (Figure 11a) that revealed an accretional trend of the shoreline, presenting a maximum of 38.4 m/year located mainly at the central and north part of the gulf, while the erosional trends appear mainly in the south part, having a maximum value of −5.7 m/year. The average rate of the LRR method was 1.22 m/year. The implementation of the EPR method provided 44 transects (Figure 11b), which specified a maximum increase in 35.5 m/year at the central and north part of the gulf and a decrease with a highest of −6.7 m/year, while the average rate was 1.41 m/year.

The two methods demonstrated similar results for the second quantitative period, revealing the overall accretion of the area. The LRR method provided lower average rates probably due to its susceptibility to outlier effects. In general, LRR tends to underestimate the rate of change relative to other statistics, such as EPR, as was noted by Salghuna and Bharathvaj [73,115].

5. Conclusions

The Sperchios River coastal area has undergone too many changes through time due to fluctuations of sea level, tectonic activity, high rates of sediment deposition, the local geomorphology and extensive human interventions. The long-term changes are related to the sea-level rise, sedimentation and tectonic movements, while the short-term changes correspond mainly to human activities.

The two methods used for the calculation of the change rates, CCA and CSA, provided similar results regarding the change rates and the spatiotemporal distribution of the changes. LRR and EPR statistical approaches of the DSAS software displayed analogous results.

The deltaic plain area increased by ~135 km² between 4500 BC and 2020 AD, corresponding to a mean annual growth of 0.02 km²/yr. The most significant shoreline change rates were recorded between 1805–1852, 1852–1908 and 1960–1984, especially at the central and south part of the gulf. During the last three decades, the coastline has remained relatively stable, manifesting a minor regression over the last decade. Engineering interventions at the coastal zone, as well as the reduction in the precipitation amounts, and the rise of irrigation and fresh-water needs (due to the expansion of water-consuming crops, population growth, tourism increase and industrial activities), have contributed significantly to the accretion/erosion balance in the deltaic area. The notable decrease in water discharge and consequently the sediments deposition in the coastal area has markedly decreased, and, along with the compaction of the sediments at the shallower areas, led to the predominance of wave corrosion over deposition and the slight receding of the coastline.

Irrespective of delivered results and analyzed parameters, the most vital aspects of a coastal changes study are the availability and proper processing of geospatial data. The accessibility to a plethora of remote sensing imagery from the second half of the previous century and the powerful characteristics of the GIS, which allow the accurate comparison of the georeferenced geospatial layers, provide an exceptional background for studies of spatiotemporal change detection. Remote sensing optical and radar data and their processing products, such as spectral indices, offer unique tools in coastal areas monitoring and management.

Conclusively, the major changes in the configuration of the study site’s coastline area in recent decades, are primarily linked to human activities, leading to the stabilization and regression of the shoreline position. This fact, along with the expected sea-level rise due to climate change, between 0.4 m and 0.63 m in 2100, according to the predictions of the Fifth Assessment Report of the Intergovernmental Panel on Climatic Change (IPCC) can infer significant threats to the coastal zone, local communities, infrastructure, and the coastal ecosystems’ equilibrium. Therefore, coastal evolution monitoring will assist policymakers and stakeholders to take action and apply appropriate measures for the protection of the inestimable coastal zone.