Evolution of the Beach–Dune Systems in Mediterranean Andalusia (Spain) Using Two Different Proxies

Abstract

Coastal environments are complex systems that are influenced by a combination of natural processes and human activities. Scientific interest in the effects of coastal erosion/accretion and climatic change-related processes has greatly increased in recent decades due to the growing human development along coastal areas. This paper investigates the state and evolution of beach–dune systems for the 1977–2001 and 2001–2019 periods of the Mediterranean coast of Andalusia (Spain) using two different proxies: the dune toe line, which was used to track foredunes evolution, and the high-water line, which was used to assess shoreline evolution. Results showed a general erosional behavior of the studied beach–dune systems and identified cases where the main trend was altered through human interventions. During the 1977–2001 period, foredunes essentially showed erosion (54%), accretion (24%), and stability (22%) and shorelines showed accretion (40%) and erosion and stability (34% each). During the 2001–2019 period, foredunes essentially showed erosion (42%), stability (30%), and accretion (28%), and shorelines showed erosion (40%), accretion (34%), and stability (26%). Combining the evolution classes of each proxy (dune toe/shoreline) allows the behavior of both shoreline proxies to be assessed together and provides insights additional to those derived from the use of a single proxy. In this regard, Erosion/erosion (EE) and Accretion/accretion (AA) were the most frequent behaviors in the first and second periods. The results obtained provide additional insights on the nature and drivers of coastal change that aid local coastal managers and administrations in understanding erosion processes. The method can be applied to other areas around the world where a similar database is available.

1. Introduction

Coastal environments are complex systems that are influenced by many interrelated physical, chemical, and biological processes [1,2,3]. In developed coastal areas, human pressure represents an added threat to valuable natural environments [4,5,6,7], and, in many places, the risk of coastal erosion/flooding represents a severe problem for human settlements [8,9].

Scientific studies on coastal erosion/accretion processes have proliferated during recent decades as a result of the increased human interest related to coastal developments and infrastructure [10,11,12,13,14,15] and the effects of climate change-related processes [16,17,18] such as sea level rise, the increased height of extreme waves, or changes in the frequency and intensity of storms [6,7,19,20,21,22,23,24,25,26,27,28].

Many different proxies are used for shoreline change analysis depending on the particular coast study and its environment, the availability of data sources, and the aim and expected outcomes of the study (e.g., high-water line, mean high-water line, wet/dry line, vegetation line, cliff top or base, etc.) [29,30], and several studies have compared different proxies in coastal evolution studies, e.g., [29,31,32,33]. The various shoreline proxies used to determine sandy beach evolution [34,35] only capture the dynamics of a certain portion of the coastal area [36]. The use of more than one indicator is recommended by many authors [30,37,38,39], as each proxy captures different processes in the different environments of a coastal system, providing a more nuanced picture of the morphodynamic behavior and processes in a given area.

An appropriate time scale approach is also necessary to avoid associated errors, e.g., errors in the definition and extraction procedures [29,30]. Depending on the purpose of the study, different time scales can be selected, e.g., if the aim of the study is the observation of interannual changes, the analysis of shorter time periods can be sufficient, but if the aim is to predict future shoreline trends, larger time periods (>60 years) are needed [40,41].

The dry beach represents a buffer zone that absorbs, reflects, and dissipates energy delivered by waves to the shore, especially during storm events, protecting in this way the areas behind it from the impact of erosion and flooding [42] and, therefore, dry beach characteristics have to be taken into consideration in coastal sensitivity determination [43]. In addition, foredunes are one of the most relevant coastal ecosystems that work as natural defenses able to reduce flood sensitivity/vulnerability [44,45,46], as they often protect large sections of low-lying coasts against flooding during extreme storms [47,48,49]. Therefore, their maintenance/emplacement has been considered as an effective coastal protection measure that is included among possible “Disaster Risk Reduction” (DRR) strategies in several European directives [45,50,51,52,53]. Many authors agree that the temporal natural variations of sediment supply and wind and wave regimes are among the most important factors controlling the natural beach–dune system relationship [54,55,56,57]. Human impacts on the beach–dune system arise from urbanization and decreases in sediment supplies (by, for example, the construction of dams, urbanization on the coast, and coastal protection structures [58,59,60,61,62,63,64]). Regarding climate change-related processes, the sea level rise affects beach–dune systems reducing the dry beach width but also increasing the frequency of storms of any given magnitude [65]. Last, the decrease in rainfall reduces river sediment contributions producing the degradation of river deltas [66,67].

