Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Citizen science for monitoring seasonal-scale beach erosion and behaviour with aerial drones

## Abstract

Sandy beaches are highly dynamic systems which provide natural protection from the impact of waves to coastal communities. With coastal erosion hazards predicted to increase globally, data to inform decision making on erosion mitigation and adaptation strategies is becoming critical. However, multi-temporal topographic data over wide geographical areas is expensive and time consuming and often requires highly trained professionals. In this study we demonstrate a novel approach combining citizen science with low-cost unmanned aerial vehicles that reliably produces survey-grade morphological data able to model sediment dynamics from event to annual scales. The high-energy wave-dominated coast of south-eastern Australia, in Victoria, is used as a field laboratory to test the reliability of our protocol and develop a set of indices to study multi-scale erosional dynamics. We found that citizen scientists provide unbiased data as accurate as professional researchers. We then observed that open-ocean beaches mobilise three times as much sediment as embayed beaches and distinguished between slowed and accelerated erosional modes. The data was also able to assess the efficiency of sand nourishment for shore protection. Our citizen science protocol provides high quality monitoring capabilities, which although subject to important legislative preconditions, it is applicable in other parts of the world and transferable to other landscape systems where the understanding of sediment dynamics is critical for management of natural or anthropogenic processes.

## Introduction

The coastal zone accommodates 40% of today’s population1,2, with densities much greater than non-coastal land especially in low income countries3. Sandy beaches are extremely dynamic morphological sub-systems of the coastal zone, offering, amongst others, coastal protection and erosion control eco-services4,5. However, global studies have reported that one-quarter of the world’s sandy beaches are eroding at rates exceeding 0.5 m/year6, contributing to a global coastal land loss of 20,000 km2 in the last 35 years7, posing an increasing threat to coastal populations and economies. Beach erosion is not uniformly distributed along the coast, but concentrated on discrete areas6,7 often referred as erosional hotspots8,9. These hotspots are spatiotemporally variable ranging from days to decades and from hundreds to thousands of meters10. As a consequence, erosion mitigation strategies (beach nourishment, rockwalls, sand fencing) tend to be highly localized11. In these areas, the measurement of high-temporal volumetric variations is critical for discerning short-term beach behaviour and recovery from erosional events12. This allows coastal planners to target intervention measures and evaluate their efficiency in protecting backshore assets from current and projected increase of erosion hazard13,14,15.

Nowadays, beach topographic data is obtainable using a variety of ground-based (graded rods16, surveyor-grade global positioning systems or total stations17, terrestrial LiDAR18) or aerial-based (airborne LiDAR19, traditional photogrammetry20,21, unmanned aerial vehicles22,23) approaches. While ground-based surveying methods are labour intensive and of limited spatial coverage, airborne LiDAR or traditional photogrammetric approaches are cost-prohibitive for monitoring purposes, especially in developing countries.

In the last decade unmanned aerial vehicles combined with structure-from-motion multi-view stereo algorithms (UAV-SfM hereinafter) have emerged in environmental research24,25,26 as the best compromise in terms of costs, precision, reproducibility and simplicity, particularly for repetitive beach topographic surveys27,28,29. In coastal areas, multi-rotor UAVs are frequently used30 and their flight time is typically between 20 and 40 min per battery, which corresponds to a flight coverage of 5–30 × 103 m2, depending on flight altitude24. As a consequence, time and budget constraints have so far restricted UAV-SfM surveys to professional researchers or commercial operators at a few representative sites, with limited revisit times27,31,32,33,34,35.

However, recent technological advances in low-cost UAVs and automated flight and positioning solutions, coupled with regulatory changes, provide new opportunities for implementing citizen science to expand the scale of monitoring programs36,37. Citizen science is the process of creating knowledge by engaging non-professional volunteers in scientific research36. Citizen scientists have allowed large-scale scientific experiments to make substantial contributions to science for hundreds of years38. Yet, errors and bias in large and longitudinal citizens science datasets are often poorly understood39. This uncertainty, being difficult to quantify, can compromise data quality or limit the integration of multiple datasets into a single coherent analysis. As a result scepticism exists among professional scientists about the quality of citizen scientists’ data40,41.

In this study, we propose a novel and cost-effective approach to monitor sandy beaches sediment dynamics at a spatiotemporal resolution previously unachievable, through citizen science with low-cost multi-rotors UAVs. We test our method on the high-energy temperate coast of Victoria, Australia, where more than 100 citizen scientists collected 83 aerial datasets in ten previously identified erosional hotspots over a 1.3 year period, approximately every 6 weeks. With this unique dataset, we first evaluate citizen scientists' data quality and bias, then we quantify and compare short-term volumetric and profile dynamics on both open-ocean and embayed beaches. We propose a novel set of indices that capture multi-scale landform dynamics, using topography timeseries alone. This allows us to evaluate beach nourishment efficiency in protecting backshore infrastructure. Lastly, we discuss the potential and limitations that our approach has not only for coastal management but also in other scientific disciplines.

## Results

### Citizen scientist’s data accuracy and bias

Accuracy and bias are objective task-independent metrics of data quality42. For an independent and realistic vertical accuracy assessment of the digital surface models (DSM), the error metrics should be calculated with independent checkpoints that haven’t been used during the digital photogrammetric procedure (Supplementary Method “Citizen Scientists, UAV surveys and photogrammetric details”) and are distributed across the landscape, along representative transects or landform elements27,43. Two benchmark surveys were performed in Warrnambool to evaluate the DSM vertical accuracy resulting from the citizen science protocol under operational (2018-11-29) and worst-case (2019-12-11) scenarios (Supplementary Method “Independent checkpoint surveys”). The checkpoints in both surveys show a very good linear match of modelled elevation (both R2 above 0.99) with slightly larger deviations observed at higher elevations (Fig. 1a).

The Q–Q plots and statistical tests indicate non-normal distribution of absolute errors (Δh, Fig. 1b), which confirms that the normalised median absolute deviation (nmad) is the most appropriate robust estimator for vertical accuracy assessment of photogrammetric datasets44. The nmad values are 0.048 and 0.054 m AHD (Australian height datum) for the operational and worst scenarios respectively. The mean errors indicate that the 2018 checkpoints survey slightly underestimated the height values (− 0.044 m AHD), the 2019 checkpoints survey overestimated (0.128 m AHD) it, which is observable from the two error distributions (Fig. 1c). Both surveys show good precisions (standard deviation) of 0.077 m AHD and 0.063 m AHD respectively. Supplementary Table S3 also reports root mean squared error (rmse) and mean absolute error for comparison purposes. Casella et al. 45 recently demonstrated that approximately 0.05 m rmse is consistently found under different surveying set-ups (varying cameras and flight altitudes), which corroborates a median rmse of 0.059 m found in the relevant literature (Supplementary Table S2). The authors state this systematic error could be due to the vertical sinking of the surveying pole on various sand types, therefore, 0.05 m rmse is a typical error in sandy beaches UAV-SfM studies.

As our operational scenario rmse (0.089 m AHD) is about 0.03 m higher than the literature median, we further explored systematic errors by mapping the spatial variability of checkpoints errors and smart ground control points (GCPs) (Supplementary Figs. S2, S3), finding that errors were generally higher with elevation, especially within foredune vegetation. We mitigated this systematic error by removing vegetation and applying specific limits of detection thresholds, which is a form of split data test46 used to obtain the expected DSMs vertical errors by computing the elevation difference between pre- and post- survey pairs over known stable areas (i.e. calibration areas22).

