Monitoring Coastal Vulnerability by Using DEMs Based on UAV Spatial Data

Abstract

The use of Unmanned Aerial Vehicles (UAVs) represents a rather innovative, quick, and low-cost methodological approach offering applications in several fields of investigation. The present study illustrates the developed method using Digital Elevation Models (DEMs) based on UAV-derived data for evaluating short-term morphological-topographic changes of the beach system and related implications for coastal vulnerability assessment. UAV surveys were performed during the summers of 2019 and 2020 along a beach stretch affected by erosion, located along the central Adriatic coast. Acquired high-resolution aerial photos were used to generate large-scale DEMs as well as orthophotos of the beach using the Structure from Motion (SfM) image processing tool. Comparison of the generated 2019 and 2020 DEMs highlighted significant morphological changes and a sediment volume loss of about 780 m³ within a surface area of about 4400 m². Based on 20 m spaced beach profiles derived from the DEMs, a coastal vulnerability assessment was performed using the CVA approach that highlighted some significant variations in the CVA index between 2019 and 2020. Results evidence that UAV surveys provide high-resolution topographic data, suitable for specific beach monitoring activities and the updating of some parameters that enter in the CVA model contributing to its correct application.

1. Introduction

The availability of Digital Surface Models (DSMs) characterized by high spatial resolution and vertical accuracy is becoming more and more important for all studies involving three-dimensional analyses and representations of the environment. In the field of coastal monitoring, for example, topographical high-resolution data are essential to perform reliable reconstructions of short-term coastal dynamics and evolution, as well as erosion and flooding risk assessments.

Most traditional methods for mapping and monitoring coastal environments, such as visual interpretations of topographic maps and aerial photos, require a lot of time and effort, and may not cover the most recent time intervals. Remote sensing techniques, such as Interferometric Synthetic Aperture Radar (InSAR), Airborne Light Detection and Ranging (LiDAR) and other geophysical techniques (gravimetric, electrical and geomagnetic methods), have meanwhile generated large amounts of high-accuracy and high-density topographical measurements over the last 40 years, allowing for regular updating of coastal environments and automated mapping [1,2,3,4,5,6]. Some of these techniques generate high-resolution Digital Elevation Models (DEMs), providing a high level of detail that allows the delineation and analysis of small-scale geomorphic elements [7,8]. However, these technologies have significant costs that make them prohibitive for several countries and institutions and difficult to be implemented on a regular basis [9].

Techniques based on Global Navigation Satellite Systems (GNSS) and total stations are among the high-quality methods available today for the reconstruction of beach and dune morphometry [10,11,12]. In particular, GNSS location (based on Network or Real-Time Kinematic) is precise and efficient, but the number of measurable points is restricted. Furthermore, Terrestrial Laser Scanner (TLS) applications might be extremely precise but need extensive survey sessions and considerable data processing times [13].

Recent advancements in Unmanned Aerial Vehicle (UAV) technologies made it possible to develop an alternative technique for coastal monitoring that efficiently captures the spatial and temporal requirements across a wide range of environmental applications, overcoming the limitations of the above-mentioned methodologies [14]. Furthermore, UAVs offer high versatility due to the significant cost reductions of aerial vehicles and sensors, as well as the recent developments in data processing software [15].

A summary with main references regarding the use of these methodologies is proposed in

Table 1. Main methodological approaches used in recent years for monitoring coastal studies. In the References column the listed papers are reported in chronological order.

Approaches	References
Airborne LIDAR	Chust et al., 2010 [1]; Stockdon et al., 2002 [4]; Coveney et al., 2010 [7]; Schmid et al., 2011 [8]
Satellite	Tralli et al., 2005 [2]
Terrestrial Laser Scanner	Rosser et al., 2005 [3]; Nield et al., 2011 [10]
GNSS	Coveney et al., 2010 [7]; Di Paola et al., 2014 [11]; Di Luccio et al., 2018 [12]
UAV	Bryson et al., 2013 [16]; Mancini et al., 2013 [13]; Turnet et al., 2016 [17]; Manfreda et al., 2018 [14]; Flores-de-Santiago et al., 2020 [15]

.

The review by Colomina and Molina [18] provides a detailed description of the evolution and state-of-the-art of the use of UAVs, country-specific restrictions, and possible applications to photogrammetry and remote sensing. UAV surveying and mapping among others include geological resource mapping [19], geomorphological process mapping [20,21,22,23], agricultural watershed studies [24], mining [25], archaeology [26], firefighting [27], forestry [28], and a variety of more traditional cadastral mapping jobs [29].

For more than half a century, stereo-photogrammetry was used to interpret aerial imagery of coastal zones. Recent and major advancements in digital photogrammetry and computer vision have propelled the possible applications at a breakneck pace [17]. Particularly, Structure from Motion (SfM) is a photogrammetric tool already in use for a couple of decades, which allows reconstructing three-dimensional surfaces from a set of overlapping images. Westoby et al. [30] describe the SfM algorithm with a technical detail approach and provide a comprehensive overview.

The use of UAVs for coastal studies has just started [13]. Harwin and Lucieer [31] compared the accuracy of point clouds generated from UAV photography of a coastal area with a discrete set of points collected by a total station. Bryson et al. [16] studied intertidal rocky landscapes with kite aerial photography and SfM, providing a quantification of the accuracy of final products.

This paper deals with the methodological approach based on UAV-derived data that we have developed to provide a valid and economically convenient alternative for coastal monitoring and the assessment of coastal vulnerability to erosion and inundation, defined as the hazard or susceptibility of a specified coastal system [32,33]. We tested this approach on the Molise coast (central Adriatic coast), a well-known coastal system under continuous study for some decades, along a small coastal stretch that we selected based on its conditions of shoreline instability and recent erosion trend.

2. The Molise Coast