Irrespective of tidal range, beach and dune responses to erosion processes are different and due to different agents [29]. Variability in position of beach proxies is linked to marine processes and variability in the position of dune proxy is linked to both marine and aeolian processes; the former determines dune erosion and the latter determines dune accretion. Therefore, dune erosion is very fast and episodic and recovery usually occurs over several years. Beach evolution is more dynamic and variable because daily, seasonal and annual cycles of erosion/accretion and recovery is fast compared with dune recovery but ranges a lot, i.e., from days, weeks to months, according to the beach morphodynamic state [68,69].

This paper investigates the evolution of the beach–dune systems of the Mediterranean coast of Andalusia (Spain) in the medium-term (between 10 and 60 years, [40]), using two different proxies: the dune toe was used to determine foredune evolution and the high-water line was used to assess the shoreline evolution. The results obtained for each proxy were then combined to obtain different beach–dune systems states, i.e., from “Erosive” and “Mixed” to “Accretion” states. The results obtained provide enhanced information for coastal managers and administrations to properly understand and deal with contemporary erosion processes but do not allow one to make future predictions on coastal behavior. The method used in this investigation can be applied to other areas around the world where a similar database is available.

2. Geographical Setting

The 564 km long Mediterranean coast of Andalusia is located in southern Spain and administratively belongs to the Cádiz, Málaga, Granada, and Almería provinces (Figure 1). It has a rectilinear E–W orientation, with two NE–SW facing sectors located at its western and eastern ends.

The coast is micro-tidal (tidal range < 0.2 m) and is mainly composed of beaches of medium to coarse dark sand and/or pebbles. Foredunes are especially well developed in Cádiz and Almería provinces and extend along ca. 76 km of coastline [59,71,72,73].

The Betic Range, a tectonically active mountain chain that, at places, reaches to >2000 m elevation close to the coast, determines coastal orography and morphology, forming cliffs, embayments, and promontories. Several small coastal plains are present at the mouth of short rivers and seasonal streams.

Large coastal towns include Málaga (>500,000 inhabitants), Marbella (150,000 inhabitants), Fuengirola (80,000), and Torremolinos (70,000). Málaga has the densest coastal occupation in Andalusia due to the development of national and international tourism [74,75]. Along the coast, there are several marinas and the main commercial port is located at Málaga [61,62,76].

The coast of Málaga province is one of the most heavily developed coastal stretches in Europe. Tourist infrastructure development during the 1960s, without any rational spatial planning, resulted in a large population increase and transformation of the coastal landscape. Although urban developments have produced economic benefits, the environmental impacts in some cases are irreversible and the coastal sediment budget has been altered in most sandy sectors [75,76,77].

Cádiz, Málaga, and Granada have a Mediterranean climate with “Humid and Subhumid” and “Tempered Dry-subhumid” areas, with an average annual temperature from 15 °C to 19 °C and annual rainfall between 600 and 1000 mm [78]. Almería Province has a semi-arid Mediterranean climate with sparse episodes of rain (<200 mm/year in some places, [78]), and average annual temperature of 21 °C, reaching 26 °C in summer [79].

The coast is generally exposed to winds blowing from E to W and from NNE to SW in the easternmost part of Andalusia, with minimum and maximum velocities ranging from 0.4 to 9.0 m/s [80]. The wave climate and storm energy are very variable as the coasts of the Málaga, Granada, and (partially) Almería provinces are exposed both to western and eastern storms, and the easternmost area of Almería province is primarily exposed to eastern storms [70,80].

Waves show a clear seasonal behavior with storm conditions being recorded during winter (November–March) [77,80,81], and significant wave heights reach 4.73 m during extreme storms [80]. A storm characterization for the study area [80], by means of the Energy Flux, classified storm events into five classes, from weak (Class I) to extreme (Class V). Formulas used to estimate the total energy (Eⁱ_tot) of each storm were:

$$P={\displaystyle \frac{\rho g^2}{64\pi }}{T}_e{H}_{m_0}^2\left[{\displaystyle \frac{W}{m}}\right]$$

(1)