Bias in citizen science projects can be introduced by allowing individuals flexibility in how, when and where to collect data47. We reduced the risk of bias prior UAV operations by (1) assigning groups to fixed locations, (2) targeting low-tide for survey, (3) providing standardised protocols and training on simple to use highly automated equipment and (4) assisting citizen scientists at each sites for the first three flights. The principle issue is the location of the portable GCPs which can be affected by citizen scientists’ spatial effort bias31, leading to positional errors due to GPS signal blockage from foredunes, tall trees or buildings. As GCPs positional errors directly impact the point cloud georeferencing process (hence, the resultant DSM), we assessed whether significant (p = 0.05) differences exist in GCPs spatial dispersion (i.e. their two-dimensional spread across the surveyed area), their positional X, Y and Z variances and the images georeferencing accuracies (X, Y and Z rmse) reported by the photogrammetric software Pix4Dmapper (V4.3.31), between locations (Fig. 2).

The Kruskal–Wallis H test49 indicated that there are statistically significant differences (p = 0.05) in all positional GCPs variances (Fig. 2a) and in image georeferencing errors (Fig. 2b), between the different locations. This can be observed in the pairwise comparison heatmaps (Fig. 2c–e), where in the great majority of cases, statistically significant differences in X,Y and Z variances are found between any one location and at least four others. This is likely to be the cause of all locations having statistically significant differences in rms georeferencing errors compared to each other (Fig. 2f). Importantly, GCPs elevation precisions were outside acceptable variance (20 mm) 7% of the time (61 out 886), with only 9 cases (1%) exceeding 64 mm (outliers). Therefore, positional precision differences between locations are largely attributed to random errors during the GCPs position recording process and no citizen scientist’s bias can be ascertained. The same applies to the images georeferencing errors. Their location-specific medians range from a minimum of 2 mm (Seaspray and Cowes) to a maximum of 15 mm (Warrnambool), which are of the same order of magnitude of the GCP positional precisions. As such, we cannot exclude the possibility that these differences are also due to random errors in the geolocation recording of the GCPs.

Regarding the GCPs two-dimensional spatial spreading (Supplementary Fig. S6), the Kruskal–Wallis H test showed that there are no statistically significant differences between locations (H = 11.12; p = 0.28, n = 78).

Additionally, to assess whether there are statistically significant differences between the end-products of our protocol across locations (i.e., the expected DSM error) and considering that more independent checkpoint surveys were unavailable, we analysed limits of detection thresholds as indicators of ‘overall’ data quality (Supplementary Method “Limit of detection analysis”). An example of limit of detection threshold derivation and error normality evaluation in Apollo Bay is shown in Supplementary Fig. S8. The Kruskal–Wallis H test found that no significant differences (H = 16.167, p = 0.063, n = 73) exist between each location specific threshold distributions at the 0.05 significance level.

### State-level volumetric monitoring of open-ocean and embayed beaches

Overall, 100,794 ± 243 m3 of sand has been transferred off the beachfaces (net erosion) during the monitored period (from the first June 2018 to the 29th August 2019) across the 10 locations (Fig. 3a, see Additional Method “Area of study” for more information). Swell-exposed open-ocean beaches displayed mean elevation changes (MEC) of greater amplitudes than embayed beaches, which are situated along fetch-limited or sheltered coastlines. All the following observations are relative to the whole aforementioned monitoring period.

On average, the absolute MEC for open-ocean beaches (0.11 ±  < 0.01 m) is almost three times larger than for embayed ones (0.04 ±  < 0.01 m). During erosional phases (negative MEC in Fig. 3b), the average MEC is − 0.13 ±  < 0.01 and − 0.04 ±  < 0.01 m in open-ocean and embayed beaches respectively. Similarly, during recovery phases (positive MEC in Fig. 3b), the average MEC is + 0.09 ±  < 0.01 and + 0.04 ±  < 0.01 m on open-ocean and embayed beaches respectively. Post-erosion recovery times are highly variable, ranging from a few days (Warrnambool) to approximately one year (Portarlington), while some locations never fully return to their initial state (Apollo Bay, Inverloch).

To model the monitored beachface dynamics and distinguish between erosional and depositional behavioural regimes, we used the residual beachface cluster dynamics index (r-BCD, Method, Fig. 3a). Therefore, a location r-BCD index is one single signed value that characterises the dominant behaviour that the beach in that location exhibited during the entire monitored period.

The open-ocean beaches of Port Fairy, Warrnambool, Marengo and Seaspray exhibit highly to slightly depositional regimes. Interestingly, Port Fairy and Seaspray r-BCDs indicated that both locations had a depositional behavioural regime during the monitoring period despite these locations experiencing net sediment losses of − 10,962 ± 66 m3 and − 7,517 ± 726 m3 respectively. In other words, their intra-annual cut and fill dynamics (Fig. 3b) are skewed toward net sediment gains. By contrast, Inverloch, Apollo Bay and Point Roadknight display highly, moderately and slightly erosional behavioural regimes which are corroborated by end-of-monitoring sediment losses of − 66,700 ± 909 m3, − 13,871 ± 2,286 m3 and − 2,159 ± 562 m3 respectively.

The embayed beaches of St. Leonards, Portarlington and Cowes are protected from the predominantly south-westerly swell being located within Port Phillip and Western Port Bays. Despite Portarlington never fully recovering from its first erosional event in the timeseries, its depositional behavioural regime indicates that it is more likely to accrete rather than erode. Conversely, Cowes displays a highly erosional behavioural regime despite never showing a negative sediment budget respective to the beginning of the monitoring.

### Sediment dynamics highlight accelerated or slowed-by-intervention erosion in critical locations

Empirical beachface cluster dynamics indices (e-BCDs, Fig. 3c) measure the importance (absolute score) and main magnitude trend (score sign) of beachface erosional, recovery, depositional and vulnerability behaviours. We use e-BCDs to detect naturally accelerating and slowed-by-intervention erosional modalities in Inverloch and Apollo Bay respectively. The e-BCD indices are representative of the whole monitoring period and are intended to depict the main trend of beachface sediment dynamics that occurred during the observation time.

From the 22th August 2018 to the 30th July 2019, the monitored area in Inverloch displayed an erosional score (+ 1.37) that is greater than its high recovery score (− 0.93), indicating that erosional hotspots mostly continued to erode the beachface rather than recover. The combination of positive erosional and negative recovery indices imply that through time, an increased amount of sediment has been eroded, while recovering areas tended to accumulate less sediment than was previously lost. The depositional score (− 0.60) is substantially lower than the vulnerability (+ 2.28), indicating that depositional areas became erosional, rather than continuing to accrete. These signs suggest that depositional hotspots gained less sediment through time and, once starting to erode, the volumes lost were typically higher than what had previously been deposited.

These dynamics indicate an accelerated beachface depletion, which is confirmed by Inverloch mean elevation change time series (Fig. 3b). An accelerating erosional phase started in mid-October 2018 and lasted until the end of the monitoring, totalling a net loss of 48 ± 0.06 m3/m, at a rate of 0.14 m3/m eroded daily, beach wide.

In Apollo Bay from the first June 2018 to the 25th July 2019, the beachface dynamics have been impacted by a weekly sand nourishment program that deposited 16,050 m3 of sand in the intertidal and foredune areas (focussing from 800 m north, zone D in Fig. 4), from the 19th June to the 4th September 2018 and from the 22nd May to the 21st June 2019.