The Molise coast is part of the central Adriatic coast and is limited to the north and south by the mouths of the Formale del Molino channel and the Saccione Stream, respectively (Figure 1A). It is approximately 36 km long and, except for its central 13 km long portion, is made of low coast. The sea cliffs in the high coast sector are mainly made of Plio-Pleistocene clayey-sandy marine successions (Montesecco Clay and Serracapriola Sands formations [34], overlain by Middle Pleistocene sandy to conglomeratic fluvial deposits (Campomarino Conglomerates, [34]). These cliffs are only locally in contact with the sea (protected cliff of Termoli promontory [35]). Moreover, they are inactive and located behind beach-dune systems between 10 s to approximately 200 m from the shoreline. Therefore, beaches are mostly continuous all along the Molise coast, and only locally interrupted by major river mouths and harbour structures. The beaches are sandy, prevailingly made of fine to medium-grained sands, and extend to variable widths between a few and 93 m [35]. Three harbour structures are present along the Molise coast (Figure 1A): the Margherita Sveva touristic port, the commercial harbour of Termoli with the connected adjoining San Pietro touristic port, and the Campomarino (Marina di San Cristina) touristic port.

The Termoli promontory separates the Molise coast in two, nearly independent sectors (physiographic sub-units, [35]). These sectors, hereinafter named the northern and southern Molise coast, are oriented in the WNW and NW-SE direction respectively. Recent investigations by Rosskopf et al. [35] and Buccino et al. [36] have confirmed a dominant longshore transport from NW to SE, consistent with previous studies.

In agreement with previous studies ([35] and references therein), the subdivision of the southern Molise coast in coastal segments S6–S9 (Figure 1A) has been maintained in the present study to better compare shoreline change rates calculated up to date.

The wave climate in the Adriatic Sea is usually mild or moderate most of the year. Analysing wave data recorded at the Ortona buoy (Sea Wave Measurement Network, RON) during the period 1990–2006, reveals annual maximum significant wave heights ranging between 3 and 6 m, approximately, approaching from the N-NE directional sector ([35,36]). In detail, the average values of H_s and H_t (related to sea states exceeding 2 m) and the corresponding periods T_s and T_t are given in

Table 2. Main wave climate features estimated for the Molise coast. Average significant wave height (H_s) and corresponding period (T_s); Average significant wave height of events exceeding 2 m (H_t) and corresponding period (T_t).

Recording Period (Ortona Buoy)	Main Wave Direction (°N)	Secondary Wave Direction (°N)	Effective Fetch (km)	Hs (m)	Ts (s)	Ht (m)	Tt (s)
1990–2006	340–10	70–100	476	0.7	3.5	3.5	6.6

.

As highlighted by several studies ([35,37] and references therein) erosion strongly controlled the long-term evolution of the Molise coast. Shoreline retreat especially affected the coastal segments that include the major river mouths (Trigno and Biferno rivers, Figure 1A) due to the decrease in fluvial sediment delivery to the coast mostly resulting from in-channel mining and hydraulic interventions at the basin and fluvial reach scale from the 1950s onwards ([36,38] and references therein). To face the ongoing coastal erosion, numerous defence structures were built from the 1980s onwards mostly in emergency conditions. In fact, hard defence structures (adherent breakwaters, emerged and submerged detached breakwaters, revetments, groins and jetties) cover about 62% of the Molise coast [35].

Regarding in particular the southern Molise coast and its coastline changes in recent years, a progressive destabilization is highlighted [35], with annual retreat rates of −0.50 m/y, −0.60 m/y and even −3.57 m/y affecting segments S7, S8 and S9 (Figure 1A), respectively, from 2011 to 2014 [35]. The comparative analysis of recent Google earth images and beach profiles surveyed respectively in 2001, 2010 and 2016 [35] also confirms this trend. Particularly it showed that erosion was very intense in the northernmost portion of segment S9, causing a shoreline retreat of approximately 50 m during the period 2001–2016, of which 30 m during the period 2010–2016 alone (profile P2 in [35]). This recent trend to increased erosion is confirmed also by Buccino et al. [36] who individuate southwards of the Campomarino harbour (Figure 1A) one of the major erosion areas, characterized by a LLR erosion rate of −4.0 m/y during period 2004–2016. Furthermore, as highlighted by the study of Aucelli et al. [37], the northernmost portion of segment S9 is characterized by the highest coastal erosion hazard level.

Starting from this situation of very recent, perhaps ongoing erosion and destabilization of the southernmost portion of the Molise coast, which also involved previously stable coastal segments (S8 and S9), to apply our UAV methodological approach we decided to locate the test area in the northernmost sector of segment S9 (Figure 1A). This sector, in fact, has given clear evidence to be in rapid evolution and affected by progressive erosion and shoreline retreat.

The selected test area (Figure 1) is represented by an approximately 200 m long coastal stretch located about 1.5 km south of the Campomarino touristic port. It is characterized by a sandy beach with a variable width of between 6 and 25 m, backed by a dune system, up to 4 m high above sea level. Due to marine erosion, the front slope of the dune has been replaced by an up to 2.5 m high, sub-vertical erosion scarp that separates the dune from the backshore (Figure 1B).

3. Materials and Methods

In this paragraph, we shortly illustrate the materials and methods used for analysing the long-term evolution of the southern Molise coast and the short-term morphological changes of the test area together with its vulnerability to beach erosion and inundation.

According to Crowell et al. [39], we refer to short-term, mid-term and long-term shoreline changes considering periods of <10 years, 10–60 years, and >60 years, respectively.

3.1. Shoreline Changes along the Southern Molise Coast from 1954 to 2019

Shoreline digitization was carried out in ArcGIS environment using aerial photos of 1954, orthophoto maps of 2004, and images taken by Google Earth in 2014, 2016 and 2019 (

Table 3. Data sources used to calculate the shoreline variations in the test area. RMSE = root mean square error.

Date	Data Source	Scale	RMSE (m)
1954	Aerial photo	1:36,000	5
2004	Orthophoto map	1:2500	3
2014	Google Earth image	1:500	1
2016	Google Earth image	1:500	1
2019	Google Earth image	1:500	1

).

As the test area is part of a microtidal environment, the shoreline position was defined as the water line at the time of the photo [40]. Moreover, as it was not possible to reconstruct the tidal conditions for each image, a maximum uncertainty of ±1.6 m was assumed for the daily water line position, considering a daily rise and fall of the sea level of ±20 cm [41] and an intertidal slope of about 12% derived from beach profiles surveyed in 2016.

To evaluate shoreline changes, we used the Digital Shoreline Analysis System (DSAS), a freely available extension to ESRI’s ArcGIS [42]. This tool allows to create automatically regularly spaced transects, orthogonally to the shoreline, along which shoreline positions are compared. In this way, it provides a distance parameter (Net Shoreline Movement—NSM) that measures for each transect the distance between the oldest and the youngest shoreline considered. It also allows getting the Linear Regression Rate (LRR), which represents the average rate of accretion (if positive) or erosion (if negative), obtained for various periods after fitting a least-square straight line to each shoreline section [42]. In this study, shoreline variations were determined for the southern Molise coast using 1451 transects placed at an equidistance of 10 m.

3.2. Short Term Changes of the Beach System in the Test Area

A UAV system was used to assess the modifications of the beach system in the short period. To date, two UAV survey campaigns have been carried out, respectively on 19 June 2019 and 20 June 2020, by using a Phantom 3 Standard (Quadcopter), which is a drone developed by Da-Jiang Innovations (DJI, China). This drone mounts a stabilized camera that compensates involuntary movements of the UAV due to wind, ensuring in this way the correct orientation of the photos with respect to the ground. Coordinates for each photo are registered by an internal GPS allowing the drone to follow the points (waypoints) previously fixed with the GPS in the flight plan. An average flight altitude of 40 m was used for both surveys allowing obtaining a ground resolution of 1.5 cm/px. Using the default camera of the Phantom 3 Standard, during each flight about 250 images with a longitudinal overlap (flight direction) of 85% and a flight strip overlap of 60% were taken.

For both flights, six targets, 40 × 40 cm in size and easily visible from above, were placed to record the position of the Ground Control Points (GCPs) and consequently orient the model in space. A Trimble R6 GNSS receiver in static nRTK mode (positioning performance horizontal 8 mm, vertical 15 mm) was used to acquire the GCPs.

To evaluate plano-altimetric variations of the study beach and related volumetric changes, 3D models were generated using Agisoft Metashape Software. To obtain these models, the (2D) photos acquired with the drone were processed using the Structure from Motion (SfM) algorithm, which internally implements the photogrammetry and computer vision methods. A functional correlation between 3D object points and 2D image points via collinearity condition is the core concept of photogrammetric image data processing [43,44]. Examination of two photos and related orientation parameters allowed identification of the common points and determining the related 3D coordinates.

To obtain the final outputs (Dense Cloud Points, DEM, Orthophoto and Volume change), we followed the work-steps illustrated in Figure 2: (1) Generation of the flight plan; (2) Fieldwork consisting in the positioning of the targets in the field, measurement of the position with the GNSS receiver and acquisition of photos with the UAV; (3) Aerial data processing (as described by [45]), error checking of the model against GCPs (

Table 4. Mean error values on GCPs.

GCPs Error	East	North	Altitude
2019	1.6 cm	1.5 cm	1.6 cm
2020	1.5 cm	1.2 cm	1.7 cm

—The errors represent the average displacement of the point cloud, orthophoto and DSM with respect to the position of the points measured with the GNSS), and export of point clouds, DEMs and orthophotos; (4) Dense cloud point processing. To obtain information about volume changes, the point clouds of 2019 and 2020 were imported into the Cloud Compare software. The volume changes were calculated with the 2.5D Volume tool.

Finally, based on DEMs 2019 and 2020 derived from the UAV surveys, for both years ten beach profiles were extracted, extending from the top of the erosion scarp limiting the dune up to the water line.

3.3. Coastal Vulnerability Assessment

The beach vulnerability to erosion and inundation was performed for the area under study using the Coastal Vulnerability Assessment (CVA) approach developed by Di Paola et al. [11], which has been applied in several other Mediterranean [11,12,46] and Atlantic [46,47] coastal contexts. The CVA method allows considering for the vulnerability assessment both the beach retreat due to storm surges and the coastal inundation due to run-up on the beach. This was possible due to the detailed collection of climatic, sedimentological and topographic data (e.g., bathymetry, beach width) performed by the authors along the investigated area since 2001.

In detail, the CVA method evaluates the coastal vulnerability for each considered period according to the following equation:

CVA = I_Ru + I_R + E + T_i

(1)

where I_Ru is the wave run-up height index, given by the run-up level divided by the beach foreshore slope [48], I_R is the short-term erosion index, which provides a measurement of the maximum beach recession due to storms normalized with the beach width [49], E is the beach erosion rate in m/y, and T_i measures the horizontal distance travelled by the tidal range.

In this study, the CVA has been determined by means of Equation (1) leaving out the contribution of the T_i index. In fact, T_i is zero because the study area is microtidal and experiences ordinary tidal excursions of 30–40 cm [35].

The I_Ru index provides the measurement of the potential inundation capacity that characterizes natural beaches with respect to wave storms. According to Stockdon et al. [48], the wave run-up height is provided by Ru_2% that is the wave run-up level exceeded by 2% of the number of incoming waves, measured vertically from the still water line. This value is projected along the beach through the calculation of X_Ru2% that corresponds to the horizontal distance traveled by the wave in the run-up process, and is obtained through the following formula:

$${\mathrm{X}}_{\mathrm{Ru}2\%}={\mathrm{Ru}}_{2\%}/\tan \left({\mathsf{\beta }}_{\mathrm{f}}\right)$$

(2)

where β_f is the foreshore beach slope and I_Ru is evaluated through the 2% exceedance level for run-up peaks (Ru_2%) on natural beaches. The latter is retrieved according to the empirical approach proposed in Stockdon et al. [48]:

$${\mathrm{Ru}}_{2\%}=1.1(0.35{\mathsf{\beta }}_{\mathrm{f}}\sqrt{\left(\mathrm{H}\cdot {\mathrm{L}}_0\right)}+\sqrt{\left[{\mathrm{HL}}_0\left(0.563{\mathsf{\beta }}_{\mathrm{f}}{}^2+0.004\right)\right]}/2)$$

(3)

where β_f is the foreshore beach slope defined over the area of significant swash activity. H can be the significant wave height (H_s) or the average significant wave height of events exceeding 2 m (H_t) depending on the type of the analysis, while L₀ is the offshore wave-length, which can be expressed in terms of the wave period by means of the linear dispersion relationship, L₀ = gT²/2π (in this case, T is the corresponding period of H_s or H_t). Equation (3) also considers the increase in water level due to wave set-up, which constitutes the main part of the increase in mean sea level, so the other terms of wind set-up and inverter barometer are properly neglected.

Therefore, I_Ru takes values that depend on the percentage associated with the maximum horizontal run-up distance of the wave on the beach (X_Ru2%) normalized with respect to the width of the emerged beach (L). In this way, the I_Ru index can be customarily clustered into four discrete levels, as shown in

Table 5. Coastal vulnerability assessment (CVA) classification scheme (after [11]).

Variable	1	2	3	4
IR (%)	≤15	16 ÷ 30	31 ÷ 50	>50
IRu (%)	≤40	41 ÷ 60	61 ÷ 80	>80
E (m/y)	≥−0.5	−0.6 ÷ −1.0	−1.1 ÷ −2.0	<−2.0
	Low	Medium	High	Very high
CVA	≤6	7 ÷ 9	10 ÷ 12	≥13

.

The I_R index provides a measurement of potential beach retreat and is used for the dynamical calculation of the shoreline retreat based on the convolution method of Kriebel and Dean [49]. I_R values depend on the percentage associated to the maximum beach retreat (R_max) normalized with the beach width L. R_max is evaluated as the maximum value of general solutions associated to the Kriebel and Dean [49] convolution method:

$${\mathrm{R}}_{\mathrm{t}}/{\mathrm{R}}_{\infty }=\left\{1-{\mathsf{\gamma }}^2/\left(1+{\mathsf{\gamma }}^2\right)\exp \left(-2\mathsf{\sigma }\mathrm{t}/\mathsf{\gamma }\right)-1/\left(1+{\mathsf{\gamma }}^2\right)\left[\cos \left(2\mathsf{\sigma }\mathrm{t}\right)+\mathsf{\gamma }\mathrm{sen}\left(2\mathsf{\sigma }\mathrm{t}\right)\right]\right\}/2$$

(4)

$${\mathrm{R}}_{\infty }=\mathrm{S}\left[\left({\mathrm{W}}_{\mathrm{b}}-{\mathrm{d}}_{\mathrm{b}}/{\mathrm{m}}_0\right)/\left(\mathrm{B}+{\mathrm{d}}_{\mathrm{b}}-\mathrm{S}/2\right)\right]$$

(5)

where $\mathsf{\gamma }=2\mathsf{\pi }{\mathrm{T}}_{\mathrm{s}}/{\mathrm{T}}_{\mathrm{D}}$, that is the ratio between the time scale of beach erosion T_s and the storm duration T_D.

In Equations (4) and (5), S is the sea level increase due to wave storm, B is the berm height, m₀ is the slope of the seabed up to the closure depth, d_b is the breaking depth of the wave storm and W_b represents the offshore breaking depth distance related to the wave storm, evaluated following the formula ${\mathrm{W}}_{\mathrm{b}}={\mathrm{y}}_0\sqrt[2]{{\left({\mathrm{d}}_{\mathrm{b}}/\mathrm{A}\right)}^3}$, where ${\mathrm{y}}_0=\left(4{\mathrm{A}}^3\right)/27{\mathrm{m}}_0^3$. The parameter A depends on the typology of the sediment forming the analyzed beach foreshore.

Therefore, I_R takes values that depend on the percentage associated with R_ma normalized with respect to the width of L. Therefore, I_R index can be customarily clustered into four discrete levels, as shown in Table 5.

Considering that the evolution of the beach is not only linked to the effects produced by coastal inundation during an extreme event (H_t = 3.5 m, Table 2), but also to events caused by ordinary wave dynamics (H_s = 0.7 m, Table 2), the parameters I_Ru and I_R were calibrated considering both conditions. In particular, a double weight was attributed to the indexes evaluated taking into account the effect of the storm surge.

$${\mathrm{I}}_{\mathrm{Ru}}=\left({\mathrm{I}}_{\mathrm{Ru}\left({\mathrm{H}}_s\right)}+2·{\mathrm{I}}_{\mathrm{Ru}\left({\mathrm{H}}_{\mathrm{t}}\right)}\right)/3$$

(6)

$${\mathrm{I}}_{\mathrm{R}}=\left({\mathrm{I}}_{\mathrm{R}\left({\mathrm{H}}_s\right)}+2·{\mathrm{I}}_{\mathrm{R}\left({\mathrm{H}}_{\mathrm{t}}\right)}\right)/3$$

(7)

Regarding the E index, shoreline changes in the medium-term (periods 1954–2004—E₁ and 2004–2016—E₂) and in the short term (periods 2016–2019—E₂₀₁₉ and 2016–2020—E₂₀₂₀) were considered. To this aim, we used the indications proposed by Crowell et al. [39], attributing to the first period a single weight, to period 2004–2016 a double weight, and to the short-term periods 2016–2019 and 2016–2020 a triple weight.

$$\mathrm{E}=\left({\mathrm{E}}_1+2·{\mathrm{E}}_2+3·{\mathrm{E}}_{2019/2020}\right)/6$$

4. Results

4.1. Shoreline Changes along the Southern Molise Coast from 1954 to 2019

The application of the DSAS method allowed the calculation of the overall shoreline changes of the southern Molise coast (segments S6–S9) which occurred during the last 65 years considering periods 1954–2019, 1954–2004, 2004–2016 and 2016–2019.

Regarding the long-term evolution of the southern Molise coast (time interval 1954–2019), results obtained (Figure 3,

Table 6. Long to short-term shoreline changes along the southern Molise coast.

Segment	1954–2019		1954–2004		2004–2016		2016–2019
	NSM (m)	LRR (m/y)	NSM (m)	LRR (m/y)	NSM (m)	LRR (m/y)	NSM (m)	LRR (m/y)
S6	76.41	1.07	49.86	1.00	9.06	0.75	17.49	5.84
S7	−166.78	−2.71	−169.80	−3.40	−14.00	−1.17	1.26	0.42
S8	36.66	0.53	26.66	0.53	3.72	0.31	6.29	2.10
S9	5.76	0.09	9.76	0.20	−5.65	−0.47	1.65	0.55

) confirm that erosion primarily affected S7, which registered an overall average shoreline retreat of −166.78 m (Table 6). Conversely, the other segments showed an overall positive balance (S6 and S8) or substantial stability (S9). Shoreline changes during period 1954–2004 largely confirm this overall trend. Considering instead the subsequent periods 2004–2016 and 2016–2019, several differences become evident. First, the shoreline changes from 2004 to 2016 highlight for S7, although still affected by erosion, a clear trend towards stabilization (−1.17 m/y. Table 6), while S6 and S8 maintain a more or less significant trend to shoreline advance. The positive shoreline changes of S6 and S8 appear clearly at least partially linked to the construction of the San Pietro (Termoli) and Santa Cristina (Campomarino) touristic ports and that of the related defence structures, as well as to the interventions of artificial nourishment. S9, instead, substantially stable in the long-term period, is affected by an overall negative shoreline change from 2004 to 2016 (−0.47 m/y, Table 6).

During the period 2016–2019, shoreline trends are positive for all segments, especially for S6, suggesting, in a first hypothesis, a general trend of amelioration of marine climate conditions, and allowing a partial recovery of the land lost during the precedent period. Nevertheless, this recovery trend did not extend to the northern part of S9 including the test area where erosion continued (Figure 3) with shoreline retreat rates locally comprised between 2 and 5.5 m/y.

4.2. Shoreline and Beach Morphology Changes in the Test Area from 2019 to 2020

Based on the data acquired during UAV flights 2019 and 2020, orthophotos (Figure 4A,B) and related DEMs (Figure 4C,D) were produced. Comparing the two orthophotos of the test area (Figure 4C,D) gives evidence about a net increase from 2019 to 2020 of the coverage in vegetation debris and other kind of beach litter in the northern and southern portions of the study beach. In the central portion, instead, the beach appears clean, and several vehicle tracks are clearly visible (Figure 4B) suggesting operations of cleaning and smoothing of the beach.

4.2.1. Verification of the Correspondence between UAV and GNSS Altimetric Data 2019 and 2020

The two DEMs were used to verify the goodness of derived z (altimetric) UAV values by comparing them with those acquired by GPS surveys (GNSS sampling points) performed at the same time. For each year, two transects of GNSS sampling points, extending from the dune to the shoreline (A1–A10 and A11–A19 in July 2019, Figure 4C; B1–B11 and B12–B22 in July 2020, Figure 4D), were compared with equivalent points extracted from the DEMs (Figure 5). In particular, the two transects correspond to profiles T1 and T7 described below and present in Figure 6, and are characterized by an equidistance of 1 m between sampling points.

In Figure 5, both types of z-data (zGNSS and zUAV) are represented and compared. In particular, the z coordinates measured by GNSS are reported on the y axis, and the z coordinates in DSMs (derived from UAV data) are reported on the x axis. The result shows an overall very good to nearly total correspondence between zGNSS data and constructed regression lines, with R² = 0.96 and 0.99, respectively. In detail, the calculated average of the differences between the two types of measurements for 2019 provide respectively a value of 0.130 m with a standard deviation of 0.156 m for 2019, and a value of 0.068 m with a standard deviation of 0.091 m for 2020. It confirms the overall very good quality of acquired zUAV data and the possibility to use them for several purposes such as the analysis of morphologic and plano-altimetric changes of the beach, and the estimation of related volumetric changes (see the following sub-sections).

4.2.2. Analysis and Differences of 3D Models 2019 and 2020

The results obtained by using the Cloud Compare’s “2.5D Volume” tool to evaluate volume changes 2019–2020 highlighted: (1) an “Added Volume” of 1713 m³; (2) a “Removed volume” of 2493 m³; (3) a total loss of 780 m³ (2493 m³–1713 m³) within a total beach area of 4400 m².

Figure 6 shows the difference between the Dense Clouds Points 2019/2020 and the location of the ten beach profiles (T1–T10) extracted from the DEMs. This difference provides important data about the areal distribution of topographic changes occurred in the test area from 2019 to 2020, helping to quantify local degrees of aggradation/lowering of the beach surface and possible retreats of the dune front.

Based on the distribution of topographic changes, it is possible to distinguish within the test area a central, substantially stable but most likely artificially modified beach portion (including profiles T4–T6), from a northern and southern portion, including respectively profiles T1–T3 and T7–T10. The latter, in fact, show some morphological changes, which are mostly due to slight beach aggradation (mainly 0.5–1 m) and a consistent topographic lowering close to the dune front (up to 2.5 m) especially in the southern beach portion.

4.2.3. Shoreline and Beach Morphology Changes 2019–2020 along Beach Profiles T1–T10

By examining in detail the beach profiles T1–T10 extracted from DEMs 2019 and 2020, respectively (Figure 7), together with their main morphometric parameters (

Table 7. Main morphometric features of beach profiles T1–T10 surveyed in 2019 and 2020.

Profiles	Backshore Width—L (m)		Backshore Slope—βb (%)		Foreshore Slope—βf (%)		Total Slope—m0 (%)		Berm—B (m)		Dune Front Retreat (m)
	2019	2020	2019	2020	2019	2020	2019	2020	2019	2020	-
T1	8.92	11.91	9.8	12.8	8.6	12.2	9.3	12.6	0.60	0.85	−0.31
T2	8.44	11.64	10.7	11.1	12.2	12.0	11.3	11.3	0.61	0.54	−0.29
T3	7.25	9.88	9.4	4.4	10.9	14.6	10.0	8.5	0.54	0.99	−0.79
T4	22.33	22.41	5.7	4.6	11.2	13.6	6.4	6.4	0.36	0.77	−0.04
T5	26.14	25.33	4.2	4.5	10.0	17.8	5.1	6.2	0.48	0.67	−0.01
T6	11.01	9.66	6.4	7.0	15.3	18.3	8.4	10.4	0.49	0.75	−0.14
T7	7.20	6.00	8.9	9.8	14.9	28.3	10.6	17.6	0.45	1.22	−0.81
T8	6.71	7.79	15.8	19.6	7.1	20.1	11.9	19.7	0.38	0.65	−0.85
T9	4.89	8.98	15.8	13.6	3.5	11.6	9.5	12.9	0.18	0.56	−1.42
T10	3.99	8.22	21.0	12.8	5.2	25.5	13.1	15.7	0.21	0.62	−2.59

), several aspects can be highlighted.

First, it is possible to note that the test area is substantially stable; positive and negative shoreline changes 2019–2020 are limited within a few meters (Figure 7) and reach maximum values of 4–5 m along the southernmost profiles T9 and T10 (Figure 7).

The front of the dune is substantially stable both in the northern and central stretches (−0.29 m ÷ −0.79 m and −0.01 m ÷ −0.14 m, Table 7), conversely is affected by retreat in the southern stretch (−0.81 m ÷ −2.59 m). Consequently, the backshore width remains rather stable along all profiles, with some major enlargement along T1, T9 and T10.

Morphological changes also find evidence in the variations of beach slopes. Backshore slopes variably increase or decrease from 2019 to 2020 but are relatively stable along T2 and in the central beach stretch (T4–T6, Table 7). Foreshore slopes, instead, with the only exception of T2, increase to various degrees, and in a more consistent way in the southern stretch (see especially T7, T8 and T10, Table 7). Consequently, total slopes also change significantly and show major increases in the southern stretch.

These data, together with those referring to the amount of dune front retreat (Table 7) allow for identifying the southern stretch as the most dynamic one, with significant retreat rates of the dune front whose sediments contribute to beach aggradation (all profiles) and progradation (T9 and T10).

The central stretch is largely affected by shoreline retreat (profiles T5 and T6), but the elevated width of the backshore together with its artificial management and the probably scarce wave attack of the dune front, have as a result its overall morphological stability.

Finally, the northern stretch appears to be in an intermediate position, with slight to nil dune front retreat and an overall positive shoreline balance (stable to slight progradation).

As regards its sedimentological features, the foreshore of the investigated coastal stretch is characterized by the presence of medium to coarse sands, with a mean size value (D50) between 0.39 and 0.78 mm

4.3. Long to Short-Term Shoreline Changes in the Test Area and Related Erosion Indexes

The long to short-term shoreline changes calculated along profiles T1–T10 for the overall time interval 1954–2020 and periods 1954–2016, 2004–2016, 2016–2019 and 2016–2020, allow the following data to be highlighted.

In the long-term (period 1954–2016,

Table 8. Shoreline changes in m (NSM) and annual shoreline change rates in m/y (LRR) calculated along profiles T1–T10 for periods 1954–2016 and 2004–2016, and evaluation of the related erosion indexes (E₁ and E₂).

Transects	1954–2016			2004–2016
	NSM (m)	LRR (m/y)	E1	NSM (m)	LRR (m/y)	E2
T1	11.42	0.24	1	−7.86	−0.65	2
T2	11.58	0.23	1	−7.67	−0.64	2
T3	14.91	0.29	1	−6.64	−0.55	2
T4	18.45	0.37	1	−5.49	−0.46	1
T5	20.46	0.41	1	−4.95	−0.41	1
T6	1.80	0.22	1	−31.39	−2.62	4
T7	0.38	0.21	1	−34.85	−2.90	4
T8	1.59	0.24	1	−36.61	−3.05	4
T9	0.12	0.22	1	−37.41	−3.12	4
T10	−2.02	0.20	1	−37.77	−3.15	4

), the test area is characterized by a substantial stability (T6–T10) to slight progradation (T1–T5), with erosion indexes (E₁, Table 8) assuming a minimum value of 1.

During the recent sub-period 2004–2016 (Table 8), the test area undergoes a general destabilization and the onset of a consistent trend to erosion is observed, especially in the southern stretch of the test area (T6–T10), where erosion indexes switch from 1 to the maximum value 4. These data are in accordance with the data previously illustrated (Section 4.1, Figure 3) about the negative shoreline trend during period 2004–2016 of the northernmost portion of coastal segment S9.

During the following period 2016–2019 (

Table 9. NSM and LLR values calculated along profiles T1–T10 for periods 2016–2019 and 2016–2020, evaluation of the related erosion indexes 2019 (E₂₀₁₉) and 2020 (E₂₀₂₀), and absolute shoreline variations (Net Shoreline Measurement) from 2019–2020 (NSM 2019–2020).

Transects	2016–2019			2016–2020			2019–2020
	NSM (m)	LRR (m/y)	E2019	NSM (m)	LRR (m/y)	E2020	NSM (m)
T1	−4.55	−1.52	3	−3.68	−0.92	2	0.87
T2	−5.49	−1.83	3	−5.14	−1.28	3	0.36
T3	−6.80	−2.27	4	−6.46	−1.61	3	0.34
T4	−8.27	−2.76	4	−7.65	−1.91	3	0.62
T5	−9.70	−3.24	4	−8.87	−2.22	4	0.84
T6	−7.36	−2.46	4	−9.81	−2.45	4	−2.45
T7	−7.50	−2.50	4	−9.88	−2.47	4	−2.38
T8	−10.12	−3.38	4	−12.02	−3.09	4	−1.90
T9	−14.01	−4.67	4	−12.52	−3.13	4	1.49
T10	−15.10	−5.04	4	−12.63	−3.16	4	2.47

), shoreline retreat rates increase on overall (with the only exception of profile T7), but more consistently in the central and southern stretches (profiles T4–T5 and T9–T10, Table 9). This acceleration in shoreline retreat causes a further increase of erosion indexes (E₂₀₁₉, Table 9) calculated along profiles T1–T5. Consequently, apart from T1 and T2, all other profiles are characterized by a maximum erosion index value 4 (T3–T10).

Comparing periods 2016–2019 and 2016–2020 (Table 9) highlights a general decrease in annual shoreline retreat rates during the second period, except for T6 and T7 that confirm previous values. Accordingly, along some beach profiles the erosion index level decreased, precisely from 3 to 2 along T1, and from 4 to 3 along T3 and T4.

An example of significant and progressive erosion that characterizes the study area and the entire southern portion of the Molise coast is provided by the beach profiles surveyed from 2001 onwards (Figure 8) in correspondence of T7.

This sequence of profiles highlights a strong shoreline retreat during the last 20 years or so, which exceeded on overall 60 m and occurred at the expense of both the beach (causing approximately a halving of the beach width) and the dune system behind.

4.4. Coastal Vulnerability Index Assessment

For the CVA assessment, besides the beach erosion indexes already examined above (Section 4.3, Figure 8), the wave run-up height index I_RU and the short-term erosion index I_R were assessed for 2019 and 2020 respectively.

I_RU2% 2019 and 2020 levels were obtained calculating run-up values R_u2% and related parameters X_Ru2% and X_Ru2%/L (width of the beach affected by run-up and relative percentage) under both average wave (H_s = 0.7 m) and average stormy wave conditions (H_t = 3.5 m) along beach profiles T1–T10 (

Table 10. Estimated run-up values Ru_2% and related parameters X_Ru2% and X_Ru2%/L under normal wave (H_s = 0.7 m) and mean storm wave conditions (H_t = 3.5 m), and resulting wave run-up height indexes (I_RU2%) for 2019 and 2020 along profiles T1–T10.

Transects	2019						2020
	Hs = 0.7 m			Ht = 3.5 m			Hs = 0.7 m			Ht = 3.5 m
	XRu2% (m)	XRu2%/L (%)	IRu2%	XRu2% (m)	XRu2%/L (%)	IRu2%	XRu2% (m)	XRu2%/L (%)	IRu2%	XRu2% (m)	XRu2%/L (%)	IRu2%
T1	3.51	39.4	1	14.82	166.1	4	3.23	27.1	1	13.61	114.3	4
T2	3.23	38.2	1	13.61	161.3	4	3.24	27.8	1	13.65	117.3	4
T3	3.30	45.6	2	13.93	192.2	4	3.13	31.7	1	13.20	133.6	4
T4	3.28	14.7	1	13.85	62.0	3	3.17	14.1	1	13.35	59.6	2
T5	3.37	12.9	1	14.22	54.4	2	3.05	12.0	1	12.85	50.7	2
T6	3.11	28.2	1	13.11	119.0	4	3.04	31.4	1	12.81	132.6	4
T7	3.12	43.3	2	13.16	182.7	4	2.90	48.4	2	12.25	204.1	4
T8	3.75	55.8	2	15.80	235.4	4	3.00	38.6	1	12.67	162.7	4
T9	5.34	109.3	4	22.53	460.8	4	3.26	36.3	1	13.75	153.1	4
T10	4.28	107.3	4	18.05	452.4	4	2.93	35.7	1	12.37	150.5	4

).

The values obtained for the horizontal distance traveled by the wave in the run-up process (X_Ru2%) are quite small for profiles T1 to T8 in 2019, and remain quite similar in 2020 for H₀ = 0.7 m. This condition changes completely when considering a H0 of 3.5 m, as the beaches are almost completely temporarily submerged by the sea advance, except for profiles T4 and T5. It is interesting to note the behaviour of profiles T9 and T10 that changes considerable from 2019 to 2020, both for H₀ = 0.7 m and for H₀ = 3.5 m: this change must be related to the morphometric changes that the profiles have undergone between 2019 and 2020.

Summarizing, results highlight some improvement of run-up height conditions from 2019 to 2020 in the southern portion of the test area (T9–T10) with I_Ru2% index levels that pass from 4 to 3 (Table 10). This result can be essentially referred to lower values obtained for the parameter X_Ru2%/L, which expresses the relationship between run-up height and width of the beach affected, i.e., the fact that run-up extends up to the inner limit of the backshore attacking the dune (X_Ru2%/L values > 100%) in average wave conditions (H_s = 0.7 m) or not.

Furthermore, short-term erosion indexes I_R2019 and I_R2020 were evaluated based on R and R/L values obtained for both average wave (H_s = 0.7 m) and average stormy wave conditions (H_t = 3.5 m) along T1–T10 (

Table 11. Evaluation of the short-term erosion indexes I_R for 2019 and 2020, along T1–T10, based on R values calculated for normal wave (H_s = 0.7 m) and average stormy wave conditions (H_t = 3.5 m), respectively.

Transects	2019						2020
	Hs = 0.7 m			Ht = 3.5 m			Hs = 0.7 m			Ht = 3.5 m
	R (m)	(R/L) %	IR	R (m)	(R/L) %	IR	R (m)	(R/L) %	IR	R (m)	(R/L) %	IR
T1	1.30	14.56	1	11.11	124.6	4	2.73	22.89	2	12.43	104.34	4
T2	3.12	36.98	3	13.72	162.5	4	3.51	30.17	3	14.04	120.64	4
T3	2.73	37.61	3	13.35	184.2	4	3.12	31.54	3	12.90	130.61	4
T4	3.59	16.06	2	14.82	66.4	4	3.39	15.12	2	13.52	60.33	4
T5	2.39	9.15	1	13.08	50.0	3	4.86	19.21	2	15.63	61.72	4
T6	4.80	43.64	3	15.94	144.8	4	4.64	47.99	3	15.33	158.71	4
T7	4.89	67.93	4	16.08	223.3	4	4.28	71.30	4	15.34	255.65	4
T8	0.55	8.20	1	10.57	157.5	4	5.42	69.60	4	16.38	210.30	4
T9	12.43	254.11	4	2.38	48.6	3	3.25	36.20	3	13.68	152.30	4
T10	5.52	138.27	4	4.13	103.4	4	6.39	77.78	4	17.68	215.02	4

).

This parameter shows generally high values, since the analysed beach is quite unstable. In particular, even low waves put in crisis several profiles in 2019 (T7, T9 and T10) confirming this trend also for 2020, showing a marked worsening for the T8 profile. Extreme wave attack causes the destabilization of the entire beach and worsening of the conditions along profiles T5 and T9. It is often observed that the R/L value exceeds 100%, with R values exceeding the width of the beach along the analysed profiles (Table 11).

Summarizing, the results obtained highlight a local increase of short-term erosion indexes in the central and southern portions of the test area with I_R values passing from 2 to 3 (T5) and from 3 to 4 (T8 and T9) (Table 11). Such higher short-term erosion levels in 2020 can be attributed to higher I_R levels obtained for average and/or stormy wave conditions.

Comparing the I_RU, I_R and E indexes and related CVA levels obtained respectively for 2019 and 2020 (Figure 9) highlights several variations. In particular, an increase in coastal vulnerability in the central and southern portions of the test-area (T5, T7 and T8), with a significant change along T5 from a medium to high CVA level. Conversely, a slight decrease in coastal vulnerability indexes occurred along T1 and T10, but not such to produce a lowering of the pre-existing CVA levels.

5. Discussion

Beaches can undergo significant changes in the short term and even from one year to the next. Therefore, short-term to annual monitoring of shoreline dynamics and morphometric changes of the beach-dune system appears important to detect in time erosion trends and increases in coastal vulnerability. However, traditional investigation approaches, based on available photogrammetric and/or satellite imagery, may not always adequately respond to the need to update previous data and to ascertain current shoreline conditions and most recent variations. In fact, available data may not have a sufficient resolution or may not cover the periods that need to be investigated. Moreover, annual or even seasonal surveys, carried out in a more traditional way (for example by using GNSS measurements, aerial photos produced ad hoc, etc.), may not provide adequate data and/or be too expensive and time-consuming.

To best overcome these limitations, we developed the methodological approach proposed in this paper. This approach, based on the use of a commercial drone and UAV-derived data, well responds to several of the limitations of major traditional investigation methods, especially to those concerning time consumption and investigation costs, as already amply demonstrated in literature [9,13,14,15,17,24,25,30].

Using the UAV methodological approach, integrated with data derived from more traditional sources (photogrammetric and satellite imagery), allowed us to investigate long to short-term shoreline changes and very recent plano-altimetric changes of the beach-dune system along the southern Molise coast.

In synthesis, the application of this methodology allowed for the following major results:

• Full verification of the correspondence between UAV and GNSS altimetric data acquired respectively in 2019 and 2020 in the test area, confirming the very good quality of acquired UAV z-data and the possibility to use them for high-resolution plano-altimetric beach change analyses.
• Rapid and precise estimation of short-term erosion trends and morphological changes of the beach-dune system in the test area.
• Estimation of volumetric changes from 2019 to 2020 related to beach morphology and shoreline changes, highlighting an overall land loss of about 780 m³ to confirm the persistence of coastal erosion in the test area.
• Verification that changes in shoreline and beach morphology from 2019 to 2020, despite the very short period, caused evident changes of some indexes that enter in the CVA approach for coastal vulnerability assessment. This highlights the importance of coastal monitoring performed at close intervals (at least annually) and carried out over several years, especially in erosion hot spot areas.

6. Conclusions

Results obtained highlight that the use of UAVs can represent a breakthrough in the study and monitoring of morphological changes induced by coastal dynamics. The standardised workflow to obtain DEMs with a high degree of vertical accuracy from UAV data suggest applications in the fields of natural hazards, disaster response, and high-resolution terrain analysis [9,13,14,15,19,20,21,22,23,24]. However, when deriving a DEM from sparse clouds, the effects introduced by the interpolation process need to be considered, and a reliable uncertainty value must be assigned to the elevation dataset.

In conclusion, the relative simplicity of using drones together with the possibility to use UAV-derived data to produce high-resolution DEMs and georeferenced images with a high temporal frequency may contribute to the improvement of methodologies and databases for the mapping of ongoing coastal changes. Moreover, the integrated use of UAV and GIS approaches has proven to be an effective instrument not only for quick spatial data analysis but also for offering an objective approach with consistent measurement and calculation processes. This instrument has furthermore demonstrated to be a valid support for the rapid definition of several parameters that are inserted within a model for coastal vulnerability assessment (in our study the CVA), allowing to keep under control some sensitive indicators of coastal erosion and beach degradation.