where P is the wave energy flux, or wave power per unit of wave-front length, ρ is water density, g is the gravity acceleration, T_e is the energy period that represents the period of the sinusoidal wave having the same energy as a real sea-state (for a JONSWAP spectrum is about 90% of the peak period T_p), and Hm₀ is the spectral significant wave height, and:

$$E_{tot}^i={{\displaystyle \int }}_0^{d_i}Pdt\left[{\displaystyle \frac{Wh}{m}}\right]$$

(2)

where Eⁱ_tot is the estimation of the total energy of each storm and d_i is the duration of i-th storm. The most energetic coast is between the Málaga and Almería provinces [80]. Shoreline orientation, predominant easterly winds, and associated storm waves give rise to a prevailing westward littoral drift [81]. An opposing drift is present in some coastal sectors and/or periods [77,82]. Sea-level rise is not relevant in the studied area since a negative trend was recorded near the Gibraltar Strait and Almeria and a slow positive trend was recorded in Málaga and Granada [83,84,85].

3. Materials and Methods

Material and methods used in this work are summarized in Figure 2.

3.1. GIS Project Creation

In this paper, aerial orthophotos from 1977, 2001, and 2019 were used to digitize and extract shoreline and foredune proxies. The aim of this paper is not to predict future coastal trends but rather to use two proxies for the characterization of recent coastal evolution trends that might prove more reliable indications on considered coastal environment (i.e., the beach–dune system) and useful tools in developing sound coastal management measures. Therefore, it is not necessary to use a long time-span interval since a medium-term interval is sufficient to compare the behavior of the two selected proxies and avoid uncertainties linked to short-term variations. The behavior of the two proxies is considered and, therefore, the great variability in the position of the shoreline proxy is partially compensated by the low variability in the position of the dune proxy. The orthophotos were obtained from the Web Map Services (WMS) (https://www.juntadeandalucia.es/medioambiente/portal/acceso-rediam/geoportal/servicios-ogc/web-map-service-wms, accessed on 1 March 2023) developed by the Regional Government according to Open Geospatial Consortium interoperability standards. All information was presented in projected coordinate system ETRS89 UTM zone 30N (Reference System EPSG: 25830). The spatial resolution of the orthophotos was 0.5 m for 1977 and 2001 (black and white) orthophotos and 0.35 m for 2019 orthophoto (color). In order to reduce uncertainty related to storm impacts and recovery sequences, which could affect the decision of using the selected orthophotos, shoreline displacement distances and rates were also calculated for intermediate periods between the selected photographs, i.e., 1956–1977, 2001–2010, and 2010–2016 using the data recorded in previous works and elaborated by means of the DSAS application [82]. The results showed a common trend confirming the validity of the data obtained. Furthermore, according to Crowell et al. [40] the influences of short-term fluctuations, e.g., the effects of high energy storm events, are minimized in the medium-term coastal studies and this is the case of this paper that investigates 24- and 18-year time intervals.

The two proxies selected in this study were the high-water line (HWL) for the shoreline and the dune toe line for foredunes [29,86]. The HWL, also called the wet/dry line, is identified by a change in the color of the sand. This line corresponds to the last tide mark that does not change very much in the micro-tidal Mediterranean Andalusia coast. The dune toe line was selected as it represents the line that separates the backshore from the shoreward limit of foredunes and can be determined in the orthophotos by a change in color between bare sand and vegetation. The dune proxy was digitized for foredunes >100 m in length.

3.2. Change Rate Calculations

Corrections of the shoreline position were carried out according to wave run-up and tidal conditions using the relation of the total uncertainty (σ_T) [87]:

$${\sigma }_{\mathrm{T}}=\pm \sqrt{{\sigma }_d^2+{\sigma }_p^2+{\sigma }_r^2+{\sigma }_{co}^2+{\sigma }_{wr}^2+{\sigma }_{td}^2}$$

(3)

Such total uncertainty combines digitization errors, and the photo characteristics (Moore, 2000), i.e., the digitalizing error (σ_d) determined by digitizing several times the same feature on the orthophoto, the accuracy linked to pixel size (σ_p) that is the pixel size of the orthophoto, the ortho-rectification error (σ_r) that is the root mean square error (RMSE) for photogrammetric blocks, the image co-registration error (σ_co) corresponds to the misalignment between single pixels from the set of orthophotos obtained by the rectification, and the onshore definition and position determination, i.e., wave run-up (σ_wr) and tidal conditions (σ_td) calculated for five areas in which the study area was divided and using the Formulas (4) and (5) given by Manno et al. [88].

$${\sigma }_{wr}={\displaystyle \frac{S_{wr}}{\tan \alpha }}$$

(4)

$${\sigma }_{td}={\displaystyle \frac{S_{td}}{\tan \alpha }}$$

(5)

Here, $\tan \alpha$ denotes the beach slope and $S_{wr}$ and $S_{td}$ denote the standard deviation of the values obtained in the run-up, using the Nielsen and Hanslow [89] equations, and tidal range, using the data from Puertos del Estado buoys [90].

The uncertainty of the position of the dune toe was calculated according to the first four parameters of Equation (3). The results are presented in

Table A1. Shoreline determination error for each document used (Equations (3)–(5)).

Year	Error Components (m)
	σd	σp	σr	σco	σwr	σtd	σT
1956	7.60	1.00	4.00	0.50	3.00	5.48	10.7
1977	2.10	0.50	3.30	0.50	3.00	5.48	7.4
2001	2.10	0.50	1.00	1.00	3.00	5.48	6.8
2010	1.90	0.50	0.50	0.00	3.00	5.48	6.6
2016	0.70	0.25	0.50	0.00	3.00	5.48	6.3
2019	0.70	0.25	0.20	0.00	3.00	5.48	6.3

and

Table A2. Dune toe line determination error for each document used (Equations (3)–(5)).

Year	Error Components (m)
	σd	σp	σr	σco	σT
1977	2.10	0.50	3.30	0.50	4.00
2001	2.10	0.50	1.00	1.00	2.60
2010	1.90	0.50	0.50	0.00	2.00
2016	0.70	0.25	0.50	0.00	0.90
2019	0.70	0.25	0.20	0.00	0.80

.

The coastline was divided into different coastal sectors, each of which corresponded to a single beach–dune system. Change rates of the two proxies were then calculated using the same baselines and transects to be able to pair the data at each point.

The Shoreline Change Envelope (SCE), the Net Shoreline Movement (NSM), and the End Point Rate (EPR) were calculated for two periods, i.e., 1977–2001 and 2001–2019, using the DSAS extension of ArcGIS 10.6, that takes into account the uncertainty of each proxy using Equation (1), [91]. The SCE method deals with shoreline position variability at each transect, taking into account the maximum spatial recorded displacement, regardless of the time span over which it was recorded; The NSM is associated with the dates of only two shorelines and it reports the distance between the oldest and youngest shorelines, i.e., 1977 and 2019, for each transect, although this movement may be not the maximum shoreline displacement; the EPR is calculated by dividing the distance of shoreline movement by the time elapsed between the oldest and the most recent shoreline. Baselines were digitized parallel to the shoreline and the spacing of transects (25 m) was determined in accordance with the regional scope of this paper. Following Anfuso et al. [92], the separation between transects strongly influences the results on shoreline evolution obtained by means of the transect-based analysis, with acceptable errors being achieved with transects spaced at <100 m. As the space between transects is generally defined according to the scale of the project, a distance of 25 m between transects was selected to avoid the smoothing of the shoreline considering the regional scale of this study.

A total of 53 beach–dune systems were observed during the 1977–2001 period and 38 during the 2001–2019 period, i.e., 15 systems disappeared during the second study period, essentially because of urban development [72]. Thus, in this paper, erosion/accretion rates were calculated for the 38 beach–dune systems that were observed during the two study periods. The evolution of these systems was not affected by beach nourishment programs that were only carried out at urban beaches and these lack dunes.

3.3. Statistical Analysis

The same categories of change were selected for both shoreline and dune proxies (

Table 1. Definition of classes based on EPR values.

Class	m/Year
Accretion	>+0.2
Stability	>−0.2; ≤+0.2
Erosion	<−0.2

). The stability class was established to comprise the uncertainty in the determination of the proxies position in aerial orthophotos. Therefore, the mean value of uncertainty obtained by Equation (3) for the shoreline position (Table A1), which was much higher than the value recorded for dunes toe position (Table A2), was divided by the largest time span considered, i.e., 24 years corresponding to the 1977–2001 period. This interval, i.e., ±0.2 m/year, was named the “Stability class”, because proxies variations recorded within this class were within the margin of error of the used methodology. Evidently, the used method contains a certain grade of subjectivity and slight differences in the results obtained may be observed if such limits are changed.

A total of 3234 transects were measured among the 38 beach–dune systems. During the study period, foredunes eroded, accreted or/and migrated laterally, and the number of transects varied in each study period. As the aim of this work is to determine the behavior of foredunes and the shoreline in front of them and associated relationships, only transects that intersected both the shoreline and the dune toe, were considered. This resulted in 2731 transects during the 1977–2001 period and 2654 transects during the 2001–2019 period.

The evolution classes obtained for each proxy were then combined, obtaining nine combinations of evolution states of the beach–dune systems:

• Accretion/accretion (AA): accretion classes were observed at the same transects for both shoreline and dune proxies.
• Accretion/erosion (AE): accretion was obtained for the dune proxy and erosion for the shoreline.
• Accretion/stability (AS): accretion was obtained for the dune proxy and stability for the shoreline.
• Erosion/erosion (EE): erosion classes were obtained for both shoreline and dune proxies.
• Erosion/accretion (EA): erosion was obtained for the dune proxy and accretion for the shoreline.
• Erosion/stability (ES): erosion was obtained for the dune proxy and stability for the shoreline.
• Stability/stability (SS): stability class was obtained for both dune and shoreline proxies.
• Stability/accretion (SA): stability was obtained for the dune proxy and accretion for the shoreline in the same transect.
• Stability/erosion (SE): stability was obtained for the dune proxy and erosion for the shoreline.

4. Results

Trends of transects used to determine the shoreline and foredune evolution are summarized in Figure 3. Regarding foredune behavior in the first period (1977–2001), 24% (644 transects) of the transects showed accretion, 22% (603 transects) showed stability, and 54% (1484 transects) showed erosion. In the second period, i.e., 2001–2019, 28% (755 transects) of transects showed accretion, 30% (786 transects) showed stability, and 42% (1113 transects) showed erosion.

Concerning the shoreline evolution during the first period (1977–2001), 40% (1082 transects) of transects showed accretion, 26% (715 transects) showed stability, and 34% (934 transects) showed erosion (Figure 3A). During the second period (i.e., 2001–2019), 34% (902 transects) of the studied transects showed accretion, 26% (702 transects) showed stability, and 40% (1050 transects) showed erosion (Figure 3B).

The EPR data obtained and classified in three evolution classes for foredunes and the shoreline were paired, and the combinations of the evolution classes of each proxy were determined to analyze their behavior in the 1977–2001 and 2001–2019 periods (Figure 4).

As shown in Figure 4, there is a weak positive correlation between data: the behavior of the two proxies was similar but the intrinsic higher variability in the position of the shoreline compared to the dune toe (Figure 3) makes the relation statistically weak. The bar graph shows that the EE combination of classes clearly dominates (Figure 4).

A general trend of the distribution of the different types of combinations of foredunes and shoreline evolution is presented in

Table 2. Brief description of the most important combinations of classes observed in the study area for the 1977–2001 and 2001–2019 periods.

Combinations of Classes	Brief Description
EE	Dominant in both periods, it occurred in both natural and urbanized areas (35% in both periods), in areas down-drift of coastal structures and at river/delta mouths (24% in 1977–2001 and 30% in 2001–2019) and natural promontories (41% in 1977–2001 and not observed in 2001–2019), and down-drift of anthropic structures as groins and ports (29% in 2001–2019 and not observed in 1977–2001).
AA	Widespread in both periods, this class was observed in natural areas that act as sinks for sediment eroded from adjacent areas (56% in 1997–2001 and 45% in 2001–2019), in areas up-drift of anthropic structures such as groins and ports (15% in 1977–2001 and 30% in 2001–2019), in artificially altered areas (19% in 1977–2001 and 20% in 2001–2019).
EA	Most common during the 1977–2001 period, this pairing was mostly observed to be interspersed with other combinations of classes such as ES and AA at human altered areas (29%), in areas up- and down-drift of anthropic structures (27%), river mouths (23%), rocks, and natural promontories (18%) and at pocket beaches (3%). This combination of classes is especially common in two large (mostly) natural areas in Almeria province (see Discussion).
ES	Relevant during the 1977–2001 period, this was mostly observed to be interspersed with other combinations of classes as EA and EE in the areas down-drift of river mouths (50%), anthropic structures (12%), and natural promontories (5%) and, in less frequently, up-drift of them (12%), at pocket beaches (12%) and at human-altered areas as described in the AA case (10%). This combination of classes is especially prevalent in two large (mostly) natural areas in Almeria province (see Discussion).
SS	Relevant during the 2001–2019 period, this pairing was mostly observed to be interspersed with other combinations of classes such as AA and SA in natural protected areas (52%), and with EE in areas down-drift of river mouths (28%) and anthropic structures (16%).
SA	Important during the 2001–2019 period, this situation was mostly observed to be interspersed with other combinations of classes such as AA and SS in natural protected areas (51%), down-drift of river mouths and structures (20%) and, less frequently, up-drift of them (17%).

.

The dominant combinations of evolution classes were EE and AA for both periods (Figure 4;

Table 3. Classification of the state of the beach–dune system considering the combination of erosion/accretion/stability classes recorded for shoreline and foredunes.

State of the System	Combination of Evolution Classes	Conditions of the Beach–Dune System
Erosion	EE	Severe degradation. Both proxies present erosion, indicating a severe deterioration of the beach–dune system.
	ES	Degradation. This is a signal of deterioration of the system as dunes experience erosion and the beach is stable.
	EA	Moderate degradation. This is a signal of moderate deterioration of the system as dunes present erosion and the beach is accreting.
Mixed	SS	Stable. Both proxies present stability indicating no significative changes in the beach–dune system.
	SE	Very variable. Changes may be expected as the dunes, that are not already receiving sediment, will probably register a future loss of sediment if beach erosion continues.
	AE	Variable. Minor changes are expected in the system as the dunes may stop growing or register a loss of sediment if beach erosion continues.
Accretion	AA	Good healthy conditions. The system presents very good health as both proxies indicate accretion.
	AS	Healthy. The system is in healthy condition as dunes are accreting and shoreline registers stability.
	SA	Moderately healthy. The system indicates moderately healthy conditions as dunes register stability and the shoreline presents accretion.

). Specifically, there were two areas in Almeria province, i.e., Punta Entinas-El Sabinar and Cabo de Gata, where this category occurred in association with other “minority” combinations (see Discussion).

Slight differences were evident between the 1977–2001 and 2001–2019 periods. The EE and AA trends were very similar for the two periods while the other combinations showed more changeability. The erosion of foredunes coincident with accretion (EA) or stability (ES) of the shoreline were quite frequent classes during the first period, but in the second period, the stability of the foredunes was more commonly associated with shoreline accretion (SA) or stability (SS) (Figure 4).

A classification of the evolution state of each beach–dune system was developed based on the combined evolution of the shoreline and dune proxies (Table 3).

Approximately half of the beach–dune systems recorded an “Erosion” state for both periods, and from the first to the second studied period, there was a slight reduction in the extent of the “Erosion” state and a slight increase in the “Mixed” and “Accretion” states (Figure 5).

5. Discussion

5.1. Spatial and Temporal Distribution of Evolution Classes

Generally, dune and shoreline proxies recorded an eroding trend (Figure 3), a behavior that was also reported by several authors and mainly attributed to the emplacement of dams that reduced sediment input to beach–dune systems, observed also by Prieto et al. [64] who recorded the erosional behavior of Mediterranean deltas, in Málaga [93,94] and Almería province [95], at other Spanish coasts [96,97], and around the world [63,98,99,100], and the effect of coastal urbanization, i.e., land occupation and the implementation of coastal protection structures, as observed in Málaga province by Gómez-Zotano et al. [59] and Malvárez et al. [61], in previous works in the Mediterranean coast of Andalusia [72,82], and in other places in Spain [101], and around the world [102,103,104]. Accretion was usually observed up-drift of structures and natural promontories, as well as at pocket beaches and in areas directly affected by coastal protection structures [103,105,106,107]. At few places, accretion was also observed at the mouth of seasonal streams. According to the definition used in this paper, Stability represented the expected error in the shoreline position that comprises the total uncertainty due to the characteristics of the orthophotos, wave run-up, tides, etc., and it was usually observed between erosion and accretion areas [108].

Dune erosion was the most frequent class during both periods but showed lower values during the second period during which accretion and stability slightly increased in frequency (Figure 3), possibly influenced by changes in management policies in Spain after the implementation of the Coastal Act in 1988. In the 60s and 70s, prior to the Spanish Coastal Act, rapid coastal development prioritized the tourism and urbanization of the coastal area, leading to the development of the back-beach [61,76] and a decrease and/or destruction of foredunes, especially in Malaga province [58,59,61,72,76,95,109,110,111].

The increase in the accretion and stability classes of the foredunes in the 2001–2019 period might reflect the implementation of the management and restoration measures [59,109,112], e.g., the establishment of protected areas such as Punta-Entinas and Laguna de Adra in Almería [113] or Cabo de Gata [114] that restrict public access to particular areas, or the development of specific management plans for dune conservation [115].

The distribution of shoreline evolution classes showed no important changes during the period investigated, although erosion was slightly more widespread during the 2001–2019 period compared to 1977–2001 (Figure 3).

5.2. Beach–Dune System Behavior

This large spatial scale study obtained a great amount of data that indicates marked spatial variability. This reflects the heterogeneity of the Mediterranean coast of Andalusia and the distribution of areas influenced by human activities. Furthermore, shoreline and foredunes have different behaviors and erosion/accretion processes affect them in different ways:

• Beach erosion is normally associated with winter storm events or groups of them and its recovery takes place over weeks to months during fair weather conditions, especially in summer [68,73].
• Dune erosion or disappearance may be linked to natural processes and/or human activities [60,72,73,116], and usually occurs very quickly. Foredune recovery requires months to years depending on sediment availability, the accommodation space, the colonization and growth of appropriate vegetation, and wind conditions [56,57,60,73,117,118,119,120].

Sea-level rise has a relevant influence in erosion/accretion processes, gradually limiting the beach space, and therefore the contribution of sand to the dunes. Another important factor that affects erosion/accretion processes is the space inland that dune systems occupy as they move forward as long as the wind pushes them: human intervention prevents this movement, promoting the degradation of the dunes and this is what is generally observed in the study area. The above-mentioned concepts are important and must be considered when analyzing the data obtained in this paper.

These results suggest that dunes are better indicators of coastal erosion/accretion trends because of their lower variability in position compared to the shoreline, as also observed by several authors, e.g., [30].

In general, the shoreline and foredunes showed the same trend. For example, the accretion of both proxies was recorded after the enlargement of the port of Algeciras (Cádiz province) up-drift of the structure (Figure 6A). The opposite behavior, i.e., erosion recorded by both proxies, was observed at the Guadalhorce river mouth (Málaga, Málaga province) (Figure 6B) mainly due to the critical reduction in river sediment supplies [64,94,110]. The areas accreted/eroded in the 2001–2019 period at the mentioned study sites is presented in Figure 7: in the first example, foredunes recorded a mean accretion rate of 3.40 m/year, creating an area of over 18,600 m² in front of the existing foredunes, and the shoreline advanced at a mean rate of 4.23 m/year, forming over 21,000 m² of new beach area in front of the remnant foredunes (Figure 6A); in the second example, a mean value of −1.30 m/year was recorded for the evolution of the dunes, i.e., the loss of over 19,600 m² in front of the previous dune toe line and the shoreline registered a mean rate of −1.78 m/year, i.e., the loss of over 26,000 m² in front of the remnant foredunes (Figure 6B).

The most frequent combination of classes were EE and AA in both periods followed by EA and ES in the 1977–2001 period and SA and SS in the 2001–2019 period (Figure 4). In places, a clear spatial trend was observed alongshore with EA and ES areas located between EE and AA areas, e.g., at Punta Entinas-El Sabinar (Figure 7). This is a large protected natural area where the system is accommodating to a new equilibrium because of the natural coastal processes that favor erosion of sandy cuspate forelands and sedimentation in adjacent areas [8]. In this case, an EE combination was observed on the shoreline salient and AA was observed at both sides of it where the eroded sediment had accumulated, i.e., the salient constitutes a divergent fixed limit [105,106]. The areas in between showed dune erosion and beach accretion (EA) or stability (ES) reflecting a changing trend.

In several places, a change was noted between the first and second study periods. In some instances, this was due to new coastal protection structures. For example, at Playa del Perdigal and Garrucha in Almería province, EA and ES combinations during the 1977–2001 period were replaced by an EE combination during the 2001–2019 period. The change was caused by the emplacement of five groins that caused down-drift erosion (Figure 8A II). During the second period, the behavior of this sector changed after the enlargement of the port of Garrucha. Up-drift of the port, EA combination recorded during the first studied period, switched to the AA combination during the second period (Figure 8B).

Changes attributed to natural processes were also recorded [4,5,19,23]. For example, the Artola foredunes (Málaga province) have a complex dynamics behavior [121]. The area is characterized by an alternation of easterly and westerly winds and high-energy storm events that approach from the east, generating an eastward directed littoral transport. The trend described by previous authors was reflected in the findings reported here: during the first period, EE and ES combinations were generally observed and AA combination was recorded close to the port, which had been enlarged in 1980 (Figure 9A). During the second period, SA and SS combinations prevailed (Figure 9B). The welding of nearshore bar at this sites [121], reflects the changes observed in this paper.

5.3. Considerations for Coastal Management

A combination of different shoreline proxies provides more information about the whole coastal system evolution as each proxy is able to capture certain specific processes [30,37,38]. The use of a combination of shoreline and dune proxies provides an opportunity for better understanding the beach–dune system behavior and the establishment of sound management measures and plans to counteract present coastal trends.

The combination of evolution classes calculated for the shoreline and foredunes reflects the state of each beach–dune system (Table 3) and prompts a range of potential management responses:

• Erosion states represent systems that present different levels of degradation and need management measures;
• Mixed states represent systems with diverse levels of changeability or stability conditions, which can shift to an erosion state in the short/medium term. These systems need to be monitored to fully comprehend their present and future behaviors;
• Accretion states represent systems in a good state of health, and they do not need management measures in the short/medium term.

The methodology developed in this paper represents a useful tool to determine areas that need more attention at both local and regional scale. For example, the dominance of the “Erosion” category in the 2001–2019 interval indicates that a great proportion of the beach–dune systems in the Mediterranean coast of Andalusia need some kind of management measures to avoid future erosion problems. It also identifies the particular areas where the problem exists.

Pranzini et al. [122] presented a review of the protection strategies carried out in Europe, concluding that there is no single solution to coastal problems but rather a range of practical possibilities. Many methods exist for beach and dunes restoration (e.g., [123,124]).

For example, the area of the Gualdalhorce river delta (Málaga province), which presents an “Erosion” state, is affected by several groins that stop longshore transport and cause up-drift accretion and down-drift erosion. A solution could be a sediment bypass from the accreting to the erosion areas and the emplacement of fences and the planting of endemic vegetation for foredune stabilization.

6. Conclusions

In this paper, the evolution of 38 beach–dune systems was investigated using two different proxies: the dune toe line for foredune evolution and the high-water line for the shoreline position changes. Evolution rates were calculated for the 1977–2001 and 2001–2019 periods. Results were compared, obtaining a general erosional behavior of both proxies along the studied coast and a relevant changeability for the shoreline versus the foredunes.

The evolution classes of foredunes and beach changes were paired, obtaining nine combinations of classes: Erosion/erosion (EE) was the most frequent and was recorded by 25% and 27% of the transects in the first and second periods, respectively, and Accretion/accretion (AA) was represented 17% of the total cases observed during both periods. Changes in the evolution trend of the beach–dune systems were also found, especially in areas where human interventions occurred between the first and second periods.

The use of different proxies to monitor coastal change is frequently recommended, as each proxy provides different information about each considered environment and captures different processes, e.g., marine and/or aeolian processes that are predominant in coastal environments. In this case, the shoreline proxy provides a very good reflection of medium-term changes due to its high variability in position, and the dune proxy better reflected long-term changes, e.g., the impact of high-energy storms or the accretion due to long-lasting fair-weather conditions.

The classification used to reflect the state of the beach–dune systems resulted in a useful application tool for coastal management purposes since it makes easy to recognize areas that need more attention; furthermore, its periodic updating could help verify whether the measures are working or not in the short- and medium-term time scales.

Future research could be devoted to fully understanding the behavior and/or evolution of the beach–dune systems investigated. The enlargement of the temporal resolution used in this paper would allow a better understanding of the beach–dune system changeability at different spatial and temporal scales, e.g., from months to seasons and to years.