The e-BCDs analysis in Apollo Bay indicates that erosional hotspots tend to turn depositional by a slight margin (− 1.10 erosional and − 1.14 recovery scores), in which case a lower amount of sediment is deposited compared to what had been previously lost (negative depositional sign). On the other hand, depositional hotspots tend to turn erosional by a greater margin (− 1.85 vulnerability and + 0.73 depositional scores), in which case the sediment lost is usually less than what was previously deposited. Overall, Apollo Bay e-BCD indices indicate a slight erosional predisposition (driven by the vulnerability index), precariously reduced by “sediment-sparing” mechanisms (negative signs of erosional and vulnerability e-BCDs) that slows the beachface sediment loss. In fact, despite being supported by sand nourishment, Apollo Bay sediment volumes never returned to its initial state at the commencement of the surveys. This indicates that the management intervention in Apollo Bay is likely to be slowing beachface erosion during the monitoring period, hence, we observed a slowed-by-intervention erosional mode.

### Site level case study: assessing beach nourishment efficiency in Apollo Bay

To demonstrate the scalability of citizen scientists’ data, we down-scaled behavioural regime, morphological and volumetric analysis from the location to the single transect level, assessing the efficiency of a beach nourishment project in protecting backshore economical assets.

The Great Ocean Road is a scenic coastal drive that contributes up to 6.1% of revenue to the regional economy51. It attracts more than 5,000,000 visitors annually, and Apollo Bay is the second most popular visitor destination51. The spatial distribution of transect-specific r-BCDs indicates three distinctive zones where behavioural regimes are mostly mixed (M), erosional (E) or depositional (D) (Fig. 4).

Zone M represents a mixed-zone extending from immediately in the lee of the harbour to 180 m alongshore. Generally, erosion during the monitoring period occurred on the lower intertidal beach and occasionally on the incipient dune, leaving the vegetated foredune intact. Swash-generated berms (Fig. 5a) occur throughout zone M and are typical landforms of low tide terrace type beaches52.

Zone E is predominantly erosional, despite this area being reported to have experienced accretion since 1956 due to the construction of the harbour and its shadowing effect from incident waves53. Here, the erosion in the intertidal zone is mostly uninterrupted throughout the surveys, leading to a gradual transgression of the erosional scarps, eventually causing up to 6 m of incipient dune recession (Fig. 5b). The most erosional transects are located 500 m alongshore, within 40 m from the main pedestrian beach access (Fig. 4 inset map c). Further to the north, the magnitude of beachface lowering gradually diminishes to finally become net accretion accompanied with foredune recession at 800–840 m alongshore, entering zone D.

Zone D has been reported to be receding since mid-198053 and is also where the majority of nourishment has occurred during the monitoring. Notwithstanding, zone D comprises two “erosional enclaves” (Fig. 4 inset maps a,b) which signal local erosional hotspots of major concern. Moreover, despite an informal rockwall effectively protects the foredune between 1,100 and 1,200 m alongshore, foredune recession has been widespread between June and July 2018 causing sections of a footpath to collapse (1,500 m alongshore, Fig. 5c).

The different behaviours in zones M, E and D can be observed by adopting a beach-wide perspective, while considering the volumes of sand introduced by the nourishment program. From the first of June 2018 to the 21th June 2018 (Supplementary Data ‘site_level_change”), particularly stormy weather eroded 9,055 ± 76 m3 of sand from the beachface, slightly alleviated by post-storm sand renourishment of 1,320 m3. From the 21th June to 26th July 2018 (Fig. 6a) the beachface recovery was assisted by 4,875 m3 of renourished sand, taking the alongshore volumetric change in this period only slightly below zero (− 0.20 ± 0.09 m3/m). The majority of deposition took place in zone D where renourishment focussed, while zone E kept losing sediment, especially from the intertidal beach (Fig. 5b). This pattern is accentuated in the next period (from 26th July 2018 to 24th September 2018, Fig. 6b) when 9,165 m3 of sand was supplied to the beach by managers mostly northward from the 800 m mark (Fig. 4).

During these renourishment periods, zones E and D behaviours are noticeably inverse to each other, whereas zone M generally follows the behaviour of the adjacent zone E. When renourishment is suspended, zones E and D are essentially undifferentiated, whereas zone M is remarkably dissimilar to the rest of the beach, especially during the spring to summer transition (from the 24th September 2018 to 11th December 2018, Fig. 6c,d).

Overall, after initial recession of the foredune, the intertidal zone showed a net accretion by the end of the monitoring, indicating that the depositional behaviour of zone D is precariously tied to the nourishment operations. In fact, during nourishment periods, the alongshore distribution of volumetric change clearly sets zone D apart, being it the only accreting zone amongst overall beach erosion (Fig. 6a,b). Conversely, when nourishment is suspended, zones E and D erode or accrete conjointly, while zone M shows a noticeably divergent behaviour (Fig. 6c,d).

Therefore, in the absence of nourishment, we would expect the whole monitored area in Apollo Bay to be erosional, corroborating our previous argument about the slowed-by-intervention erosional mode of Apollo Bay. We can conclude that the current sand nourishment strategy is likely to be locally decelerating a recent beach wide erosive trend, but it is not sufficient to maintain or fully recover subaerial beach sediment budget in the long-term. Additionally, mostly depositional areas seem to be accreting the intertidal beach, while the foredune is still receding.

## Discussion

The vertical accuracy of the citizen scientists’ protocol (rmse = 0.089 m) obtained in an operational scenario is 0.03 m less accurate than the median accuracy reported across the relevant UAV-SfM beach monitoring literature (median rmse = 0.059 m). However, using robust statistical estimators which put less weights on error outliers (focussed around foredune vegetation) led to an error estimation of 0.048 m, slightly better than what others obtained in similar works27,54 (Supplementary Discussion “Citizen scientists’ accuracy”). Our bias analysis indicates that citizen scientists’ ground control points (GCPs) spatial arrangement was consistent through time and space. GCPs positional precision and image georeferencing errors have significant inter-groups differences which cannot be attributed to citizen scientists as errors are below the reported GCPs precision. As no significant differences of limit of detection distributions between locations have been found, we can assert that limited bias impacted the quality of data during the full data creation pipeline, with our protocol. Therefore, as citizen scientists' data is of comparable accuracy and bias to professionally acquired UAV-SfM datasets, it is considered of high-quality36.

Besides monitoring beachfaces volumetric changes, our approach temporal resolution allowed us to conceive the residual beachface cluster dynamics index (r-BCD) to model hotspots behavioural regimes in the absence of precise wave data, but based on topography alone. The fundamental assumption of r-BCDs is tied to the classic morphodynamic feedback loops of process and morphology55. As r-BCDs are derived from high frequency elevation changes, they incorporate the morphological variability that landforms experienced due to the locally active geomorphic processes. Citizen scientists allow r-BCDs to be used over multiple locations at site and transect scales, providing government authorities a rapid and powerful beach dynamics assessment that managers can use to prioritise erosion hazard mitigation expenditure. This approach becomes especially favourable in situations where high quality nearshore waves and currents data is lacking and can complement other cost-effective remote sensing approaches, such as satellite-based coastal monitoring.

Multispectral imagery acquired by optical satellites, such as Sentinel-2 or Landsat, allows the use of two-dimensional shoreline landward/seaward shifts as erosion/deposition proxy for large scale erosion monitoring. However, shorelines derived from space are difficult to validate in-situ, consequently, only a few studies used ancillary beach topographic data or coincident shoreline GPS surveys to test the accuracy of the extracted shorelines56,57,58,59. Moreover, sandy beaches topographic features such as beachface slope60, intertidal extent61 and elevation62,63 at continental-scale have also started to be derived from space-based observations. These methods need beach topographic data at higher resolution to be validated against. Our study demonstrates that citizen science UAV-SfM is in an advantageous position to provide the high-quality and reliable 3D data that is essential for ground truthing space-based observations.

Alternatively, the closest comparable topographic monitoring methods are repetitive airborne light detection and ranging (LiDAR) surveys or fixed coastal imaging stations (ARGUS). Coastal LiDAR is often used for quantifying wide-scale shoreline spatial variability of single-event post-storm erosion19,64,65, however, high operational costs have traditionally limited its temporal resolution within geoscientific research66, often missing the fine seasonal dynamics that our approach provides. ARGUS stations on the other hand can typically provide hourly shoreline images of key sites, which are then processed to obtain shoreline elevations and later used to monitor erosion/accretion patterns in a relatively cost effective way. However, ARGUS-derived shoreline elevations are computed with a variety of empirical, semi-empirical and complete numerical models, depending on the quality and amount of site-specific field data (topography, tide and offshore wave measurements) used to calibrate the shorelines. As a consequence, vertical accuracies (in the range of 0.1 to 0.4 m) and methodologies are often site-specific17, possibly limiting the integration of multi-site observations into one integrative study. With our approach, 10 UAVs and 100 smart ground control points were used to repeatedly survey (7 to 10 times) 10 key locations at a total cost of US$250,590 (including equipment, online data processing and hosting and the time to train and support a group of 3–4 citizen scientists per location). This equals to US$3,020 per survey (n = 83), which is expected to be approximately US\$1,290 per survey in the next years, as groups become established and less training and support hours are budgeted (Supplementary Discussion “Monetary costs''). Therefore, we consider our approach to be cost-effective as high-quality, high temporal and very high spatial resolution beach topographic monitoring is consistently achieved at 10 high-priority locations (site scale), which are representative of the Victorian coast (regional scale).

The cost-effectiveness and geographical coverage offered by our approach has important implications in contexts of erosion mitigation in vulnerable communities, such as those in small island developing states (SIDS). In SIDS, a lack of quality coastal data at relevant spatiotemporal scale has been recognised67,68 and their economic dependence to tourism and coastal eco-services within a changing climate context put their subsistence at risk69. Our approach would allow SIDS communities to reliably monitor coastal erosion in sensitive areas cost-effectively, producing continuous survey-grade topographic data over wide and often dispersed geographic areas, with some limitations.

Without environmental context information (especially wave climate), our methods reliably depict seasonal beach change trends but cannot explain the geomorphic forcings responsible for those changes. In fact, the magnitude of changes could differ depending on how far in time from major storms (not investigated here) the surveys are. Although we demonstrated the reliability and usefulness of our approach, coastal processes data is needed to formulate more causative sand dynamics and geomorphological interpretations, which would allow more accurate sand dynamics to be evaluated. This is the next logical step and focus of future studies.

Additionally, as opposed to other coastal citizen science projects where anyone can take part as a citizen scientist70,71, in our approach, legislation plays a central role in defining who can be a pilot, directly challenging its feasibility and reproducibility around the world (Supplementary Discussion “Role of legislation”). In general, the key legislation requirements that need to be met in order to perform and replicate our protocol are that (1) UAV flight are permitted over interest areas, (2) automatic (waypoint-based) UAV flight mode is allowed and (3) no mandatory UAV flight license is required for research applications using sub-2 kg airframes.

In Australia, scientific operations with UAVs below 2 kg of weight are part of the “sub-2 kg excluded category”, which allows individuals older than 16 years of age to fly small UAVs (1) within standard operating conditions, (2) without the need of a remote pilot licence and (3) without a mandatory risk assessment to be approved. Elsewhere, various UAV laws can challenge, impede or allow UAV citizen science applications. In SIDS, UAV regulations range from total ban (Cuba and Barbados) to “Australia-like” approach (Papua New Guinea), with important variants that limit the feasibility (heavy air-traffic and multiple authorities permit systems in Maldives) or replicability (automatic flight procedures not allowed in the Dominican Republic) of our protocol. On the other hand, from the 31st December 2020 the European Union Aviation Safety Agency (EASA) UAVs regulation will allow UAV citizen science projects to take place in 32 European countries (25 of which are coastal), thanks to a regulation very similar to the one adopted in Australia. Despite the combination of a greatest number of coastal no-fly zones (due to controlled aerospace from the military, proximity to aerodromes or natural protected areas) with crowded beaches can discourage UAV citizen science coastal applications in EASA countries, our protocol can also be applied to other environments or scientific disciplines where accurate cost-effective topographic monitoring is needed. In fact, multi-temporal UAV-SfM has already been used by professional researchers around the world for monitoring erosion in mudflats72, badlands73, agricultural watersheds74, rivers75 and open-pit mines76. Additionally, UAV-SfM topographic data has also been used for non-erosion purposes to monitor both natural processes, such as landslide dynamics77, sediment retention dams filling78, crops growth variability79, forest trees growth80, snow depth81 and glaciers melting dynamics82, or anthropogenic processes, such as landfills growth rates83, environmental contamination84 or hiking trails conditions85.

In conclusion, our results not only demonstrated the value of citizen scientist’s high-quality and unbiased data for multi-scale sandy beaches sediment dynamics monitoring, but can also encourage further application of citizen science with drones into all the aforementioned scientific applications. This not only would greatly expand the spatiotemporal scale of scientific experiments, but also democratise scientific engagement access, enhance global environmental awareness and transform citizen scientists into key stakeholders within an adaptive environmental management system.

## Methods

In this work, we used machine learning and limit of detection analysis to detect and reliably quantify subaerial beachface dynamics due to sand-only and above detection limits changes. Spatial autocorrelation analysis is also employed at the site-level to remove spatial outliers and detect statistically significant (p = 0.05) hotspot of changes to better characterise sediment dynamics through time. This information is then explored with novel indices (empirical and residual beachface cluster dynamics indices, e-BCD and r-BCD respectively) based on discrete Markov chain models, at the site and transect scales.

### Virtual network of elevation profiles

A virtual network of digital transects was created for every site and kept fixed during the analysis. Transects are uniformly distributed alongshore with a spacing of 20 m, normal to the shoreline, with an across-shore length of 80–150 m, depending on the beach width discernible from the earliest orthophoto available (Supplementary Fig. S1). We extracted synchronous elevation and colour information along each transect, with a 0.1 m sampling step, in all 10 sites and surveying dates (8 to 10 surveys per site), resulting in a total of 6,809,610 points on 999 transects.

Swash zones were excluded from the analysis. We only retained points which were classified as sand and within the subaerial beachface, which is defined as the area from the upper swash to 2–3 m landward of the vegetation line or where anthropogenic barriers are present.

### Machine learning sand classification

We used machine learning to restrict the analyses to those extracted points that are sand, removing the ones representing beach wrack or coastal vegetation that would otherwise skew our volumetric and behavioural computations. For each survey, we performed the Silhouette Analysis86 to find a sub-optimal number of clusters (k) to partition the points with KMeans clustering algorithm87, using spectral (red, green, and blue bands) and topographic (slope, curvature and distance from the transect seaward origin) features. By iteratively run KMeans (parameters: initial cluster centres selected with KMeans ++, 300 iterations per run, inertial tolerance of 0.001, pre-computed distances) increasing k by 1 at every iteration (up to k = 20), we were able to compute the overall silhouette coefficient associated at every k. We chose as sub-optimal k the value after which a greater k would not substantially reduce the overall clustering performance. Once the sub-optimal k has been found for every dataset and KMeans run using it, clusters were displayed in QGIS (version 3.2.3) and visually labelled as sand or no-sand (Supplementary Fig. S1). No-sand points across each transects have been replaced by an interpolated value using a linear model. Minor manual editing was required to correct sand points erroneously assigned to a non-sand cluster. This mainly occurred in very dark shadows cast by tall trees along coastal walking paths or right below near vertical foredune.

### Limit of detection thresholds

The morphological method88 has been used in a variety of environments, including sandy beaches23,29,31,89. It involves the subtraction of two digital surface models (DSMs) of the same location at different times to obtain a DSM of Difference (DoD), which represents surface change over a shared elevation datum for each period. Our transects are subject to the same error estimation techniques used for DoD analysis. Therefore, in order to account for areas of apparent elevation change (Δh) due to inherent DSM uncertainties, we computed limit of detection thresholds for each DoD (n = 78). Changes within these thresholds are considered uncertain and their contribution is expressed as error margins when altimetric or volumetric change is reported.

To obtain limit of detections, we firstly identified invariant objects as close as possible to the beachface in each site (paths, roofs, rockwall, anthropogenic structures), which we used as calibration zones. Then, we manually digitised vector lines across calibration zones and created a checkpoint every 0.1 m along those lines, resulting in 23,660 check points. Finally, we measured the checkpoints Δh in each DoD and obtained the threshold computing the normalised median absolute deviation45 as robust statistical estimator.

### Volumetric computations and mean elevation change

The exclusion of points within the swash extent and variations in UAV survey coverage resulted in irregularities in the number of transects and valid points per survey (see Additional Data). Therefore, we compared subaerial changes across sites and time using the mean elevation change (MEC), as follows:

$$Mean\,elevation\,change \left( {MEC} \right) = \frac{1}{n}\mathop \sum \limits_{z = 0}^{n} (z_{post} - z_{pre} ),$$

where n is the total number of valid elevation points, $$z_{pre}$$ and $$z_{post}$$ are the elevation above Australian height datum (equivalent to mean sea level) values occurring at the same location in both pre and post surveys.

Additionally, when no inter-site comparisons were involved, we approximated the alongshore volumetric change (in m3/m) as:

$$Along. beachface \,change = \mathop \smallint \limits_{{x_{swash} }}^{{x_{limit} }} \left( {z_{post} - z_{pre} } \right)dx,$$

where $$x_{swash}$$ and $$x_{limit}$$ are the upper swash and landward limit respectively. Plus or minus ( ±) error intervals for both MEC and volumetric change estimates represent the uncertainty related to changes within the limit of detection thresholds.

### Hotspot analysis: local Moran’s I

In order to obtain spatially explicit and statistically significant hotspots of erosion or deposition at the site level, the Local Moran-I90 ($$I_{i}$$) statistics with false discovery rate correction was computed for every inter-surveys elevation difference (Δh) points. The $$I_{i}$$ statistics is defined as:

$$Local Moran^{\prime}s Index \left( { I_{i} } \right) = \frac{{z_{i} - \mu_{z} }}{{\sigma^{2} }}\mathop \sum \limits_{j = 1, j \ne i }^{n} \left[ {w_{ij} \left( {z_{j} - \mu_{z} } \right)} \right],$$

where $$z_{i}$$ is the value of the variable at location i with $$\mu_{z}$$ and $$\sigma^{2}$$ the respective mean and variance, as calculated on the $$n$$ number of observations; $$w_{ij}$$ is the spatial weight between the observation at location i and j and $$z_{j}$$ is the value of the variable at all locations different than i90,91.

A binary row-standardised spatial weight matrix conceptualises spatial regions (neighbourhoods) within a 35 m of radial Euclidian distance from each focal point, obtaining the weights $$w_{ij}$$ for all points in relation to each other. The 35 m distance band has been chosen to include in the neighbourhood sand points from the two adjacent transects. We used 999 random permutations to compute the reference distribution and obtain a minimum pseudo p-value of 0.05 (95% level of confidence), which we used to discard non-significant $$I_{i}$$.

In this analysis, only significant High-High (areas where high values are surrounded by high values) and Low-Low (areas where low values are surrounded by low values) hotspots have been retained, discarding spatial outliers.

Significant hotspots of Δh have been classified into five magnitude classes, using the Jenks-Caspall optimised-natural breaks method92, based on the totality of absolute Δh values (Table 1). Weights of each magnitude class are used to represent the severity of change during the e-BCD sign computations and are obtained deriving the median absolute values for each Δh magnitude class.

### Empirical and residual beachface cluster dynamic indices

The empirical and residual beachface cluster dynamics indices (e-BCD and r-BCD respectively) are purposefully designed metrics to leverage the very high spatiotemporal resolutions and three-dimensionality of our data for studying subaerial beach landform dynamics (morphodynamics). With elevation change (Δh) magnitude classes as transient states (Table 1), we used finite discrete Markov chain models to compute first-order stochastic transition matrices and steady-state probability vectors, used to derive e- and r-BCD respectively.

Following Lambin (1994), a discrete Markov process can be represented as:

$$s_{t + 1} = Ms_{t} ,$$

where $$s_{t}$$ is a column vector, $$s = \left( {s_{1} , \ldots ,s_{m} } \right)$$ having as elements the valid points (within the beachface and beyond LoD sand-only Δh points) in one of the $$m$$ states (i.e. Δh magnitude classes) at time $$t$$. M is a $$m \times m$$ matrix holding the first-order (from $$t$$ to $$t + 1$$) transition probabilities $$p_{ij}$$, derived as:

$$p_{ij} = \frac{{n_{ij} }}{{\mathop \sum \nolimits_{k = 1}^{m} n_{ik} }},$$

where $$n_{ij}$$ is the number of transitions from an initial state $$i$$ to state $$j$$ and $$m$$ is the number of states (i.e. elevation change classes in Table 1) in which each observation can be. The matrix M is row-standardised, so that the sum of transition probabilities from a given state is always equal to one.

The e-BCDs divide the first-order transition matrix M into four sub-matrices (Supplementary Fig. S7), each representing site-level erosional, depositional, recovery and vulnerability behaviours of the subaerial beach over the monitoring period. The e-BCDs indices are computed for every sub-matrices ($$sub$$) as follows:

$$Empirical\, BCD_{sub} = \mathop \sum \limits_{i,j = 1}^{n} \left[ { ws^{\prime}_{ij} } \right]\times p_{ij} ,$$

where $$ws_{ij} ^{\prime}$$ is a transformed weight that reflects the magnitude trend of such transition, which is defined as:

$$ws^{\prime}_{ij} = \left\{ {\begin{array}{*{20}l} { ws_{i} \times\left( { - ws_{j} } \right) \leftrightarrow i > j } \\ {ws_{i} \times ws_{j} \leftrightarrow i \le j } \\ \end{array} } \right\},$$

where $$ws_{i}$$ and $$ws_{j}$$ are the weights (Table 1) related to the initial $$i$$ and transitioning $$j$$ states respectively. The brackets “[ ]” indicate that the $$ws_{ij} ^{\prime}$$ transformation is implemented separately to determine the e-BCD sign only. The e-BCD absolute score computation does not implement this multiplication, capturing the process importance only. Any state to either no-hotspot or no-data transition probabilities are not included in the e-BCD interpretation.

The r-BCDs are computed from the steady-state probability vector.

The steady-state of a Markov chain returns a unique probability vector representing the states limiting probability distribution, which, once attained, one additional (or more) time steps will return the exact same initial states probabilities, signalling a situation of dynamic equilibrium has been achieved. This is represented as:

$$\pi M = \pi ,$$

where $$\pi$$  is the vector containing the limiting probabilities $$\pi_{j}$$ for each $$j$$ state in $$s$$. This vector $$\pi$$ is derived by solving a system of $$m$$ equations with $$m$$ unknowns, each equation represented as:

$$\pi_{j} = \mathop \sum \limits_{k = 1}^{m} \pi_{k} p_{kj} ,$$

given that:

$$\mathop \sum \limits_{j = 1}^{m} \pi_{j} = 1.$$

The steady-state can be seen in a descriptive way as representing the states hierarchy, which is unique to the system being modelled93, from which we derive the stochastic tendency the system had towards depositional or erosional states at the end of monitoring. We interpret this tendency as the most likely behavioural regime the system was subjected to, given the drivers of change and boundary conditions that influenced its evolution during the monitoring period.

The computation of the r-BCDs is as follows:

$$Residual \,BCD_{SS} = 100 \times \mathop \sum \limits_{i = 1}^{m} (\pi_{ie} - \pi_{id} ),$$

where $$ss^{ }$$ is the steady-state probability distribution of one location, $$\pi_{ie}$$ are the limiting probabilities of the erosional classes (ue, se, me, he, ee) and $$\pi_{id}$$ the limiting probabilities of the depositional classes (ud, sd, md, hd and ed) (Table 1). The r-BCDs are not signed as no transitions are represented in the resultant vector. Any state to either no-hotspot or no-data transition probabilities are not included in the r-BCD interpretation. The multiplication by 100 is performed for index readability purposes.

We computed r-BCDs using erosional/depositional hotspots at the site level only (Fig. 3a), while at the transect level (Fig. 4), r-BCDs are computed with the full Δh cross-shore profiles. This has been done due to the relatively narrow beach width which resulted in a limited number of valid sand-only observations for hotspot analysis.

## Data availability

More than 220 3D datasets are already freely accessible to anyone via a user friendly web-platform to share and communicate information, promote coastal awareness, build knowledge and further increase the impact of our efforts. Link https://www.propelleraero.com/. Credentials email: vcmp@deakin.edu.au; Password: propellervcmp.

## Code availability

The code used throughout this study is freely available here: Draft code: https://github.com/npucino/sandpiper/tree/master. Except for digital photogrammetric procedure, all the analysis has been automated using open-source Python 3 packages, especially PySAL 2.0.093 for LISA and Markov Chain computations, Scikit-Learn94 for machine learning classification and Seaborn and Matplotlib95 for quantitative data visualisation. Maps and some minor geoprocessing tasks have been accomplished using QGIS versions 3.0.2 (Madeira) and 3.2.3 (Bonn).

## References

1. Adger, W. N., Hughes, T. P., Folke, C., Carpenter, S. R. & Rockström, J. Social-ecological resilience to coastal disasters. Science 309, 1036–1039 (2005).

2. Martínez, M. L. et al. The coasts of our world: Ecological, economic and social importance. Ecol. Econ. 63, 254–272 (2007).

3. Neumann, B., Vafeidis, A. T., Zimmermann, J. & Nicholls, R. J. Future coastal population growth and exposure to sea-level rise and coastal flooding—A global assessment. PLoS ONE 10, e0118571 (2015).

4. Barbier, E. B. et al. The value of estuarine and coastal ecosystem services. Ecol. Monogr. 81, 169–193 (2011).

5. Defeo, O. et al. Threats to sandy beach ecosystems: A review. Estuar. Coast. Shelf Sci. 81, 1–12 (2009).

6. Luijendijk, A. et al. The state of the World’s beaches. Sci. Rep. 8, 6641 (2018).

7. Mentaschi, L., Vousdoukas, M. I., Pekel, J.-F., Voukouvalas, E. & Feyen, L. Global long-term observations of coastal erosion and accretion. Sci. Rep. 8, 12876 (2018).

8. McNinch, J. E. Geologic control in the nearshore: Shore-oblique sandbars and shoreline erosional hotspots, Mid-Atlantic Bight, USA. Mar. Geol. 211, 121–141 (2004).

9. List, J. H., Farris, A. S. & Sullivan, C. Reversing storm hotspots on sandy beaches: Spatial and temporal characteristics. Mar. Geol. 226, 261–279 (2006).

10. List, J. H., Farris, A. S. & Sullivan, C. Evaluating the persistence of shoreline change hotspots, Northern North Carolina. In American Geophysical Union Fall Meeting, San Francisco. CA, EOS Transactions vol. 83 (2002).

11. Lazarus, E. D., Ellis, M. A., Murray, A. B. & Hall, D. M. An evolving research agenda for human-coastal systems. Geomorphology 256, 81–90 (2016).

12. Burningham, H. & French, J. Understanding coastal change using shoreline trend analysis supported by cluster-based segmentation. Geomorphology 282, 131–149 (2017).

13. Nicholls, R. J. & Cazenave, A. Sea-level rise and its impact on coastal zones. Science 328, 1517–1520 (2010).

14. Vousdoukas, M. I. et al. Sandy coastlines under threat of erosion. Nat. Clim. Chang. 10, 260–263 (2020).

15. Masselink, G. & Lazarus, E. D. Defining coastal resilience. Water 11, 2587 (2019).

16. Andrade, F. & Ferreira, M. A. A simple method of measuring beach profiles. J. Coastal Res. 224, 995–999 (2006).

17. Harley, M. D., Turner, I. L., Short, A. D. & Ranasinghe, R. Assessment and integration of conventional, RTK-GPS and image-derived beach survey methods for daily to decadal coastal monitoring. Coast. Eng. 58, 194–205 (2011).

18. Meng, X., Zhang, X., Silva, R., Li, C. & Wang, L. Impact of high-resolution topographic mapping on beach morphological analyses based on terrestrial LiDAR and object-oriented beach evolution. ISPRS Int. J. Geo-Inf. 6, 147 (2017).

19. Harley, M. D. et al. Extreme coastal erosion enhanced by anomalous extratropical storm wave direction. Sci. Rep. 7, 6033 (2017).

20. Nikolakopoulos, K., Kyriou, A., Koukouvelas, I., Zygouri, V. & Apostolopoulos, D. Combination of aerial, satellite, and UAV photogrammetry for mapping the diachronic coastline evolution: The case of Lefkada island. IJGI 8, 489 (2019).

21. Lentz, E. E., Hapke, C. J., Stockdon, H. F. & Hehre, R. E. Improving understanding of near-term barrier island evolution through multi-decadal assessment of morphologic change. Mar. Geol. 337, 125–139 (2013).

22. Ierodiaconou, D., Schimel, A. C. G. & Kennedy, D. M. A new perspective of storm bite on sandy beaches using Unmanned Aerial Vehicles. Z. Geomorphol. Suppl. 60, 123–137 (2016).

23. Turner, I. L., Harley, M. D. & Drummond, C. D. UAVs for coastal surveying. Coast. Eng. 114, 19–24 (2016).

24. Tmušić, G. et al. Current practices in UAS-based environmental monitoring. Remote Sens. 12, 1001 (2020).

25. Navarro, A. et al. The application of unmanned aerial vehicles (UAVs) to estimate above-ground biomass of mangrove ecosystems. Remote Sens. Environ. 242, 111747 (2020).

26. Carrera-Hernández, J. J., Levresse, G. & Lacan, P. Is UAV-SfM surveying ready to replace traditional surveying techniques?. Int. J. Remote Sens. 41, 4820–4837 (2020).

27. Brunier, G., Fleury, J., Anthony, E. J., Gardel, A. & Dussouillez, P. Close-range airborne Structure-from-motion photogrammetry for high-resolution beach morphometric surveys: Examples from an embayed rotating beach. Geomorphology 261, 76–88 (2016).

28. Cunliffe, A. M. et al. Rapid retreat of permafrost coastline observed with aerial drone photogrammetry. The Cryosphere 13, 1513–1528 (2019).

29. Gonçalves, J. A. & Henriques, R. UAV photogrammetry for topographic monitoring of coastal areas. ISPRS J. Photogramm. Remote Sens. 104, 101–111 (2015).

30. Klemas, V. V. Coastal and environmental remote sensing from unmanned aerial vehicles: An overview. J. Coast. Res. 315, 1260–1267 (2015).

31. Gómez-Pazo, A., Pérez-Alberti, A. & Trenhaile, A. Recording inter-annual changes on a boulder beach in Galicia, NW Spain using an unmanned aerial vehicle. Earth Surf. Process. Landforms 44, 1004–1014 (2019).

32. Long, N., Millescamps, B., Guillot, B., Pouget, F. & Bertin, X. Monitoring the topography of a dynamic tidal inlet using UAV imagery. Remote Sens. 8, 387 (2016).

33. Pagán, J. I., Bañón, L., López, I., Bañón, C. & Aragonés, L. Monitoring the dune-beach system of Guardamar del Segura (Spain) using UAV, SfM and GIS techniques. Sci. Total Environ. 687, 1034–1045 (2019).

34. Ruessink, B. G., Arens, S. M., Kuipers, M. & Donker, J. J. A. Coastal dune dynamics in response to excavated foredune notches. Aeol. Res. 31, 3–17 (2018).

35. Westoby, M. J., Brasington, J., Glasser, N. F., Hambrey, M. J. & Reynolds, J. M. ‘Structure-from-motion’photogrammetry: A low-cost, effective tool for geoscience applications. Geomorphology 179, 300–314 (2012).

36. Bonney, R. et al. Citizen science. Next steps for citizen science. Science 343, 1436–1437 (2014).

37. Newman, G. et al. The future of citizen science: Emerging technologies and shifting paradigms. Front. Ecol. Environ. 10, 298–304 (2012).

38. McKinley, D. C. et al. Investing in citizen science can improve natural resource management and environmental protection. Issues Ecol. 2015, 1–27 (2015).

39. Dickinson, J. L. et al. The current state of citizen science as a tool for ecological research and public engagement. Front. Ecol. Environ. 10, 291–297 (2012).

40. Cohn, J. P. Citizen science: Can volunteers do real research?. Bioscience 58, 192–197 (2008).

41. Silvertown, J. A new dawn for citizen science. Trends Ecol. Evol. 24, 467–471 (2009).

42. Kosmala, M., Wiggins, A., Swanson, A. & Simmons, B. Assessing data quality in citizen science. Front. Ecol. Environ. 14, 551–560 (2016).

43. James, M. R., Robson, S., d’Oleire-Oltmanns, S. & Niethammer, U. Optimising UAV topographic surveys processed with structure-from-motion: Ground control quality, quantity and bundle adjustment. Geomorphology 280, 51–66 (2017).

44. Höhle, J. & Höhle, M. Accuracy assessment of digital elevation models by means of robust statistical methods. ISPRS J. Photogramm. Remote Sens. 64, 398–406 (2009).

45. Casella, E, Drechsel J, Winter, C, Benninghoff, M, Rovere, A.: Accuracy of sand beach topography surveying by drones and photogrammetry. Geo-Marine Letters 40(2), 255–268 (2020).

46. James, M. R. et al. Guidelines on the use of structure-from-motion photogrammetry in geomorphic research. Earth Surf. Processes Landforms 44, 2081–2084 (2019).

47. August, T., Fox, R., Roy, D. B. & Pocock, M. J. O. Data-derived metrics describing the behaviour of field-based citizen scientists provide insights for project design and modelling bias. Sci. Rep. 10, 11009 (2020).

48. Dunn, O. J. Multiple comparisons among means. J. Am. Stat. Assoc. 56, 52–64 (1961).

49. Kruskal, W. H. & Allen Wallis, W. Errata: Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 48, 907 (1953).

50. Jenks, G. F. The data model concept in statistical mapping. Int. Yearb. Cartogr. 7, 186–190 (1967).

51. Dean, D. L., Novianti, S. & Noor, A. A. An assessment of the international and domestic tourists behavior in Australia. Int. J. Appl. Bus. Res. 2, 46–57 (2020).

52. Wright, L. D. & Short, A. D. Morphodynamic variability of surf zones and beaches: A synthesis. Mar. Geol. 56, 93–118 (1984).

53. Miner, T. & Rosengren, N. Apollo Bay Coastal Geomorphology Study. 1–116 (The State of Victoria Department of Environment, Land, Water and Planning, 2019). https://s3.ap-southeast-2.amazonaws.com/hdp.au.prod.app.vic-engage.files/3315/7438/2555/APOLLO_BAY_COASTAL_GEOMORPHOLOGY_STUDY_NOV_2019.pdf, [Accessed the 5th March 2020].

54. Gonçalves, G. R., Pérez, J. A. & Duarte, J. Accuracy and effectiveness of low cost UASs and open source photogrammetric software for foredunes mapping. Int. J. Remote Sens. 39, 5059–5077 (2018).

55. Wright, L. D. & Thom, B. G. Coastal depositional landforms: A morphodynamic approach. Progr. Phys. Geogr. Earth Environ. 1, 412–459 (1977).

56. Chen, W.-W. & Chang, H.-K. Estimation of shoreline position and change from satellite images considering tidal variation. Estuar. Coast. Shelf Sci. 84, 54–60 (2009).

57. Hagenaars, G., de Vries, S., Luijendijk, A. P., de Boer, W. P. & Reniers, A. J. H. M. On the accuracy of automated shoreline detection derived from satellite imagery: A case study of the sand motor mega-scale nourishment. Coast. Eng. 133, 113–125 (2018).

58. Selvan, S. C., Kankara, R. S., Prabhu, K. & Rajan, B. Shoreline change along Kerala, south-west coast of India, using geo-spatial techniques and field measurement. Nat. Hazards 100, 17–38 (2020).

59. Kelly, J. T. & Gontz, A. M. Using GPS-surveyed intertidal zones to determine the validity of shorelines automatically mapped by Landsat water indices. Int. J. Appl. Earth Obs. Geoinf. 65, 92–104 (2018).

60. Vos, K., Harley, M. D., Splinter, K. D., Walker, A. & Turner, I. L. Beach slopes from satellite-derived shorelines. Geophys. Res. Lett. 47, 502 (2020).

61. Sagar, S., Roberts, D., Bala, B. & Lymburner, L. Extracting the intertidal extent and topography of the Australian coastline from a 28 year time series of Landsat observations. Remote Sens. Environ. 195, 153–169 (2017).

62. Bishop-Taylor, R., Sagar, S., Lymburner, L. & Beaman, R. J. Between the tides: Modelling the elevation of Australia’s exposed intertidal zone at continental scale. Estuar. Coast. Shelf Sci. 223, 115–128 (2019).

63. Khan, M. J. U. et al. High-resolution intertidal topography from sentinel-2 multi-spectral imagery: synergy between remote sensing and numerical modeling. Remote Sens. 11, 2888 (2019).

64. Beuzen, T., Harley, M. D., Splinter, K. D. & Turner, I. L. Controls of variability in berm and dune storm erosion. J. Geophys. Res. Earth Surf. 124, 2647–2665 (2019).

65. Splinter, K. D., Kearney, E. T. & Turner, I. L. Drivers of alongshore variable dune erosion during a storm event: Observations and modelling. Coast. Eng. 131, 31–41 (2018).

66. Shrestha, R. L., Carter, W. E., Sartori, M., Luzum, B. J. & Slatton, K. C. Airborne laser swath mapping: Quantifying changes in sandy beaches over time scales of weeks to years. ISPRS J. Photogramm. Remote Sens. 59, 222–232 (2005).

67. Faivre, G. et al. Coastal processes within a coral reef Lagoon system: Erakor Lagoon, Efate island. Vanuatu. J. Coast. Res. 95, 1427–1432 (2020).

68. Bove, G., Becker, A., Sweeney, B., Vousdoukas, M. & Kulp, S. A method for regional estimation of climate change exposure of coastal infrastructure: Case of USVI and the influence of digital elevation models on assessments. Sci. Total Environ. 710, 136162 (2020).

69. Powers, M., Begg, Z., Smith, G. & Miles, E. Lessons from the pacific ocean portal: Building pacific island capacity to interpret, apply, and communicate ocean information. Front. Mar. Sci. https://doi.org/10.3389/fmars.2019.00476 (2019).

70. Hart, J. & Blenkinsopp, C. Using citizen science to collect coastal monitoring data. J. Coast. Res. 95, 824–828 (2020).

71. Harley, M. D., Kinsela, M. A., Sánchez-García, E. & Vos, K. Shoreline change mapping using crowd-sourced smartphone images. Coast. Eng. 150, 175–189 (2019).

72. Dai, W. et al. UAV photogrammetry for elevation monitoring of intertidal mudflats. J. Coast. Res. 85, 236–240 (2018).

73. Neugirg, F. et al. Erosion processes in calanchi in the Upper Orcia Valley, Southern Tuscany, Italy based on multitemporal high-resolution terrestrial LiDAR and UAV surveys. Geomorphology 269, 8–22 (2016).

74. Pineux, N. et al. Can DEM time series produced by UAV be used to quantify diffuse erosion in an agricultural watershed?. Geomorphology 280, 122–136 (2017).

75. Duró, G., Crosato, A., Kleinhans, M. G. & Uijttewaal, W. S. J. Bank erosion processes measured with UAV-SfM along complex banklines of a straight mid-sized river reach. Earth Surf. Dyn. 6, 933 (2018).

76. Esposito, G., Mastrorocco, G., Salvini, R., Oliveti, M. & Starita, P. Application of UAV photogrammetry for the multi-temporal estimation of surface extent and volumetric excavation in the Sa Pigada Bianca open-pit mine, Sardinia, Italy. Environ. Earth Sci. 76, 103 (2017).

77. Rossi, G. et al. Multitemporal UAV surveys for landslide mapping and characterization. Landslides 15, 1045–1052 (2018).

78. Cucchiaro, S. et al. Geomorphic effectiveness of check dams in a debris-flow catchment using multi-temporal topographic surveys. CATENA 174, 73–83 (2019).

79. Bendig, J., Bolten, A. & Bareth, G. UAV-based imaging for multi-temporal, very high resolution crop surface models to monitor crop growth variability. Photogramm.—Fernerkundung—Geoinfor. 2013, 551–562 (2013).

80. Guerra-Hernández, J. et al. Use of multi-temporal UAV-derived imagery for estimating individual tree growth in Pinus pinea Stands. For. Trees Livelihoods 8, 300 (2017).

81. Fernandes, R. et al. Monitoring snow depth change across a range of landscapes with ephemeral snowpacks using structure from motion applied to lightweight unmanned aerial vehicle videos. Cryosphere 12, 3535–3550 (2018).

82. Rossini, M. et al. Rapid melting dynamics of an alpine glacier with repeated UAV photogrammetry. Geomorphology 304, 159–172 (2018).

83. Incekara, A. H., Delen, A., Seker, D. Z. & Goksel, C. Investigating the utility potential of low-cost unmanned aerial vehicles in the temporal monitoring of a landfill. ISPRS Int. J. Geo-Inf. 8, 22 (2019).

84. Messinger, M. & Silman, M. Unmanned aerial vehicles for the assessment and monitoring of environmental contamination: An example from coal ash spills. Environ. Pollut. 218, 889–894 (2016).

85. Ćwiąkała, P. et al. Assessment of the Possibility of using unmanned aerial vehicles (UAVs) for the documentation of hiking trails in alpine areas. Sensors 18, 81 (2017).

86. Rousseeuw, P. J. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 20, 53–65 (1987).

87. Lloyd, S. Least squares quantization in PCM. IEEE Trans. Inf. Theory 28, 129–137 (1982).

88. Lane, S. N., Westaway, R. M. & Murray Hicks, D. Estimation of erosion and deposition volumes in a large, gravel-bed, braided river using synoptic remote sensing. Earth Surf. Process. Landforms 28, 249–271 (2003).

89. Mancini, F. et al. Using unmanned aerial vehicles (UAV) for high-resolution reconstruction of topography: The structure from motion approach on coastal environments. Remote Sens. 5, 6880–6898 (2013).

90. Anselin, L. Local indicators of spatial association—LISA. Geogr. Anal. 27, 93–115 (1995).

91. Yuan, Y., Cave, M. & Zhang, C. Using Local Moran’s I to identify contamination hotspots of rare earth elements in urban soils of London. Appl. Geochem. 88, 167–178 (2018).

92. Jenks, G. F. & Caspall, F. C. Error on choroplethic maps: Definition, measurement, reduction. Ann. Assoc. Am. Geogr. 61, 217–244 (1971).

93. Brown, L. A. On the use of Markov chains in movement research. Econ. Geogr. 46, 393–403 (1970).

94. Rey, S. J. & Anselin, L. PySAL: A Python Library of Spatial Analytical Methods 175–193 (Springer, Berlin, 2010).

95. Pedregosa, F. et al. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).

96. Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 9, 90–95 (2007).

Download references

## Acknowledgements

Funding for this project was provided by the State of Victoria Department of Environment, Land, Water and Planning (DELWP) under the Victorian Coastal Monitoring Program, Deakin University and The University of Melbourne. Apollo Bay nourishment data has been provided by Hannah Fallon (DELWP). We are grateful to Karina Sorrell and all the citizen scientists involved in the collection of the data.

## Author information

Authors

### Contributions

N.P.: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing—Original Draft, Visualization. D.K.: Writing—Review & Editing, Supervision, Project administration, Funding acquisition. R.C.: Writing—Review & Editing, Supervision. B.A.: Writing—Review & Editing, Data Curation. D.I.: Writing—Review & Editing, Supervision, Project administration, Funding acquisition.

### Corresponding author

Correspondence to Nicolas Pucino.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

### Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

## About this article

### Cite this article

Pucino, N., Kennedy, D.M., Carvalho, R.C. et al. Citizen science for monitoring seasonal-scale beach erosion and behaviour with aerial drones. Sci Rep 11, 3935 (2021). https://doi.org/10.1038/s41598-021-83477-6

Download citation

• Received:

• Accepted:

• Published:

• DOI: https://doi.org/10.1038/s41598-021-83477-6

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

## Search

### Quick links

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing